© ESO, 2014

1. Introduction

In cold clouds (~10 K), water is predominantly found on dust grains in the form of water ice. Its main formation pathway is in-situ: atomic oxygen is accreted from the gas phase onto the grain surface and is successively hydrogenated to form water ice (e.g. Tielens & Hagen 1982; Hiraoka et al. 1998; Miyauchi et al. 2008; Ioppolo et al. 2008, 2010; Mokrane et al. 2009; Oba et al. 2009; Cuppen et al. 2010; Dulieu et al. 2010). Water vapour can form in the gas phase through ion-molecule chemistry at low temperatures and through neutral-neutral reactions at high temperatures (for recent review, see van Dishoeck et al. 2013). The gas and ice phases are linked through freeze-out from the gas phase onto dust grains, and through thermal and non-thermal desorption of ice back into the gas phase. Because water ice formation is efficient and starts already in molecular clouds prior to collapse, most models of pre- and protostellar evolution turn the bulk of the available oxygen into water ice in the cold parts of the cores (e.g., Aikawa et al. 2008; Hollenbach et al. 2009; Caselli et al. 2012).

Water ice can be observed through infrared absorption of vibrational bands superposed on the continuum of an embedded young stellar object or a background star. Surveys of large samples of low- and high-mass protostars as well as background sources reveal typical water ice column density ratios N_{s - H2O}/N_H ~ 5 × 10^-5 (e.g., Smith et al. 1993; Gibb et al. 2004; Pontoppidan et al. 2004; Boogert et al. 2004, 2008, 2011; Whittet et al. 2001, 2007, 2013; Öberg et al. 2011). Here N_H is the column density of hydrogen nuclei, inferred either from the silicate optical depth or the colour excess toward the star. Although these values can vary by up to a factor of 2, the implication is that water ice contains only a small fraction, <20%, of the overall interstellar oxygen abundance with respect to hydrogen atoms of 5.75 × 10^-4 (Przybilla et al. 2008). Water gas can lock up a large fraction of oxygen, but only in hot gas where high-temperature neutral-neutral reactions drive most of the oxygen not contained in refractory grains into water (e.g., van Dishoeck et al. 2011, and references therein). Indeed, in cold regions the water gas abundance has been found to be very low, 10^-10−10^-8, as inferred from observations with the Submillimeter Wave Astronomy Satellite (SWAS) and subsequent missions (Snell et al. 2000; Bergin et al. 2000; Caselli et al. 2012).

While the overall picture of high ice and low gas-phase water abundances in cold clouds appears well established, there are only a few sources for which both ice and gas have been observed along the same line of sight. Using the Infrared Space Observatory, the column densities of ice and warm water vapour have been determined through mid-infrared absorption lines toward a dozen high-mass infrared-bright sources (van Dishoeck & Helmich 1996; Boonman & van Dishoeck 2003). The warm water absorption lines originate in the inner envelope where the dust temperature is above the water ice sublimation temperature, which is of the order of ~100 K (e.g., Fraser et al. 2001; Burke & Brown 2010), and column densities comparable to those of water ice have been found. The infrared absorption data of warm water have subsequently been combined with SWAS submillimetre emission lines of cold water of the same high-mass sources to infer the gas-phase water abundance profile in both the cold and the warm gas (Boonman et al. 2003). Using different trial abundance structures applied to physical models of the sources constrained from continuum data, a jump in the gas-phase water abundance of ~4 orders of magnitude from cold to warm regions was established. Using a standard ice abundance of 10^-4, the ice column in these models was found to be a factor of 3–6 above the observed values. Although no full gas-grain model was adopted, even this simple empirical analysis showed difficulties in getting the water gas and ice chemistry to be consistent.

With the increased sensitivity of ground- and space-based infrared and submillimetre instrumentation, the combined study of water gas and ice can now be extended to low-mass protostars. In particular, ESO-VLT 3 μm and Spitzer mid-infrared water ice spectra exist for about 50 infrared-bright low-mass protostars (e.g., van Broekhuizen et al. 2005; Boogert et al. 2008), and the Herschel Space Observatory (Pilbratt et al. 2010) has observed submillimetre lines of water vapour of a comparable sized sample. Here we investigate the overlapping set of eleven low-mass protostars.

The goal of this paper is to constrain the relative importance of the main processes that shape the water vapour and ice abundance profiles in the cold parts of protostellar envelopes. To this end, a Simplified Water Network (SWaN) is developed, which is a dedicated tool that only incorporates the key processes that control the water gas and ice abundances in regions with temperatures ≲100 K. This model is then combined with physical envelope models from Kristensen et al. (2012). Together these data provide insight into the cold water chemistry as well as the puzzle as to why water ice occupies only a minor fraction of the available oxygen (Whittet 2010).

The paper is structured as follows. In Sect. 2 we introduce the Herschel HIFI observations of our sample, and give a brief overview of the observations and existing physical models of the sources. In Sect. 3 we introduce our simplified network (with the benchmarking results against full networks shown in the Appendix B) and analyse its sensitivity to key parameters. In Sect. 4 we compare the models with the observations of water ice column densities and Herschel spectra, which is followed by a discussion on the important parameters in Sect. 5 and the conclusions in Sect. 6.

2. Observations and physical models

2.1. Observations

Our observations of water vapour lines were obtained within the framework of the Herschel key project Water in Star-forming Regions with Herschel (WISH; van Dishoeck et al. 2011). The sub-sample analysed here was specifically selected to contain observations of both water vapour and water ice column densities. Thus, our target list is limited to eleven infrared-bright low-mass protostars (two Class 0, nine Class I; Table 1). The observed column densities of water ice range from 2.5 × 10¹⁸ to 1.5 × 10¹⁹ cm^-2 (Boogert et al. 2008; Zasowski et al. 2009; Aikawa et al. 2012). Since the water ice column density is estimated by observing water ice absorption in MIR spectra towards the embedded protostar, these column densities are always upper limits for the ice content of the envelope. Other contributions can be provided by, e.g., foreground clouds or disks.

Herschel observations are performed with the Heterodyne Instrument for the Far-Infrared (HIFI de Graauw et al. 2010), which provides line profiles of cold water vapour down to an unprecedented sensitivity and spatial resolution and opens up the spectral window to observe emission lines from higher excited levels. For all eleven protostars, observations have been performed in the ortho-ground state transition H₂O 1₁₀−1₀₁ at 557 GHz. For most sources there are also data of the para ground-state H₂O 1₁₁−0₀₀ at 1113 GHz, plus two additional excited transitions – H₂O 2₀₂−1₁₁ at 988 GHz, and H₂O 2₁₁−2₀₂ at 752 GHz. A list of observing dates and IDs is found in Table D.1 in the Appendix. It should be noted that there are two more sources with observations of both the water ice and water vapour, but they are excluded from our sample due to contamination with foreground absorption (Elias 29) and outflow emission (GSS30-IRS), which do not allow to draw any conclusions from their Herschel observations.

For the data reduction, the two polarisations (H, V) are combined and corrected for main beam efficiency as described in Kristensen et al. (2012). For low-mass sources, the HIFI lines are dominated by broad features due to the outflow, which is not of interest here. A Gaussian decomposition following Mottram et al. (2014) has been be used to subtract the contribution of the water emission from the outflows and spot shocks. The continuum-subtracted spectra alongside the best-fits for the non-envelope contribution are depicted in Fig. 1.

2.2. Physical models

To analyse the water gas and ice data, a physical model of the protostellar envelope is needed which specifies the density and temperature profiles of the protostellar envelope. Spherically symmetric model fits for each source have been made by Kristensen et al. (2012). Following the procedure by Jørgensen et al. (2002) the density profile is characterised by a power-law density structure of the form (1)The dust temperature has been determined self-consistently by performing a full continuum radiative transfer calculation assuming a central source with the observed luminosity as input. The best fitting values of α, envelope mass, and envelope extent have been obtained by comparison with the spectral energy distribution and sub-millimetre continuum images of the sources. The envelope is truncated at the point where either the dust temperature is <10 K or the hydrogen number density is <2 × 10⁴ cm^-3, whichever comes first. These points mark the transition between the envelope and the ambient cloud, which in our simulations is assumed to be chemically inert and devoid of water. The gas temperature is mostly coupled to the dust temperature, but following the approach by Bruderer et al. (2012) we take into account that in regions with elevated UV radiation field T_gas>T_dust.

