© ESO 2021

1 Introduction

Understanding the formation and evolution of water during star and planet formation is essential to our understanding of the conditions for life in other planetary systems. Water is a prerequisite for life as we know it and, furthermore, also an important molecule for the planet formation process: it contributes significantly to the solid mass reservoirs outside the ice line and impacts the thermal evolution of the gas and the coagulation of dust particles (see, e.g., van Dishoeck et al. 2014).

It is well established that water predominantly forms on dust grain surfaces during the molecular cloud phase, where water constitutes the bulk of the ice (e.g., van Dishoeck et al. 2014). The evolution from the onset of star formation, through the protostellar collapse and protoplanetary disk phases, and finally onto planetary bodies is however still uncertain. Key questions include to what degree water is processed during star and planet formation, namely, whether planets accrete pristine water inherited from the molecular cloud, or if the water is a product of local processes within the envelope and disk. Crucially, it remains unclear to what extent variations in the local cloud environment impact the water chemistry of the final planetary system, and hence ultimately influence the conditions for biology in extrasolar systems.

Water deuterium fractionation serves as a powerful tracer of the chemical and physical evolution of water during the star and planet formation process, since the degree of deuterium fractionation depends sensitively on the formation environment(e.g., Caselli et al. 2012). Water formed under molecular cloud conditions, where the mean temperature is T_gas ~ 20 K and the visual extinction is low (~1–3 mag), shows a moderate degree of deuterium fractionation, with D/H ratios around 10⁻⁴ –10⁻³, compared tothe canonical value of 1.5 × 10⁻⁵ in the local interstellar medium (ISM) (Linsky 2003)¹. Conversely, water formed in prestellar cores, where temperatures are lower (~ 10 K) and the visual extinction higher (A_v > 5 mag), is highly enriched in deuterium, with D/H ratios around 10⁻² –10⁻¹. Hence, the D/H ratio of water can record information about the formation environment. Deuterium fractionation during star formation is primarily driven by the gas-phase exchange reaction H${}_3^{+}$ + HD ⇌ H₂ D⁺ + H₂ + ΔE, where the exact value of ΔE depends on the spin state of the involved reactants (Pagani et al. 1992; Hugo et al. 2009). This reaction is exothermic and the backward reaction is inhibited at low temperatures (T ≲ 50 K), leading to an enrichment of H₂ D⁺ relative to H${}_3^{+}$. H₂ D⁺ dissociatively recombines with free electrons to form atomic D, thus increasing the local atomic D/H ratio in the gas phase and ultimately on dust grain surfaces where water and other molecules are formed through hydrogenation. Detections of gas-phase H₂ D⁺ confirm this enrichment occurs in prestellar cores (Caselli et al. 2003; Vastel et al. 2004).

Observations toward young embedded protostars have revealed a complex evolution of water deuteration during the earliest stages of star formation. Single-dish observations of the cold extended envelope show HDO/H₂O ratios of the order of 10⁻² (Parise et al. 2005), while interferometric observations of the hot corino regions show lower ratios in the range 10⁻⁴ −10⁻³ (Coutens et al. 2013; Persson et al. 2014; Jensen et al. 2019). Meanwhile, the gaseous D₂ O/HDO ratio in hot horinos appear to be of the order of 10⁻² which is an order of magnitude higher than the HDO/H₂O ratio in the same region (Coutens et al. 2014). This is schematized in Fig. 1. A priori, a higher D₂ O/HDO ratio than HDO/H₂O is unexpected: if H₂O, HDO, and D₂O were formed at the same time through surface reactions the [D₂ O/HDO]/[HDO/H₂O] ratio would be the statistical value of 0.25 (Rodgers & Charnley 2002). The observed fractionation ratios have been explained by detailed chemical models invoking a layered ice structure with at least two notable components: a water-rich ice mantle with low D/H ratio and a surface ice component enriched in deuterium (e.g., Taquet et al. 2013; Furuya et al. 2016, and Fig. 1). In the cold outer envelope of protostars, the gas-phase HDO/H₂O ratio reflects the D/H ratio of the surface ice chemistry, namely, water formed during the prestellar core phase with efficient deuterium fractionation, and direct gas phase water formation in this region. Conversely, the HDO/H₂O and D₂ O/HDO ratios in hot corinos reflect the bulk ice reservoir, since the entire ice mantle is thermally desorbed in this region. The complete desorption of the ice layers lowers the HDO/H₂O ratio in the gas-phase, since the bulk water ice formed during the molecular cloud phase with a lower degree of deuterium fraction. Meanwhile, the D₂ O/HDO ratio remains roughly constant, as both D₂O and HDO primarily formed during similar physical conditions, in the prestellar core phase. This scenario can explain the observed D/H ratios toward young embedded protostars and indicates that the water chemistry observed in hot corinos predominantly originates from the molecular cloud, in other words, an inheritance scenario in which the water is largely inherited from the molecular cloud. Physicochemical models of protoplanetary disks generally support this notion: the D/H ratios measured in comets cannot alone be explained by local processing within the disk, and hence some degree of inheritance from the prestellar phase is required to reproduce the high D/H ratios of water observed in pristine Solar System material (see, e.g., Cleeves et al. 2014; Furuya et al. 2017; Altwegg et al. 2017).

As the water chemistry appears to be mainly inherited from the environment (i.e., the molecular cloud), the question becomes whether and how the large-scale environment can impact the water chemistry of the final planetary system. In a recent work, evidence for a correlation between the HDO/H₂O ratios and the local cloud environment is observed. Jensen et al. (2019) find a higher HDO/H₂O ratio towardisolated protostars than what has previously been detected toward sources in clustered environments. The authors proposed two explanations for such a correlation: (1) either temporal differences between isolated and clustered cores: a slower collapse of an isolated prestellar core could prolong the prestellar core phase leading to a higher D/H ratio in the water, or (2) higher temperatures or a stronger radiation field in clustered cores, where nearby (proto)stars and turbulent cloud dynamics may heat the gas compared to isolated counterparts. Chemical diversity induced by variations in the natal cloud environment have scarcely been studied theoretically, as the majority of physicochemical models of protostellar collapse have utilized one-dimensional models that fail to capture the complex nature of star formation indynamic molecular clouds. If local cloud variations drive a difference in the D/H ratio of water in the hot corino phase, such a difference may also impact the complex organic molecules (COMs) which are characteristic of hot corinos and hot cores. In a recent study, Aikawa et al. (2020) simulate the impact of various prestellar core conditions on the abundance of COMs and warm carbon chain chemistry (WCCC) in hot corinos. These authors findthat the WCCC is more pronounced in cores with lower initial temperature, lower extinction, or a longer prestellar core phase. No clear pattern emerged for the hot corino chemistry; some molecules, such as CH₃ OH, are less abundant when the temperature is higher, while other molecules show no impact of variation in prestellar core conditions.

In this work, we combine a detailed 3D MHD model of a star-forming region with radiative transfer and a three-phase chemical model to explore how the local environment may change the chemical evolution through protostellar collapse and determine which conditions can reproduce the observed variation in HDO/H₂O ratios toward low-mass protostars. Throughout the paper, we denote HDO/H₂O as f_D1 and D₂ O/HDO as f_D2. Finally, we define the ratio f_D2 /f_D1 as α.

Fig. 1 Schematics of the predicted ice structure in the prestellar core stage (right) and the observed water deuterium fractionation toward embedded low-mass protostars (left). The outer layer of the ice mantle, which has high water D/H ratio, could sublimate in the cold envelope, while the whole mantle is sublimated in the hot corino. The schematical structure is proposed by Furuya et al. (2016). Observational values are from Persson et al. (2014); Coutens et al. (2013, 2014); Jensen et al. (2019).

Fig. 2 Projected densities in the molecular cloud simulation. A total of 233 protostars, denoted with white crosses, form by the end of the simulation (t = 1.65 Myr). Most protostars form in clusters along filaments. The protostars studied in this work are denoted with red circles.

2 Methods

2.1 Physical collapse models

2.1.1 Molecular cloud models

A simulation of a 4 pc³ region of a molecular cloud is carried out with the adaptive mesh refinement (AMR) code RAMSES (Teyssier 2002). The entire box is shown in Fig. 2 toward the end of the simulation, where 233 protostars formed in the box. Initial conditions for the molecular cloud, including turbulence driving the star formation process, is identical to the models presented in Haugbølle et al. (2018). The RAMSES model presented in this work was modified to include sink particles and tracer particles. Sink particles represent protostars in the model and record the mass accretion onto the protostars. Tracer particles are passive particles injected into the simulation where they are advected along the gas flow, recording the dynamical evolution for chemical post-processing. A detailed description of the model, including these modifications, is available in Haugbølle et al. (2018). The molecular cloud simulation reproduces a number of significant empirical relations for star-forming clouds, such as the Larson scaling relations, the star formation rate (SFR), the core mass function (CMF) and the initial mass function (IMF). The model is therefore well suited to study star formation as a heterogeneous process, where the physical evolution on larger scales may impact the physical and chemical evolution on smaller scales. Compared to Haugbølle et al. (2018), the physical model used in this work has a factor of two higher resolution with a root grid of 512³. With six levels of refinement this results in a minimum cell-size of 25 au. This resolution is sufficient to study the dynamics from the larger molecular cloud scales down to the collapse of individual protostellar cores, but is not suitable to resolve any emerging accretion disk. The accretion rate onto sink particles within the simulation was recorded at a cadence of ~100 yr, while the physical structure of the simulation was recorded in ‘snapshots’ at a lower cadence, once every 5 kyr. We calculated the positions of the tracer particles between each snapshot from the position, velocity, and acceleration vectors at each snapshot. For each sink particle, we followed tracer particles for a duration of ~350 kyr, starting with a pre-collapse phase 100 kyr before the formation of the sink particle (i.e., the onset of protostellar collapse) and continuing through the Class 0 and Class I phases. Including a longer pre-collapse phase was not possible, as several of the studied sink particles form early in the simulation. An overview of the different phases of the simulation is shown in Fig. 3. During the pre-collapse phase, we averaged the density in 10 kyr windows to reduce computational time in this phase, where the temperature remains fixed.

