EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 708 (April 2026) A&A, 708 (2026) A279 Full HTML
Open Access
Issue
A&A
Volume 708, April 2026
Article Number A279
Number of page(s) 17
Section Interstellar and circumstellar matter
DOI https://doi.org/10.1051/0004-6361/202659027
Published online 16 April 2026

© The Authors 2026

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model.

Open Access funding provided by Max Planck Society.

1 Introduction

Pre-stellar cores are dense, gravitationally bound fragments of molecular clouds that are dynamically unstable and on their way to form stars or brown dwarfs. The study of the physical and chemical structure of pre-stellar cores therefore provides important constraints on the initial conditions of the star and planet formation process.

Previous studies of pre-stellar cores have revealed a rich chemistry, with a broad range of molecules detected, ranging from simple diatomic molecules and key ions (CO and HCO+) to complex organic molecules (COMs) such as CH3OH, CH3OCH3, and CH3CHO (e.g., Jiménez-Serra et al. 2016). These studies indicated that a significant fraction of the molecules available during planet formation may already be formed in the molecular cloud phase or during the pre-stellar core phase (Caselli & Ceccarelli 2012; Booth et al. 2021; Ceccarelli et al. 2023). This is also supported by isotopic fractionation studies, which generally support a high degree of chemical inheritance (Jensen et al. 2021b; Drozdovskaya et al. 2021; Nomura et al. 2023; Leemker et al. 2025; Scibelli et al. 2025).

Recent studies have revealed a clear spatial segregation between different classes of molecules in pre-stellar cores (Spezzano et al. 2017, 2020). The origin of the observed chemical structure is not yet understood. One hypothesis is that the observed segregation arises from the local cloud structure and environment. Star formation occurs across a wide range of different environments ranging from isolated Bok globules to dense star-forming regions with tens or hundreds of stars forming within parsec-scale regions (1-10 parsec). These diverse environmental conditions may affect the star formation process and the chemistry that evolves through a number of mechanisms: dynamical effects on the accretion process (density perturbations and differing timescales), and uneven irradiation through either structural differences in the core or cloud environment or because of local irradiation by nearby OB stars. Another effect that might contribute to the chemical structure of prestellar cores is the dynamic nature of the accretion process itself. Recently, a growing number of accretion streamers have been detected toward pre-stellar cores (Choudhury et al. 2025) and protostellar sources (Pineda et al. 2020; Gieser et al. 2024). Streamers may add fresh material with a distinct chemical composition, having undergone different physical and chemical histories (Kuffmeier et al. 2016). To which extent these effects affect the chemical composition and structure of pre-stellar cores and the later stages of star formation is still an open question. In a recent study, Giers et al. (2025) applied a density-based clustering analysis to molecular emission maps toward three young cores. They reported a clear link between the H2 column density N(H2 ), the gradient in N(H2), and the chemical morphology and associate this with environmental effects around the cores that shape the chemical morphology.

To assess the effects of the local environment, models need to account for realistic cloud environments, which resemble the stellar nurseries we observe in the local galaxy. Furthermore, models need to account for a broad range of physics and chemistry, including magnetohydrodynamics (MHD), grain-surface chemistry, and photochemistry. Because chemical inheritance likely plays a key role in the chemistry during star and planet formation, it is crucial to assess to which degree the local cloud environment might affect the early chemistry, and hence, ultimately the chemistry of the later planetary system.

While astrochemical models serve as an important tool for understanding observations, they regularly struggle to provide accurate predictions across multiple different species (e.g., Ruaud et al. 2015; Vasyunin et al. 2017). In these cases, the physical and chemical parameters may have to be fine-tuned to specific species or sub-categories of molecules to provide a good agreement with observed line intensities (e.g., Maret et al. 2013; Giers et al. 2023; Harju et al. 2025). These fundamental challenges may be a result of the many simplifications in astrochemical models. For example, rate equations for grainsurfaces reactions are crude approximations of the stochastic processes involved, and many parameters in these models are highly uncertain (Caselli et al. 1998; Cuppen et al. 2017). However, the apparent discrepancies in astrochemical models might also arise from limitations in the physical modeling. Often, 0D or 1D models are preferred to 3D models because they are simpler to develop and run and are often more straightforward to interpret (e.g., Garrod 2013; Borshcheva et al. 2025). It is essential to identify the primary source of limitations in current astrochemical models, regardless of whether they stem from the choice of physical model or from challenges such as poorly constrained chemical parameters and oversimplifications in chemical modeling. This understanding will guide the prioritization of the most critical problems to solve to improve the model.

We present a 3D MHD simulation of a pre-stellar core formed in a larger molecular cloud environment. The physical simulation is post-processed to compute the local temperature, radiation field, and the cosmic-ray ionization rates self-consistently around the core. The chemical evolution is then computed along tracer particle trajectories that accrete onto the core. This allows us to model the dynamical evolution of the matter that is accreted into the pre-stellar core in a realistic cloud environment. We compare the model with a previous static model of the same 3D core and to the numerous observations of the pre-stellar core L1544, located in the Taurus molecular cloud.

The paper is organized as follows. In Sect. 2, we introduce the model. In Sect. 3 we present the results of the model and compare them with the previous static model of the well-studied pre-stellar core L1544 and with observations of the same core. In Sect. 4, we discuss the implications of this work for future models, the comparison with L1544, and the next steps. In Sect. 5, we summarize our results.

Thumbnail: Fig. 1 Refer to the following caption and surrounding text. Fig. 1

Projected H2 column densities in the molecular cloud simulation. The white crosses show the location of protostars formed in this snapshot from the simulation (t = 305 kyr). The red diamond marks the protostar that is studied here during the pre-stellar phase. The pre-stellar core studied in this work is centered in the figure for clarity because the box has periodic boundary conditions.

2 Model description

The overall framework of the model we used was previously introduced for the static model in Jensen et al. (2023). In this section, we provide a brief overview of the central elements of the model and a more in-depth description of the changes implemented in this work.

2.1 Physical model

2.1.1 Molecular cloud simulation

The underlying physical simulation is similar to the models presented in Haugbølle et al. (2018), but the linear resolution is twice higher, and there are eight times more cells overall to resolve the stellar environment and turbulent cascade better. It is a 3D MHD simulation of a 4 pc3 star-forming molecular cloud using the code RAMSES with a number of modifications (Teyssier 2002). Star formation in the cloud is driven by an interplay between supersonic turbulence and gravitational collapse. The supersonic turbulence is maintained by solenoidal driving on the largest scale, resulting in an RMS velocity of 2 km/s in the cloud. The simulation reproduces the empirically derived core mass function and initial mass function to a high degree (Haugbølle et al. 2018; Pelkonen et al. 2021). The overall characteristic of the cloud simulation is therefore comparable to the observed star formation process. To create realistic initial conditions and erase any trace of the homogeneous initial conditions, the box was first evolved with pure MHD for 20 turbulent turnover times (20 Myr), and only then were self-gravity and the star formation recipe turned on. The simulation uses adaptive mesh refinement, and the smallest cells are 25 au across. The entire model domain is shown in Fig. 1. The RAMSES simulation has already been used in numerous studies to investigate late-stage infall (Kuffmeier et al. 2023), chemistry in pre-stellar cores and embedded protostars (Jensen et al. 2021a, 2023), angular momentum in Class II systems (Pelkonen et al. 2025), early Bondi-Hoyle accretion (Padoan et al. 2025), and RAMSESbased zoom-in studies within molecular clouds (Jørgensen et al. 2022; Tuhtan et al. 2024; Al-Belmpeisi et al. 2024).

The model is isothermal and does not consider radiation processes. These are included by post-processing snapshots from the simulation. Each RAMSES snapshot was post-processed using RADMC3D (Dullemond et al. 2012) to estimate the local radiation field within the model. This step was made on a subset of the full simulation domain with a radius of 50 000 au, centered on the pre-stellar core, for computational efficiency. From the RADMC-3D simulations, we derived the dust temperature and the local radiation field, which we used to derive the local visual extinction Av,eff. We denote this as the effective visual extinction to distinguish it from the line-of-sight extinction as seen by observers.

The radiation source in the RADMC-3D simulation was the interstellar radiation field (ISRF) surrounding the box. The ISRF was parameterized following Hocuk et al. (2017). No internal radiation source was used, as the core is pre-stellar and does not yet harbor a hydrostatic core. The ISRF was attenuated by a visual extinction of the ambient cloud of either Av,effambMathematical equation: $A_\mathrm{v,eff}^\mathrm{amb}$ = 1 mag or 2 mag, consistent with simulations presented in Sipilä & Caselli (2018) and Sipilä et al. (2022). This attenuation was intended to mimic the effects of the core being embedded in a larger ambient molecular cloud that provides additional shielding against the ISRF.