Fractional abundances with respect to hydrogen nuclei of various species X_i ≡ n_i/n_H are specified at each radius, r, measured from the centre of the core. Alternatively, the visual extinction, A_V, measured from the outer edge of the envelope, r_env, can be used to describe the chemistry, acknowledging the fact that substantial changes in abundance profiles are introduced with the attenuation of the interstellar radiation field (ISRF). The extinction in the radial direction is obtained from the models through where the conversion factor in the denominator is taken from the empirical determination by Bohlin et al. (1978) and Rieke & Lebofsky (1985).

3. Model abundances in protostellar envelopes

3.1. Simplified Water Network (SWaN)

With the aid of observations of both water vapour and ice and a prescription for the temperature and density structure of the sources, we attempt to understand the key processes that shape the water abundance profiles in the cold parts of the protostellar envelopes. For this purpose a Simplified Water Network (SWaN) is developed, which is reduced to a minimum set of ingredients and reaction channels needed to reliably determine the abundance structure of water vapour and water ice. A comparison of SWaN to more sophisticated chemical networks (Visser et al. 2011, V11; Albertsson et al. 2013, A13; Walsh et al. 2013, W13) can be found in Appendix B.

Previous studies used step- or drop-abundance profiles for water vapour (Herpin et al. 2012; Coutens et al. 2012, 2013). These profiles are characterised by distinct regions of constant abundance. In the case of drop-abundance profiles, these are the inner region (T ≳ 100 K) with X_H2O ~ 10^-5, the outer region with X_H2O ~ 10^-8, and a photodesorption layer at the envelope edge with X_H2O ~ 10^-7 (Coutens et al. 2012). The location of the transition between the inner and outer region is placed at T ~ 100 K, but the extent of the photodesorption layer is a-priori unknown and has to be estimated by other means. The advantage of SWaN over these phenomenological profiles is that we have full control over the parameters that shape the abundance profile, and the photodesorption layer and its extent comes as a natural consequence of the interplay of the physical and chemical processes. The strength of these processes and their relative importance can be assessed, which allows us to study the link between water gas and ice. Our approach is similar to the work by Caselli et al. (2012) and Keto et al. (2014) on water vapour abundance profiles in prestellar cores, except our study with SWaN focuses on the understanding of the connection between water gas and ice in the cold, outer envelope of protostellar cores, which are characterised by a temperature increase towards the centre.

Our simplified network, SWaN, consists of three species (water ice on the grain surfaces, water vapour, and atomic oxygen), which are connected by four reaction channels (Fig. 2). In the following description of the chemical network, all reaction rates R_i correspond to the number of atoms/molecules that are transformed through a particular reaction channel per unit volume and unit time.

Water ice can desorb through far-ultraviolet (FUV) photodesorption at a rate (2)where σ_H is the grain cross section per hydrogen nuclei, and is given by σ_H = σ_grn_gr/n_H (with the grain cross section , the grain radius a_gr, and the grain volume density n_gr), F_FUV is the flux of FUV photons at the grain surface¹, and Y_pd is the photodesorption yield. Following recent lab results (Öberg et al. 2009; Bertin et al. 2012), FUV photodesorption is treated as a zeroth-order process (i.e., molecules can only desorb from the top few layers). Öberg et al. (2009) give a photodesorption yield of (3)where Y_pd,0 = 10^-3 is the photodesorption yield for thick ice, and θ_M is the monolayer coverage factor (4)where l = 0.6 is the diffusion length. M is the number of monolayers, which is given by (5)where N_s = 1.5 × 10¹⁵ cm^-2 is the density of sites on the grain surface (Hasegawa et al. 1992). This approach is supported by molecular dynamics simulations (Andersson & van Dishoeck 2008; Arasa et al. 2010). In reality, FUV absorption of water ice results not only in H₂O but also OH molecules escaping from the ice (Andersson & van Dishoeck 2008; Öberg et al. 2009), but the latter channel is not taken into account explicitly.

In our model, the flux of FUV photons consists of two components F_FUV = F_FUV,isrf + F_FUV,cr, where the first term is the contribution of FUV photons from the ISRF, and the second term reflects the secondary UV field caused by cosmic rays interacting with molecular hydrogen (Prasad & Tarafdar 1983). The flux of ISRF FUV photons penetrating the ice mantle at a certain depth into the cloud can be calculated via (6)with G_isrf as the scaling factor for the standard ISRF photon flux F₀, and as the spherically averaged extinction for photodesorption (Appendix A). The flux of secondary FUV photons induced by cosmic rays is independent of extinction, and results in an isotropic, constant FUV photon flux (7)with G_cr depending on the assumed energy distribution of the cosmic rays (Shen et al. 2004). Other mechanisms such as direct cosmic-ray desorption have no significant impact (Hollenbach et al. 2009), and are therefore not included in SWaN. Chemical desorption using the excess energy produced by water ice formation has also been proposed as a desorption mechanism (Dulieu et al. 2013) but is still poorly understood qualitatively and quantitatively and is therefore neglected.

Desorption of water ice from the grain surface can also occur through thermal desorption (e.g., Visser et al. 2011) at a rate (8)where ν_H2O is the lattice vibrational frequency of a water molecule in its binding site, T_b,H2O is the binding energy of water expressed as temperature, and T_dust is the dust temperature. For reasons of simplicity, the same monolayer coverage θ_M as for the photodesorption is also adopted for thermal desorption. The lattice vibrational frequency of water ν_H2O = 2.8 × 10¹² Hz is calculated through the harmonic oscillator approach of Hasegawa et al. (1992).

Water vapour can be photodissociated through FUV photons, again taking into account contributions from both the ISRF and the CR-induced secondary field, at a rate (9)where k_phdis is the unshielded photodissociation rate of water in a FUV field with G_isrf = 1, and is the spherically averaged mean extinction through the envelope for photodissociation (see Appendix A). It is assumed that all the water, which is photodissociated, is turned into atomic oxygen, essentially leaving out the intermediate product OH. The slightly different exponential dependence on extinction in Eq. (6) vs. (9) arises in the FUV absorption spectrum of water ice, which is shifted to higher energies by about 1 eV compared with water vapour (Andersson & van Dishoeck 2008).

Water ice is formed directly through the freeze-out of water vapour, or indirectly through freeze-out of atomic oxygen. In this second step, intermediate steps through O₂ are ignored, but atomic oxygen is instantaneously hydrogenised to form water ice. The benchmarking (Appendix B) reveals that the formation timescale of water ice can be longer than the freeze-out timescale of atomic oxygen, but also that the abundance profiles of water vapour and water ice are only marginally affected. For species i, freeze-out occurs at a rate of (10)where n_i is the number density, v_i is the thermal velocity, and S_i the sticking probability. For water vapour a sticking probability of unity is assumed. For atomic oxygen we introduce an effective sticking probability S_O, which is determined by the balance between freeze-out and thermal desorption. We take the relative reaction rates k_fr,O = σ_Hn_Hv_i and k_thdes,O = ν_O e^{−T_b,O/T_dust} to determine k_thdes,O/k_fr,O. For the protostellar cores in our sample this generally means that in regions with T_dust ≳ 15 K the thermal desorption rate is higher than the freeze-out rate, i.e., oxygen atoms cannot freeze out, and the formation of water ice through the atomic oxygen route is inhibited.

The number densities n_O(t), n_H2O(t), and n_{s - H2O}(t) are determined by a set of three differential equations: These are solved with the aid of the Python function odeint, which is part of the scipy.integrate² package and makes use of the Fortran library odepack. The standard model parameters are summarized in Table 2.

All models use the same parameters as in Table 2, but nevertheless show significant differences, both among each other and with SWaN. In spite of these uncertainties, the derived abundance structures for water vapour and water ice are robust in all models. In contrast, atomic oxygen in our simple network is very different from the detailed networks and only serves as a proxy for other oxygen-bearing species within the full water-chemistry network (e.g., OH, H₂O₂).

3.2. Water ice and gas in a representative model