As no disks form within the simulation, the evolution toward the end of the simulated time range is less accurate as the formation of an accretion disk changes the dynamics and radiation field in the inner part of the core. We note that higher-resolution simulations performed within the same framework, such as the models presented inKuffmeier et al. (2018), show that accretion disks form when the resolution is increased. We focus on the earlier stages of star formation in this work, namely the protostellar collapse down to hot corino scales (~ 100 au). A sample oflow-mass protostars, both clustered and isolated, were selected from the simulation. Throughout this work, we define isolated protostars as sink particles where no other protostars enter within 20 000 au during the simulation. For the clustered sources, we require at least two neighboring protostars within 10 000 au at some point in the simulated evolution. The characteristics of the selected sample is summarized in Table 1.

Fig. 3 Overview of the different phases of the model. In the initial prestellar (static) phase, the physical conditions are inferred previous simulations and observations (e.g., Keto & Caselli 2008). In the pre-collapse phase, which has a fixed duration of 100 kyr, the physical evolution is derived from RAMSES, but since the protostar is not yet formed, no radiative transfer is performed. During the collapse, the physical evolution is derived from RAMSES and post-processed using RADMC-3D. The duration of the collapse varies between individual tracer particles, with an upper duration of ~ 250 kyr.

2.1.2 Radiative transfer models

The dust temperatures and radiation field around each sink particle was post-processed using the Monte Carlo radiative transfer code RADMC-3D (Dullemond et al. 2012). For each sink particle in the RAMSES simulation, a spherical cutout of 30 000 au around the sink particle was converted to a RADMC-3D AMR grid. Creating smaller cutouts of the total grid is necessary toreduce the memory load, which can be substantial for large photon packages. The dust temperature is calculated via the MCTHERM routine; we injected 15 million photons in each radiative transfer model to sample the grid adequately. Models with a higher number of photons were tested and the results of the radiative transfer converge to the same result. We utilized the dust opacities from Semenov et al. (2003) in the radiative transfer modeling.

The protostars, represented by the sink particles, are included in the radiative transfer models as point sources. In addition to the central protostar in each cutout, we included neighboring protostars entering within 10 000 au of the protostar in focus. The protostellar luminosity is interpolated from the pre-main sequence (PMS) models from D’Antona & Mazzitelli (1997). We included a 100 kyr offset between the onset of core collapse and the PMS age, similar to Young & Evans (2005), corresponding to approximately one free-fall timescale. In the first 100 kyr of the collapse, the surface temperature of the protostar is fixed at 2000 K. The radius is set to be 2.5 R_⊙ for the calculation of the accretion luminosity. This implementation is similar to that presented in Frimann et al. (2016).

In addition to the protostellar luminosity derived from the PMS models, we included the accretion luminosity of the central protostar from Eq. (1), $$L_{\text{acc}}=G\frac{M\dot{M}}{2R}.$$(1)

The accretion rates onto the sink particles were recorded with a higher cadence than the RAMSES output. To account for this difference in sampling, we ran a total of five RADMC-3D models for each RAMSES output with varying accretion luminosities: we created a linear grid of four models, ranging from the lowest accretion rate to the highest accretion rate associated with any given RAMSES output. We calculated the fifth model using the mean accretion rate during the period between the snapshots. When estimating the dust temperature for a tracer particle at a given position and time, we calculated the total luminosity (L_acc + L_protostar) of the protostar at that time, and interpolated the temperature from the values provided by the five RADMC-3D models. We assume that the dust and gas is thermally coupled, such that the gas temperatures and dust temperatures are equal. This approach is reasonable at higher densities present during the protostellar collapse (n ≳ 10⁴, Goldsmith 2001), but may not be valid in the prestellar phase where densities are lower. The limitations of the model are discussed in Sect. 4.2. A lower limit on the temperature of 10 K is set throughout the chemical modeling. While lower temperatures can occur in the center of dense cores, this is a reasonably average temperature for dense cores (e.g., Pagani et al. 2007). Furthermore, many reactions in the chemical network are validated to a lower limit of 10 K.

Along with the thermal dust calculation, the mean UV radiation field is calculated in the RADMC-3D grid with the MCMONO routine. In addition to the protostellar sources, the models include an external radiation field in the form of the interstellar radiation field (ISRF) from Draine (1978), denoted as G₀, which irradiates the 30 000 au RADMC-3D domain around each sink particle. For the MCMONO calculations, we opted to run just one model using the mean accretion luminosity for each RAMSES snapshot. As the UV flux emerging from the central protostar is heavily attenuated, the calculated UV field is dominated by the external radiation field, except for the innermost 25–50 au around a protostar. Running five models with varying accretion luminosity would hence not have a notable impact on the surrounding envelope and the chemical evolution herein. Figure 4 shows an example of the thermal structure calculated around one of the sinks studied in this work.

We note that the earliest protostellar evolution may not be well explained by the PMS models, as highlighted by the luminosity problem (see, e.g., Jensen & Haugbølle 2018, and references therein). We do not consider this a major issue for this work for two reason. First, the accretion luminosity, not the protostellar luminosity, is the primary luminosity source at the very early stages modeled here. Second, the purpose of this work is to study how variations in environment influences the chemistry during star formation and does not depend on the exact choice of PMS model. Additionally, we note that the PMS phase occurs toward the end of the range of the time-domain studied in this work, namely, the ~ 250 kyr of the protostellar collapse after the initial collapse of the core.

Fig. 4 Structure around the binary sink M048-C with a final mass of M = 0.48 M⊙. The color map shows the density and reveals that the sink, along with several siblings, is born along a filament, extending from northwest to southeast in the imaged plane. Contours show the temperature, calculated via RADMC-3D, in levels of [25.0, 50.0, 100.0] K. Right panel: yellow 25 K contour is incomplete. Left panel: full structure of this region. The image shows a slice through the grid, not a projection.

2.2 Calculating extinction in the cores

For photodissociation and photodesorption, rate coefficients are given as a function of the visual extinction A_v. From the calculated mean field intensities, we follow the approach of Visser et al. (2011) and derive the optical depth τ_UV at UV wavelengths from: $${\tau }_{\text{UV}}=-\ln \frac{F_{\text{UV},\text{RADMC-3D}}}{F_{\text{UV},\text{0}}},$$(2)

where F_UV,RADMC-3D is the integrated UV flux from the radiativetransfer calculation, and F_UV,0 is either theintegrated UV flux from the ISRF or from the central protostar, whichever is higher in each grid point. Initially, the UV flux from the protostar is quite low owing to the low surface temperature of the protostar (2000 K) and during this period, the ISRF dominates the UV flux. Later on, the UV flux from the protostar may increase, however, the flux is not able to penetrate the envelope, because no outflow cavities are resolved in the simulation.

From the optical depth, the visual extinction is calculated as $A_{\text{v}}^{\text{RADMC-3D}}={\tau }_{\text{UV}}/3.02$. In addition to the extinction calculated with RADMC-3D, we assume an ambient extinction everywhere in the molecular cloud, (i.e., surrounding the 30 000 au radiative transfer domain). The total extinction during the collapse is then given as $A_{\text{v}}=A_{\text{v}}^{\text{cloud}}+A_{\text{v}}^{\text{RADMC-3D}}$, where we assume $A_{\text{v}}^{\text{cloud}}=5$ mag in our fiducial model. Models with a lower ambient extinction, $A_{\text{v}}^{\text{cloud}}=2$ mag, are presented in Appendix D. The column densities are estimated using the empirical relation N_H = A_V × 1.59 × 10²¹ cm⁻² (Bohlin et al. 1978; Diplas & Savage 1994). These parameters are necessary to derive the self-shielding functions and photochemical rates for the chemical models.

As mentioned in the previous section, each tracer particle trajectory includes 100 kyr of evolution prior to the formation of the sink particle. Since the protostar is not formed in the pre-collapse phase, no radiative transfer is calculated, and we use a fixed temperature of 10 K and a fixed visual extinction of A_v = 5 mag during this period.

2.3 Chemical model