We assumed perfect thermal coupling between the gas and dust particles. This assumption is valid at higher densities (≥104−5 cm−3), which characterize the core itself, but not accurate in the outer diffuse boundary layers (Goldsmith 2001). Nevertheless, the differences between the dust and gas temperatures are expected to remain modest (≤2 K) at the edge of cores as well (Sipilä et al. 2017).

Jensen et al. (2023) selected a snapshot that matched the conditions of L1544 best as derived from observations: a central density of 106 cm−3 within the high-density kernel with a radius of 1800 au (Caselli et al. 2019). Because this snapshot occurred 300 kyr after tracer particles were injected and star formation enabled in the simulation, we were able to study the evolution of each tracer particle for a total of 300 kyr before the snapshot in question. In Appendix A we visualize the physical structure of the anisotropic core using eight radial profiles.

2.1.2 Dynamical evolution tracks using tracer particles

In the static model of the core, the chemistry was evolved for a given period of time at fixed physical conditions. Such models are broadly used to interpret observations of cores (e.g., Jiménez-Serra et al. 2016; Lu et al. 2018; Redaelli et al. 2021; Riedel et al. 2025; Borshcheva et al. 2025). However, given the highly dynamical nature of star-forming regions and the nonspherical accretion processes, it is important to improve our understanding of how the accretion dynamics affect the chemical evolution and the resulting chemical structure. We used tracer particles to capture the physical evolution of the gas that is accreted onto the core (Coutens et al. 2020; Ferrada-Chamorro et al. 2021; Jensen et al. 2021a; Panessa et al. 2023; Priestley et al. 2023; Navarro-Almaida et al. 2024). Tracer particles are passive, massless particles in the MHD simulations that are advected along with the local gas velocity field. During the simulation, the tracer particles record (trace) the local physical conditions, namely the local density and velocity field. The local temperature and visual extinction were derived from the RADMC-3D post-processing as described in the previous section and applied to each tracer particle trajectory. Hence, the trajectories have the necessary information for computing the chemical history: density, temperature, and irradiation. Figure 2 shows the tracer particles (red dots) accreted onto the pre-stellar core in four different snapshots. Examples of the physical evolution for five tracer particles are shown in Appendix B.

2.2 Chemical model

We employed a two-phase gas-grain chemical model based on the framework presented in Jensen et al. (2023) and initially introduced in Jensen et al. (2021a). Most of the chemical network and parameters remained unchanged from previous works, and we here limit ourselves to a brief introduction of the model and the updates applied in this work. The gas-phase chemical network was based on the network initially presented in Majumdar et al. (2017), based on KIDA2014 (Wakelam et al. 2015), with the addition of deuterated species and spin chemistry for light hydrogen-bearing species (H2, H2+, H3+, and deuterated isotopologs). The gas-phase network was updated to use corrected rates for the H3+ + H2 reactions following Sipilä et al. (2017). Furthermore, photodissociation and photoionization rates were updated based on Heays et al. (2017). The full network and chemical model is available for the community on GitHub1. The grain-surface chemical network was derived from Furuya et al. (2015). We updated the grain-surface formation routes for methanol (CH3OH) and deuterated isotopologs to follow the activation barriers presented in Riedel et al. (2023). The model considers grain-surface reactions via the Langmuir-Hinshelwood mechanism and desorption through chemical reactions, grainheating by cosmic rays, and UV photons. Cosmic-ray-induced grain heating is implemented following Leger et al. (1985). Chemical desorption follows the description from Garrod et al. (2007) with an efficiency factor of 0.01. For CH3OH, nonthermal desorption primarily proceeds through chemical desorption. The assumed photodesorption yield for CH3OH was γPD = 10−6, following the results from Bertin et al. (2016) for ice mixtures with CH3OH and CO. The gas-phase network included the same reaction types as in the standard KIDA network (three-body reactions were excluded). The model included self-shielding of H2, CO, HD, and N2 (Visser et al. 2009; Wolcott-Green & Haiman 2011; Wolcott-Green et al. 2011; Heays et al. 2014). The column densities used to calculate the self-shielding factors were not computed self-consistently due to computational constraints. Instead, N(H) was derived from the observational relation N(H) = Av/1.59 × 1021 cm−2 (Diplas & Savage 1994), where we assumed Av ≈ Av,eff. To estimate the N2 and CO selfshielding factors, the local gas-phase ratio with respect to H2 was used to estimate the column densities. The complete network contained ~950 species and ~50 000 reactions.

To improve the self-consistency of the model, we adopted the parametric description of the cosmic-ray attenuation as a function of H2 column density presented in Padovani et al. (2018). Specifically, we opted for the low model. The H2 column density was computed from the local visual extinction, Av,eff, similarly to the approach described for the self-shielding factors. In the previous model, presented in Jensen et al. (2023), a fixed cosmic-ray ionization rate of ζ = 1.3 × 10−17 s−1 was assumed throughout the core.

We ran models with two different sets of initial elemental abundances, as listed in Table 1. One set of abundances was similar to the abundances used in Jensen et al. (2023) and is labeled LM for “low metal.” The second set contained the high-metal (HM) abundances from Wakelam & Herbst (labeled EA2, 2008). These initial abundances are closer to the cosmic abundances and represent the abundances expected at the earliest stages of molecular clouds, prior to any chemical evolution, which may cause depletion of heavier elements (S, Si, Mg, Fe, and Cl) in the gas phase due to accretion onto dust grains or incorporation into larger molecules.

For the chemical simulations, several different initial conditions were tested. Our fiducial model included a molecular cloud phase with static physical conditions for 1 Myr with nH = 1000 cm−3, Tgas= Tdust= 14 K, and Av,effambMathematical equation: $A_\mathrm{v,eff}^\mathrm{amb}$ = 2 mag. Other models were tested with a lower ambient extinction of 1 mag, higher Tgas/dust = 15 K, and different durations of this static phase, including models without any static phase before the tracer particle evolution. After the initial static phase, the tracer particle phase was computed for 300 kyr of the tracer trajectories. The ambient extinction assumed during the initial phase was also applied during the tracer particle phase, that is, Av,eff=Av,effamb+Av,effMathematical equation: $A_\mathrm{v,eff}$=$A_\mathrm{v,eff}^\mathrm{amb}$+$A_\mathrm{v,eff}$. One set of models was also run with a higher density of nH = 104 cm−3, a higher ambient extinction of Av,effambMathematical equation: $A_\mathrm{v,eff}^\mathrm{amb}$ = 4 mag, and a lower temperature of 10 K. For this model, the ambient extinction for the tracer particle trajectories remained at Av,effambMathematical equation: $A_\mathrm{v,eff}^\mathrm{amb}$ = 2 mag, which provides a better fit with the predicted temperature profile of L1544 (see Fig. E.1). An overview of the different initial conditions is presented in Table 2.

Thumbnail: Fig. 2 Refer to the following caption and surrounding text. Fig. 2

Four snapshots showing the accretion onto the pre-stellar core at different stages of the simulation. The red dots mark the position of the tracer particle accreted onto the pre-stellar core. Only 20% of the tracer particles are plotted for clarity. The color map shows the projected column density, similar to Fig. 1. Zero-age is defined as the onset of star formation in the simulation.

2.3 Radiative transfer modeling: Synthetic images and spectral cubes

The final step of the model presented here involved the nonlocal thermodynamic equilibrium radiative transfer simulation using the code LIME (Brinch & Hogerheijde 2010). The computation of the radiative transfer for selected molecular transitions allowed us to directly compare the model with observations of L1544. The LIME models solve the radiative transfer equation in a domain with a radius of ~50 000 au around the core. The physical structure of the core was based on the RAMSES and RADMC-3D simulations, which provided the density, velocity, and temperature fields. Thermal turbulence was included self-consistently for each molecule based on the local gas temperature. As in the chemical modeling, we assumed Tgas = Tdust. The collisional rates for the molecules presented here were accessed through the LAMDA (Schöier et al. 2005)2 and EMEAA3 databases. A complete list of the references for the collisional rates is provided in Appendix D. The output of LIME models are spectral cubes with a predefined bandwidth and spectral resolution. The distance to the modeled core was set to 170 pc, similar to the distance to L1544 (Galli et al. 2019). The spectral cube was convolved with a Gaussian beam with a beam width similar to that of the IRAM 30 m telescope (~30″) that was used to observe L1544 in Spezzano et al. (2016) and Ferrer Asensio et al. (2025). Hence, a direct comparison between the convolved spectral cube and observations of L1544 is possible.

Table 1

Initial elemental abundances for the LM and HM models.

Table 2

Initial conditions for the chemical simulation before the calculation of the evolution during the tracer particle phase.

3 Results

3.1 Studying the chemical morphology of c-C3H2 and CH3OH