Using SWaN as described in the previous section, the abundance profiles of water vapour and water ice can be determined. The overall volatile oxygen abundance, i.e., the oxygen not contained in silicate grains, is taken to be X_O,ISM = 3.2 × 10^-4, or 320 ppm³, as determined for diffuse clouds (Meyer et al. 1998). This is the amount of oxygen that can cycle between the various forms of oxygen-containing gas-phase and ice species. Of particular interest for this work is the cold, outer part of the envelope where water is frozen out on the dust grains, henceforth denoted as the water freeze-out zone.

Our model consists of two stages:

• The pre-collapse or prestellar phase, where a full chemical network under the assumption of dark cloud conditions (W13 with T_dust = 10 K, A_V = 10 mag, n_H = 2 × 10⁴ cm^-3) sets the initial molecular abundances for the second step. Table 3 lists abundances of water vapour, water ice, and atomic oxygen for different pre-collapse times. The freeze-out of atomic oxygen is closely connected to the steady rise of the water ice abundance. At times 0.1 − 1.0 Myr, a considerable amount of oxygen is also found in other oxygen bearing species (mainly CO), reducing the amount of oxygen within the water chemistry network. At t_pre> 1 Myr, oxygen returns into the water network and water ice then becomes the only considerable oxygen reservoir.

• The second stage is the so-called post-collapse phase, where the abundance structure is modelled with SWaN on a static envelope with the temperature and density structure from Kristensen et al. (2012).

In our representative model a pre-collapse time of t_pre = 0.1 Myr is chosen, which is motivated by our finding that a rather short pre-collapse time is required to match the observed low water ice abundances (see Sect. 4.2). The initial abundance of oxygen within SWaN at that timestep is X_O,SWaN = X_H2O + X_{s - H2O} + X_O = 270 ppm. The remaining X_O,ISM − X_O,SWaN = 50 ppm can be attributed to oxygen locked up in CO, s-CO and other oxygen species not contained in the simple network. The post-collapse timescale in this representative model is chosen to be t_post = 1.0 Myr, which is about the time when equilibrium is reached at all radii. As fiducial values for the FUV fluxes in this simulations, G_isrf = 1 and G_cr = 10^-4 are used.

The abundance profiles for water vapour, water ice and atomic oxygen vs. the radial extinction A_V are shown in Fig. 3 for L1551-IRS5, which is chosen as a representative core from our sample. The overall structure of the water vapour profile can be roughly separated into four different regions, which are denoted as Regions A to D. Depending on the choice of parameters and the envelope temperature and density structure the extent of these regions can vary.

The outermost layer (Region A) of the protostellar envelope is dominated by atomic oxygen. Water ice and water vapour are only present in trace amounts, since the impinging ISRF efficiently photodissociates water molecules, but due to the long freeze-out timescale in the tenuous envelope, the regeneration of water ice is significantly slower. The attenuation of the ISRF deeper into the core leads to a build up of water ice and water vapour. When the grains are fully covered with water molecules (in this model run at around A_V ~ 0.7 mag), the water vapour abundance reaches a plateau (Region B). The extinction value at which this plateau is reached strongly depends on G_isrf (cf. Sect. 3.3.1) The reaction network in these two outermost regions is characterised by the cycle O → s - H₂O → H₂O → O.

That changes deeper into the cloud (Region C). Due to the increased number densities in concert with decreased FUV fluxes, the freeze-out of water vapour starts to outrun the photodissociation. The reaction network reduces to H₂O ⇆ s - H₂O, where the water abundances are determined by the balance between freeze-out of water vapour and photodesorption of water ice. At A_V ≳ 4−5 mag (assuming standard values for G_isrf and G_cr), the CR-induced FUV field is the main contributor of FUV photons, which maintain an approximately constant number density of gas phase water. Since the hydrogen volume density in this region spans almost three orders of magnitude, the water vapour abundance drops considerably from the outer towards the inner edge of Region C. Water is predominantly in the form of water ice. For most of our sources, 50% of the pencil beam water ice column density can be found within the central 20−80 au, which is typically around 1% of the envelope extent.

The water ice abundance profile shows a characteristic hump at A_V ~ 2−3 mag, which is the result of the freeze-out threshold for atomic oxygen. Only in the oxygen freeze-out zone (T_dust ≲ 15 K), can atomic oxygen be transformed further into water ice. Typically, ~5−15% of the column density in the water freeze-out zone is also part of the oxygen freeze-out zone. However, deep in the envelope (T_dust ≳ 15 K) the water ice abundance remains constant at its initial value. No water ice can form, but there are also no mechanisms which can considerably decrease its abundance.

Even deeper in the envelope, the dust temperature reaches the sublimation temperature of water ice, which consequently marks the boundary of the water freeze-out zone. With SWaN we do not attempt to model any abundances beyond this point (Region D). Thermal desorption only plays a role in Region D. Nevertheless, this reaction channel is included in our simplified network, since it allows us to determine the extent of the water freeze-out-zone. In our simulations, the limit of the water freeze-out zone is defined as the radius at which the water ice abundance drops by three orders of magnitude from its plateau abundance. This is typically at temperatures of T_dust = 110−140 K, and its location can – depending on the envelope density and temperature structure – vary by up to a few au compared to the 100 K-radius, which is generally the first order assumption for the sublimation radius.

The water vapour abundance profiles have mostly reached an equilibrium stage at t_post = 0.01 Myr. In the outer envelope (Region A) this is established through the high photodesorption and photodissociation rates, and deeper into the envelope (Region C) through the high freeze-out rates. In the transition region between Regions B and C the situation is a bit different. The hydrogen number densities are ~10⁵ cm^-3 (which results in low freeze-out rate), but the ISRF FUV photon flux is already considerably attenuated (which results in low photo-rates). Therefore, in this region it takes much longer (~0.1–1.0 Myr) to reach an equilibrium stage. However, the discussed timescales can change upon the choice of initial conditions. It should be noted that this region does not affect the total ice or gas column density (Sect. 5.2), but it does affect the water emission line profiles (Sect. 3.3.2).

3.3. Dependence on model parameters

3.3.1. Water abundances and column densities

The focus of this work is to determine the key factors that shape the water gas and ice abundance profiles and regulate their column densities, with particular attention to its dependence on the FUV photon fluxes and initial abundances. In the following set of simulations, the initial abundances of a pre-collapse time of 0.1 Myr are used (Table 3). The initial overall oxygen abundance in the system is X_O,SWaN = 270 ppm. About 7% of the oxygen is found in the form of water ice, but for demonstrative purposes we use a fiducial value for the initial fraction of water ice of ξ_{s - H2O} = 10% to allow easier scaling in order of magnitudes. The fiducial values for the FUV fields are G_isrf = 1 and G_cr = 10^-4. In three different runs, each parameter is scaled by a factor of f = 10^-2−10¹, leaving the other two parameters at their fiducial values.

An enhanced G_isrf can be due to UV from a nearby bright star but potentially also from from fast shocks related to the protostellar jet impinging on the envelope. An enhanced ISRF can lead to an increased gas temperature (e.g., Hollenbach & Tielens 1997). A lower than standard value of G_isrf could arise from shielding by low density gas from the surrounding molecular cloud. The value of G_cr is linked to the spectrum of cosmic rays and the cosmic ray ionization rate (Webber & Yushak 1983; Shen et al. 2004). Variations in the initial fraction of water ice can be caused by a different pre-collapse time: the longer the pre-collapse time, the higher the fraction of water ice (Table 3).

Figure 4 shows the resulting abundance profiles (left panels) and column densities (right panel) for water ice and water vapour at a post-collapse time of t_post = 1 Myr. Because the focus of this paper is on the cold chemistry, the column densities are not computed throughout the entire envelope, but rather in the water freezeout-zone. This particularly affects the column density of water vapour, which is abundant in the inner hot core (Boonman et al. 2003).

Variations of the ISRF shape the abundance profiles in the outer envelope (Fig. 4, left panel). The extinction threshold for the appearance of the water vapour abundance plateau, A_V,f, which is a result of the formation of a first monolayer of ice around the grains, strongly depends on G_isrf. This behaviour, which has already been described by Hollenbach et al. (2009), is a result of the attenuation of the ISRF, which is accompanied by decreased photodesorption and photodissociation rates of water ice and water vapour, respectively. At A_V ≳ A_V,f the desorption turns effectively from a first order (all water molecules on the grain surface can desorb) into a zeroth order process (where only the top layers contribute to desorption).