The chemical evolution was simulated using the rate-equation approach (see, e.g., Cuppen et al. 2017). The model adopts the three-phase model introduced by Hasegawa & Herbst (1993), with a number of modifications detailed below. Three-phase models consider the ice mantle and ice surface as distinct phases, as opposed to two-phase models where the ice species all reside in the same (bulk) phase (Hasegawa et al. 1992). As the chemistry evolves, two-phase models may produce inaccurate results, since these models cannot distinguish between species present in the reactive surface layers and species buried deep in the ice mantle where the molecules are inert or less likely to react. Three-phase models overcome some of these issues by separating ice species in surface and mantle phases, however other caveats are introduced with these models, which we will discuss later. We do not consider bulk diffusion in our model; the effectiveness of bulk reactivity is still a matter of debate (see, e.g., Shingledecker et al. 2019).

Experiments have shown that surface chemistry likely occurs in the top two to four monolayers². Motivated by these results, we adopted a similar approach to that of Furuya et al. (2015), in which we extended the surface layer to thickness of four monolayers. Other three-phase models, such as the NAUTILUS, model have adopted a similar approach, albeit with two monolayers instead of four (Ruaud et al. 2016). Henceforth, we refer to the surface phase of the model as the surface layer, although the thickness of the phase is four monolayers.

2.3.1 Chemical network

For the gas-phase chemistry, we adopted the chemical network presented by Majumdar et al. (2017), which extends the KIDA³ network to include spin-state chemistry of H₂, H${}_3^{+}$ and other key molecules necessary to accurately model deuterium fractionation processes in the cold ISM. The basis of the spin-state chemistry was presented in Sipilä et al. (2015). The gas-phase network was initially benchmarked in Majumdar et al. (2017) for dark cloud conditions and found to agree with observations. In this work, the network was updated to include more recent photodissociation rates from the Leiden photodissociation database (Heays et al. 2017), a different treatment of grain recombination reaction, and other minor adjustments. We refer to the papers of Wakelam et al. (2015) and Majumdar et al. (2017) for a detailed description of the network and reactions. We excluded three-body reactions from these calculations, as they are inefficient under typical ISM conditions (Cuppen et al. 2017).

For the surface network, we used a network derived from Taquet et al. (2013), which was further extended to include reactions from Garrod (2013). The network has been extended to include single and doubly deuterated isotopologues for molecules with up to four carbon atoms. For a number of key species such as methanol, we included all deuterated isotopologues, such as CD₃OD.

Binding energies from Garrod & Herbst (2006) were used for the majority of the species, with updated estimates from Wakelam et al. (2017) used where applicable. For deuterated isotopologs, we used the same binding energy as the main species, except for atomic D, which is set 21 K higher than H atoms (Caselli et al. 2002). For the extensive network used in this work, a few species do not have any binding energy estimates. For these, we use the binding energy of H₂ O as a “default” value.

The full chemical network and the chemical code is publicly available for reference and provided for the community for future work⁴.

2.3.2 Adsorption and desorption

For the adsorption onto the grain surface and thermal desorption of surface species, we followed the method of Hasegawa et al. (1992). The adsorption rate coefficient is given by: $$k_{\text{ads}}(i)=S\langle v(i)\rangle {\sigma }_{\text{dust}}{n}_{\text{dust}}({\text{s}}^{-1}),$$(3)

where ${\sigma }_{\text{dust}}=\pi r_{\text{dust}}^2$ is the cross-section of dust grains, n_dust is the number density of dust grains, and $\langle v(i)\rangle =\sqrt{8k_{\text{b}}T/\pi m_i}$ is the mean thermal velocity of the species i. There parameter S is the sticking coefficient, which is estimated from the experimentally derived prescription presented by He et al. (2016).

The thermal desorption is calculated from the first-order expression as follows: $$k_{\text{des}}(i)={\nu }_i\exp \left(\frac{-E_{\text{bind}}}{k_{\text{b}}T}\right)({\text{s}}^{-1}),$$(4)

where E_bind is the binding energy of the species on the surface. ${\nu }_i=\sqrt{(2N_{\text{s}}{E}_{\text{bind}})/({\pi }^2{m}_i)}$ is the characteristic frequency for the species, calculated assuming a harmonic potential for the surface binding sites (Tielens & Allamandola 1987). The quantity N_s is the number of binding sitesper dust grain and m_i is the mass of the species.

Over time, incident cosmic-ray particles heat dust grains (e.g., Hasegawa & Herbst 1993). This leads to a cosmic-ray induced desorption rate, which we calculated following Leger et al. (1985): $$k_{\text{des},\text{CR-induced}}=f(70\text{\hspace{0.17em}K})k_{\text{des}}(70\text{\hspace{0.17em}K})({\text{s}}^{-1}),$$(5)

where the factor f determines the fraction of time the dust grains spend at 70 K. This value is determined through the balance between heating rate of dust grains by cosmic-ray particles and the cooling rate through radiation as follows: $f(70\text{\hspace{0.17em}K})=\frac{{\zeta }_{C}R}{1.3\times {10}^{-17}}\times 3.16\times {10}^{-19}$.

2.3.3 Surface reactions

We modeled surface chemistry reactions using the Langmuir-Hinshelwood mechanism, which assumes that both reactants adsorb onto the grain surface before any reaction occurs (Cuppen et al. 2017). After adsorption, the reactants diffuse thermally across the grain surface and upon encounter they may react, depending on the reaction probability. The rate coefficient k_ij for this mechanism is given as $$k_{\text{LH}}(i,j)={\kappa }_{ij}\left(k_{\text{hop}}(i)+k_{\text{hop}}(j)\right)\frac{1}{N_{\text{layer}}{N}_{\text{s}}{n}_{\text{dust}}}({\text{cm}}^3 {\text{s}}^{-1}),$$(6)

where N_layer = n_s∕(n_dustn_s) is the current number of surface layers in monolayers. In this equation, n_s denotes the sum of the number densities of all surface species, n_dust denotes the number density of dust grains, κ_ij is the reaction probability, and k_hop is the thermal hopping rate given by k_hop = ν_i exp(−E_diff∕k_bT). We assume E_diff = 0.5 × E_bind in this work. The reaction probability κ_ij is unity for barrierless reaction. For reactions with barriers, the barrier is overcome either through quantum tunneling or thermal hopping, whichever is the fastest. We included reaction-diffusion competition as introduced by Garrod & Pauly (2011). This mechanism considers that two reactants encountering one-another in a binding site may have multiple chances to react before one of them diffuses out of the site, hence increasing the reaction probability of reactions with barriers. The inclusion of reaction-diffusion competition is pronounced. For instance, the primary formation pathway of H₂ O changes from the barrierless reaction H + OH → H₂ O to the reaction H₂ + OH → H₂ O + H, which has a barrier of 2100 K, when the mechanism is included (Furuya et al. 2015). The reaction probability when assuming reaction-diffusion competition is calculated from: $${\kappa }_{ij}=\frac{{\kappa }_{ij}^0}{{\kappa }_{ij}^0+k_{\text{hop}}(i)+k_{\text{hop}}(j)},$$(7)

where ${\kappa }_{ij}^0$ is the probability to overcome the reaction barrier. The probability of crossing an activation barrier E_A through thermal hopping is given by ${\kappa }_{ij}^0=\exp (-E_{\text{A}}/k_{\text{b}}T)$, while the expression for quantum tunneling depends on the choice of potential barrier. For a rectangular barrier of width a, the expression is ${\kappa }_{ij}^0=\exp [-2(a/\hslash ){(2\mu E_{\text{A}})}^{1/2}]$. In this expression, μ denotes the reduced mass of the reactants. We adopted a rectangular barrier with a width of 1 Å, except for a number of reactions listed in Appendix A. For exothermic reactions, products may desorp from the grain surface owing to the excess energy released during the reaction. This process is known as chemical desorption and efficiency of this process is poorly constrained, but could play an important role for the interplay between gas-phase and grain-surface chemistry in star-forming regions. We included chemical desorption using the prescription of Garrod et al. (2007) for the probability of desorption for reactions with one product, with an efficiency parameter of a = 0.01 (see, e.g., Cuppen et al. 2017). Reactions with multiple products are not allowed to chemically desorp.

2.3.4 Photochemistry

The chemical model includes photodissociation and photodesorption of species in the surface layer. Data on photodissociation cross-sections for surface species are sparse and we adopted a similar approach as in Furuya et al. (2015) for the surface photochemistry. In this approach, gas-phase photodissociation rates are rescaled for the corresponding surface reaction. The scaling factor is determined from the ratio of the gas-phase to surface rates for water: $k_{ij}^{\text{surface}}=k_{ij}^{\text{gas}}\times k_{{\text{H}}_{\text{2}}\text{O}}^{\text{surface}}/k_{{\text{H}}_{\text{2}}\text{O}}^{\text{gas}}$. The rate coefficient for the photodissociation is given as: $$k_{\text{photodiss}}(i)={\theta }_i{P}_{\text{abs}}(i){\sigma }_{\text{dust}}{n}_{\text{dust}}{F}_{\text{UV}}(i)\times \min (N_{\text{layer}},4)({\text{s}}^{-1}),$$(8)