Multiple studies have revealed distinct morphological patterns among different molecular groups in young cores (Spezzano et al. 2017, 2020). These spatial differences appear to be linked to the formation or destruction mechanism of the molecules. A key example is the contrast between methanol and carbon-chain molecules, which follow fundamentally different formation pathways: methanol forms on dust grain surfaces via hydrogenation of CO ice, while carbon-chain molecules are primarily synthesized in the gas phase from atomic carbon and are typically suppressed when the majority of carbon is converted into CO. Following Spezzano et al. (2017), we relied on c-C3 H2 and CH3OH as tracers of carbon-chain and COMs, respectively.

Figure 3 shows a comparison between the molecular emission of c-C3H2 and CH3OH and the physical structure of the core for model #3 with higher metallicities. The figure divides the model into three principal planes (x-y, x-z, and z-y). The transitions modeled in this figure are identical to those observed in Spezzano et al. (2016, 2017). We can therefore directly compare the observed morphologies and intensities. In the x-y plane and x-z planes, the peak emission and contours of CH3OH are offset from the dust-continuum peak position, while the c-C3H2 emission peak coincides with the dust peak. The offsets of CH3OH coincide with the more shielded regions of the core, as shown by the Av,eff and Tdust slice plots. This is consistent with the morphology observed toward L1544 for these two species in L1544. In the third principal plane, no clear offset is seen, and the contours of both molecules trace the denser and more shielded region of the core. The static model presented in Jensen et al. (2023) was able to qualitatively reproduce the observed morphology in L1544 for c-C3H2 and CH3OH. Given the updated chemical model in this work, we reran the static model and compared the new results with the dynamical model presented in this work. Figure 4 shows the morphology of the static model with the updated chemical model in the top panel. The static model has changed morphology since the one presented in Jensen et al. (2023): the methanol emission is now centered in the dust continuum peak in all planes. This change is primarily due to the updated methanol formation pathways, which include new activation energies from Riedel et al. (2023), one additional abstraction reaction, and the correction of branching ratios, which have altered the efficiency of methanol formation in the high-density kernel of the core. The bottom panel shows the results for dynamical model #3 with standard metallicities. Similarly to the model presented in Fig. 3, the CH3OH emission is offset from the dust peak toward the regions with higher extinction.

A direct comparison between static and dynamical models under identical conditions is not directly possible: the static model was evolved for a fixed time at fixed physical conditions, while the dynamical model considered an initial static phase followed by ~300 kyr of dynamical evolution as the gas accretes toward the core. In Fig. 4 we chose to compare the fiducial model presented in Jensen et al. (2023, tstatic = 106 yr, Av,effambMathematical equation: $A_\mathrm{v,eff}^\mathrm{amb}$) with model #3 LM in this work (tinit = 106 yr, Av,effambMathematical equation: $A_\mathrm{v,eff}^\mathrm{amb}$ = 2 mag), but we note that the two models are fundamentally different in terms of physical and chemical evolution. In the static model, tstatic denotes the duration of the chemical simulation. Hence, in the fiducial model from Jensen et al. (2023), the chemical evolution was modeled for a duration of 1 Myr.

An overview of the morphologies in all 18 dynamical models is presented in Fig. 5. Only two models show a significant difference for the morphologies: model #9 LM and HM show a different morphology for c-C3H2, which under these conditions traces less dense clumps of gas that are offset from the core itself.