Variations in the CR-induced FUV field lead to changes in the abundance deep in the envelope (Fig. 4, centre panel). Due to the reduction of the network to H₂O ⇆ s - H₂O, the water vapour abundance scales practically 1:1 with the G_cr. However, the abundance is too low to considerably affect the water ice abundance.

Scaling the initial fraction of water ice has mostly an effect on the abundance of water ice deep in the envelope (Fig. 4, right panel). The temperatures of T_dust ≳ 15 K inhibit the formation of water ice through freeze-out of atomic oxygen. This has as a result that the abundance, and thus the water ice column density, is already imprinted by the initial conditions (cf. Sect. 3.2). The effect of a higher binding energy of O is discussed in Sect. 5.2.

The column densities of water vapour are only mildly affected by parameter variations within our parameter space. It should be noted that a combination of low FUV fluxes from both the ISRF and the CR-induced FUV field can indeed cause a reduction of this column density by a factor of 2−3. However, the ice column density remains hardly affected by the choice of FUV photon fluxes, and envelope-averaged gas-to-ice ratios of ~10^-4 within the water freeze-out zone are found. Variations of the initial fraction of water ice ξ_{s - H2O} only affect the water ice. Since ξ_{s - H2O} varies with the length of the pre-collapse time, the observed water ice column densities strongly depend on its initial abundance at the beginning of the post-collapse phase.

3.3.2. Water emission/absorption profiles

In contrast to the observations of water ice, which are derived from the attenuation of the light of a background source (and thus, only dependent on the column density along the line-of-sight), the spectrally resolved water vapour lines allow us to model the underlying number density and velocity profile of the envelope. Due to the large beam size of Herschel, the observed spectra do not show the pencil beam spectrum towards the core centre, but the major contribution actually originates in lines-of-sight that intersect the envelope at impact parameters of up to a few thousand au. Therefore, the spectrum has to be modelled through radiative transfer tools to take into account the complex interplay of excitation conditions, abundances, and velocities within the sampled area on the sky. A simple absorption study will at best only yield a lower limit on the water vapour column of typically 10¹³ cm^-2, even for lines such as the p-H₂O 1₁₁−0₀₀ line at 1113 GHz which are primarily in absorption (Kristensen et al. 2010).

To demonstrate the influence of various parameters on the emission profiles, we model the o-H₂O 1₁₀−1₀₁ ground state line of L1551-IRS5 and explore its dependence on FUV fluxes, post-collapse time, and velocity profile with the aid of Ratran (Hogerheijde & van der Tak 2000), following a similar approach as Mottram et al. (2013). The initial abundances are determined from our dark cloud model with a standard pre-collapse time of t_pre = 0.1 Myr. Figure 5 shows synthetic spectra for four different parameter combinations. For Model A, the post-collapse time is set to t_post = 1.0 Myr, and it assumes high FUV fluxes of G_isrf = 10² and G_cr = 10^-3, a static envelope with zero radial velocity, and a constant turbulent broadening, which is characterised by Doppler-β = 0.4 km s^-1. The synthetic spectrum shows poor agreement with the observations. In Model B, the FUV fluxes are reduced to G_isrf = 10⁰ and G_cr = 10^-4. This results in a shift of the abundance peak to lower extinctions. In contrast to the line flux, which is in good agreement with the observations, the skewed line profile cannot be reproduced. Including expansion motions in Model C results in a skewed line and a good match of the absorption feature, but a poor match to the peak emission. Only a further decrease of the FUV photon fluxes result in good agreement with the the observations (Model D). This example highlights the complex interplay of abundance structure and velocity profile to shape the emission and absorption lines as observed with Herschel.

4. Comparison with observations

4.1. Water vapour

The observed line profiles exhibit a remarkable diversity (Fig. 1), even though they all share the same overall structure of the water vapour abundance profiles. As shown in Sect. 3.3.2 (Fig. 5), the interplay of the abundance and velocity structure has a significant impact on shaping the line profiles.

In the following parameter study, we vary parameters that actively influence the shape of the water vapour abundance profile (the post-collapse time t_post, and the FUV fluxes G_isrf and G_cr). In addition, the radial velocity distributions, namely the velocity centroid and the Doppler broadening, are investigated. For the sake of simplicity, only Hubble-like radial infall/expansion velocities are considered, which are of the form (14)where r_env is the envelope radius, and v_r,env the velocity at this point. Negative velocities mark infall, positive values represent expansion. To account for possible discontinuities in the Doppler-β distribution, for which hints have been discovered in both low- and high-mass cores (Herpin et al. 2012; Mottram et al. 2013), various cases of Doppler-β distributions are tested. Firstly, a constant Doppler-β for the whole envelope is tested. Secondly, an extinction threshold A_V,j, which marks a jump in the Doppler-β distribution, is introduced. Doppler-β is only varied in regions A_V ≥ A_V,j, whereas in the outer regions A_V<A_V,j a constant β = 0.2 km s^-1 is chosen, which is motivated by the presence of narrow absorption features in some sources (e.g., RCrA-IRS5). Our parameter space is therefore defined by .

The initial abundances are fixed after a pre-collapse time of t_pre = 0.1 Myr, since it has been shown in Sect. 3.3.1 that the choice of the pre-collapse time (i.e., the inital abundance of water ice ξ_{s - H2O}) does not considerably affect the water vapour abundance profiles in this study.

Owing to the line shape complexity and the plethora of parameters, we do not attempt to find a singular “best fit” point in our parameter grid. We rather aim to find general trends in the solutions, analysing the influence of each parameter individually on the overall fit quality through averaging over the values (i.e. marginalisation) of all other parameters. This is pursued by adopting a Bayesian approach, which is summarised in more detail in Appendix C. A parameter p can take m values V = {V₁,...,V_m} (Table 4). The overall fit quality of our model with a particular value V_i is quantified by its evidence E_p(V_i). This number becomes only meaningful when the evidences for all possible parameter values V are compared to each other. For each value V_i we therefore normalize the evidence by defining the Bayes Factor as B_p(V_i) = E_p(V_i) / max [E_p(V)]. The parameter value with maximised likelihood has therefore a Bayes Factor of unity. Generally, it is assumed that a parameter value can be rejected if B_p(V_i) ≪ 0.1. It should be noted that this approach is only a relative comparison of the grid points in the parameter space, but does not quantify the absolute quality of the fit.

The envelope models of Kristensen et al. (2012) do not take into account any deviation from spherical symmetry (e.g., due to the presence of a disk) at scales of ≲300−500 au. Therefore, we focus on the outer regions of the envelope, and fit only ortho and para ground state lines H₂O 1₁₀−1₀₁ and H₂O 1₁₁−0₀₀.

A summary of the parameters that yield the best representation of the observed line profiles for all sources is presented in Table 5, and the spectra are shown in Fig. 6. In some cases the overplotted best-fit spectra seem to be poor fits, but this is owed due to the fact that we have chosen a coarse grid in our Bayesian analysis plus an overall complexity in the source structure, which is heavily simplified by our assumption of a radially symmetric envelope. Instead of fitting every detail of the spectral lines, the goal of our analysis is to find global trends, i.e., if a certain region in the parameter space yields significantly better fits. In particular the infall profiles of the Class 0 sources lack the deep absorption feature in the H₂O 1₁₀−1₀₁ line. Mottram et al. (2013) showed that a detailed modelling of the influence of the absorption against the outflow is able to recover the absorption depth, but this kind of modelling is beyond the scope of this paper. Moreover, L 1551-IRS5 and HH 46 show emission peaks in the higher excited lines, which we attribute to the presence of structure in the inner regions that is not part of our radially symmetric envelope model.

A few general trends for our sample of protostars can be derived. All sources (except RCrA-IRS5) are characterised by a chemical age for water of around 1 Myr, but in some, an age of 0.1 Myr cannot be ruled out. Typically, our sources are exposed to low CR-induced FUV fields (G_cr ≲ 10^-4). Mostly, we also find weak ISRF (G_isrf ≲ 1). In terms of velocity, the two Class 0 sources are characterised by infall motions, and Class I sources show expansion motions of the envelope. However, due to weak emission, the velocity field is hard to constrain in some sources. There seems to be no general trends for the turbulent broadening (Doppler-β), but a few sources show better results when including a discontinuity in the Doppler-β-distribution.