where θ_i is the surface coverage factor of the species, θ_i = n(i)∕n_s, F_UV is the UV flux in photons cm⁻² s⁻¹, and we allowed photodissociation in the top four monolayers, that is, the surface layer in the three-phase model. In our model, it is implicitly assumed that the photo-fragments in the ice mantle recombine immediately. While this may lead to an increase in the abundance of radicals and trigger the formation of COMs, the effect on the water chemistry studied in this work is expected to be limited as a result of the high extinction in these environments. The photodesorption rates are given by $$k_{\text{photodes}}(i)=Y_i{\theta }_i{\sigma }_{\text{dust}}{n}_{\text{dust}}{F}_{\text{UV}}(i)\times \min (N_{\text{layer}}/4,1)({\text{s}}^{-1}),$$(9)

where Y_i is the effective yield per photon for thick ice. For photodesorption, we likewise allowed reactions in up to four monolayers. Photodesorption yields for CO, CO₂, and N₂ are derived from experiments (Oberg et al. 2007; Öberg et al. 2009; Fayolle et al. 2011, 2013). The CO yield for thick ice is 2.7 × 10⁻³, while for CO₂ the total yield is 10⁻³, with an equal branching ratio between the outcomes CO^gas + O^gas and CO${}_2^{\text{gas}}$. The desorption yield for N₂ is set to 5 × 10⁻⁴. For H₂ O, HDO, and D₂O, we adopted the branching ratios for photodissociation and photodesorption presented by Arasa et al. (2015). The absorption probability for H₂ O per incident UV photon is set to 5 × 10⁻³. For the remaining species we used a generic yield of 10⁻³. While the network was developed to include complete photochemistry, we note that in the models presented in this work the photon flux is heavily attenuated for the majority of the time, and therefore the effect of photochemistry is limited. The model included self-shielding for H₂, HD, CO, and N₂. For CO and N₂, we interpolated the shielding factor from the tables of Visser et al. (2009) and Heays et al. (2014). For H₂ and HD we used the parametric approximations derived by Draine & Bertoldi (1996); Richings et al. (2014) and Wolcott-Green & Haiman (2011), respectively.

2.4 Benchmarking the network

The chemical model was benchmarked against the gas-grain model presented in Furuya et al. (2015) and Furuya et al. (2016), which successfully reproduced the observed HDO/H₂O and D₂ O/HDO ratios toward protostellar cores. We refer to this model as the reference model throughout this section. Initially, a 10 Myr steady-state simulation with n_H = 2 × 10⁴ cm⁻³, T = 10 K, and A_v = 5 mag was compared, shown in Fig. 5. The cosmic-ray ionization rate for H₂ is set to ${\xi }_{{\text{H}}_2}=1.3\times {10}^{-17} {\text{s}}^{-1}$. The abundances of major species, such as H₂O, CH₃OH, and CO₂ show good agreement, while the CO abundance differ by a factor of ~2. This is due to the inclusion of the hydrogen abstraction reaction: $$\text{HCO}+\text{H}\to \text{CO}+{\text{H}}_2.$$(10)

This reaction is not included in the reference model, but was recently demonstrated as a pathway for CO formation (Chuang et al. 2016). As in the reference model, H₂ O and HDO are primarily formed in the surface chemistry, while D₂O is dominantly formed in the gas phase. The principal formation pathway of H₂ O is the surface reaction H₂ + OH →H₂O + H, with a barrier of 2100 K. This pathway dominates when reaction-diffusion competition is included in the model. Without reaction-diffusion competition, the dominant pathway to H₂ O formation would be the surface reaction OH + H →H₂O, as in Furuya et al. (2015). The deuterated variant H₂ + OD →HDO + H can form HDO efficiently, but D₂O formation is inefficient owing to the lower tunneling probability for the heavier HD and D₂ molecules. Gas-phase formation of water proceeds through ion-molecule reactions, namely, D₂ O is formed in the reaction HD₂O⁺ + e⁻ $\underset{}{\overset{}{\to }}$ D₂ O + H (see, e.g., van Dishoeck et al. 2014). Overall, the reactions leading to the formation of water and deuterated isotopologs are similar to the reference model.

The chemical models were subsequently compared using the 1D radiative hydrodynamics (RHD) model utilized in Furuya et al. (2016). The physical model follows a tracer particle during the collapse phase of a 3.9 M_⊙ core that forms a solar mass star; the simulation is explained in detail in Masunaga & Inutsuka (2000). The initial conditions prior to the protostellar collapse are a 3 Myr static phase with n_H = 2 × 10⁴ cm⁻³, T = 10 K, and an visual extinction ofA_v = 5 mag. The model is started from the atomic abundances detailed in Table 2. We adopted the ISRF of Draine (1978), corresponding to F_UV ~ 2 × 10⁸ photons cm⁻² s⁻¹. The cosmic-ray induced UV flux is set to 3 × 10³ photons cm⁻² s⁻¹. We define these initial conditions as our fiducial model; we return to the impact of the conditions during the static phase in subsequent sections. Figure 6 shows a comparison between the results obtained at the end of the simulation for a particle trajectory that reaches the protostar (4.9 au) after ~ 350 kyr. The final gas-phase ratios are comparable, with f_D1 of 1.3 × 10⁻³ in the reference model, compared to 1.6 × 10⁻³ the model used here. For f_D2, the results are 10⁻² and 9.5 × 10⁻³ for the reference model and this model, respectively. The minor difference in ratios are attributed to differences in the gas-phase network, where there is less formation of HDO and D₂O compared with the reference model during the final stages of the collapse. The model presented here agrees with the chemical evolution of water proposed by Furuya et al. (2016), with lower D/H ratios in the mantle and higher D/H ratios in the surface ice (see Fig. 6).

Fig. 5 Steady-state surface chemistry at 10 K, nH = 2 × 104, and Av = 5 mag. All species are bulk (mantle + surface) ice, expect when specified otherwise. Left panel: model used here with the inclusion of CO formation through reaction (10), while the middle panels show the network with this reaction switched off. Right panel: reference model from Furuya et al. (2015, 2016).

3 Results

3.1 Comparing 1D and 3D models: D/H ratios of water

To establish how the chemical evolution of water differs between 1D and 3D models of the protostellar collapse, we compare the results of the 1D RHD model presented in Masunaga & Inutsuka (2000) with the 3D MHD model presented in this work. These models share the same static phase of 3 Myr as the fiducial model presented in the previous section. In Fig. 7, the evolution of gas-phase f_D1 and f_D2 during protostellar collapse are shown for a solar-mass protostar. The left panel shows the evolution for the 1D RHD model. The tracer particle reaches the hot corino (T ≳ 100 K) at τ~ 10², where τ = t_final − t. Once in thehot corino, the complete desorption of the water ice lowers the D/H ratios to f_D1~ 10⁻³ and f_D2~ 10⁻², respectively,as the surface and mantle ice D/H ratios are released and mixed with the gas-phase ratios. The right panel shows ten tracer particle trajectories from the RAMSES models accreted onto sink for the sink M101-C with a final mass of 1 M_⊙. In the figure, the shaded regions show the spread among the individual trajectories. In the 3D model, the gas-phase ratios show significant scatter during protostellar collapse and, generally, f_D1 and f_D2 are higher than in the 1D RHD model. The higher D/H ratios are due to higher densities, which leads to more efficient gas-phase water deuteration. Once the particles reach the hot corino (τ ~ 10³−10²), the complete desorption of the water-ice reservoir reduces the variation drastically and the final ratios are similar for the ten tracer particles. The final D/H ratio in the gas phase is primarily set by the D/H ratio in the ice, which is similar among trajectories. To illustrate this, Fig. 8 shows the change in the ice composition for the same ten particle trajectories during the pre-collapse phase (100 kyr prior to formation of the protostar) and through the collapse. The selected tracer particles all reach the hot corino at ~ 200−250 kyr after the onset of the collapse. Initially, at τ ≳ 2.5 × 10⁵ kyr, the D/H ratios in both surface and mantle ice show significant scatter. This scatter is caused by variation in the densities in the early evolution of the trajectories. Once the trajectories reach the protostellar core, however, the D/H ratios in the ice converge to similar values.

This change is driven by the convergence toward higher densities for the trajectories once they reach the dense core. This behavior is illustrated in Fig. 9, where the physical evolution of the ten particle trajectories is shown. The densities for the trajectories mainly range between 10⁵ –10⁶ cm⁻³ in the first 100 kyr, that is, around an order of magnitude higher than in the prestellar phase. At higher densities the thickness of the ice is increased. In the fiducial model this occurs at an epoch when water deuterium fractionation is highly efficient, which increases the D/H ratio in the water ice. Generally, the densities within the protostellar cores are similar and this behavior is seen for all cores in the simulation, independent of mass and environment (see Appendix C for the density profiles of the cores).