To quantify the observed patterns, we computed the column density maps for the two molecules in each plane and compared them with the integrated column of H2, Av,eff, and Tgas4. While the synthetic maps from the LIME models allow a direct comparison between observations, studying the column densities provides an alternative view of the underlying chemical structure along the line of sight, independent of excitation effects. To determine to which degree the different molecules correlate with the physical conditions along the line of sight, we computed the Pearson correlation coefficients for the molecular column densities and the visual extinction, local gas density, and the gas temperature. The results for two models are reported in Tables 3 and 4. As expected, the results show a strong correlation between NH2 and the integrated column of the visual extinction N̂(Av,eff), while an anticorrelation is seen for the integrated column of the gas temperature. In general, N(CH3OH) and N̂(c-C3H2) show positive correlations with N(H2) and N(Av,eff). For CH3OH, this is expected because the formation pathway depends on the freeze-out of CO, which requires higher densities and higher extinction. For c-C3H2, the results are more surprising because the formation of c-C3H2 also occurs efficiently at lower densities and extinction. The correlation coefficients between N(CH3OH) and N(c-C3H2) with N(Av,eff) are shown in Fig. 6 as a function of tinit. A general trend is present in the correlation coefficients between c-C3 H2, CH3OH, and Av,eff. The models with lower Av,effambMathematical equation: $A_\mathrm{v,eff}^\mathrm{amb}$ = 1 mag (#5-#8) show a stronger correlation between methanol and the integrated extinction (N̂(Av,eff)), while the correlation is weaker for models with Av,effambMathematical equation: $A_\mathrm{v,eff}^\mathrm{amb}$=2 mag (#1-#4). Additionally, the correlation is more pronounced in models with a shorter initial phase, representing a younger chemistry. The trend is visualized in Fig. 6. We interpret this trend as a result of the minimum requirements for CH3 OH formation. Tracer particles need sufficient time at higher density and Av,eff (n(H) ≳ 104, Av,eff ≳ 3) for CH3 OH formation to proceed. In models with lower Av,effambMathematical equation: $A_\mathrm{v,eff}^\mathrm{amb}$, there is a stronger dependence on the individual tracer particle trajectories. Similarly, in models with a shorter tstatic, the dependence on the individual tracer particle trajectories is enhanced.

It should be noted that the abundance of c-C3H2 is very low in models with a more evolved chemistry (#2-#4). Hence, the correlation with high N̂(Av,eff) and NH2 is a result of the very weak emission that is confined to the region with the highest densities. For methanol, on the other hand, this is not the case. Methanol is weaker in the models with a younger chemistry (#5-#6), which shows the strongest correlation with the integrated H2 column density and N̂(Av,eff).

Thumbnail: Fig. 3 Refer to the following caption and surrounding text. Fig. 3

LM model with 106 yr static phase and Av,effambMathematical equation: $A_\mathrm{v,eff}^\mathrm{amb}$ = 2 mag (#3 LM). Top : integrated intensity levels for CH3 OH (white contours) and c-C3 H2 (red contours) along the three principal axes of the simulation. The color map shows the continuum emission at 1.1 mm. Middle : slice plot showing the local log Av,eff for each principal axes. The upper range of the color map is limited to 10 mag for clarity. Bottom : slice plot showing the local dust temperature for each principal axes. The red and black contours show the 10 K and 7 K limits, respectively.
Thumbnail: Fig. 4 Refer to the following caption and surrounding text. Fig. 4

Comparison of the fiducial static model and the fiducial dynamic model, both with the updated network from this work. Top : static model with LM abundances. The contours show the integrated-intensity levels for CH3OH (white) and c-C3H2 (red) along the three principal axes of the simulation. The color map shows the continuum emission at 1.1 mm. Bottom : similar to the top panel but for dynamical model #7 LM.
Thumbnail: Fig. 5 Refer to the following caption and surrounding text. Fig. 5

Comparison of the 18 dynamical models. The red and white contours indicate the integrated emission contours for c-C3H2 and CH3OH, respectively. The background shows the continuum map of the core at 1.1 mm. The c-C3H2 emission generally peaks at the dust continuum peak, while CH3OH is offset toward the more shielded part of the core.
Table 3

Pearson correlation coefficients for model #3 LM in (x-y, x-z, z-y) planes.

Table 4

Pearson correlation coefficients for model #6 LM in (x-y, x-z, z-y) planes.

Thumbnail: Fig. 6 Refer to the following caption and surrounding text. Fig. 6

Correlation between the column densities of c-C3 H2, CH3 OH, and Av,eff for models #1-#8 for standard and higher metallicities. With Av,effambMathematical equation: $A_\mathrm{v,eff}^\mathrm{amb}$=1 mag, CH3OH shows a stronger correlation, especially at lower tinit. For Av,effambMathematical equation: $A_\mathrm{v,eff}^\mathrm{amb}$=2 mag, the correlation is lower for CH3OH and stronger for c-C3H2
Thumbnail: Fig. 7 Refer to the following caption and surrounding text. Fig. 7

Line profiles as predicted by the fiducial model with LM abundances, 106 yr static phase, and Av,effambMathematical equation: $A_\mathrm{v,eff}^\mathrm{amb}$ = 2 mag. The spectra were extracted toward the continuum peak in the x-y plane. A fixed 16O/18O ratio of 560 was used for the C18O abundance profiles.

3.2 Comparing the molecular intensities of L1544 and the model

Figures 7 and 8 shows synthetic spectra for selected molecules: c-C3H2, CH3OH, SO, HCN, HCO+, CS, C18O, and N2H+. All spectra were extracted toward the peak in the dust continuum in the x-y plane. The behavior of the lines was qualitatively the same in all planes.

To evaluate how well the model reproduces the observed molecular intensities toward the dust peak of L1544, we focused not only on absolute values, but primarily on the relative differences between species. Achieving exact agreement in absolute intensities, even with a chemically perfect model, would require the physical model to precisely replicate the density structure, temperature, irradiation environment, and evolutionary stage of L1544. Because the modeled core differs from the real L1544 in these aspects, a one-to-one comparison of absolute intensities is not meaningful scientifically. Instead, by examining the relative differences in intensities for the different molecules, we identified the species that are simultaneously well reproduced by the model and those that require further refinement, offering a more robust and insightful assessment of the model performance. For L1544, we compared our results with the observations of CS, SO, and HCO+ reported in Ferrer Asensio et al. (2025), c-C3H2 and CH3OH reported in Spezzano et al. (2016), N2H+ reported in Redaelli et al. (2019), and HCN (priv. comm.). All observations were carried out with the IRAM 30m telescope at ~3 mm.

Figures 9 and 10 show the ratio of the peak emission observed toward the continuum emission peak of L1544 and the values extracted from the model for models #1-#8 with LM and HM abundances. Several trends are notable. The intensity of the HCO+ 1-0 transitions are well reproduced by models with standard metallicities. For the HM models, on the other hand, the intensities are too weak by factors of 4-10. The intensity of the central component of the N2H+ 1-0 hyperfine transition is weaker by factors of ≥10 for the LM and HM models. The intensity of the central component of the HCN 1-0 hyperfine transition is underproduced by factors of 2-10 in the LM and HM models. Sulfur-bearing species are underproduced in the LM models, but perform reasonably well in the HM models, which can be directly understood as a consequence of the higher initial abundance of sulfur (Laas & Caselli 2019). For methanol, the differences between the two metallicities are minimal, which is due to the similar abundances of carbon and oxygen in the LM and HM models. However, the carbon-chain molecule c-C3H2 is better reproduced in models #6-#8 by the HM models. Nonetheless, c-C3H2 remains underproduced in all models. Models #1-#8 all agree poorly with the observed intensities. To assess whether this is linked to a too diffuse initial phase, we also ran an additional set of models with a denser initial phase and higher shielding (#9) to determine whether this would lead to better agreement. The results are shown in Fig. 11. Overall, a denser initial phase does lead to better agreement, but N2H+ and c-C3 H2 remain weaker than expected. HCO+ and the sulfur-bearing species also differ. The former is best reproduced with LM abundances, while SO and CS are significantly closer to the observed intensities with the HM conditions.

Figure 12 shows the comparison of the peak intensities in the fiducial static model for both LM and HM models. In this case, a better agreement with the observed ratios toward L1544 is observed for all species. In particular, the higher-metal model is within a factor of ~2 in intensity for all molecules.

Several of the targeted transitions are optically thick, with τ > 0.5. The comparable peak intensities therefore do not directly relate to the total column density, but instead indicate that the combination of the density, abundance, and temperature profile of the modeled core leads to excitation conditions that are qualitatively similar to L1544.

Thumbnail: Fig. 8 Refer to the following caption and surrounding text. Fig. 8

Line profiles as predicted by the fiducial model with HM abundances, 106 yr static phase, and Av,effambMathematical equation: $A_\mathrm{v,eff}^\mathrm{amb}$ = 2 mag. The spectra were extracted toward the continuum peak in the x-y plane. A fixed 16O/18O ratio of 560 was used for the C18O abundance profiles.
Thumbnail: Fig. 9 Refer to the following caption and surrounding text. Fig. 9

Comparison table for the peak intensity of the modeled emission lines toward L1544 (dust peak position) and toward the dust peak in the x-y plane. This figure shows the comparison for models with standard metallicities.
Thumbnail: Fig. 10 Refer to the following caption and surrounding text. Fig. 10

Similar to Fig. 9 but for higher metallicities.
Thumbnail: Fig. 11 Refer to the following caption and surrounding text. Fig. 11

Comparison table for the peak intensity of the modeled emission lines toward L1544 (dust peak position) and toward the dust peak in the x-y plane. This figure shows the comparison for the model #9 with either standard or higher metallicities.
Thumbnail: Fig. 12 Refer to the following caption and surrounding text. Fig. 12

Comparison table for the peak intensity of the modeled emission lines toward L1544 (dust peak position) and toward the dust peak in the x-y plane. This figure shows the comparison for the static model with either standard or higher metallicities.

4 Discussion

4.1 Effects of dynamic versus static modeling on intensity predictions

One of the principal motivations for us to include the dynamic tracer particle trajectories instead of relying on a static approach was to assess whether this time-dependent treatment leads to an improved agreement with observational constraints, particularly in terms of the spatial morphology of the emission regions and the relative intensities of key molecular species.

Our results do not lend themselves to a straightforward conclusion. With the updated chemical network, the predicted emission morphology of CH3OH and c-C3H2 is qualitatively similar for the dynamical model to the observed distribution toward L1544, whereas the predicted morphology for the static model is notably different. This is different for the predicted intensities. The dynamical models (#1-#8) predict lower intensities for all molecules when compared with the static model for similar initial elemental abundances. The lower abundances and intensities predicted by the dynamical model can be attributed to a central limitation in the current setup: the dynamical tracer particles follow trajectories over a total duration of 300 kyr, resulting in shorter residence times at high densities than in the static model, where chemical evolution is computed over 1 Myr across the entire parameter space. While a shorter residence time agrees better with our current understanding of a dynamic star formation process, it explains why the static model produces more N2H+, which requires substantial time to form at higher densities (Tafalla et al. 2021; Priestley et al. 2023). For dynamical models with a longer and/or a denser initial phase (#4 and #9), the N2H+ emission is brighter, although still not as bright as in the static model. An even denser initial phase is likely to bring the abundances and intensities of the dynamical model closer to that of the static model, but this comes with additional caveats.

In the current setup, this challenge cannot be easily circumvented. The core we studied here reaches central densities comparable to those of L1544 after 300 kyr of evolution and evolves to form a protostar shortly thereafter. Adopting a denser initial phase improves the agreement with observations, but at the cost of compromising the self-consistency of the model. The temperature profile of L1544 is well reproduced using an ambient extinction of Av,effamb≃2Mathematical equation: $A_\mathrm{v,eff}^\mathrm{amb} \simeq2$ mag, which is consistent with a lower background density of nH⪅103Mathematical equation: $n_\mathrm{H}\lessapprox$10$^{3}$ cm−3. Increasing Av,effambMathematical equation: $A_\mathrm{v,eff}^\mathrm{amb}$ to ≃ 4 mag, consistent with densities nH ≳ 104 cm−3, leads to temperatures that are too low at the edge of the core (see Fig. E.1), which contradicts the previous models of L1544 and the temperature maps derived from Herschel/SPIRE continuum maps (Keto & Caselli 2010; Bianchi et al. 2023). Given this inconsistency, we are cautious to suggest a denser ambient medium with higher shielding as a realistic option. An alternative approach might be to extend the duration tinit of the initial static phase in the dynamical models beyond 2 Myr. However, increasing tinit from 1 Myr to 2 Myr between models #3-#4 leads to a rather modest increase in the intensity of N2H+ by 25%. This suggests that increasing the duration alone does not lead to significantly higher abundances of N2H+ ; higher densities are necessary as well.

In contrast, the static model predicts stronger intensities for all molecules. The predicted intensities are notably closer to what is observed toward L1544, except for CH3 OH, which is generally reproduced well by both approaches. The HM static model predicts intensities that are all within a factor of ~2 from the observed intensities (Fig. 12). This level of agreement is quite satisfactory given that the modeled core is not a perfect replica of L1544. This suggests that, despite their simplified assumptions about physical evolution, static models can still serve as reliable tools for studying pre-stellar cores in three dimensions, especially when the primary goal is to reproduce observed emission intensities.

For c-C3H2, the intensity is underproduced by a factor of >10 in the dynamical models. For the static models, the agreement is within a factor of ~4, but the intensity and abundance are still underproduced. This behavior is a consequence of a rather low steady-state abundance of c-C3H2 in the chemical network. At higher densities and longer chemical timescales, carbon and hydrogen are locked in COMs and the abundance of c-C3H2 is reduced. The effect is less pronounced in the static model because the outer less dense regions of the core remain chemically young, meaning that a smaller fraction of the carbon is locked in CO. In the dynamical models, all trajectories experience periods of higher densities (nH ≳104 cm−3), which leads to a lower c-C3H2 abundance, and the effect is not reversed when the particles enter more diffuse regions later on.

For N2H+ and c-C3H2, the reaction rates used here from KIDA2014 were recently updated in KIDA2024 (Wakelam et al. 2024). The update contains revised rates that contribute to the formation and destruction of both these species and may consequently affect the results presented here. Updating the chemical network used here to the latest rates and benchmarking the updates is the subject of future work.

Future studies could focus on another of the 200 protostars that form in the simulation. However, when longer tracer particle tracks are used, another limitation might play a role. As the 4 pc3 RAMSES domain has periodic boundary conditions, the tracer particles will pass the domain multiple times when longer simulation timescales are considered. This may not pose a problem compared to very dense star-forming clusters, but it may not represent a realistic scenario if the simulation is compared to star-forming clouds with lower densities.

We caution that these results do not necessarily reflect that dynamical models in general are worse than static models, but rather that the current setup presented here lacks a long enough dynamical timescale to reproduce the physico-chemical history of the gas well enough to produce more realistic results than the static model. Sipilä et al. (2022) studied the effects of dynamical and static 1D models for L1544 and found generally similar results for the static and dynamical models, although the dynamic models did produce a slightly better agreement with observed line intensities.

We did not include a comparison of the line profiles between the model and L1544. The specific shape of the emission line profiles computed for the 3D model are highly dependent on the viewing angle, although the peak intensities remain comparable. A deeper analysis of the effect of the viewing angle on the modeled spectra is left for future work.

4.2 Interpretation of the morphologies

While the static model reproduces the intensities toward L1544 well, it fails to capture the observed characteristic structure for CH3OH and c-C3H2 toward L1544 and other young cores. In contrast, the dynamical models presented here reproduce the observed dichotomy in the molecular emission maps better. Because both approaches employ the same underlying physical structure and chemical network, the differences in morphology are attributed to the distinct evolutionary histories encoded in the tracer particle trajectories as opposed to the chemical evolution of the static model at fixed physical conditions. The insensitivity of the modeled morphology to the initial chemical conditions further supports the suggestion that the spatial distributions are primarily governed by the trajectories of individual tracers and not by the initial chemical state. Correlation analyses of column densities and integrated physical parameters reveal that CH3OH and c-C3H2 correlate with higher visual extinction and H2 column densities. This correlation is stronger for methanol in models with shorter initial chemical phases, suggesting that the effect of individual tracer trajectories diminishes when the initial chemistry is allowed to evolve over longer timescales. The correlation for methanol is stronger in models with lower ambient shielding, indicating that the observed offset in methanol emission (particularly its concentration in more shielded regions) is strongly affected by radiative processes. Higher ambient extinction leads to lower temperatures and diminished differences in photochemical processing along different particle paths. These conditions promote efficient CO ice formation and subsequent hydrogenation, enhancing methanol production, which contributes to the offset seen in the maps (Fig. 3). This finding agrees with previous studies by Spezzano et al. (2016, 2020), who concluded that methanol preferentially traces the most shielded parts of the core. It is worth noting that the correlation between the local physical parameters (N(H2) and N(Av,eff)) and the molecular emission does not indicate a reduction of the dependence on the initial conditions and tracer particle history. Instead, the observed correlation suggests that gas in regions with a higher density and visual extinction have, as part of their tracer particle histories, undergone a physical evolution promoting formation of, for example, CH3OH. Conversely, regions with lower density and visual extinction in the core studied here tend to have physical histories that are less suitable for CH3OH formation. This is different in the static model, where the abundance of CH3OH in any given point of the model only depends on the local physical conditions at that specific point. Compared with the intensity of c-C3H2 observed toward L1544, the static and dynamic models are too faint, which suggests that some aspects of both models fail to capture the formation of this molecule.

4.3 Comparison of the metallicities of sulfur-bearing species

Throughout this work, we have compared two sets of initial abundances. The standard metallicities presented here are commonly used in a broad range of astrochemical studies of prestellar cores (Aikawa et al. 2001; Vasyunin et al. 2017; Redaelli et al. 2021). They are depleted in sulfur compared to the solar values (Federman et al. 1993; Daflon et al. 2009). In our simulations, the modeled intensities of CS and SO are comparable (within factors of ≲2-3) to observed values for L1544 using the un-depleted (higher) metallicities, while the depleted abundances underproduce the intensities by a factor of ≳5. A similar result was also reported in Priestley et al. (2024) for a 3D smoothed-particle hydrodynamic simulation of a star-forming cloud and by Laas & Caselli (2019), who studied the depletion of sulfur from diffuse to dense clouds. Likewise, Hily-Blant et al. (2022) reported little evidence of sulfur depletion across a sample of young cores, with L1544 showing slightly more depletion than a number of younger cores. These results suggest that the high sulfur depletion (a factor of ~200) introduced in the standard metallicities may be a poor choice for modeling starless and pre-stellar cores. The initial motivation behind the depleted sulfur abundance was a series of studies that showed a clear deficit in sulfur-bearing molecules compared to cosmic abundances (Woods et al. 2015). This led to the suggestion that sulfur might be locked in a unidentified solid phase (e.g., Jiménez-Escobar et al. 2012; Martín-Doménech et al. 2016; Fuente et al. 2025). However, based on current astrochemical models, a high sulfur depletion appears unnecessary in L1544 and several other cores. This agrees with the recent study by Fuente et al. (2023), who compared observations and astrochemical models of CS across various star-forming regions. They found that the extent of sulfur depletion likely varies between star-forming clouds due to environmental effects. The regions studied in that work required a more modest depletion, with depletion factors of ~10, rather than ~100 as in the LM models presented here.

4.4 Caveats of the model

As discussed in earlier sections, the underproduction of the species in question within the dynamical model is likely due to the specific initial chemical conditions, such as the density, temperature, and the duration of the initial phase. Although the current model incorporates a more realistic representation of core formation by accounting for the physical evolution of the gas within the core, it remains sensitive to the initially selected conditions. This dependence on initial conditions is a common challenge in astrochemical modeling (e.g., Vidal et al. 2019; Aikawa et al. 2020). Overcoming this limitation would necessitate longer 3D cloud simulations, capable of fully capturing the transition from the diffuse interstellar medium through molecular cloud formation to the development of dense cores.

In the present model, we assumed a perfect coupling between gas and dust temperatures. This assumption is not valid at the lower densities present in the outer parts of the core and might specifically affect the conditions in the lower-density regions of the core, where carbon-chain molecules such as c-C3H2 are more abundant. However, this assumption applies to the static and dynamical models and therefore does not explain the poorer agreement for the modeled intensities in the dynamical model.

One possible problem of the tracer particle approach we used is that the tracer particles may not adequately sample all regions of parameter space. The current setup includes ~20000 tracer particles. These predominantly trace the denser regions of the gas, which can contribute to the underproduction of c-C3H2 in the model as c-C3H2 is destroyed in higher-density regions. It is computationally challenging to circumvent this behavior. The underlying RAMSES model already includes a high number of tracer particles (> 109), but not all particles will trace any given core in the molecular cloud simulation. To assess to which extent insufficient sampling may affect the results, we compared the densities traced by the tracer particles with the densities of the underlying grid. They agree well overall, but reveal slight differences at lower densities that might contribute to a poorer sampling of the chemistry in these regions of the dynamical models (see Figs. C.1 and C.2).

5 Summary

We presented a dynamical 3D physico-chemical model of a prestellar core resembling the pre-stellar core L1544. We compared the chemical morphology and modeled molecular intensities with L1544 and a static physico-chemical model of the same core. Our results are listed below:

  • Qualitatively, the dynamical model reproduced the chemical structure of carbon chains and COMs (traced by c-C3H2 and CH3OH, respectively) as observed in L1544. This was not the case for the static model. This effect was stronger when a weaker ambient extinction was applied around the core, indicating that it is directly linked to the structure of the core as traced by individual tracer particles;

  • The dynamical model was unable to reproduce the observed intensities of the modeled molecules simultaneously. Instead, the model had to be fine-tuned for individual molecules to produce a good agreement with observations of L1544. This was not the case for the static model, where the modeled species all agreed simultaneously within a factor of —2;

  • The weak molecular intensities predicted by the dynamical model are primarily due to the relatively short evolutionary timescale of the tracer particle trajectories. In particular, the limited duration prevents a sufficient formation of N2H+ (a key tracer of high-density regions), leading to underproduction compared to observations. One potential remedy is to pre-evolve the initial chemical composition of the tracer particles at higher densities, which helps to reduce the disagreement with the observed L1544 data and the static model. However, this approach is inconsistent with the inferred temperature profile of L1544, thereby compromising the physical self-consistency of the model;

  • The intensities derived from the dynamical model are highly sensitive to the conditions and duration of the initial static phase, reflecting a high degree of chemical inheritance. While this underscores the importance of the initial state in shaping the final chemical composition, it also highlights the need for careful and physically motivated choices when setting the initial conditions in astrochemical modelling;

  • Models with cosmic sulfur abundance agree well with observed intensities of SO and CS, while depleted abundances lead to underestimated intensities of the two molecules.

The static model we presented agrees well with the observed intensities toward L1544. However, the intensities predicted by the dynamical model using tracer particles are too low for most molecules, suggesting that some critical aspects are not adequately accounted for in the current framework. Future work should explore alternative approaches to improve the agreement for all molecules in the dynamical model. This might involve adjusting the initial phase or studying longer trajectories of the tracer particles to better capture the chemical evolution.

Acknowledgements

We wish to thank the anonymous referee for helpful comments which improved the manuscript. S.S.J. and S.S. wish to thank the Max Planck Society for the Max Planck Research Group funding. All other authors affiliated to the MPE wish to thank the Max Planck Society for financial support. The research leading to these results has received funding from the Independent Research Fund Denmark through grant No. DFF 10.46540/4283-00305B (TH). The Tycho supercomputer hosted at the SCIENCE HPC center at the University of Copenhagen, was used for carrying out the simulations, the analysis, and long-term storage of the results. This research has made use of spectroscopic and collisional data from the EMAA database (https://emaa.osug.fr and https://dx.doi.org/10.17178/ENAA). EMAA is supported by the Observatoire des Sciences de l’Univers de Grenoble (OSUG). This paper makes use of MATPLOTLIB (Hunter 2007) and SCIPY (Virtanen et al. 2020).

References

  1. Ahrens, V., Lewen, F., Takano, S., et al. 2002, Z. Naturforsch. A, 57, 669 [NASA ADS] [CrossRef] [Google Scholar]
  2. Aikawa, Y., Ohashi, N., Inutsuka, S.-i., Herbst, E., & Takakuwa, S. 2001, ApJ, 552, 639 [Google Scholar]
  3. Aikawa, Y., Furuya, K., Yamamoto, S., & Sakai, N. 2020, ApJ, 897, 110 [Google Scholar]
  4. Al-Belmpeisi, R., Tuhtan, V., Christensen, M. B., Kuruwita, R., & Haugbølle, T. 2024, MNRAS, 534, 3194 [Google Scholar]
  5. Bertin, M., Romanzin, C., Doronin, M., et al. 2016, ApJ, 817, L12 [NASA ADS] [CrossRef] [Google Scholar]
  6. Bianchi, E., Remijan, A., Codella, C., et al. 2023, ApJ, 944, 208 [NASA ADS] [CrossRef] [Google Scholar]
  7. Booth, A. S., Walsh, C., Terwisscha van Scheltinga, J., et al. 2021, Nat. Astron., 5, 684 [NASA ADS] [CrossRef] [Google Scholar]
  8. Borshcheva, K., Fedoseev, G., Punanova, A. F., et al. 2025, ApJ, 990, 163 [Google Scholar]
  9. Brinch, C., & Hogerheijde, M. R. 2010, A&A, 523, A25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  10. Caselli, P., & Ceccarelli, C. 2012, A&A Rev., 20, 56 [NASA ADS] [Google Scholar]
  11. Caselli, P., Hasegawa, T. I., & Herbst, E. 1998, ApJ, 495, 309 [Google Scholar]
  12. Caselli, P., Pineda, J. E., Zhao, B., et al. 2019, ApJ, 874, 89 [NASA ADS] [CrossRef] [Google Scholar]
  13. Cazzoli, G., Cludi, L., Buffa, G., & Puzzarini, C. 2012, ApJS, 203, 11 [NASA ADS] [CrossRef] [Google Scholar]
  14. Ceccarelli, C., Codella, C., Balucani, N., et al. 2023, in Astronomical Society of the Pacific Conference Series, 534, Protostars and Planets VII, eds. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, 379 [Google Scholar]
  15. Chandra, S., & Kegel, W. H. 2000, A&AS, 142, 113 [Google Scholar]
  16. Choudhury, S., Kim, J., Caselli, P., Lee, C. W., & Pineda, J. E. 2025, A&A, 696, A190 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  17. Clark, W. W., & De Lucia, F. C. 1976, J. Mol. Spectrosc., 60, 332 [CrossRef] [Google Scholar]
  18. Coutens, A., Commerçon, B., & Wakelam, V. 2020, A&A, 643, A108 [EDP Sciences] [Google Scholar]
  19. Cuppen, H. M., Walsh, C., Lamberts, T., et al. 2017, Space Sci. Rev., 212, 1 [Google Scholar]
  20. Daflon, S., Cunha, K., de la Reza, R., Holtzman, J., & Chiappini, C. 2009, AJ, 138, 1577 [Google Scholar]
  21. Denis-Alpizar, O., Stoecklin, T., Guilloteau, S., & Dutrey, A. 2018, MNRAS, 478, 1811 [Google Scholar]
  22. Denis-Alpizar, O., Stoecklin, T., Dutrey, A., & Guilloteau, S. 2020, MNRAS, 497, 4276 [Google Scholar]
  23. Diplas, A., & Savage, B. D. 1994, ApJ, 427, 274 [Google Scholar]
  24. Drozdovskaya, M. N., Schroeder I, I. R. H. G., Rubin, M., et al. 2021, MNRAS, 500, 4901 [Google Scholar]
  25. Dullemond, C. P., Juhasz, A., Pohl, A., et al. 2012, RADMC-3D: A multipurpose radiative transfer tool, Astrophysics Source Code Library [record ascl:1202.015] [Google Scholar]
  26. Endres, C. P., Schlemmer, S., Schilke, P., Stutzki, J., & Müller, H. S. P. 2016, J. Mol. Spectrosc., 327, 95 [NASA ADS] [CrossRef] [Google Scholar]
  27. Federman, S. R., Sheffer, Y., Lambert, D. L., & Gilliland, R. L. 1993, ApJ, 413, L51 [Google Scholar]
  28. Ferrada-Chamorro, S., Lupi, A., & Bovino, S. 2021, MNRAS, 505, 3442 [NASA ADS] [CrossRef] [Google Scholar]
  29. Ferrer Asensio, J., Jensen, S. S., Spezzano, S., et al. 2025, ApJ, 992, 124 [Google Scholar]
  30. Fuente, A., Rivière-Marichalar, P., Beitia-Antero, L., et al. 2023, A&A, 670, A114 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  31. Fuente, A., Esplugues, G., Rivière-Marichalar, P., et al. 2025, ApJ, 986, L17 [Google Scholar]
  32. Furuya, K., Aikawa, Y., Hincelin, U., et al. 2015, A&A, 584, A124 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  33. Galli, P. A. B., Loinard, L., Bouy, H., et al. 2019, A&A, 630, A137 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  34. Garrod, R. T. 2013, ApJ, 765, 60 [Google Scholar]
  35. Garrod, R. T., Wakelam, V., & Herbst, E. 2007, A&A, 467, 1103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  36. Giers, K., Spezzano, S., Caselli, P., et al. 2023, A&A, 676, A78 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  37. Giers, K., Spezzano, S., Lin, Y., et al. 2025, A&A, 699, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  38. Gieser, C., Pineda, J. E., Segura-Cox, D. M., et al. 2024, A&A, 692, A55 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  39. Goldsmith, P. F. 2001, ApJ, 557, 736 [Google Scholar]
  40. Harju, J., Caselli, P., Sipilä, O., et al. 2025, A&A, 700, A141 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  41. Haugbølle, T., Padoan, P., & Nordlund, Å. 2018, ApJ, 854, 35 [CrossRef] [Google Scholar]
  42. Heays, A. N., Visser, R., Gredel, R., et al. 2014, A&A, 562, A61 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  43. Heays, A. N., Bosman, A. D., & van Dishoeck, E. F. 2017, A&A, 602, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  44. Hernández Vera, M., Lique, F., Dumouchel, F., Hily-Blant, P., & Faure, A. 2017, MNRAS, 468, 1084 [Google Scholar]
  45. Hily-Blant, P., Pineau des Forêts, G., Faure, A., & Lique, F. 2022, A&A, 658, A168 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  46. Hocuk, S., Szúcs, L., Caselli, P., et al. 2017, A&A, 604, A58 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  47. Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [NASA ADS] [CrossRef] [Google Scholar]
  48. Jensen, S. S., Jørgensen, J. K., Furuya, K., Haugbølle, T., & Aikawa, Y. 2021a, A&A, 649, A66 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  49. Jensen, S. S., Jørgensen, J. K., Kristensen, L. E., et al. 2021b, A&A, 650, A172 [EDP Sciences] [Google Scholar]
  50. Jensen, S. S., Spezzano, S., Caselli, P., Grassi, T., & Haugbølle, T. 2023, A&A, 675, A34 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  51. Jiménez-Escobar, A., Muñoz Caro, G. M., Ciaravella, A., et al. 2012, ApJ, 751, L40 [Google Scholar]
  52. Jiménez-Serra, I., Vasyunin, A. I., Caselli, P., et al. 2016, ApJ, 830, L6 [Google Scholar]
  53. Jørgensen, J. K., Kuruwita, R. L., Harsono, D., et al. 2022, Nature, 606, 272 [CrossRef] [Google Scholar]
  54. Keto, E., & Caselli, P. 2010, MNRAS, 402, 1625 [Google Scholar]
  55. Keto, E., Caselli, P., & Rawlings, J. 2015, MNRAS, 446, 3731 [NASA ADS] [CrossRef] [Google Scholar]
  56. Kuffmeier, M., Frostholm Mogensen, T., Haugbølle, T., Bizzarro, M., & Nordlund, Å. 2016, ApJ, 826, 22 [NASA ADS] [CrossRef] [Google Scholar]
  57. Kuffmeier, M., Jensen, S. S., & Haugbølle, T. 2023, Eur. Phys. J. Plus, 138, 272 [NASA ADS] [CrossRef] [Google Scholar]
  58. Laas, J. C., & Caselli, P. 2019, A&A, 624, A108 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  59. Leemker, M., Tobin, J. J., Facchini, S., et al. 2025, Nat. Astron., 9, 1486 [Google Scholar]
  60. Leger, A., Jura, M., & Omont, A. 1985, A&A, 144, 147 [Google Scholar]
  61. Lique, F., Senent, M.-L., Spielfiedel, A., & Feautrier, N. 2007, J. Chem. Phys., 126, 164312 [NASA ADS] [CrossRef] [Google Scholar]
  62. Lique, F., Daniel, F., Pagani, L., & Feautrier, N. 2015, MNRAS, 446, 1245 [NASA ADS] [CrossRef] [Google Scholar]
  63. Lu, Y., Chang, Q., & Aikawa, Y. 2018, ApJ, 869, 165 [Google Scholar]
  64. Majumdar, L., Gratier, P., Ruaud, M., et al. 2017, MNRAS, 466, 4470 [Google Scholar]
  65. Maret, S., Bergin, E. A., & Tafalla, M. 2013, A&A, 559, A53 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  66. Martín-Doménech, R., Jiménez-Serra, I., Muñoz Caro, G. M., et al. 2016, A&A, 585, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  67. Müller, H. S. P., Schlöder, F., Stutzki, J., & Winnewisser, G. 2005, J. Mol. Struct., 742, 215 [Google Scholar]
  68. Navarro-Almaida, D., Lebreuilly, U., Hennebelle, P., et al. 2024, A&A, 685, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  69. Nomura, H., Furuya, K., Cordiner, M. A., et al. 2023, in Astronomical Society of the Pacific Conference Series, 534, Protostars and Planets VII, eds. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, 1075 [Google Scholar]
  70. Padoan, P., Pan, L., Pelkonen, V.-M., Haugbølle, T., & Nordlund, Å. 2025, Nat. Astron., 9, 862 [Google Scholar]
  71. Padovani, M., Ivlev, A. V., Galli, D., & Caselli, P. 2018, A&A, 614, A111 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  72. Panessa, M., Seifried, D., Walch, S., et al. 2023, MNRAS, 523, 6138 [NASA ADS] [CrossRef] [Google Scholar]
  73. Pelkonen, V.-M., Padoan, P., Haugbølle, T., & Nordlund, Å. 2021, MNRAS, 504, 1219 [Google Scholar]
  74. Pelkonen, V.-M., Padoan, P., Juvela, M., Haugbølle, T., & Nordlund, Å. 2025, A&A, 694, A327 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  75. Pickett, H. M., Poynter, R. L., Cohen, E. A., et al. 1998, J. Quant. Spec. Radiat. Transf., 60, 883 [Google Scholar]
  76. Pineda, J. E., Segura-Cox, D., Caselli, P., et al. 2020, Nat. Astron., 4, 1158 [NASA ADS] [CrossRef] [Google Scholar]
  77. Priestley, F. D., Clark, P. C., Glover, S. C. O., et al. 2023, MNRAS, 524, 5971 [CrossRef] [Google Scholar]
  78. Priestley, F. D., Clark, P. C., Glover, S. C. O., et al. 2024, MNRAS, 531, 4408 [NASA ADS] [CrossRef] [Google Scholar]
  79. Rabli, D., & Flower, D. R. 2010, MNRAS, 406, 95 [Google Scholar]
  80. Redaelli, E., Bizzocchi, L., Caselli, P., et al. 2019, A&A, 629, A15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  81. Redaelli, E., Sipilä, O., Padovani, M., et al. 2021, A&A, 656, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  82. Riedel, W., Sipilä, O., Redaelli, E., et al. 2023, A&A, 680, A87 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  83. Riedel, W., Sipilä, O., Redaelli, E., et al. 2025, A&A, 701, A291 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  84. Ruaud, M., Loison, J. C., Hickson, K. M., et al. 2015, MNRAS, 447, 4004 [Google Scholar]
  85. Schöier, F. L., van der Tak, F. F. S., van Dishoeck, E. F., & Black, J. H. 2005, A&A, 432, 369 [Google Scholar]
  86. Scibelli, S., Drozdovskaya, M. N., Caselli, P., et al. 2025, A&A, 702, A127 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  87. Sipilä, O., & Caselli, P. 2018, A&A, 615, A15 [Google Scholar]
  88. Sipilä, O., Harju, J., & Caselli, P. 2017, A&A, 607, A26 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  89. Sipilä, O., Caselli, P., Redaelli, E., & Spezzano, S. 2022, A&A, 668, A131 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  90. Spezzano, S., Bizzocchi, L., Caselli, P., Harju, J., & Brünken, S. 2016, A&A, 592, L11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  91. Spezzano, S., Caselli, P., Bizzocchi, L., Giuliano, B. M., & Lattanzi, V. 2017, A&A, 606, A82 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  92. Spezzano, S., Caselli, P., Pineda, J. E., et al. 2020, A&A, 643, A60 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  93. Tafalla, M., Usero, A., & Hacar, A. 2021, A&A, 646, A97 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  94. Teyssier, R. 2002, A&A, 385, 337 [CrossRef] [EDP Sciences] [Google Scholar]
  95. Tinti, F., Bizzocchi, L., Degli Esposti, C., & Dore, L. 2007, ApJ, 669, L113 [NASA ADS] [CrossRef] [Google Scholar]
  96. Tuhtan, V., Al-Belmpeisi, R., Christensen, M. B., Kuruwita, R., & Haugbølle, T. 2024, MNRAS, 534, 3176 [Google Scholar]
  97. Vasyunin, A. I., Caselli, P., Dulieu, F., & Jiménez-Serra, I. 2017, ApJ, 842, 33 [Google Scholar]
  98. Vidal, T. H. G., Gratier, P., Vaytet, N., Coutens, A., & Wakelam, V. 2019, MNRAS, 486, 5197 [NASA ADS] [Google Scholar]
  99. Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261 [Google Scholar]
  100. Visser, R., van Dishoeck, E. F., & Black, J. H. 2009, A&A, 503, 323 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  101. Wakelam, V., & Herbst, E. 2008, ApJ, 680, 371 [NASA ADS] [CrossRef] [Google Scholar]
  102. Wakelam, V., Loison, J. C., Herbst, E., et al. 2015, ApJS, 217, 20 [NASA ADS] [CrossRef] [Google Scholar]
  103. Wakelam, V., Gratier, P., Loison, J.-C., et al. 2024, A&A, 689, A63 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  104. Wolcott-Green, J., & Haiman, Z. 2011, MNRAS, 412, 2603 [Google Scholar]
  105. Wolcott-Green, J., Haiman, Z., & Bryan, G. L. 2011, MNRAS, 418, 838 [NASA ADS] [CrossRef] [Google Scholar]
  106. Woods, P. M., Occhiogrosso, A., Viti, S., et al. 2015, MNRAS, 450, 1256 [Google Scholar]
  107. Yang, B., Stancil, P. C., Balakrishnan, N., & Forrey, R. C. 2010, ApJ, 718, 1062 [Google Scholar]
4

The integrated columns for the visual extinction Av,eff and Tgas do not have a direct physical meaning here, but represent the structure of the core along the line of sight.

Appendix A The physical structure of the core

Figure A.1 shows radial profiles along 8 directions for the simulated core marked with different colors. The map in the lower right corner shows the continuum emission at 1.1 mm and the 8 directions in which the profiles were extracted. In this figure, an external extinction of (Av,amb) = 2 mag was adopted.

Thumbnail: Fig. A.1 Refer to the following caption and surrounding text. Fig. A.1

Physical structure of the core in the x -y plane. Left and top right : Radial profiles for the dust temperature, number density and visual extinction are shown along 8 different radial cuts on the x-y plane. Lower right : Continuum image of the core at 1.1 mm (in grayscale) and 30%, 45%, 60%, and 75% of the peak intensity (contours). Colored lines indicate the direction along which the radial profiles in the remaining panels were extracted.

Appendix B Example of tracer particle trajectories

In Fig. B.1 we show an example of five tracer particle trajectories. The corresponding gas-phase abundance of CH3OH and c-C3H2 are shown in the lower-right panel. The tracer particle tracjectories start at different position, with different physical conditions. All tracer particles except one (orange color) finish at higher densities than they started, with varying degree of freezeout of both gas-phase species.

Thumbnail: Fig. B.1 Refer to the following caption and surrounding text. Fig. B.1

Trajectories for five tracer particles. Abundances are shown for model #3 LM.

Appendix C Testing the sampling of the tracer particles

To visualize how well the tracer particles sample the actual core, we compared the structure of the 3D core when interpolated from tracer particle positions with the actual computational grid. Figure C.1 visualizes a kernel density estimation (KDE) of the density structure based on the actual model (top panel) and the KDE derived after interpolating the tracer particle densities onto the grid (bottom panel). The figures shows minor deviations along the edges of the KDE, while the overall density distributions are almost identical. Figure C.2 shows the deviation between the KDE and the underlying grid in percent.

Thumbnail: Fig. C.1 Refer to the following caption and surrounding text. Fig. C.1

Top: KDE for the density structure of the core as calculated using the true grid. Bottom: KDE for the density structure of the core when interpolated from the density of the tracer particles.
Thumbnail: Fig. C.2 Refer to the following caption and surrounding text. Fig. C.2

Deviation between the density structure of computational RAMSES grid and the density structure when interpolated from the ~20 000 tracer particles used in this study. The figure shows the deviation between the kernel density estimates for the density structure derived from either grid or tracer particles. Positive value indicate that the density derived from the tracer particles exceed the local grid value.

Appendix D References for spectroscopic data and collisional rates for radiative transfer calculations

An overview of the transitions studied using the nonlocal thermodynamic equilibrium line radiative transfer modeling is included in Table D.1. The table also list references for the collisional rates used for the modeling.

Table D.1

Transitions used for the line radiative transfer modeling.

Appendix E Temperature profiles as a function of ISRF extinction

Figure E.1 shows a comparison of the dust temperatures computed for different external shielding: Av,effambMathematical equation: $A_\mathrm{v,eff}^\mathrm{amb}$ = 1 mag, Av,effambMathematical equation: $A_\mathrm{v,eff}^\mathrm{amb}$ = 2 mag, and Av,effambMathematical equation: $A_\mathrm{v,eff}^\mathrm{amb}$ = 4 mag. Also included are the gas and dust temperature profiles from the 1D model of Keto et al. (2015).

Thumbnail: Fig. E.1 Refer to the following caption and surrounding text. Fig. E.1

Mean temperature profiles for different degrees of ISRF attenuation. Black lines show the gas and dust temperature profiles from the 1D L1544 presented in Keto et al. (2015).

All Tables

Table 1

Initial elemental abundances for the LM and HM models.

In the text
Table 2

Initial conditions for the chemical simulation before the calculation of the evolution during the tracer particle phase.

In the text
Table 3

Pearson correlation coefficients for model #3 LM in (x-y, x-z, z-y) planes.

In the text
Table 4

Pearson correlation coefficients for model #6 LM in (x-y, x-z, z-y) planes.

In the text
Table D.1

Transitions used for the line radiative transfer modeling.

In the text

All Figures

Thumbnail: Fig. 1 Refer to the following caption and surrounding text. Fig. 1

Projected H2 column densities in the molecular cloud simulation. The white crosses show the location of protostars formed in this snapshot from the simulation (t = 305 kyr). The red diamond marks the protostar that is studied here during the pre-stellar phase. The pre-stellar core studied in this work is centered in the figure for clarity because the box has periodic boundary conditions.
In the text
Thumbnail: Fig. 2 Refer to the following caption and surrounding text. Fig. 2

Four snapshots showing the accretion onto the pre-stellar core at different stages of the simulation. The red dots mark the position of the tracer particle accreted onto the pre-stellar core. Only 20% of the tracer particles are plotted for clarity. The color map shows the projected column density, similar to Fig. 1. Zero-age is defined as the onset of star formation in the simulation.
In the text
Thumbnail: Fig. 3 Refer to the following caption and surrounding text. Fig. 3

LM model with 106 yr static phase and Av,effambMathematical equation: $A_\mathrm{v,eff}^\mathrm{amb}$ = 2 mag (#3 LM). Top : integrated intensity levels for CH3 OH (white contours) and c-C3 H2 (red contours) along the three principal axes of the simulation. The color map shows the continuum emission at 1.1 mm. Middle : slice plot showing the local log Av,eff for each principal axes. The upper range of the color map is limited to 10 mag for clarity. Bottom : slice plot showing the local dust temperature for each principal axes. The red and black contours show the 10 K and 7 K limits, respectively.
In the text
Thumbnail: Fig. 4 Refer to the following caption and surrounding text. Fig. 4

Comparison of the fiducial static model and the fiducial dynamic model, both with the updated network from this work. Top : static model with LM abundances. The contours show the integrated-intensity levels for CH3OH (white) and c-C3H2 (red) along the three principal axes of the simulation. The color map shows the continuum emission at 1.1 mm. Bottom : similar to the top panel but for dynamical model #7 LM.
In the text
Thumbnail: Fig. 5 Refer to the following caption and surrounding text. Fig. 5

Comparison of the 18 dynamical models. The red and white contours indicate the integrated emission contours for c-C3H2 and CH3OH, respectively. The background shows the continuum map of the core at 1.1 mm. The c-C3H2 emission generally peaks at the dust continuum peak, while CH3OH is offset toward the more shielded part of the core.
In the text
Thumbnail: Fig. 6 Refer to the following caption and surrounding text. Fig. 6

Correlation between the column densities of c-C3 H2, CH3 OH, and Av,eff for models #1-#8 for standard and higher metallicities. With Av,effambMathematical equation: $A_\mathrm{v,eff}^\mathrm{amb}$=1 mag, CH3OH shows a stronger correlation, especially at lower tinit. For Av,effambMathematical equation: $A_\mathrm{v,eff}^\mathrm{amb}$=2 mag, the correlation is lower for CH3OH and stronger for c-C3H2
In the text
Thumbnail: Fig. 7 Refer to the following caption and surrounding text. Fig. 7

Line profiles as predicted by the fiducial model with LM abundances, 106 yr static phase, and Av,effambMathematical equation: $A_\mathrm{v,eff}^\mathrm{amb}$ = 2 mag. The spectra were extracted toward the continuum peak in the x-y plane. A fixed 16O/18O ratio of 560 was used for the C18O abundance profiles.
In the text
Thumbnail: Fig. 8 Refer to the following caption and surrounding text. Fig. 8

Line profiles as predicted by the fiducial model with HM abundances, 106 yr static phase, and Av,effambMathematical equation: $A_\mathrm{v,eff}^\mathrm{amb}$ = 2 mag. The spectra were extracted toward the continuum peak in the x-y plane. A fixed 16O/18O ratio of 560 was used for the C18O abundance profiles.
In the text
Thumbnail: Fig. 9 Refer to the following caption and surrounding text. Fig. 9

Comparison table for the peak intensity of the modeled emission lines toward L1544 (dust peak position) and toward the dust peak in the x-y plane. This figure shows the comparison for models with standard metallicities.
In the text
Thumbnail: Fig. 10 Refer to the following caption and surrounding text. Fig. 10

Similar to Fig. 9 but for higher metallicities.
In the text
Thumbnail: Fig. 11 Refer to the following caption and surrounding text. Fig. 11

Comparison table for the peak intensity of the modeled emission lines toward L1544 (dust peak position) and toward the dust peak in the x-y plane. This figure shows the comparison for the model #9 with either standard or higher metallicities.
In the text
Thumbnail: Fig. 12 Refer to the following caption and surrounding text. Fig. 12

Comparison table for the peak intensity of the modeled emission lines toward L1544 (dust peak position) and toward the dust peak in the x-y plane. This figure shows the comparison for the static model with either standard or higher metallicities.
In the text
Thumbnail: Fig. A.1 Refer to the following caption and surrounding text. Fig. A.1

Physical structure of the core in the x -y plane. Left and top right : Radial profiles for the dust temperature, number density and visual extinction are shown along 8 different radial cuts on the x-y plane. Lower right : Continuum image of the core at 1.1 mm (in grayscale) and 30%, 45%, 60%, and 75% of the peak intensity (contours). Colored lines indicate the direction along which the radial profiles in the remaining panels were extracted.
In the text
Thumbnail: Fig. B.1 Refer to the following caption and surrounding text. Fig. B.1

Trajectories for five tracer particles. Abundances are shown for model #3 LM.
In the text
Thumbnail: Fig. C.1 Refer to the following caption and surrounding text. Fig. C.1

Top: KDE for the density structure of the core as calculated using the true grid. Bottom: KDE for the density structure of the core when interpolated from the density of the tracer particles.
In the text
Thumbnail: Fig. C.2 Refer to the following caption and surrounding text. Fig. C.2

Deviation between the density structure of computational RAMSES grid and the density structure when interpolated from the ~20 000 tracer particles used in this study. The figure shows the deviation between the kernel density estimates for the density structure derived from either grid or tracer particles. Positive value indicate that the density derived from the tracer particles exceed the local grid value.
In the text
Thumbnail: Fig. E.1 Refer to the following caption and surrounding text. Fig. E.1

Mean temperature profiles for different degrees of ISRF attenuation. Black lines show the gas and dust temperature profiles from the 1D L1544 presented in Keto et al. (2015).
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.