4.2. Water ice

As discussed in Sect. 3.3.1, the final abundance of water ice in our models is practically independent of the FUV photon fluxes over the parameter space considered, but strongly depends on the initial oxygen abundances (namely, the fraction of water ice ξ_{s - H2O}, but also the total abundance of oxygen X_O,SWaN in the chemical network). In our subsequent modelling, this parameter is taken from a dark cloud model with different pre-collapse timescales t_pre = 0.01, 0.1 and 1.0 Myr (Table 3). Then, our simplified chemistry is run for post-collapse time of t_post = 0.1 and 1.0 Myr. To determine the abundance from the column densities, both the observed and simulated water ice column densities are divided by the column density of hydrogen atoms in the water freeze-out zone, .

We note that the abundances strongly vary with pre-collapse time, but varying the post-collapse from 0.1 to 1.0 Myr makes only a marginal difference. A comparison of the observations and the modelled abundances is found in Fig. 7. In general, water ice abundances range from X_{s - H2O} = 30−80 ppm, except for TMC1 and RCrA-IRS5. It should be noted that the observed water ice column densities are always upper limits for the column density that is part of the envelope. In some cases, a considerable amount of water ice column density could originate from elsewhere along the line-of-sight (e.g., foreground clouds, disks). The exceptionally high water ice abundance in TMC1, e.g., might be explained by the presence of a large disk (Harsono et al. 2014), which is intersected by the line-of-sight towards the central protostar. In that case, the major fraction of the absorbing water ice column in TMC1 would originate in the disk rather than the envelope. The disks around TMC1A and TMR1, which are also in the sample of Harsono et al. (2014), seem to not affect the water ice abundances, which lets us suggest that the line-of-sights towards these protostars are pristine, and truly only through the envelope.

The modelled abundances show strong source-to-source variations in the abundance of water ice for a short pre-collapse time, 0.01 Myr. These variations become smaller with increasing pre-collapse time, and almost disappear for t_pre = 1.0 Myr. To understand this behaviour, one has to understand that in the prestellar core phase, water ice slowly builds up through freeze-out of atomic oxygen. This build-up is stopped in the post-collapse phase, where parts of the envelope are too warm to allow freeze-out of atomic oxygen. The observed differences in the abundance structure on short pre-collapse timescales therefore reflect the volume density and temperature structure of the individual protostellar envelopes. In the envelopes of L1551-IRS5, HH 46, TMC1A, TMR1, and L 1489, the bulk of material in the post-collapse stage is found at T_dust ≳ 15 K where water formation through freeze-out of oxygen is inhibited. The ice column density in these sources is, therefore, almost completely imprinted during the pre-collapse phase. Other sources (L 1527, IRAS 15398, TMC1, RCrA-IRS5, RNO91, and HH 100) have a post-collapse envelope with a large atomic oxygen freeze-out zone, and the transformation of oxygen to water ice can continue throughout the post-collapse stage.

In contrast to older sources with abundances ≳60 ppm, young sources (with the exception of TMC1) are generally characterised by abundances 20−50 ppm. To speak of an evolutionary effect would be an overinterpretation of the data, but it is clear that the observed low abundances require a rather short pre-collapse phase of ≲0.1 Myr. This is an apparent contradiction to the estimated lifetime of prestellar cores of ~0.5 Myr (Enoch et al. 2008), which will be discussed in Sect. 5.2.

4.3. The connection between water gas and ice

From the radiative transfer modelling in Sect. 4.1 the likelihood-averaged column densities are determined. Since the scatter in the derived values is mostly low, the error bar is defined by the separation of the two column density values in our discrete grid of parameters, which are closest to the best-fit value. The water vapour column densities are on the order of 10¹⁴−10¹⁵ cm^-2. These are significantly lower than the water ice column densities of 10¹⁸−10¹⁹ cm^-2, and the resulting gas-to-ice ratios are on the order of 10^-5−10^-4. Correlations of water vapour and ice with N_H and T_bol are depicted in Fig. 8. The only significant correlation is found for column densities of water ice and hydrogen. Water vapour, on the other hand, is largely uncorrelated with any of these parameters. This only reflects the nature of water ice as a bulk tracer, whereas water vapour rather traces the surface layers of the envelope. That obviously also means that the gas-to-ice ratios, which vary by up to an order of magnitude from source to source, are an intrinsic property of each individual source rather than a global indicator.

5. Discussion

5.1. Water abundances with SWaN

In Sect. 3.1, our simplified water network SWaN is introduced. Through benchmarking with other chemical codes, the chemical network can be limited to a small number of species and reaction channels to understand the connection between water gas and ice, and to reliably predict the abundance profiles for both these species in the cold regions of the protostellar envelope. The only species that needs to be added is atomic oxygen, which takes over the role as a proxy for other oxygen bearing species in the full water-chemistry network.

Simplifying the network to O → s - H₂O → H₂O → O in the outer envelope and to H₂O ⇆ s - H₂O further in, results in equilibrium abundance profiles that are in good agreement with sophisticated chemical networks. Photodesorption, photodissociation and freeze-out are sufficient to explain the abundance structure of water vapour and water ice in the water freeze-out zone.

5.2. Water ice

A key ingredient in understanding the observed water ice column densities is a low initial abundance of water ice after the pre-collapse phase in concert with the presence of a freeze-out barrier for atomic oxygen, caused by a relatively low binding energy.

The initial abundances for the post-collapse phase are determined through a dark-cloud model (T_dust = 10 K, A_V = 10 mag, n_H = 2 × 10⁴ cm^-3; W13). The formation of water ice is controlled by the freeze-out of atomic oxygen. This freeze-out timescale, and thus the water ice formation time scale, depends strongly on the product of the hydrogen number density, n_H, and the grain cross section per hydrogen atom, σ_H (Eq. (10)), which in turn equals n_grσ_gr. Through a reduction of either n_gr or σ_gr, or both, the formation of water ice can be slowed down significantly. This could eradicate the apparent need for a short pre-collapse time of ~ 0.1 Myr that was found in our model, which conflicts with the observed prestellar core lifetime of ~0.5 Myr (Enoch et al. 2008).

Our adopted grain abundance of X_gr = 6.5 × 10^-13 combined with a grain size of 0.1 μm corresponds to σ_H = 2.0 × 10^-22 cm² in our standard model. This is almost an order of magnitude lower than the canonical value of σ_H = 1.0 × 10^-21 cm² derived from diffuse cloud observations (e.g. Prasad & Tarafdar 1983), but consistent with dense cores in which grains have grown to somewhat larger sizes as predicted theoretically and found observationally (e.g., Pagani et al. 2010; Steinacker et al. 2010). Growth to even bigger sizes than assumed here would be needed to bring the two timescales into agreement.

The lifetime estimate of Enoch et al. (2008) refers to the dense prestellar cores with mean densities above 2 × 10⁴ cm^-3, to which their millimetre continuum data were sensitive. Water ice formation can already start to form in lower density, translucent clouds with densities of a few 10³ cm^-3 (Whittet et al. 2001; Cuppen & Herbst 2007). The time that the cloud spends in this low density phase is unknown, but it would only add to the amount of water ice after the pre-collapse phase, exacerbating the need for a short dense, prestellar phase. Our short inferred pre-stellar phase is in contrast with Yıldız et al. (2013b) who argued for a long prestellar phase of at least 1 Myr at n_H = 10⁵ cm^-3 to explain the absence of gas-phase O₂ toward the NGC 1333 IRAS4A protostellar core. Since water ice is not observed directly toward this core, it is not clear whether there is a similar discrepancy.