Toward the end of the collapse, the surface ice f_D1 varies up to a factor of ~4 between tracer trajectories, while f_D2 varies by less than a factor of 2. Meanwhile, in the mantle reservoirs, both f_D1 and f_D2 are almost constant between the tracer particle trajectories. We note that while the initial variations among trajectories are lost in Fig. 8, variations in the physical conditions in the static phase, still impact the final D/H ratio as shown in the subsequent sections.

Overall,since the bulk ice show little variation among tracer particle trajectories in the fiducial model, a tracer particle accreted during the main accretion phase (e.g., during the Class 0 stage) shows little variation in the final D/H ratio, if we assume the same initial conditions for all tracer particles.

Fig. 6 Comparison between the chemical models presented in Furuya et al. (2016) (top panels) and the model utilized here (bottom panels). The 1D RHD simulations of a protostellar collapse are from Masunaga & Inutsuka (2000). The solid lines indicate the ice mantle, the dashed lines present the ice surface, and the dotted lines show the gas-phase ratios in the three-phase models. The log-scale on the x-axis is reversed on the second row.

Fig. 7 Gas-phase evolution of fD1 and fD2 during protostellar collapse. Left panel: evolution of a particle accreted ~350 kyr after core collapse toward a 1 M⊙ protostar. Right panel: ten particle trajectories toward the 1 M⊙ protostars M101-C in the 3D simulations.

Fig. 8 Ice surface and mantle evolution for fD1 and fD2 during the protostellar collapse phase. The variation in densities for different trajectories in the pre-collapse phase induces large shifts in the ratios in the very beginning of the collapse. Once the tracer particles reach the protostellar core, the D/H ratios in the ice mantle converge to similar values, as the trajectories reach similar densities. The y-axes differ between the left and right panels.

Fig. 9 Evolution in density, Av, and temperature for the ten tracer particles shown in Fig. 8. The shaded region indicates the 100 kyr pre-collapse phase, before the protostar is formed. In the pre-collapse phase, the densities are averaged in 10 kyr windows, to reduce simulation time.

Fig. 10 HDO/H2O (fD1) and D2 O/HDO (fD2) ratios toward three protostars in the simulation, with similar final mass of ~ 0.5 M⊙. The tracer particles are binned according to the time at which they reach the hot corino, in bins with a width of 50 kyr. The time corresponds to the time after the onset of collapse, t0. The third protostar, sink M048-I, is classified as isolated since no protostar enters within 20 000 au during the simulation, while sink M048-C and sink M049-C are clustered. Each bar shows the median values of tracer particles accreted in the time interval, and the error shows the [15.9, 84.1] percentiles. The number of tracer particles within each bin is denoted in blue above the first row.

3.2 Comparing isolated and clusters protostars

In Figs. 10 and 11, we present comparisons between f_D1 for protostars with a similar final mass. The protostars labeled M048–I and M072–I are isolated, meaning that no neighboring protostars enter within 20 000 au of the protostar for the duration of the simulation. The remaining protostars are born in clusters; there are at least two neighboring protostars within 10 000 au of the protostars at some point during the simulation. Further characteristics of the protostars are presented in Table 1.

The top panels show f_D1 in the gas phase as the tracer particles reach the hot corino region, that is, the gas-phase HDO/H₂O ratio once the water ice has been sublimated off of the dust grains. We focus on this region to provide ratios directly comparable to the ratiosobserved toward hot corinos (Jensen et al. 2019). Tracer particles are binned according to the time at which they reach the hot corino. A minor temporal evolution is evident; the tracer particles that accrete earlier showing slightly lower D/H ratios than the tracer particles accreted toward the end of the simulation, but not to a degree where the effect is observationally resolvable. As mentioned in the previous section, the bulk of the ice reservoir remains similar for all tracer particles, limiting the degree to which the D/H ratios can vary. The bottom panels of the figures show the f_D2, with an identical binning on the x-axis. f_D2 shows less trend with time. This is attributed to the contemporary formation of the deuterated water isotopologs. As a longer accretion time leads to a higher overall D/H ratio, the ratio of deuterated to non-deuterated water is increased, while the ratio between the two deuterated isotopologs remain roughly constant.

No significant differences in the D/H ratios for the isolated and clustered protostars are seen in the simulations, in contrast with the observations presented in Jensen et al. (2019).

3.3 Varying the physical conditions in the prestellar phase

The results presented in Sect. 3.2 assume the same initial conditions at the onset of protostellar collapse. This assumption is unlikely to be fulfilled in dynamic star-forming regions, where temperature and radiation fields vary significantly depending on, for example, nearby massive stars and turbulence in the parental cloud. To assess the impact of the initial conditions, several combinations of initial conditions during the static phase were tested and are reported in Table 3.

3.3.1 Duration of prestellar phase

In Fig. 12, we show the f_D1 and f_D2 for static phases with durations of 1, 2, and 3 Myr. A shorter initial phase leads to a lower degree of deuterium fractionation. For a duration of 1 Myr, f_D1 ≈ 4 × 10⁻⁴ and f_D2 ≈ 2 × 10⁻³. Meanwhile, for a duration of 2 Myr, the D/H ratios are increased and the ratios are closer to the 3 Myr fiducial case, with f_D1 ≈ 1 × 10⁻³ and f_D2 ≈ 2 × 10⁻². The shorter static phase of 1 Myr introduces two additional effects: an increased scatter among different particle trajectories, as indicated bythe error bars in the figure, and a larger temporal variation among the tracer particles accreted early and late in the protostellar collapse. The latter effect is attributed to the shorter time to build the water ice, which has not reached a steady state at 1 Myr, making the chemical evolution more susceptible to variations in, for example, density, among the individual trajectories.

Fig. 11 Similar to Fig. 10, but for protostars with final masses of ~ 0.72 M⊙. The third protostar, sink M072-I, is classified as isolated since no protostar enters within 20 000 au during the simulation, while the remaining protostars are clustered.

3.3.2 Varying temperature

Since the efficiency of deuterium fractionation depends on the temperature, several temperatures were tested in the prestellar phase. In Fig. 13, we show the results of increasing the temperature during the static phase. The radiation field and cosmic-ray flux are unaltered with respect to the fiducial model, and the prestellar duration is kept at 3 Myr. Increasing the temperature in the initial phase has only limited impact on the water D/H ratio in the hot corino for a prestellar phase. While a higher temperature in the prestellar phase reduces the f_D1 and f_D2 in the ice prior to collapse, the 100 kyr pre-collapse phase where the temperature is fixed at 10 K followed by low temperatures during the early phases of the protostellar collapse are sufficient to increase the f_D1 and f_D2 to ≳ 10⁻³. The efficient deuterium fractionation on these short timescales is a consequence of the high densities in the cores (n_H ≳ 10⁵ cm⁻³). For models #13 and #14 listed in Table 3, minimum temperatures of 15 and 20 K, respectively, were enforced during the protostellar collapse. In these cases, the final D/H ratios in the hot corino remain lower, because the deuterium fractionation processes are less efficient throughout the entire chemical evolution. A stronger impact of higher temperatures is expected at temperatures >30 K, or at lower densities, where CO freeze-out is less efficient.

Fig. 12 fD1 and fD2 ratios toward the isolated protostar sink M048-I for varying prestellar phase durations. The tracer particles are binned according to the time at which they reach the hot corino, in bins with a width of 50 kyr. The time corresponds to the time after the onset of collapse, t0. Each bar shows the median values of tracer particles accreted in the time interval, and the error shows the [15.9, 84.1] percentiles.

3.3.3 Varying ISRF and cosmis-ray ionization rate

Increasing the cosmic-ray ionization rate leads to a higher f_D1, from ~ 10⁻³ at ξ_{H 2} = 10⁻¹⁷ s⁻¹, to f_D1 > 10⁻² at ξ_{H 2} = 10⁻¹⁵ s⁻¹. The impact on f_D2 is more complex. At ξ_H2 = 10⁻¹⁶ s⁻¹, f_D2 is slightly lower than the fiducial case, while at ξ_H2 = 10⁻¹⁵ s⁻¹, f_D2 is higher than the fiducial model early in the collapse and similar toward the end. Higher cosmic-ray ionization rates lead to a larger spread in the D₂ O/HDO ratios between particles trajectories, as seen in Fig. 14. The increased deuterium fractionation at stronger cosmic-ray fluxes is driven by the formation of ${\text{H}}_3^{+}$ and deuterated isotopologs, which are central to the fractionation processes in the cold ISM. The formation rates of ${\text{H}}_3^{+}$ and isotopologs are regulated by the cosmic-ray ionization of H₂ and consequently an increased cosmic-ray flux can increase rate of the fractionation process. Furthermore, higher ionization rates reduce the ortho-to-para conversion timescale for H₂, which further enhances deuteration. The different formation pathways for HDO and D₂ O (i.e., HDO formation through surface reactions and D₂O formation in the gas-phase) leads to the discrepancy between f_D1 and f_D2, as the efficiencies of the formation pathways differ.