A key parameter is the freeze-out barrier for atomic oxygen, which helps to maintain the low water ice abundance even after the pre-collapse phase. Owing to the lack of laboratory measurements for the binding energy of atomic oxygen on water ice or other surfaces, a binding energy of T_b,O = 800 K (Tielens & Hagen 1982) was assumed. This results in a freeze-out temperature of T_dust ~ 15 K. Recent theoretical work and laboratory data point to binding energies that could be more of the order of T_b,O = 1500 K on amorphous silicate or graphite surfaces (Bergeron et al. 2008; He et al. 2014), for which the freeze-out threshold temperature would be raised to around T_dust ~ 35 K. If such a high value would also apply to the binding of O on water ice, about 2−5 times more material (15−35% of the column density within the water freeze-out zone) would be found in a region where atomic oxygen can still be converted into water ice, effectively increasing the overall water abundance and thus requiring an even shorter pre-stellar phase.

Efficient cosmic ray desorption in the chemical model of Caselli et al. (2002) helped to retain a considerable amount of atomic oxygen in the gas phase. Their assumption of a low binding energy (T_b,O = 600 K) decreased the cosmic ray desorption timescale of atomic oxygen and resulted in desorption from the grain before further processing. However, this is an effect that we have not seen in our full chemical model benchmarking which include cosmic ray desorption. Even for T_b,O ≳ 800 K, the residence time on the surface is sufficient for hydrogenation reactions to convert atomic oxygen into more tightly bound molecules such as OH or H₂O.

As a result of this qualitative analysis, the apparent contradiction between the modelled and observed prestellar core lifetime remains. A possible way out of this is to assume that after the prestellar core phase of ~0.5 Myr there are mechanisms at play, which help to generate a water ice abundance that resembles the situation after a pre-collapse time of t_pre ~ 0.1 Myr. Our simple two-stage approach with constant density and temperature profiles in each of the stages can only be seen as a first approach, but should be followed by more detailed modelling of a self-consistent protostellar core collapse, taking into account the change in density, temperature and velocity structure. A key player in influencing the abundance profiles of water in a more realistic collapse scenario could be episodic accretion. During bursts in accretion the source luminosity can increase by ~2 orders of magnitude (Vorobyov et al. 2013), which heats the envelope and can increase the radius of the hot core by up to an order of magnitude (Johnstone et al. 2013). Our simplified network is not designed to simulate such a situation. Nevertheless, the temperature increase during an accretion burst could trigger chemistry which drives oxygen into other, more tightly bound species that are not included in our simple chemistry network.

5.3. Water vapour chemistry

The water vapour emission lines as observed with Herschel turn out to be an invaluable tracer for the FUV field and envelope kinematics (cf. Mottram et al. 2013). Amongst our sample, 8/11 sources are characterised by low ISRF-FUV fluxes G_isrf ~ 10^-2. Assuming that the host star-forming region of these sources is embedded in a standard ISRF with G_isrf = 1, a water-free cloud with an extinction of A_V ~ 2−3 mag would be needed to provide the required attenuation. For example, all sources in the Taurus Molecular Cloud (TMC1, TMC1A, L1551-IRS5, TMR1) show a trend towards G_isrf = 10^-2. This result is in agreement with the findings of an extinction threshold for water ice in Taurus at A_V ~ 3 mag (Smith et al. 1993; Teixeira & Emerson 1999; Whittet et al. 1988, 2001). Following Eq. (17) in Hollenbach et al. (2009), the formation threshold for water ice can be stretched to A_V = 2−3 mag for a tenuous ambient cloud with G_isrf = 1 and n_H ~ 10³ cm^-3, which can remain almost devoid of water due to the interplay of FUV field and low freeze-out rates. It would have an extent of around 1 pc, and extinction maps of the Taurus Molecular Cloud suggest that this is not unreasonable (Kainulainen et al. 2009).

In some sources (HH46, RCrA-IRS5, HH100-IRS – and with less significance TMR1), evidence for elevated ISRF FUV fields G_isrf ≳ 1 is found, which can originate in interaction of the outflow with the envelope or in an photodissociation region (PDR). Particularly the latter seems to be the case for HH100-IRS and RCrA-IRS5, which are known to be cocooned by the extended strong radiation field of the PDR in the RCrA star-forming region.

When it comes to assessing the CR-induced secondary FUV field, the observations are generally best represented by G_cr ≲ 10^-4, with a trend towards 10^-5. However, recalling Eq. (2), the photodesorption rate depends, amongst other things, also on the product σ_grn_gr. This product is proportional to the inverse grain radius a^-1, i.e., grain growth will result in a decreasing surface area per unit volume. Our measurements probe – assuming that all other factors are well known – the product σ_grn_grG_cr rather than G_cr alone. Therefore, low values of the CR-induced FUV field in our models could indeed be the result of a reduced cosmic ray flux deep inside dense cores (e.g., Padovani et al. 2013). Alternatively, it could also be due to dust growth to micron-sized particles deep in the envelope (e.g., Pagani et al. 2010; Steinacker et al. 2010). Unfortunately, these two effects cannot be distinguished by our analysis.

In terms of the post-collapse timescale, all sources show significant evidence for t_post = 1.0 Myr. In some sources, t_post = 0.1 Myr can be ruled out, but this does not correlate with evolutionary stage. In particular, the Class 0 source L 1527 is best represented by a post-collapse time of 1.0 Myr, which is an order of magnitude higher than the generally estimated lifetime of sources at that evolutionary stage of ~0.1 Myr (Evans et al. 2009). But that is only a contradiction at first glance. As seen earlier, the initial molecular abundances that match the observations of water ice are best represented by a pre-collapse time of 0.1 Myr (Sect. 5.2), and thus this was our choice for the determination of the initial abundances. But we have argued that some mechanism, maybe episodic accretion, helps to reset the chemical clock in the inner regions of the envelope to a situation that resembles a prestellar lifetime of ≲0.1 Myr. But the outer envelope would not be affected by this, and its chemical age would represent its true age. The water vapour profile in that region is therefore governed by the cumulative timespan of the prestellar stage (~0.5 Myr) plus the time the source has already spent in Class 0 phase. Therefore, even the water vapour profiles of the youngest Class 0 sources are characterised by a chemical age of ~1.0 Myr, rather than 0.1 Myr.

5.4. Kinematics

Our analysis shows that water emission lines are not only a good tracer for the FUV fluxes, but they also prove to be a sensitive and invaluable tracer of the kinematics of protostellar envelopes. In the two Class 0 sources we find strong evidence for infall motion. The best fit spectrum seems to be a poor representation of the observed spectrum (Fig. 6). However, the focus of this paper is on finding general trends (i.e., expansion vs. collapse) rather than determining the true nature of the collapse. A free-fall collapse model, which has been used by Mottram et al. (2013) on sources with inverse P-Cygni profiles, does indeed show better overall agreement with the observations. But quantifying the nature of the infall is beyond the scope of this paper.

In our sample of Class I sources, 7/9 sources are characterised by expansion motion of the outer envelope. In 2/9 the kinematics remain mostly unconstrained due to weak emission. This fits the picture of envelope dispersal during this evolutionary stage (e.g., Kristensen et al. 2012).

Our inferred radial motions raise the question of how reasonable the assumption of a static envelope is. Obviously, the determined density, temperature, and velocity profiles are only snapshots of the current situation, but have no information about their respective history. At expansion velocities of the order of km s^-1 as seen in L1551-IRS5 and HH100-IRS, an envelope with an extent of ~10⁴ au would have started from a singular point less than 0.1 Myr ago. Therefore, this expansion must have started fairly recent in time. However, our analysis shows that most other Class I sources show support for expansion, but it is impossible to quantify its rate. Therefore, the assumption of a static envelope in the post-collapse phase is our best approach, but it can only be seen as a first step towards a self-consistent evolution model.

5.5. The oxygen budget of protostars

The benchmarking of our simplified chemistry against full gas-grain models also provides insight into the total oxygen budget of protostellar envelopes. In Appendix B.2 the predicted reservoirs of oxygen in three other chemical networks are analysed. A snapshot of the abundances after t_post = 0.1 Myr (Fig. 9) shows that the two networks with complex grain-surface chemistry (A13, W13) efficiently drain atomic oxygen out of the system, and convert it into other species. This has consequences for the oxygen budget puzzle and the nature of the unidentified depleted oxygen (UDO; Whittet 2010). Following the trail of the oxygen in these complex chemical networks leads us to its possible hiding place. In W13, ~30% of the oxygen is in species which are not listed in Whittet (2010). It predicts s-CO₂, s-H₂CO, s-H₂O₂, O₂, and s-O₂ at abundances ≳10^-5, and NO, s-NO, s-HNO and s-HCOOH and more complex CHON species contributing X ≳ 10^-6. Instead of having a single, large oxygen reservoir, the UDO could actually be distributed amongst many molecules – both in icy and gaseous form.

6. Conclusions and summary

In this paper, Herschel-HIFI observations of water vapour and previously published observations of water ice column densities are used to understand the connection between the water gas and ice in the cold environment of protostellar envelopes.

• 1.

  We develop a simple chemistry network, SWaN. The water vapour and water ice abundances in the cold regions of pre- and protostellar cores can be reliably determined by only including freeze-out, photodesorption and photodissociation for water vapour, water ice, and atomic oxygen (as a proxy for other oxygen-bearing species). In the outer layers of protostellar envelopes (A_V ≲ 3 mag) the water abundance structure can be determined by only considering the cycle O → s - H₂O → H₂O → O. At higher extinctions (and also higher densities), where the freeze-out rather than photodissociation dominates the destruction of water vapour, the network reduces to H₂O ⇆ s - H₂O.

• 2.

  In cold, prestellar cores, the bulk of atomic oxygen is mostly converted into water ice on a timescale of 1−10 Myr. In our sample of protostellar cores we find that only 10−30% of the oxygen is found in the form of observed water ice. The key to model such a low abundance is a short prestellar core lifetime of ≲0.1 Myr together with the fact that the formation of water ice is inhibited in protostellar envelopes in regions with T_dust ≳ 15 K as a result of the assumed oxygen binding energy of T_b,O = 800 K. Use of a higher binding energy reinforces our conclusion of a short pre-collapse timescale. The apparent contradiction of such short pre-collapse phases to observed lifetimes of the dense phase of prestellar cores (~0.5 Myr) can be circumvented by introducing a mechanism that efficiently reduces the water ice abundance in the transition from a pre- to protostellar core. One hypothesis is that of episodic accretion and its effects on the abundance should be critically analysed.

• 3.

  We find infall motions of the envelope for our two Class 0 sources, and expansion motions in the majority of the Class I sources. This is consistent with infall during early stages, and envelope dispersal during later stages. In addition, the water vapour emission lines prove to be an excellent tracer for the FUV field strengths. We find CR-induced FUV field strengths of G_cr ≲ 10^-4 for all sources, and mostly low ISRF of G_isrf ~ 10^-2−10⁰ – amongst them three sources from the Taurus star-forming region. This attenuated ISRF supports the finding of previous authors of an ice formation threshold at A_V ~ 3 mag.

• 4.

  The finding of a low water ice abundance consistent with the observations sheds light on the question of the question of the unidentified depleted oxygen (UDO). Chemical modelling shows that upon introduction of an extended grains surface chemistry network, oxygen can be distributed over many molecules (mostly ices). Instead of a single, large oxygen reservoir, the UDO would consist of a plethora of various species, which makes it hard for observers to track down the complete oxygen budget in these sources.

---

Acknowledgments

M.S. would like to thank Ted Bergin, Eric Keto and Paola Caselli for a useful discussions and help with the development of SWaN, and Coryn Bailer-Jones for assistance with the statistical analysis. M.S. acknowledges support from NOVA, the Netherlands Research School for Astronomy. M.S. also acknowledges the use of astropy (Astropy Collaboration 2013), NumPy, SciPy, and matplotlib (Hunter 2007). This research has made use of NASA’s Astrophysics Data System Bibliographic Services (ADS). R.V. is supported by NASA through an award issued by JPL/Caltech and by the National Science Foundation under grant 1008800. C.W. acknowledges support from the European Union A-ERC grant 291141 CHEMPLAN and financial support (via a Veni award) from the Netherlands Organisation for Scientific Research (NWO). Astrochemistry in Leiden is supported by the Netherlands Research School for Astronomy (NOVA), by a Royal Netherlands Academy of Arts and Sciences (KNAW) professor prize, by a Spinoza grant and grant 614.001.008 from the Netherlands Organisation for Scientific Research (NWO), and by the European Community’s Seventh Framework Programme FP7/2007−2013 under grant agreement 238258 (LASSIE). HIFI has been designed and built by a consortium of institutes and university departments from across Europe, Canada and the United States under the leadership of SRON Netherlands Institute for Space Research, Groningen, The Netherlands and with major contributions from Germany, France and the US. Consortium members are: Canada: CSA, U. Waterloo; France: CESR, LAB, LERMA, IRAM; Germany: KOSMA, MPIfR, MPS; Ireland, NUI Maynooth; Italy: ASI, IFSI-INAF, Osservatorio Astrofisico di Arcetri-INAF; Netherlands: SRON, TUD; Poland: CAMK, CBK; Spain: Observatorio Astronómico Nacional (IGN), Centro de Astrobiología (CSIC-INTA). Sweden: Chalmers University of Technology − MC2, RSS & GARD; Onsala Space Observatory; Swedish National Space Board, Stockholm University − Stockholm Observatory; Switzerland: ETH Zurich, FHNW; USA: Caltech, JPL, NHSC.

References

Appendix A: Photodesorption and spherically averaged extinction

The photodesorption rate is closely connected to the rate at which FUV photons impinge on the grain surface. An isotropic ISRF is characterised by a flux of photons F_phot,0 through a unit surface per unit time, which can be directly related to the isotropic specific intensity I_phot,0 of FUV photons (A.1)Deeper into the envelope, the FUV field is attenuated according to exp(−γA_V), where γ is a pre-factor that depends on the dominant FUV range for the process in question. At an arbitrary distance r from the core centre, the specific intensity of FUV photons is isotropic, but rather depends on the azimuth angle θ (Fig. A.1). The FUV photon density can be calculated by In the last step, the directionally dependent extinction is replaced with the spherically averaged extinction (A.4)The rate at which FUV photons impinge on the surface of a single grain is given by (A.5)Since photodesorption of ice and photodissociation of water vapour are dominated by slightly different wavelength ranges, their respective γ-factors are 1.8 and 1.7. This does produce slightly different average extinctions. It should also be noted that differs substantially from A_V, which denotes the extinction from the envelope edge inwards in radial direction and is defined in Sect. 2.2. For example, at the envelope edge A_V = 0.0 mag we find a , at A_V = 1.0 mag we get , and at A_V = 10 mag we derive . These two values converge towards the center.

Appendix B: Chemical network benchmarking

Appendix B.1: Other chemical networks

Section 3.1 introduced our simplified water network SWaN. To assess its reliability, it is compared to three full chemical networks with different databases, reaction channels, reaction mechanisms, and computational details. The protostellar core NGC 1333-IRAS4A is selected as a benchmark object, since this source has been tested in detail against observations (Mottram et al. 2013).

• The network V11 is based on Visser et al. (2011), but fully updated to the Rate12 release of the UMIST Database for Astrochemistry (UDfA⁴; McElroy et al. 2013). Its grain-surface chemistry is limited to the formation of H₂, H₂O, NH₃, CH₄, H₂S; the latter four are formed simply through successive hydrogenation of the respective heavy elements.

• The network W13 is based on Walsh et al. (2013) with additions from Walsh et al. (2014), which also employs the full Rate12 gas-phase chemistry. The network is supplemented with grain-surface reactions from the Ohio State University network (OSU; Garrod et al. 2008), which is much more extensive than the simple surface chemistry in V11. The grain-surface network in W13 forms both “simple” and “complex” ices, such as H₂O, CH₄, CH₃OH, HCOOCH₃, and CH₃OCH₃. The ice chemistry is treated as a single phase, i.e., the bulk and surface ice are not treated separately. Grain-surface reactions occur via the Langmuir-Hinshelwood mechanism only and reactive (or chemical) desorption is included with branching ratio of 1% (Garrod et al. 2007).

• The network A13 is based on the deuterium chemistry model of Albertsson et al. (2013), which has been extended to include ortho-para chemistry (Albertsson et al. 2014a) and high-temperature reactions (Albertsson et al. 2014b). It emerged from the protoplanetary disk model of Semenov et al. (2010), who adapted the gas-grain model of Garrod et al. (2008). A13 has a similar grain-surface network as W13.