Increasing the ISRF to 10 G₀ in the static phase leads to a slight increase of f_D2, while f_D1 remains roughly constant (see Fig. 15). At G = 100 G₀, f_D1 and f_D2 are increased compared with the fiducial model. These effects are due to different branching ratios for photodissociation of HDO, where the branch producing H + OD is favored over the D + OH branch. This leads to more efficient reformation of deuterated water, increasing deuterium fractionation through photodissociation in both gas-phase and grain-surface chemistry (see, e.g., Arasa et al. 2015, and references therein). We note that a visual extinction of A_v = 5 mag was kept during the testing of initial conditions, and the ISRF hence remained heavily attenuated. At a lower visual extinction, stronger radiation fields lower the D/H ratio, as CO freeze-out is reduced by photodesorption.

Fig. 13 fD1 and fD2 ratios toward the isolated protostar sink M048-I for varying temperatures during the prestellar phase. The tracer particles are binned according to the time at which they reach the hot corino, in bins with a width of 50 kyr. The time corresponds to the time after the onset of collapse, t0. Each bar shows the median values of tracer particles accreted in the time interval, and the error shows the [15.9, 84.1] percentiles.

Fig. 14 fD1 and fD2 ratios toward the isolated protostar sink M048-I for cosmic-ray ionization rates ξH 2. The tracer particles are binned according to the time at which they reach the hot corino, in bins with a width of 50 kyr. The time corresponds to the time after the onset of collapse, t0.

3.3.4 Reproducing observed HDO/H₂O and D₂O/HDO ratios in hotcorinos

Current estimates of f_D1 and f_D2 in hot corinos of low-mass protostars are limited to a few sources. Coutens et al. (2014) detect D₂ O toward IRAS 2A and derive α = f_D2∕f_D1 ≳ 7. Similar high values of α are also seen for other protostars (Jensen et al. 2021). We require α ≳ 5 and f_D1 in the range 10⁻⁴−10⁻³ to reproduce current observations of the water D/H ratios in hot corinos (Persson et al. 2014; Jensen et al. 2019). The fiducial model reproduced the observed f_D1 ~ 10⁻³ toward isolated protostars presented in Jensen et al. (2019), while the lower D/H ratio toward the clustered protostars require different initial conditions. With the model presented in this work, reproducing 10⁻⁴ < f_D1 < 10⁻³ requires a shorter prestellar core phase or higher temperatures during the collapse. Model #12, featuring a 1.5 Myr prestellar phase, yields f_D1~ 7 × 10⁻⁴ and α ≈13, within the observed range for clustered sources. Similarly, models #13 and #14, where minimum temperatures of 15 and 20 K are enforced during the entire collapse, yield f_D1 ~ 7 × 10⁻⁴ and f_D1 ~ 3 × 10⁻⁴, respectively, which is well within the observed ranges. These results are in line with the predictions from Jensen et al. (2019), who proposed that shorter timescale or higher temperatures could drive the observed differentiation between isolated and clustered protostars. Increasing the cosmic-ray flux changes α substantially, as shown in Table 3, and the observed D/H ratios toward low-mass Class 0 protostars appear to be incompatible with higher cosmic-ray fluxes (${\xi }_{{\text{H}}_2}={10}^{-15}$). The impact of a stronger ISRF on the D/H ratios of water is limited within the model, where the visual extinction is relatively high during the entire simulation (A_v ≳ 5).

Fig. 15 fD1 and fD2 ratios toward the isolated protostar sink M048-I for varying ISRF, 1 G0 corresponds to the ISRF of Draine (1978). Tracer particles are binned according to the time at which they reach the hot corino, in bins with a width of 50 kyr. Time corresponds to the time after the onset of collapse, t0. Each bar shows the median values of tracer particles accreted in the time interval, and the error shows the [15.9, 84.1] percentiles.

4 Discussion

4.1 Impact of the environment during protostellar collapse

The comparison of water D/H ratios in realistic 3D protostellar collapse models show no significant variation between clustered and isolated protostars. The mass of the protostar does not appear to influence the water chemistry notably either. The homogenity of the D/H ratios results from several factors. First, the physical structure of the low-mass protostellar cores presented in this work are similar; the density profiles resemble Bonnor-Ebert spheres at the onset of collapse⁵, and once a tracer particle reaches the core, the D/H ratios in the ice converge, independent of the surrounding environment. Second, the identical initial conditions in the model prior to the protostellar collapse limit the possible range of D/H in the hot corino. Our fiducial model includes a 3 Myr static core phase, during which the water D/H ratio is established in the ice. Throughout the protostellar collapse, variations in tracer trajectories (e.g., n(t) and T(t)) can induce variations in the gas-phase D/H ratios, but the bulk water reservoir shows limited variation between trajectories and once the thick ice mantle is evaporated, the gas-phase D/H ratio is largely similar to the ice D/H ratio which is set by the physical conditions during ice formation (i.e., the evolution in the molecular cloud and static core). While the observed differentiation in water deuteration toward isolated and clustered low-mass protostars is not explained by chemical processing during the protostellar collapse, the observed feature can be explained by variations in the physical evolution prior to protostellar collapse. A shorter static phase of ≲ 1.5 Myr results in f_D1 ≲ 7 × 10⁻⁴, which agrees with the results toward clustered sources (Persson et al. 2014). Furthermore, these conditions produce f_D2~ 10⁻², yielding α ≳5 as observed toward IRAS 2A (Coutens et al. 2014). Increasing the temperature to 20 K during the entire collapse leads to a lower f_D1~ 3 × 10⁻⁴, while α ~7 is reproduced for a prephase duration of >1.5 Myr. Overall, temperatures in the range T ≲ 20 K can reproduce the observed D/H ratios.

These models do not self-consistently predict the physical and chemical evolution during the molecular cloud and static core phase, however, both a shorter static phase and higher temperatures are feasible variations in dynamic molecular clouds. A shorter collapse timescale for clustered sources can be driven by higher densities or through external forces acceleration the collapse (e.g., Ward-Thompson et al. 2007; Enoch et al. 2008). Similarly, higher temperatures in clustered regions can occur through irradiation from nearby massive protostars or through shock heating in turbulent cloud environments (e.g., Krumholz 2014).

4.2 Caveats of the physical models

Variations in cloud conditions prior to the protostellar collapse are not included in the models presented, except for a 100 kyr period prior to the protostellar collapse where the temperature is fixed at 10 K. Variations in the initial conditions are necessary to reproduce the observed variation in D/H ratios among clustered and isolated protostars. To self-consistently model the conditions in the molecular cloud and static core phase, several enhancements to the model setup are required. Currently, radiative transfer is calculated locally around protostars once the protostellar collapse is initiated. The temperature in the gas prior to collapse is therefore unknown. Implementing global radiative transfer coupled with a thermochemical code, could provide better constraints on the temperature and radiation field during the molecular cloud and static phases (e.g., Glover et al. 2010; Grassi et al. 2014). To assess whether the fixed temperature of 10 K during the 100 kyr pre-collapse phase influenced the results, we reran the fiducial models with variable temperatures and extinctions derived from parametric relations during this phase. The extinction was derived from the total number density through the relation: $$A_{\text{v}}={\left(\frac{n_{\text{H}}}{{10}^3{\text{cm}}^{-3}}\right)}^{2/3} \text{mag}.$$(11)

The above relation is derived in Grassi et al. (2017) by fitting a power law to the density-extinction relation in the RHD models of Glover et al. (2010). A semi-analytical relation between extinction and temperature is derived and tested in Hocuk et al. (2017) as follows: $$T_{\text{dust}}=11+5.7\times \tan \text{h}[0.61-{\log }_{10}(A_{\text{v}})] \text{K}.$$(12)

Implementing these relations during the pre-collapse phase had minimal effect on the final D/H ratios and did not influence the conclusions presented in this work. This emphasizes that the D/H ratio is predominantly a product of the initial conditions, with limited variation among different tracer particle trajectories.

The RAMSES models are limited to by the grid resolution of 25 au. The resolution limit means that neither circumstellar disks nor outflowsare sufficiently resolved. Resolving the physical structure of both outflow cavities and the disk is necessary to accurately predict the radiation and temperature field in the later stages of the protostellar formation. We are therefore limited to studying the initial collapse phase, where most of the matter is accreted. Increasing the resolution allows the disk and outflow cavities to be resolved and the results presented in Kuffmeier et al. (2018) demonstrate that disks are formed early on around the Class 0/I stage.

The three-phase model presented in this work does not consider bulk diffusion, that is, diffusion between surface and mantle or reactions within the ice mantle. Including bulk diffusion may impact the chemical evolution of three-phase models, especially during a protostellar collapse (Garrod 2013). Without bulk diffusion, molecules can remain locked in the ice mantle at temperatures exceeding their sublimation temperature, since only the surface molecules are free to sublimate. For instance, CO remains in the ice mantle at temperatures above ~30 K, since only the surface CO can desorb to the gas phase. This can alter the evolution of COMs and the change in CO gas-phase abundance; however, this does not impact the gas-phase deuterium fractionation process notable (Furuya et al. 2015). Furuya et al. (2017) study the evolution of HDO/H₂O and D₂ O/HDO during protostellar collapse and find that the inclusion of bulk-diffusion does not impact the D/H ratio of water post-collapse. Owing to the high binding energy of water, bulk diffusion of water is likely inefficient and the effect is unlikely to alter the strong inheritance shown in the models.