All three networks contain several hundred species and several thousand reactions. In addition to standard gas-phase chemistry, they all allow freeze-out of neutral molecules onto cold dust grains. Desorption can occur thermally, through direct cosmic ray heating, or through absorption of UV photons. Also included in all three networks are photo-ionisation and photodissociation, along with grain-surface formation of H₂.

All networks use the same set of parameters listed in Table 2, together with a primary cosmic-ray ionisation rate of ζ_cr = 5.0 × 10^-17 s^-1. Desorption of ices is treated as a zeroth-order process from only the top two monolayers (Sect. 3). The abundances relative to hydrogen nuclei X_i ≡ n_i/n_H (Table B.1) are the same in all models, and result from a dark cloud model (T = 10 K, n_H = 2 × 10⁴ cm^-3, A_V = 10 mag) at t_pre = 0.1 Myr (Walsh et al. 2013) using the Rate06 coefficients with the OSU grain surface network. Table B.1 highlights in bold the three species included in SWaN: water vapour, water ice and atomic oxygen. Their cumulative abundance, defined as X_O,SWaN ≡ X_H2O + X_{s - H2O} + X_O, is initially 212 ppm⁵ out of the available X_O,ISM = 320 ppm (Meyer et al. 1998). Other oxygen-bearing species like carbon monoxide or methanol are not considered, which limits the amount of oxygen in SWaN to ~66% of the abundance of volatile oxygen in the ISM.

The abundance profiles for water vapour, water ice, and atomic oxygen (plus the sum of these three species) at t_post = 0.1 Myr are presented in Fig. B.2. At all radii, the three full networks agree to within a factor of 2 on the water ice abundance. V11 and W13 show equally good agreement for water vapour, but A13 differs by up to two orders of magnitude. This originates in the different treatment of the desorption: In contrast to V11 and W13, where only molecules from the top layers can desorb, A13 allows desorption from the full ice mantle. SWaN recovers the basic abundance structure for water ice and vapour to within the uncertainties from the full networks, except for an underestimate of the ice abundance at the very outer edge of the envelope.

The situation for atomic oxygen is more complicated, and the predicted abundance structure of SWaN is less reliable. All networks agree on a high O abundance at low A_V due to photodesorption and photodissociation of water and other oxygen-bearing species, but in the shielded regions with higher densities, the oxygen abundance in SWaN strongly deviates from the other networks. The reason for this lies in the fact that the sophisticated chemical networks include certain channels that allow atomic oxygen to react further to form other species like O₂, CO, which are not included in SWaN. The bottom panel of Fig. B.2 shows the cumulative abundance of all species present in SWaN. At high extinctions (and high densities), the complex networks feature a net flow of oxygen into species that are not part of SWaN. Hence, the abundance profile of atomic oxygen in SWaN does not reflect the true profile in protostellar envelopes, but rather acts as a oxygen reservoir that represents all other oxygen-bearing species that are present in in the environment. Overall, considering the limitations of the network, the agreement of both the water vapour and water ice abundances with the complex chemical networks is excellent.

Appendix B.2: Implications for observations

The modelled abundance profiles can be turned into observable column densities by integrating over radius. Comparison to the total hydrogen column density then gives the average line-of-sight column density ratios, or “abundances”, which can be compared to observations. As discussed in Sect. 2.2, our analysis is limited to the water freeze-out zone, i.e., radii where T_dust ≲ 100 K. Practically all water ice is found in this region, which makes this point a convenient choice – in particular when comparing to observations of water ice.

In Table B.2 the column density ratios are summarised for all three species in the water freeze-out zone for the benchmark NGC 1333 IRAS4A model. The agreement between the chemical networks is good for water vapour (all are within a factor of 2.5) and water ice (factor of 1.5). This results in column-density-averaged gas-to-ice ratios of (3−11) × 10^-5. On the other hand, atomic oxygen shows large discrepancies of more than an order of magnitude.

In Fig. 9 in the main text, the column density ratios for all chemical networks alongside the initial abundances from Table B.1 are shown. This figure illustrates how the oxygen budget evolves in the protostellar phase. We distinguish between the chemical species present in SWaN (H₂O, s - H₂O, O), CO gas, carbon-oxygen-ices (s - CO, s - CH₃OH, s - H₂CO, which are the most abundant molecules from that group after the pre-collapse phase), and other oxygen bearing species. Water vapour is only present in trace amounts.

All chemical networks have only a small fraction of the oxygen in water ice, about 50−100 ppm, consistent with observations (e.g. Pontoppidan et al. 2004; Whittet et al. 2007). One of the key findings of Fig. 9 is the apparent drain of oxygen in the sophisticated chemical networks towards species that are not part of SWaN during the denser and warmer protostellar phase. Obviously, in the simplified chemical network the abundance of oxygen has to be conserved, since no other molecules are part of the network, and atomic oxygen acts as the oxygen reservoir. However, in particular the chemical networks A13 and W13 actively transform oxygen into other species. The abundances of the CO gas and carbon-ices increase marginally, but the atomic oxygen abundance drops from its initial value of ~157 ppm to 5−15 ppm. Other oxygen bearing molecules , which in the beginning were only present in trace amounts, make up ~106 ppm after only 10⁵ yr.

Appendix C: Bayesian analysis of the water emission profiles

The likelihood that the observed data D (a list of N data points {y_k}) can be explained by model M with input parameters θ is, assuming a Gaussian error distribution, given by with Y_k is a synthetic spectrum, which is the outcome of model M; in our case a radiative transfer model with an abundance and velocity profile determined through the parameters . The noise term σ has contributions from the noise of the observations σ_obs (typically 10−30 mK), and an uncertainty in the synthetic RATRAN spectra of ~40 mK. Different weighting schemes (e.g., increasing uncertainties in the RATRAN spectra, but also including velocity dependent terms) showed no influence on the general results. Each grid point in the parameter space is assigned the minimum χ² that is reached when shifting the spectrum by ± 0.4 km s^-1 around the individual systemic velocity, v_LSR (van Dishoeck et al. 2011; Yıldız et al. 2013a) to account for uncertainties in the estimate of v_LSR.

In our analysis (Sect. 4.1) we aim to find general trends rather than a singular best fit point in our parameter grid ϕ with N_ϕ grid points. To achieve this we follow a Bayesian approach (Bailer-Jones 2011), where the overall fit quality is assessed by the so-called evidence, which is in the case of a uniform prior distribution (i.e., all grid points are treated equally) given by (C.1)The evidence is equivalent to the probability of the data D given the model M. The dependence on parameter θ has vanished due to marginalisation of the full parameter grid.

In our analysis, however, we are interested in variations of the fit quality within the model upon variation of the value for particular parameter P. Therefore, we modify our definition of the evidence in Eq. (C.1). We define a new parameter space φ ≡ ϕ_{P = Vi} with N_φ grid points, where parameter P takes value V_i. All other parameters remain free. We then calculate the evidence by marginalising all free parameters, which is given by (C.2)These evidences give an overall estimate of the fit quality, but their absolute values are irrelevant. They only gain relevance when comparing to each other. We therefore define the Bayes Factor for each value V_i of parameter P as (C.3)The parameter value with maximum likelihood has, therefore, a Bayes Factor of unity. Other parameter values can also exhibit high likelihoods, but a value is assumed to yield a poor representation if its Bayes Factor is ≪0.1.

This approach is illustrated by the following example: We analyse the influence of the post-collapse timescale t_post ≡ p with two values V = {0.1 Myr,1.0 Myr}. The full parameter space is split up into two sub-spaces with . Marginalisation of all parameters except t_post (Equation C.2) leads to evidences E_p(0.1 Myr) and E_p(1.0 Myr). Let us now assume, that maximum evidence is found for a post-collapse time of 0.1 Myr. We then get B_p(0.1 Myr) = 1 and B_p(1.0 Myr) < 1, i.e., the most likely post-collapse time is 0.1 Myr. If we find B_p(1.0 Myr) ≪ 0.1, we can even assume that this is the only reasonable solution, whereas we can reject t_post = 1.0 Myr.

The summary for all parameters and all sources is shown in Table 5.

Appendix D: Observation

The observating dates and IDs for each transition and source are listed in Table D.1.

All Tables

All Figures