4.3 Inheritance of water D/H ratio

Determining whether the D/H ratio is inherited or reset during star and planet formation is critical for our understanding of the evolution and delivery of water to young planets (e.g., van Dishoeck et al. 2014). Through simulations of water deuterium fractionation in a protoplanetary disk, Cleeves et al. (2014) showed that cometary D/H ratios are too high to originate from local processing within the disk. Hence, some degree of inheritance is necessary to reproduce the deuterium fractionation in the Solar System, including the Earth’s oceans. In a later work, Furuya et al. (2017) study the evolution of the D/H ratio of water, combining 1D protostellar collapse models with a 2D axisymmetric protoplanetary disk model. These authors find that the majority of the water in the protoplanetary disk remained unprocessed, again favoring inheritance for the majority of the water present in the disk. The results presented in this work favor chemical inheritance for the D/H ratio of water in the hot corino as opposed to a “reset” chemistry post-collapse. Variation in physical evolution during the protostellar collapse has little impact on the degree of water deuteration in the hot corino, which is dominated by the conditions in the static phase. However, the D/H ratio set in the static phase is not fixed and may evolve as the gas and dust enters the protostellar core, where densities increase. While this process adjusts the D/H ratio in the ice, it is not a source of differentiation in the final D/H ratios for water in the hot corinos since the structures of low-mass protostellar cores are largely similar between isolated and clustered protostars, owing to the physical processes that dominate the stability of cloud cores (e.g., Kuffmeier et al. 2017).

Furuya et al. (2017) suggest that α is a better tool for distinguishing between inheritance and resetting of the deuterium fractionation, since their simulations yield α > 1 under prestellar conditions, while resetting of the chemistry in the protoplanetary disk leads to α < 1. Our results confirm that α > 1 in the case of inherited D/H from the static phase when standard ISRF and cosmic-ray ionization rates are considered. This further motivates the study of D₂O/HDO ratios in hot corino to provide further constrains on the chemical evolution of water from cloud to core.

If the deuterium fractionation of water in the hot corino is inherited, the D/H ratio of water at later stages in the star and planet formation process could be linked to the initial conditions in the cloud environment prior to the onset of star formation. However, the nature of hot corinos, and the chemical link from hot corino to protoplanetary disk is not well established at this stage (e.g., Belloche et al. 2020; Jørgensen et al. 2020). In a recent work, Coutens et al. (2020) study the chemical evolution from the collapse of a prestellar core to the formation of a protoplanetary disk. As in this work, the physical evolution during the collapse is tracked using tracer particles. Their model couples a 3D non-ideal MHD RAMSES simulation of a collapsing core with the NAUTILUS chemical model and focuses on a broad range of molecules. These authors find a clear chemical link between the chemical composition in the prestellar core and in the protoplanetary disk. Generally, the principal elemental carriers remain the same throughout the collapse, with the exception of sulfur and phosphorus, where the main reservoirs evolve. Furthermore, the abundances of the most abundant molecules (e.g., H₂ O, CH₄, and CH₃OH) remain largely unaltered through the collapse, although the abundance of “large” COMs (e.g., CH₃CHO, CH₃OCH₃) increases during the warm-up phase. Overall, these conclusions agree with the results presented in this work for water, namely, that the initial conditions have a high impact on the final chemical composition; this supports a high degree of inheritance during the protostellar collapse.

In recent years, rotational signatures have been observed in a limited sample of young embedded Class 0 sources (e.g., Tobin et al. 2012; Murillo et al. 2013), including one source hosting a hot corino (Bjerkeli et al. 2019; Imai et al. 2019). This could suggest thatthe chemistry observed in hot corinos reflects the same chemical reservoir that constitutes the initial disk structure. Still, detections of the warm gas component in sources with clear disk structures remains elusive (e.g., Artur de la Villarmois et al. 2018). Furthermore, once the disk is formed, in situ chemical processes may locally alter the chemical composition throughout the disk, introducing chemical gradients in the radial and vertical directions (e.g., Henning & Semenov 2013). Recently, a correlation between COMs in the hot corino of IRAS16293–2422 and the comet 67P were established, supporting a potential link between the chemical composition of the hot corino and protoplanetary disk, where comets are formed (Drozdovskaya et al. 2019). This could imply limited chemical processing in the protosolar nebula at larger orbital distances where comets were formed. Furthermore, a recent study of COMs in the disk around V883 Ori found similarities with the chemistry in hot corinos and comets (Lee et al. 2019).

5 Summary

In this work we present physicochemical models of protostellar collapse to study how variations in protostellar cloud environments may impact the chemistry, specifically the D/H ratio of water. The main results are as follows:

• The HDO/H₂O and D₂O/HDO ratios show no clear correlation with the cloud environment. Physical variations during the protostellar collapse are not sufficient to drive a notable difference in the deuterium fractionation between different protostars as the D/H ratio is governed by the initial cloud conditions prior to protostellar collapse.

• To reproduce the observed differentiation in HDO/H₂O ratios between isolated and clustered protostar, variations in the duration or temperature of the prestellar phase are necessary. As such, the D/H ratio is determined by these conditions, favoring inheritance of the water D/H ratio observed in hot corinos. Static phases with durations of 1–3 Myr and temperatures 10–20 K can reproduce the observed HDO/H₂O and D₂O/HDO ratios.

• Increasing the cosmic-ray ionization rates to ξ_H2 ≈ 10⁻¹⁵ s⁻¹ leads to D/H ratios that are incompatible with the observed values in hot corinos toward low-mass Class 0 sources.

• Global chemical models of star-forming regions are needed to establish the origin of chemical differentiation in star-forming regions. Tracing the chemical evolution from molecular cloud to hot corino is a necessary step to self-consistently determine which processes drive chemical differentiation of water and more complex molecules.

This work illustrates how isotope fractionation can trace the chemical evolution through the star and planet formation process, by combining detailed models with the latest observational advances. Furthermore, present results highlight the need for physicochemical models spanning the entire star formation process, to accurately model the origin and evolution of chemical complexity during the formation of stars and planets.

Acknowledgements

The authors are grateful to Tommaso Grassi, who provided the underlying framework for the chemical model. The authors would like to thank the referee for constructive comments which helped improve the manuscript. The group of J.K.J. acknowledges support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 646908) through ERC Consolidator Grant “S4F”. TH was supported by the Independent Research Fund Denmark through grant No. DFF 8021–00350B. We acknowledge PRACE for awarding us access to Curie at GENCI@CEA, France, and to Marconi at CINECA, Italy used to carry out the simulation. Resources at the University of Copenhagen HPC centre, funded in part by the Carlsberg, Novo, and Villum foundations, were used for the data analysis. This work made use of MATPLOTLIB (Hunter 2007) and NUMPY (Harris et al. 2020). Y.A. acknowledges the support by Grant-in-Aid for Scientific Research (S) 18H05222, and Grant-in-Aid for Transformative Research Areas (A) 20H05844 and 20H05847.

Appendix A Surface reactions with custom barriers

Appendix B Comparison of D/H ratios toward 0.22 M_⊙ protostars

Fig. B.1 HDO/H2O (fD1) and D2 O/HDO (fD2) ratios toward three protostars in the simulation, with similar final mass of ~ 0.22 M⊙. Tracer particles are binned according to the time at which they reach the hot corino, in bins with a width of 50 kyr. Time corresponds to the time after the onset of collapse, t0. Each bar shows the median values of tracer particles accreted in the time interval, and the error shows the [15.9, 84.1] percentiles. The number of tracer particles within each bin is denoted in blue above the first row.

Appendix C Density structure of the cores at onset of collapse

Figure C.1 shows the radial density profile for the 9 protostars studied in this work at the onset ofcollapse. The three isolated sink particles show a smooth Bonner-Ebert-like profile with ρ ∝ r⁻². The clustered sources also resemble a Bonner-Ebert profile, however several sinks show indications of additional over-densities within 10 000 au, e.g., sink M048-C and sink M101-C. Overall, the profiles of the isolated cores also exhibit slightly less scatter around the medians.

Fig. C.1 Radial density profile for 40 000 tracer particle for the sinks studied in this work. The red line shows the median in each radial bin, while the errorbars show the 25 and 75 percentiles for each bin. Blue circles indicates individual tracer particles.

Appendix D Loweringthe ambient extinction A${}_{\text{v}}^0$

In the fiducial model, the ambient cloud is assumed to shield the cores with an extinction of 5 mag. To determine the impact of a lower ambient extinction, models with an ambient extinction of A${}_{\text{v}}^0=2$ mag are presented here for 6 of the 9 protostars. Physical and chemical parameters are similar to the fiducial model, i.e., T = 10 K, n_tot= 2 × 10⁴ cm⁻³, and a duration of 3 Myr. In this case, the D/H ratios of water in the hot corino are slightly higher, and a larger spread in D/H ratios is observed. Figures D.1 and D.2 show the comparison between isolated and clustered protostars with similar final masses. No apparent correlation with the environment in present.

Fig. D.1 Similar to Fig. 10, with an ambient cloud extinction of 2 mag. HDO/H2O (fD1) and D2 O/HDO (fD2) ratios toward three protostars in the simulation, with similar final mass of ~ 0.5 M⊙. Tracer particles are binned according to the time at which they reach the hot corino, in bins with a width of 50 kyr. Time corresponds to the time after the onset of collapse, t0. The third protostar, sink M048-I, is classified as isolated since no protostar enter within 20 000 au during the simulation, while sink M048-C and sink M049-C are clustered. Each bar shows the median values of tracer particles accreted in the time interval, and the error shows the [15.9, 84.1] percentiles. The number of tracer particles within each bin is denoted in blue above the first row.

Fig. D.2 Similar to Fig. 11, with an ambient cloud extinction of 2 mag. HDO/H2O (fD1) and D2 O/HDO (fD2) ratios toward three protostars in the simulation, with similar final mass of ~ 0.7 M⊙. Tracer particles are binned according to the time at which they reach the hot corino, in bins with a width of 50 kyr. Time corresponds to the time after the onset of collapse, t0. Each bar shows the median values of tracer particles accreted in the time interval, and the error shows the [15.9, 84.1] percentiles. The number of tracer particles within each bin is denoted in blue above the first row.

References

All Tables

All Figures

Fig. 1 Schematics of the predicted ice structure in the prestellar core stage (right) and the observed water deuterium fractionation toward embedded low-mass protostars (left). The outer layer of the ice mantle, which has high water D/H ratio, could sublimate in the cold envelope, while the whole mantle is sublimated in the hot corino. The schematical structure is proposed by Furuya et al. (2016). Observational values are from Persson et al. (2014); Coutens et al. (2013, 2014); Jensen et al. (2019).

In the text

	Fig. 2 Projected densities in the molecular cloud simulation. A total of 233 protostars, denoted with white crosses, form by the end of the simulation (t = 1.65 Myr). Most protostars form in clusters along filaments. The protostars studied in this work are denoted with red circles.
In the text

Fig. 3 Overview of the different phases of the model. In the initial prestellar (static) phase, the physical conditions are inferred previous simulations and observations (e.g., Keto & Caselli 2008). In the pre-collapse phase, which has a fixed duration of 100 kyr, the physical evolution is derived from RAMSES, but since the protostar is not yet formed, no radiative transfer is performed. During the collapse, the physical evolution is derived from RAMSES and post-processed using RADMC-3D. The duration of the collapse varies between individual tracer particles, with an upper duration of ~ 250 kyr.

In the text

Fig. 4 Structure around the binary sink M048-C with a final mass of M = 0.48 M⊙. The color map shows the density and reveals that the sink, along with several siblings, is born along a filament, extending from northwest to southeast in the imaged plane. Contours show the temperature, calculated via RADMC-3D, in levels of [25.0, 50.0, 100.0] K. Right panel: yellow 25 K contour is incomplete. Left panel: full structure of this region. The image shows a slice through the grid, not a projection.

In the text

	Fig. 5 Steady-state surface chemistry at 10 K, nH = 2 × 104, and Av = 5 mag. All species are bulk (mantle + surface) ice, expect when specified otherwise. Left panel: model used here with the inclusion of CO formation through reaction (10), while the middle panels show the network with this reaction switched off. Right panel: reference model from Furuya et al. (2015, 2016).
In the text

Fig. 6 Comparison between the chemical models presented in Furuya et al. (2016) (top panels) and the model utilized here (bottom panels). The 1D RHD simulations of a protostellar collapse are from Masunaga & Inutsuka (2000). The solid lines indicate the ice mantle, the dashed lines present the ice surface, and the dotted lines show the gas-phase ratios in the three-phase models. The log-scale on the x-axis is reversed on the second row.

In the text

	Fig. 7 Gas-phase evolution of fD1 and fD2 during protostellar collapse. Left panel: evolution of a particle accreted ~350 kyr after core collapse toward a 1 M⊙ protostar. Right panel: ten particle trajectories toward the 1 M⊙ protostars M101-C in the 3D simulations.
In the text

Fig. 8 Ice surface and mantle evolution for fD1 and fD2 during the protostellar collapse phase. The variation in densities for different trajectories in the pre-collapse phase induces large shifts in the ratios in the very beginning of the collapse. Once the tracer particles reach the protostellar core, the D/H ratios in the ice mantle converge to similar values, as the trajectories reach similar densities. The y-axes differ between the left and right panels.

In the text

	Fig. 9 Evolution in density, Av, and temperature for the ten tracer particles shown in Fig. 8. The shaded region indicates the 100 kyr pre-collapse phase, before the protostar is formed. In the pre-collapse phase, the densities are averaged in 10 kyr windows, to reduce simulation time.
In the text

Fig. 10 HDO/H2O (fD1) and D2 O/HDO (fD2) ratios toward three protostars in the simulation, with similar final mass of ~ 0.5 M⊙. The tracer particles are binned according to the time at which they reach the hot corino, in bins with a width of 50 kyr. The time corresponds to the time after the onset of collapse, t0. The third protostar, sink M048-I, is classified as isolated since no protostar enters within 20 000 au during the simulation, while sink M048-C and sink M049-C are clustered. Each bar shows the median values of tracer particles accreted in the time interval, and the error shows the [15.9, 84.1] percentiles. The number of tracer particles within each bin is denoted in blue above the first row.

In the text

	Fig. 11 Similar to Fig. 10, but for protostars with final masses of ~ 0.72 M⊙. The third protostar, sink M072-I, is classified as isolated since no protostar enters within 20 000 au during the simulation, while the remaining protostars are clustered.
In the text

Fig. 12 fD1 and fD2 ratios toward the isolated protostar sink M048-I for varying prestellar phase durations. The tracer particles are binned according to the time at which they reach the hot corino, in bins with a width of 50 kyr. The time corresponds to the time after the onset of collapse, t0. Each bar shows the median values of tracer particles accreted in the time interval, and the error shows the [15.9, 84.1] percentiles.

In the text

Fig. 13 fD1 and fD2 ratios toward the isolated protostar sink M048-I for varying temperatures during the prestellar phase. The tracer particles are binned according to the time at which they reach the hot corino, in bins with a width of 50 kyr. The time corresponds to the time after the onset of collapse, t0. Each bar shows the median values of tracer particles accreted in the time interval, and the error shows the [15.9, 84.1] percentiles.

In the text

	Fig. 14 fD1 and fD2 ratios toward the isolated protostar sink M048-I for cosmic-ray ionization rates ξH 2. The tracer particles are binned according to the time at which they reach the hot corino, in bins with a width of 50 kyr. The time corresponds to the time after the onset of collapse, t0.
In the text

Fig. 15 fD1 and fD2 ratios toward the isolated protostar sink M048-I for varying ISRF, 1 G0 corresponds to the ISRF of Draine (1978). Tracer particles are binned according to the time at which they reach the hot corino, in bins with a width of 50 kyr. Time corresponds to the time after the onset of collapse, t0. Each bar shows the median values of tracer particles accreted in the time interval, and the error shows the [15.9, 84.1] percentiles.

In the text

Fig. B.1 HDO/H2O (fD1) and D2 O/HDO (fD2) ratios toward three protostars in the simulation, with similar final mass of ~ 0.22 M⊙. Tracer particles are binned according to the time at which they reach the hot corino, in bins with a width of 50 kyr. Time corresponds to the time after the onset of collapse, t0. Each bar shows the median values of tracer particles accreted in the time interval, and the error shows the [15.9, 84.1] percentiles. The number of tracer particles within each bin is denoted in blue above the first row.

In the text

	Fig. C.1 Radial density profile for 40 000 tracer particle for the sinks studied in this work. The red line shows the median in each radial bin, while the errorbars show the 25 and 75 percentiles for each bin. Blue circles indicates individual tracer particles.
In the text

Fig. D.1 Similar to Fig. 10, with an ambient cloud extinction of 2 mag. HDO/H2O (fD1) and D2 O/HDO (fD2) ratios toward three protostars in the simulation, with similar final mass of ~ 0.5 M⊙. Tracer particles are binned according to the time at which they reach the hot corino, in bins with a width of 50 kyr. Time corresponds to the time after the onset of collapse, t0. The third protostar, sink M048-I, is classified as isolated since no protostar enter within 20 000 au during the simulation, while sink M048-C and sink M049-C are clustered. Each bar shows the median values of tracer particles accreted in the time interval, and the error shows the [15.9, 84.1] percentiles. The number of tracer particles within each bin is denoted in blue above the first row.

In the text

Fig. D.2 Similar to Fig. 11, with an ambient cloud extinction of 2 mag. HDO/H2O (fD1) and D2 O/HDO (fD2) ratios toward three protostars in the simulation, with similar final mass of ~ 0.7 M⊙. Tracer particles are binned according to the time at which they reach the hot corino, in bins with a width of 50 kyr. Time corresponds to the time after the onset of collapse, t0. Each bar shows the median values of tracer particles accreted in the time interval, and the error shows the [15.9, 84.1] percentiles. The number of tracer particles within each bin is denoted in blue above the first row.

In the text