© ESO, 2013

1. Introduction

For stars to form from the interstellar medium the density must increase from ~10 cm^-3 to ~10²⁴ cm^-3. While some of this process is achieved through the formation of molecular clouds (n ~ 10⁴ cm^-3) and then dense filaments and cores (n ~ 10⁶ cm^-3) within these clouds, the remaining 18 orders of magnitude increase must be accomplished through infall of gravitationally bound and unstable envelope material onto the central protostar. This infall may proceed on an inside-out (Shu 1977) or outside in basis (Foster & Chevalier 1993). When angular momentum is taken into account (e.g. Ulrich 1976; Cassen & Moosman 1981; Terebey et al. 1984), this leads to a flattened inner structure which may eventually form a rotationally supported disk, with infall proceeding from the envelope to the disk, then through the disk onto the star.

Observational attempts to test such models and understand how infall varies both spatially and in time have proven more difficult. The observational signature associated with infalling gas is an asymmetric line profile with more blue than redshifted emission (see e.g. Evans 1999; Myers et al. 2000). In some tracers, strong absorption by the outer parts of the envelope results in an inverse P-Cygni line profile where the absorption is red-shifted with respect to the source velocity, usually accompanied by blue-shifted emission. Various authors have used dense gas tracers such as CS, HCO⁺, HCN, N₂H⁺ and H₂CO to probe the infall signatures of cores and protostars using both single-dish and interferometric observations (e.g. Zhou 1992; Zhou et al. 1993; Gregersen et al. 1997; Mardones et al. 1997; Di Francesco et al. 2001; Fuller et al. 2005; Attard et al. 2009). However counter-claims that some profiles can be reproduced with a combination of rotation and foreground layers (Menten et al. 1987; Choi 2001) have made it difficult to conclusively identify infall without the need for sophisticated modelling (e.g. Brinch et al. 2009).

In order to trace infall, one needs a good tracer of motion in envelope material. The molecules previously used were chosen because they have critical densities and abundances which make them reasonable probes of the regions where infall is expected to take place. Water, on the other hand, is particularly sensitive to motion: even though sub-millimetre water transitions have higher critical densities of order 10⁶ − 10⁸ cm^-3, they are easily sub-thermally excited due to large Einstein-A coefficients. In addition, the freeze-out of water onto dust grains at temperatures below ~100 K and the photodesorption or photodissociation of water by both interstellar and cosmic-ray induced UV radiation make the chemistry of water particularly sensitive to temperature and the local radiation field (e.g. Hollenbach et al. 2009). Indeed, the details of water line profile shapes actually allow this chemistry to be constrained on much smaller spatial scales than the observing beam (e.g. Caselli et al. 2012).

The properties that make H₂O such a good tracer of motion also make it impossible to observe the vast majority of transitions from the ground due to atmospheric absorption, so space satellites are required to study water comprehensively. The “Water in star-forming regions with Herschel” (WISH; van Dishoeck et al. 2011) survey has obtained velocity-resolved spectra of multiple water lines towards a sample of ~80 targets which cover the luminosity range from 1 to 10⁶L_⊙ and vary in evolutionary stage from pre-stellar cores to gas-rich disks. While most of the water emission detected towards the sub-sample of 29 low-mass Class 0/I protostars is related to outflow motions (e.g. Bergin et al. 2003; Kristensen et al. 2010, 2012), seven sources were identified by Kristensen et al. (2012) which show inverse P-Cygni profiles in the ground-state H₂O 1₁₀ − 1₀₁ line at 557 GHz.

The analysis method used by many previous authors, including Kristensen et al. (2012), assumed two infalling slabs of material, sometimes with the addition of a third central static slab in order to reproduce the “true” inverse P-Cygni profile. While such an approach is a good first step, it only provides information about the characteristic infall velocity at a characteristic radius from the central protostar. The first goal of this paper is to constrain the spatial scales on which infall takes place within these seven WISH sources. This will allow accurate measurement of the mass infall rate (Ṁ_inf) of material from the envelope towards radii where a flattened disk-like structure may be present. This is distinct from the mass accretion rate (Ṁ_acc) of material from any central disk-like structure onto the central protostar. The mass infall and accretion rates need not be directly related, particularly if the disk stability varies as a function of time. Understanding what the mass infall rate is, and over what spatial scales infall takes place, will allow us to predict whether a central flattened disk-like structure is stable during the early phases of star formation. The second goal of this paper is to use the intensity and shape of gas-phase water lines to constrain water chemistry in the cold regions of protostellar envelopes, a topic more extensively discussed in Schmalzl et al. (in prep.).

Both these goals require more detailed and complex modelling than previously employed. We therefore perform full 1–D spherically symmetric non-LTE radiative transfer modelling of water line profiles, taking into account an envelope with varying density, temperature, velocity and water abundance. The paper is structured as follows. In Sect. 2, we present the observations with which we will constrain our models, as well as removal of outflow-related emission using Gaussian fitting. The modelling approach, including the simple water chemistry required to calculate how the water abundance varies within the envelope, is discussed in Sect. 3. The results and analysis relating to infall are presented in Sect. 4, while those relating to water chemistry in protostellar envelopes are presented in Sect. 5. We then discuss the wider implications of these results in Sect. 6 before reaching our conclusions in Sect. 7.

2. Observations

The seven sources discussed in this paper are part of the WISH sub-sample of embedded low-mass Class 0/I protostars and were selected because they have inverse P-Cygni line profiles in the H₂O 1₁₀ − 1₀₁ line (E_up/k_b = 61 K), as identified by Kristensen et al. (2012). Their general properties are given in Table 1; six are Class 0 and one (GSS30-IRS1) is a Class I protostar. All sources were observed in two ortho-H₂O and three para-H₂O lines with the Heterodyne Instrument for the Far-Infrared (HIFI; de Graauw et al. 2010) on the Herschel Space Observatory (Pilbratt et al. 2010) as part of the WISH survey. These lines have upper-level energies ranging from 50 to 250 K. In addition, the 3₁₂ − 3₀₃ ortho-H₂O line was observed for two sources (NGC 1333-IRAS4A, hereafter referred to as IRAS4A, and Ser-SMM4) and the 2₁₂ − 1₀₁ ortho-H₂O line observed towards IRAS4A only. The observation identification numbers are given in Table A.1, which also serves to summarise which lines were observed towards each source. All data were obtained using both the wide-band spectrometer (WBS) and high-resolution spectrometer (HRS) in both horizontal and vertical polarisation. Table 2 presents details of the observed lines, including the line frequency, main-beam efficiency, spatial and spectral resolutions and the upper level energy of the transition.

The data were reduced using hipe (Ott 2010) version 8.2 with all further processing performed using python. A low-order (i.e. ≤2) polynomial was used to fit and subtract the baseline from the two WBS polarisations for each line separately, after which the data were co-added. HIFI is a dual-sideband receiver, so each spectrum is a combination of simultaneous observations in both the upper and lower sidebands. As such, where a constant baseline can be fit to a sub-band this gives twice the continuum brightness temperature as it includes continuum photons from both sidebands. Some emission lines are broader than the more limited velocity range covered by the HRS observations (typically 50−100 km s^-1), so only a constant baseline was used for these data, but otherwise the process was the same. All co-added data were then converted to the main-beam temperature scale using the beam efficiencies from Roelfsema et al. (2012). A more detailed discussion of the data processing and quality for all WISH low-mass protostars can be found in Kristensen et al. (2012) for the H₂O 1₁₀ − 1₀₁ line, and will be presented in Mottram et al. (in prep.) for all other water observations of low-mass YSOs.

Due to the higher spectral resolution available with the HRS, and the fact that envelope emission is always within ν_LSR   ±   8 km s^-1, only these data will be used in the remainder of this paper with the exception of the 2₁₂ − 1₀₁ transition at 1670 GHz. For this high-frequency line the WBS spectral resolution is comparable to the HRS resolution in the other lines and so the improved noise of the WBS over the HRS (a factor of $\sqrt{\mathrm{2}}$) makes these data a better choice. All spectra shown in this paper have been shifted so that the systemic velocity of the source is at 0 km s^-1.

Figure 1 provides an example of the WISH H₂O observations for IRAS4A. The H₂O emission in many of the lines is complex, with multiple components within the beam. As discussed in Kristensen et al. (2010) for IRAS4A, the broad (FWHM ≳ 20 km s^-1) Gaussian component is associated with emission from a molecular outflow while additional components which are offset from the source velocity and have more moderate linewidths (FWHM ~ 5 − 10 km s^-1) likely originate from currently shocked gas (Kristensen et al. 2013). These components are distinct from the narrow envelope emission near the source velocity, and result from very different physical processes which are not directly included in our model.

In order to isolate only the emission relating to the protostellar envelope, Gaussian fits to the line profiles are used to remove the broader and shifted components (see also Kristensen et al. 2012). A combination of up to three Gaussian profiles is required to remove the additional distinct physical components. We assume that the line width and central velocity of each component is the same for all lines of a given source. Therefore, after masking emission and absorption associated with the inverse P-Cygni profile, all lines and components are fitted simultaneously with common central velocities and linewidths for each component. An example fit is shown in Fig. 1, with the residuals for two lines presented in Fig. 2. All three components identified in this source are observed in multiple lines, though each component is not necessarily detected in all lines. In particular, for IRAS4A the fainter component near the source velocity, which was not included in the fit to the 1₁₀ − 1₀₁ line by Kristensen et al. (2012), is required to fit the higher-excited lines. Removing this component leads to a demonstrably worse fit to the data. This difference in the data and components considered leads to a slight difference in our fit results compared to those presented in Kristensen et al. (2012) but are consistent.

Inverse P-Cygni profiles are observed in the ground-state 1₁₀ −1₀₁ and 1₁₁ − 0₀₀ lines for all sources, with only IRAS4A showing such a profile in the excited 2₀₂ − 1₁₁ line. There are no detections of narrow envelope emission or absorption in the 2₁₁ − 2₀₂, 3₁₂ − 2₂₁ or 3₁₂ − 3₀₃ lines unlike for high-mass protostars (van der Tak et al. 2013). While some Gaussian emission is sometimes observed near the source velocity, this is too broad (FWHM ~ 10 km s^-1) to be consistent with envelope emission. In many sources the absorption in the ground-state water lines due to the envelope is saturated, i.e. the absorption extends below the continuum and is flattened (e.g. see Fig. 1) indicating that all photons, both from the continuum and the outflow component, have been absorbed by the envelope. The difference between the baseline and the saturated absorption then gives a measure of the continuum level at this frequency.

3. Models

In order to constrain and understand the velocity field within the sample of sources discussed above we simulate the water line profiles using 1-dimensional spherically symmetric ratran models (Hogerheijde & van der Tak 2000). ratran uses a Monte-Carlo method with accelerated lambda iteration, solving for the level populations including radiative pumping by dust continuum, then performing ray-tracing to create synthetic images. Finally, the synthetic images are convolved to the appropriate Herschel beam and the spectrum towards the centre of the model is extracted.

The details and assumptions of these models are discussed in the following subsections, but first it is important to discuss the impact of assuming spherical symmetry. These sources are all known to have prominent outflows which we do not include in the model grid. While we have removed components associated with the outflow from the observations to account for this, the outflow cavity will still lead to the removal of dense material from some sight-lines in the envelope. In addition, the cores around these sources may be flattened, particularly if they reside in filaments (e.g. L1157; Tobin et al. 2012a). Both these effects could lead to 2/3–D variations in the density, and thus temperature, compared to the spherically averaged case. However, using 2–D or 3–D models would require many more free geometrical parameters which the data presented here cannot constrain due to the lack of spatial information, so such models would only add degeneracy. In addition, such models are more computationally expensive, particularly if non-LTE line radiative transfer is included, which for H₂O is important. A 1–D non-LTE model is therefore the best approach with the data in hand.

On small (≲500 AU) scales, the distribution of material around a protostar is unlikely to be spherical due to rotational flattening, and so the 1–D approximation breaks down (Jørgensen et al. 2005). While we include these regions in our models, most of the observed emission comes from regions well outside such radii because this is inside the τ = 1 surface for central velocities and beam dilution minimises the contribution of emission at velocities where the optical depth is lower. Nevertheless, care must be taken not to extrapolate these results to conditions in the inner envelope where 2/3–D modelling is necessary and where the physical conditions are less constrained.

3.1. Density and temperature

The density is assumed to follow a power-law (i.e. n    ∝    r^− p). The exponent (p), dust temperature profile and inner and outer radii were determined by Kristensen et al. (2012) using dusty 1–D continuum radiative transfer models (Ivezić & Elitzur 1997) to fit the spectral energy distributions (SEDs) from 70 μm to 1.3 mm, scaled to the source bolometric luminosity, and the radial intensity profiles from SCUBA 450 μm and 850 μm images. For more details, see Appendix C of Kristensen et al. (2012). We do not fit the wavelength dependance at shorter wavelengths as this primarily probes the small-scale geometry, which would require a 2/3–D continuum model to reproduce correctly. Such a model would have additional free parameters that would be poorly constrained at best, thus adding degeneracy. The DUSTY results used are given in Table 3. The gas temperature is assumed to be the same as the dust temperature and interpolated onto the cell grid (see Sect. 3.3). The ortho and para H₂ densities are calculated from the density power law assuming a thermal ortho-to-para ratio, but with a lower limit of 10^-3 (Pagani et al. 2009). The OH5 dust opacities for grains with thin ice mantles from Ossenkopf & Henning (1994) are used when calculating the continuum for both the dust and line radiative transfer. An example of the density and temperature profile for IRAS4A is shown in Fig. 3, along with the radius of the continuum τ = 1 surface as a function of frequency. Continuum optical depth is not an issue for our observations, as the line optical depth is in all cases much higher than the continuum.

3.2. H₂O abundance structure

The models also require a profile for how the water abundance with respect to H₂ varies with radius. As will be shown in Sect. 5, water line profiles are particularly sensitive to the shape of this profile. Previous water modelling studies have used parameterised water abundance profiles which either have two or three constant abundances within different zones defined based on temperature and density limits (e.g. van Kempen et al. 2008; Coutens et al. 2012; Herpin et al. 2012), similar to those often used for other simple molecules such as CO (e.g. Jørgensen et al. 2002; Yıldız et al. 2012). However, it was impossible to convincingly reproduce the observed inverse P-Cygni profiles using such a scheme. This will be discussed further in Sect. 5.

The SWaN (Simplified WAter Network) chemical model of Schmalzl et al. (in prep.), based on the chemical description used by Hollenbach et al. (2009), is therefore used to calculate the water abundance profile as a function of the density, temperature and extinction in the envelope. This takes into account the formation of water on dust grain surfaces through the freeze-out of atomic oxygen and subsequent hydrogenation, the freeze-out of water vapour onto dust grains, the thermal and photodesorption of water ice into the gas phase and the destruction of water vapour through photodissociation. All equations are solved in a time-dependent manner until steady-state is achieved (usually <10⁵ yr). The approach used in SWaN is similar to that used in Caselli et al. (2012), which followed the earlier work on CO of Keto & Caselli (2008). Use of SWaN naturally leads to a photodesorption layer near the outer edge of the model where the A_V is low enough for the interstellar radiation field to penetrate and lower densities lead to less efficient freeze-out, rather than having to add an external (absorbing) layer as implemented by Coutens et al. (2012). The abundance profiles produced by SWaN are consistent with those produced by full chemical networks, unlike the drop abundance profile (Schmalzl et al., in prep.) The best-fit abundance profile for IRAS4A is shown in Fig. 4. Separate profiles are calculated for each source based on the source physical structure.

The main free parameters which can alter the water abundance profile are the interstellar UV radiation field (ISRF, G_isrf) and the cosmic-ray induced UV field (G_cr). Both are scaling factors relative to the flux of interstellar photons, which we take to be F_isrf = 10⁸ photons cm^-2 s^-1 (from an energy flux of 1.6 × 10^-3 erg cm^-2 s^-1 and an energy range of 6 − 13.6 eV; Hollenbach et al. 2009). All volatile oxygen (i.e. all oxygen not locked up in refractory grains) is assumed to be in the form of water, leading to a total water abundance relative to H₂ of 6 × 10^-4 (from O/H = 3 × 10^-4; Meyer et al. 1998). This does not take into account the fraction of oxygen in CO, and is therefore an upper limit. However, the abundance of gas-phase water is relatively insensitive to this quantity when T ≲ 100 K, the primary region probed by our observations. We also assume a photo-desorption yield of 10^-3 (Öberg et al. 2009; Arasa et al. 2010) and grains with a radius of 0.1 μm (Bergin et al. 1995). A more detailed discussion and comparison with more complex chemical networks will be presented in Schmalzl et al. (in prep.).

The water abundance profile and water lines are sensitive to the density at the outer edge, particularly in the 1₁₀ − 1₀₁ line which shows significant emission at central velocities from the photodissociation layer if the density at the outer edge is above ~10⁵ cm^-3. As this is not seen in the observations, we set G_isrf to 0 if the outer density is ≥10⁵ cm^-3 (see Table 3), as it will almost certainly be attenuated by cloud material at lower densities. In all other cases we set G_isrf to 1.

A constant ortho-to-para ratio for water of 3 is also assumed. A lower ortho-to-para ratio (e.g. 1) would lead to an increase (decrease) in the para (ortho) water abundance and thus an increase (decrease) in the intensity of emission and absorption features of all para (ortho) lines. However, most measurements are consistent with a value of 3 (e.g. Emprechtinger et al. 2013). When performing the radiative transfer calculations, we use the latest collisional rate coefficients from Daniel et al. (2011) as downloaded from the Leiden Atomic and Molecular Database (LAMDA; Schöier et al. 2005).

3.3. Grid

As with all grid-based simulations, it is important that the grid of cells defined in ratran resolves the key regions of the envelope. To do this, we begin by defining a grid of 23 cells which are logarithmically spaced in density. Due to the abundance structure of the outer part of the model being so important to the lines we are interested in and because the abundance profile can change significantly over small changes in radii in that region, the outer two cells are then replaced with six cells which have a power-law spacing i.e. the number of cells increases towards the outer edge. Finally, the two cells which cover where the temperature crosses 100 K are replaced by four cells to better resolve this jump (see Fig. 3).

The inner radius for the grid is taken from the dusty model. This parameter is not particularly well constrained, but because it is inside the τ = 1 surface for all H₂O lines, it does not affect our results. The outer radius of the model is restricted to either the radius at which the temperature drops below 8 K or the density drops below 10⁴ cm^-3, whichever is the smaller. This ensures that the model envelopes do not extend to unphysically large radii. Material from the cloud outside this radius along the line of sight has a negligible impact on the line profiles unless a significant foreground cloud layer is present, in which case this layer is not added to the grid but is included later during the creation of synthetic images (see Sect. 3.5 for more details).

For Ser-SMM4 and L1157 the full DUSTY model, the SCUBA 850 μm images and low-J CO maps Yildizip et al. (in prep.) all suggest that the envelope extends out to radii ≳30′′ but the cuts discussed above would lead to a smaller outer radius. In these cases we therefore set the outer radius of the model to be 30′′ and limit the temperature to be no lower than 8 K, which corresponds to thermal balance between CO line-cooling and cosmic-ray heating (Evans et al. 2001; Galli et al. 2002; Bergin et al. 2006), transferred to the dust via gas-grain collisions.

In the case of IRAS4A the model is large enough that the outer layers of the model would encompass IRAS4B and almost touch IRAS4C if these sources are all at a common distance. However, the largest Herschel beam is smaller than the model and the separation between the sources, so the properties of the model at radii much larger than the beam size only contribute to the model spectra at the front and back of the model along the line of sight, as is also the case for the observed spectra.

At densities below ~10⁵ cm^-3 the dust temperature can rise above the gas temperature (Doty & Neufeld 1997), contrary to our assumption that they are the same. External heating in the outer envelope (e.g. Jørgensen et al. 2006) can also have a similar effect which was not considered in the dusty models. However, any underestimate in the dust continuum caused by our possible underestimate of the dust temperature at the outer edge of the model will have a negligible effect on the total continuum flux, which is dominated by the inner envelope.

3.4. Velocity and non-thermal motion profiles

While ratran takes thermal broadening into account automatically, two other contributions to the line-width must be defined – systematic motions (i.e. infall or expansion) and non-thermal turbulent motions. Due to the observed inverse P-Cygni profiles in the observations, the material within the model is assumed to be in free-fall towards the protostar at the centre of the envelope/disk. The infall velocity (ν_inf) is scaled to the velocity at 1000 AU such that: $\begin{array}{ccc}{\nu }_{\inf }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\nu }_{\mathrm{1000}}{\left(\frac{\mathit{r}}{\mathrm{1000}\hspace{0.17em}\mathrm{AU}}\right)}^{-0.5}\mathrm{·}& & \end{array}$(1)A completely static envelope leads to a symmetric line profile, while increasing infall velocity leads to increasingly red-shifted absorption and blue-shifted emission. Whether this infall is throughout the whole envelope, or just on certain radii will be explored in Sect. 4.2. Rotation is not considered, but would take the form of additional line broadening which increased with decreasing radius. There is little sign of this in the observations and it would impossible to distinguish from turbulent broadening using 1–D models.

In order to take account of turbulence within the envelope, we also define the doppler b within the model, which is the 1/e half-width (0.601 × FWHM). Increasing b broadens both the absorption and emission components of the line profile while decreasing the peak intensity. We set b to be constant as a function of radius, though in some cases variation with radius may play a secondary role (see Sect. 4.1).

3.5. Additional emission/absorption components

For absorption to take place, a significant population of molecules must be in the lower level of a transition and a background emission source must provide photons of the correct frequency. In protostellar envelopes the dust continuum emission usually provides these photons, but absorption of line photons from deeper into the envelope or from other emission mechanisms is also possible. In addition, any material along the line of sight such as intervening clouds can also contribute to the absorption, though these clouds must be relatively cold and have low densities as they are not seen in emission and so only contribute to the ground-state absorption.

Figure 5 shows a cartoon of this process. In a simple infalling envelope the emission from the far side of the envelope is blue-shifted while the absorption and emission from the near side of the envelope is red shifted. If the absorption is deep enough, all line and continuum photons are absorbed and so the absorption is saturated. Absorption against the continuum emission is included automatically by ratran (cf. case 1 in Fig. 5).

For foreground clouds, the level populations of the emitting (i.e. envelope) and absorbing (i.e. foreground cloud) material are physically unrelated and so the absorbing layer can be added as an additional source of opacity along the line of sight during the production of synthetic images (cf. case 2 in Fig. 5). In ratran, the required variables are the velocity offset and line-width common to all lines along with the brightness temperature and line-centre optical depth of this component for each transition.

As noted in Sect. 2, the absorption in the ground-state lines is saturated at an intensity below the continuum and so must be absorbing both the continuum and the broad outflow component (cf. case 3 in Fig. 5). We do not know the detailed properties of the outflow material on sub-beam scales, so we do not include the outflow in the physical model. Instead, in order to reproduce the observed line profiles in the models we use a similar approach as for foreground layers, by including the observed emission from the outflow after the excitation calculation but during the calculation of model images, with the difference that the emission is added as each ray passes from the back half to the front half of the envelope. This is illustrated in the left-hand cartoon in case 3 of Fig. 5.

During the synthetic image creation, i.e. after the level population calculation, each ray starts at the back (left-hand side) and passes horizontally to the front (right-hand side). The intensity in each velocity channel is updated as the ray passes through each cell. To include the outflow, we add a broad Gaussian velocity component to the intensity as the ray passes the mid-plane of the model. The emission is then absorbed at some velocities by the front-half of the envelope, in a similar way to the continuum. To create this broad outflow component, we use the line properties (i.e. line-width, brightness temperatures etc.) obtained from the Gaussian fit to the data, assuming that the optical depth in this component is ≪1. This additional emission is then attenuated by the front half of the envelope. We subsequently remove the broad Gaussian emission from after ray-tracing but before the image is convolved to the observing beam, leaving only the absorption.

The main purpose of adding such a simple “outflow” component is to create a broad emission feature against which absorption by the front half of the envelope can occur. By subsequently subtracting the same broad emission feature so only the absorption remains, we ensure that the details of the broad emission feature are unimportant as long as the overall level of emission was correct so that the opacity of the front half of the envelope results in the right absorption depth in terms of antenna temperature. Including the outflow emission in this way assumes that the emission fills the beam and is unrelated to the excitation of the envelope. While this implementation is artificial and both of these assumptions may not strictly be true, this is probably a good approximation without including a full outflow – something which is a field of study in its own right.

While we assume that the outflow layer has a low optical depth, if it instead has a high optical depth then no emission from the back half of the envelope would be visible. In this case we would see all the red-shifted absorption but very little of the blue-shifted emission. Our assumption should not change the best-fit infall velocities, as these are derived primarily by concentrating on the absorption part of the profile. However, higher optical depth would hide emission and thus we would underestimate the abundance profile which is constrained primarily based on the emission. In reality, the outflow will have high-tau but only in some fraction of the beam while the rest of the emission is unhindered, resulting in something of a compromise situation.

4. Infall

ratran models as described above were run simultaneously for ortho and para H₂O for a range of values of ν₁₀₀₀, b and G_cr. Using a χ² statistic to find the best-fit model is not straightforward because of the variation in S/N between the lines. We therefore find the best fit by varying the different parameters and using by-eye comparison to decide which model is the most consistent with the data, within the uncertainties of the spectra. The properties of the best-fit model for all sources are presented in Table 3, along with approximate uncertainties obtained by changing each parameter individually until a noticeably worse fit was achieved. The columns give the density power law (p), inner and outer model radii and densities (r₀, r_out, n₀ and n_out), the infall velocity at 1000 AU (ν₁₀₀₀), the doppler b and G_cr. Note that though we extrapolate the velocity field to the inner radius of the model, in the case of Ser-SMM4 this results in unfeasibly large velocities. This will be discussed more in Sect. 4.3.

4.1. IRAS4A: a representative model

Let us now examine in detail the best-fit model for IRAS4A, and variations in parameter space around it, as this source shows an inverse P-Cygni profile in the excited 2₀₂ −1₁₁ line. This rules out foreground absorption as the explanation for the inverse P-Cygni profiles because absorption in this line in typical cloud conditions (T ~ 10 K, n ~ 10⁴ cm^-3, X_H2O ~ 10^-7) would require water column densities of order 8 × 10¹⁵ cm^-2, and therefore path-lengths greater than 2.5 pc, so there must be infall in the envelope. In addition to absorption against both the outflow and continuum (see Fig. 1), IRAS4A also shows absorption in the ground-state lines due to a low density foreground cloud which is redshifted by 0.8 km s^-1 with respect to the systemic velocity of IRAS4A (i.e. at ν_LSR = 8.0 km s^-1). This cloud component was identified in CO observations by Liseau et al. (1988) and has recently also been tentatively detected in O₂ emission by Yıldız et al. (2013a) with Herschel. They find a line FWHM of 1.3 km s^-1 and total column density of N(H₂) = 10²² cm^-2 with temperature in the range 20−70 K and density in the range 7 − 2 × 10³ cm^-3. Using the average density and temperature, we used the non-LTE slab-modelling code radex (van der Tak et al. 2007) to calculate optical depths and peak line intensities to be input into ratran as a foreground layer. We find that the total water abundance in this layer must be of order 10^-8 in order to reproduce the width of the observed absorption. Such an abundance is typical of gas-phase chemistry (e.g. McElroy et al. 2013). Figure 6 shows a comparison between the best-fit model (parameters given in Table 3) and the observations both with and without including the foreground and outflow absorption to illustrate the impact this has on the model line-profiles.

The best model fits the observed absorption and emission in the 2₀₂ −1₁₁ and 2₁₂ −1₀₁ lines very well, and reproduces the observed non-detections in the 2₁₁ −2₀₂, 3₁₂ −2₂₁ and 3₁₂ −3₀₃ lines. The absorption in the 1₁₁ −0₀₀ line is well reproduced though there is more emission in the model than in the observations. The emission in the model 1₁₀ −1₀₁ line is more intense and narrower than the observations by approximately a factor of two, while the absorption is still underproduced. The latter is most likely due to small differences between the strength of the continuum and outflow component in the models with respect to the observations. The absorption observed in the 2₀₂ − 1₁₁ is offset from the source velocity by 0.3 km s^-1, which is inconsistent with it being due to the foreground cloud at 0.8 km s^-1, and so the line profile can only be reproduced with an infalling envelope.

To explore what parameters the different lines are sensitive to, we show in Fig. 7 how the model line-profiles vary with different values for b and ν₁₀₀₀. For the remainder of the paper we only show lines which are sensitive to the parameter under consideration. The 1₁₁ −0₀₀ and 1₁₀ −1₀₁ lines are sensitive to b but not strongly sensitive to ν₁₀₀₀ while the 2₀₂ −1₁₁ and 2₁₂ −1₀₁ lines are sensitive to both parameters. While increasing b would improve the fit to the emission component of the 1₁₀ −1₀₁ line, this broadens the absorption in the other lines too much.

The sensitivity of each line to infall and/or turbulence can be understood if we consider which regions within the envelope these lines probe. The top panel of Fig. 8 shows the fraction of water molecules at a given radius in the first three energy levels of ortho and para water. The majority of water molecules in the outer envelope are in the ground-state, meaning that the excited lines are emitted in a much smaller region and so are subject to stronger beam dilution, explaining why no emission or absorption is detected in the 2₁₁ −2₀₂, 3₁₂ −2₂₁ and 3₁₂ −3₀₃ lines. The lower panel of Fig. 8 shows the cumulative integrated intensity of the four detected lines within a given radius in the model (I(r)) relative to the maximum absolute integrated intensity in each line (| I | _max). This is calculated for each radius by taking the level populations from the full model and re-running the ray-tracing and beam convolution with all cells outside that radius removed. The figure therefore gives a measure of how different radii within the model contribute to the emission and absorption for each line. The foreground layer was only included for the full model, hence the drop at the outer edge in the three affected lines. The emission from the 1₁₀ −1₀₁ line peaks around ~7000 AU, outside of which absorption dominates, while the emission for the other lines peaks between 700 and 2000 AU.

The 1₁₀ − 1₀₁ line is therefore much less sensitive to ν₁₀₀₀ because at the radii it probes, ν≲ b. The best-fit model does not reproduce the line-width of the emission from the 1₁₀ − 1₀₁ line well but the b = 0.6 km s^-1 model has a broader emission component. A higher turbulent line-width is therefore required in the region where the line emission comes from. This zone of higher turbulence cannot extend all the way to the outer edge as this would still lead to much broader absorption than is observed, which comes from the very outer layers. It also cannot extend to the inner edge of the model because the other lines would then be broader (cf. Fig. 7).

Figure 9 shows a model with a jump in b from 0.2 to 0.6 km s^-1 between 5000 and 8000 AU, confirming this reasoning. This also solves the overproduction of emission in the 1₁₀ − 1₀₁ and 1₁₁ − 0₀₀ lines. A decrease in turbulence to smaller envelope radii might be expected due to the increasing density allowing turbulence to be damped more efficiently, and this has been observed before towards the high-mass YSO W43-MM1 (Herpin et al. 2012) and in the B5 region of the Perseus molecular cloud (Pineda et al. 2010). However, it is unclear what physical mechanism would cause such a region of higher turbulence inside the envelope without causing a shock, which would cause variations in the infall velocity profile which we do not see in this region (see Sect. 4.2). The most likely candidate is that this is caused by the outflow driving turbulence in the wall of the outflow cavity, which is at ~45–60° to the line of sight (Yıldız et al. 2012). However, a detailed investigation would require 2/3–D modelling including a physical model for the outflow which is beyond the scope of this paper.

For consistency of analysis and because this step model results in a worse fit to the other lines we do not consider this model in the following sections, reverting to the best-fit model with a constant b. Since we are primarily interested in the infall velocity throughout the envelope, and the 1₁₀ − 1₀₁ is primarily sensitive to motion on large scales due to a combination of beam dilution and optical depth effects, this will not impact that analysis.

As discussed in Sect. 3, our physical structure assumes spherical symmetry but some of these sources have elongated or flattened envelopes which could lead to a lower or higher density power-law index (p) than is in fact the case. We show the difference caused by changing p to 1.5 and 2.0, as opposed to the value of 1.8 determined by Kristensen et al. (2012), in Fig. 10. The inner radius and density were recalculated for these models but the optical depth at 100 μm and ratio of inner to outer radius were kept the same as in the best-fit case. As can be seen, while the change in density profile does indeed lead to changes in the model water line profiles, these are similar to or smaller than the changes caused by varying other free parameters. It is therefore unlikely that non-spherical envelope structures would cause the best-fit values to be outside the estimated uncertainties, or that infall is not observed in these sources.

4.2. Infall vs. radius

Let us now consider the envelope velocity profile more closely, using IRAS4A as a working example. The optical depth in a given line decreases with offset from the source velocity. Thus even those lines which become optically thick quickly at the line centre probe a range of radii within the model in the line wings. However, there are limits to how far into the envelope our observations can probe the infalling material (cf. Figs. 7 and 8). Figure 11 shows the velocity profiles and the resulting line profiles for models where the infall velocity inside 700, 1000, 3000 and 10 000 AU was set to zero. For the model where infall stops at 10 000 AU, the absorption in the 2₀₂ −1₁₁ line is too deep and shifted too far towards the systemic velocity to be consistent with the observations. The emission in both the 2₀₂ − 1₁₁ and 2₁₂ − 1₀₁ lines is also too weak for the models where infall is stopped at 10 000 or 3000 AU. There is little difference within the uncertainties in the data between the best-fit model and those where infall stops at r ≤ 1000 AU, though the model where infall stops at 700 AU is a slightly better fit to the absorption for the 2₀₂ − 1₁₁ line. Thus, infall in IRAS4A must continue to radii of at least 1000 AU, consistent with the interferometric observations of Di Francesco et al. (2001). The observations are not able to constrain whether infall really stops at or inside this radius, but they can constrain the velocity profile outside this radius.

If infall proceeds in an inside-out fashion (Shu 1977), the outer-most radius at which infall still takes place, the free-fall radius (r_ff), is also of interest. Figure 12 shows a comparison of the velocity and line profiles for the best-fit model with models where the infall velocity outside 7000 and 10 000 AU was set to zero. Both models have too much absorption near the source velocity compared to the observations, so in IRAS4A infall must extend at least to the outer radius considered by the model. Assuming that infall proceeds in gravitational free-fall at the sound speed from the inside outwards, the free-fall timescale since infall began is given by: $\begin{array}{ccc}{\mathit{t}}_{\mathrm{ff}}\mathrm{=}\mathrm{2435}\left(\frac{\mathit{p}}{\mathrm{2.0}}\right){\left(\frac{{\mathit{r}}_{\mathrm{0}}}{\mathrm{10}\hspace{0.17em}\mathrm{AU}}\right)}^{-1}{\left(\frac{{\mathit{n}}_{\mathrm{0}}}{{\mathrm{10}}^{\mathrm{9}}\hspace{0.17em}{\mathrm{cm}}^{-3}}\right)}^{-1}{\mathrm{\int }}_{{\mathit{r}}_{\mathrm{0}}}^{{\mathit{r}}_{\mathrm{ff}}}{\left(\frac{{\mathit{r}}_{\mathrm{ff}}}{{\mathit{r}}_{\mathrm{0}}}\right)}^{\frac{\mathit{p}\mathrm{-}\mathrm{2}}{\mathrm{2}}}\mathrm{d}\mathit{r}\mathrm{yr}\mathit{,}& & \end{array}$(2)taking into account the variation of density through the envelope. r_ff is given in AU and p is again the slope of the density power-law which is given for our sources in Table 3. For IRAS4A this gives a free-fall timescale of ≳3 × 10⁵ yr.

Therefore, for IRAS4A infall must be taking place at least from 1000 AU to the outer edge of the model (11 200 AU). We define the minimum radius outside of which infall must take place as the minimum detectible infall radius (r_mdi), for which we can then calculate the corresponding velocity (ν_mdi) using Eq. (1) and ν₁₀₀₀. Having obtained the infall velocity, this can be used to calculate the central gravitational mass (M_g) from: $\begin{array}{ccc}{\mathit{M}}_{\mathrm{g}}\mathrm{=}\frac{{\nu }_{\mathrm{mdi}}^{\mathrm{2}}\hspace{0.17em}{\mathit{r}}_{\mathrm{mdi}}}{\mathrm{2}\mathit{G}}\mathrm{=}\mathrm{0.56}\hspace{0.17em}{\left(\frac{{\nu }_{\mathrm{1000}}}{\mathrm{1}\hspace{0.17em}{\mathrm{km}\hspace{0.17em}\mathrm{s}}^{-1}}\right)}^{\mathrm{2}}{\mathit{M}}_{\mathrm{\odot }}\mathit{,}& & \end{array}$(3)and the instantaneous mass infall rate at r_mdi: $\begin{array}{ccc}\mathit{Ṁ}\inf \mathrm{=}\mathrm{4}\mathit{\pi }\hspace{0.17em}\mathit{\mu }\hspace{0.17em}{\mathit{m}}_{\mathrm{H}}\hspace{0.17em}{\mathit{r}}_{\mathrm{mdi}}^{\mathrm{2}}\hspace{0.17em}{\mathit{n}}_{\mathrm{mdi}}\hspace{0.17em}{\nu }_{\mathrm{mdi}}\mathit{,}& & \end{array}$(4)where μ is the mean molecular weight, for which we assume 2.8 from Kauffmann et al. (2008). When Eq. (1) and the power-law density dependence are taken into account, this becomes: $\begin{array}{ccc}\mathit{Ṁ}\inf & \mathrm{=}& \mathrm{2.08}\mathrm{\times }{\mathrm{10}}^{-8}\left(\frac{{\mathit{n}}_{\mathrm{0}}}{{\mathrm{10}}^{\mathrm{9}}\hspace{0.17em}{\mathrm{cm}}^{-3}}\right)\\ & & \end{array}$(5)Due to the assumption of spherical symmetry, Ṁ_inf is primarily a measure of the accretion rate of material from the envelope onto any disk-like structure rather than necessarily the accretion rate of material directly onto the central object (Ṁ_acc). Therefore, Ṁ_inf is not expected to be directly related to the stellar luminosity, though higher Ṁ_inf probably leads to higher Ṁ_acc. Dividing M_g by Ṁ_inf gives the infall timescale (t_inf), assuming that the accretion is steady over time at the rate measured at r_mdi. If p ≠ 1.5 the accretion will, by definition, not be constant over time, but this assumption enables a simple analytic estimate of the infall timescale, such that: $\begin{array}{ccc}{\mathit{t}}_{\inf }& \mathrm{=}& \mathrm{2.7}\mathrm{\times }{\mathrm{10}}^{\mathrm{7}}{\left(\frac{{\mathit{n}}_{\mathrm{0}}}{{\mathrm{10}}^{\mathrm{9}}\hspace{0.17em}{\mathrm{cm}}^{-3}}\right)}^{-1}\\ & & \end{array}$(6)The values of M_g, r_mdi, Ṁ_inf, t_ff and t_inf are given in Table 4. The free-fall and infall timescales give characteristic timescales for collapse of the envelope under different sets of assumptions, which probably lead to them being accurate to a factor of a few at best. That r_ff is always the outer edge of the model might be an indication that the collapse is global, as in e.g. the Larson-Penston solution (Larson 1969; Penston 1969), or outside-in rather than inside-out. Whether constant mass infall and accretion is a good assumption will be discussed in Sect. 6, as well as the wider implications of our results and comparison to other studies.

4.3. Whole sample

The results for the other sources which show inverse P-Cygni profiles can be divided into two groups, with comparison between the best-fit models and the data shown in Figs. A.1 to A.6. IRAS 15398, L1527, L1157 and BHR71 can all be fit with models similar to IRAS4A with similarly reasonable infall velocities and b values. Though the model for BHR71 produces too much emission in the 1₁₁ −0₀₀ line, the absorption component of this line is well fit, as is the 1₁₀ − 1₀₁ line. Non-detections are confirmed in the higher-excited lines. The reasons for this will be discussed in Sect. 5, however this results in the infall velocities being less well constrained because the ground-state lines are less sensitive to infall. It also means that our observations only probe infall down to ~3000 AU (i.e. core to envelope) scales in these sources.

Inspection of archival single-dish HCO⁺ (J = 4−3) observations for all these sources, as well as IRAS4A, show asymmetric blue-shifted profiles (Gregersen et al. 1997; Hogerheijde & Sandell 2000; van Kempen et al. 2009b), suggesting that infall likely continues beyond this radius in the dense inner envelope. For L1157 the HCO⁺ observations present a more complex picture, as the J = 3−2 spectrum presented by Gregersen et al. (1997) shows an asymmetric blue-shifted profile while the J = 4−3 is also asymmetric but with slightly stronger emission in the red than the blue peak. Tobin et al. (2012a) concluded that the L1157 envelope is elongated, with the velocity field consistent with material flowing along the filament onto the envelope on similar scales to our model. Material at the outer edge of the envelope must be infalling in all sources.

On the other hand, Ser-SMM4 requires a much higher infall velocity which seems inconsistent with infall on protostellar scales, particularly given the fact that it would result in a gravitational mass more than an order of magnitude larger than the envelope mass for this heavily embedded source. Indeed, the velocity for this model at the outer radius is large enough that even the ground-state lines are dominated by bulk motions rather than turbulence. The single-dish HCO⁺ (J = 4−3) observations of Ser-SMM4 are inconclusive, with a very faint second peak redshifted by ~1.5 km s^-1 from the source velocity (Ramchandani 2012).

Duarte-Cabral et al. (2010) found two emission components separated by ~1.5 km s^-1 in C¹⁸O (J = 1−0) observations towards the Serpens Main cloud, which Duarte-Cabral et al. (2011) successfully modelled as a collision of two clouds. SMM4 is situated on the boundary between these two velocity components which correspond to the absorption and emission components observed in the water ground-state lines. In addition, the velocity at the outer edge of the best-fit model for SMM4 is ~1.3 km s^-1. Our observations of Ser-SMM4 would therefore seem more consistent with tracing this cloud-scale flow or a foreground cloud rather than infall on envelope scales.

Judging by the fit values alone, GSS30 would seem to lie in the former group. However, this source actually lies near the prominent outflow VLA1623 and a diffuse layer has been detected in front of the Ophiuchus cloud where both sources are situated (van Kempen et al. 2009c). Bjerkeli et al. (2012) observe inverse P-Cygni profiles across the whole of their H₂O 1₁₀ − 1₀₁ map towards VLA1623 which are very similar to those observed towards GSS30 (see Fig. 13). Large-scale flows or foreground clouds would therefore seem to be the cause of the observed line profiles, which is also more consistent with GSS30’s older evolutionary stage.

5. Chemistry

As discussed in Sect. 3.2, most previous attempts to model water lines have used drop or jump abundance profiles, rather than the simplified chemical model employed here. Figure 14 shows a comparison between line profiles for the best-fit model with the best possible fit using a drop abundance profile. The adopted abundance profiles are shown in Fig. 4. For the drop profile, the inner abundance is held constant at 6 × 10^-4 where T ≥ 100 K, consistent with the simple water profile. The photodesorption layer extends up to A_V = 5 but is still within the grid rather than as a separate layer as in Coutens et al. (2012). The abundance in the outer envelope and photodesorption layer are free parameters. The best fit was obtained for X_out = 3 × 10^-10 and X_pdl = 3 × 10^-7. The difference between the two models is relatively small for the lower excited line, but the step profile results in unobserved emission in the 2₁₁ − 2₀₂ line and has a much too wide emission component and too narrow absorption component in the 2₀₂ − 1₁₁ line. Indeed, only a very small set of profiles produce absorption below the continuum in the 2₀₂ − 1₁₁ line and none reproduce the equal intensity and width in the absorption and emission components observed.

This comparison also serves to show that if the gas-phase water abundance were higher in the inner, warmer part of the envelope then the 2₀₂ −1₁₁ and 2₁₁ −2₀₂ lines would be detected. This is true for the other sources as well, so that our non-detection of these lines is confirmed by the models also confirms that freeze-out of water is strong in the inner envelope outside the T = 100 K radius. That we detect the 2₀₂ −1₁₁ line in IRAS4A is probably due to a unique combination of the density and temperature profiles leading to just the right abundance and excitation structure. The non-detection of the 3₁₂ − 2₂₁ and 3₁₂ −3₀₃ lines in all sources is because these lines are primarily sensitive to the T ≥ 100 K region, which for IRAS4A has a diameter of ~1′′ and so is too heavily beam-diluted.

This sensitivity of the water observations to the water abundance profile means that at least simple water chemistry including density-dependent freezeout must be taken into account when performing modelling of these lines. These observations do not probe the hot core (i.e. T ≥ 100 K) region, for which we refer the reader to Persson et al. (2012) and Visser et al. (2013).

Let us now consider the impact of changing the key variables involved in the simple water chemistry model. Figure 15 shows how the abundance profile varies with varying interstellar radiation field and cosmic-ray UV rate. As might be expected, increasing G_isrf from the standard value increases the destruction of water in the very outer layer of the model and thus decreases the abundance. Increasing or decreasing G_cr leads to more or less water in the gas phase respectively due to increasing or decreasing the rate of cosmic-ray induced photodesorption.

The resulting model line profiles are shown in Fig. 16. Increasing G_cr leads to a detectable emission feature in the 2₁₁ −2₀₂ line, which is not observed, and too broad and intense emission and absorption in the 2₀₂ −1₁₁ line. Decreasing G_cr leads to too little emission in the 2₀₂ −1₁₁ and 2₁₂ −1₀₁ lines. Increasing G_isrf has a minimal effect on the line profiles. This is because the water affected by this change is closest to the source velocity and predominantly in the ground state where the absorption is already saturated. For those sources where the outer envelope density is above 10⁵ cm^-3 there is a significant difference between G_isrf = 0 and G_isrf = 1 (see lower plot of Fig. 16). This is because the photodesorption layer caused by G_isrf = 1 results in a significant jump in water abundance (see right-hand panel of Fig. 15) and the density is high enough that photodissociation only occurs in a very thin shell at the outer edge.

The sensitivity of the water line profiles not only to the radial abundance profile but also to its absolute scaling makes it possible to constrain the cosmic-ray induced UV field. As discussed in Sect. 3.2, the cosmic-ray induced UV flux is scaled by G_cr to the interstellar photon flux, which leads to photodesorption of H₂O with an efficiency (η) which we assume to be 10^-3. We find that the canonical value for the cosmic-ray induced UV flux (i.e. G_cr = 10^-4, Prasad & Tarafdar 1983; Gredel et al. 1989; Shen et al. 2004) gives the best fit, in contrast to the 1 dex enhancement (i.e. G_cr = 10^-3) suggested by Hollenbach et al. (2009) and used for the prestellar core L1544 by Caselli et al. (2012). Indeed, even an enhancement of 0.25 dex is already a noticeably poorer fit to our observations. Formally, we are only able to constrain the product of G_cr and η. However, given that η is well constrained from laboratory experiments and theoretical calculations (Öberg et al. 2009; Arasa et al. 2010), this parameter is already well constrained.

Most sources require a cosmic-ray induced UV field in the range (0.2−1) × 10⁴ photons cm^-2 s^-1, which is consistent with estimates for standard molecular cloud conditions (Shen et al. 2004). The exception is GSS30 which requires G_cr = 3 × 10^-4, but as discussed above the observed line profiles likely trace a foreground cloud or cloud-scale flow. In addition, the region is disrupted by the outflow from VLA1623 and is embedded in a weak PDR powered by S1 and/or HD147889. Both of these factors could lead to abundance enhancements which are not accounted for in our model.

For those sources which require lower values of G_cr, it should be pointed out that this does not necessarily mean that they actually have lower cosmic-ray rates, but rather this is the main variable within our model which varies the abundance and thus the emission intensity. It could be instead that the outflow blocks some fraction of the emission from the back half of the envelope (see Sect. 3.5). Indeed, for BHR71 there is still a significant excess of emission in the 1₁₁ − 0₀₀ line which cannot be improved by varying any of the model parameters. The outflow of BHR71 has an hourglass shape at the base in CO observations (Parise et al. 2006; van Kempen et al. 2009a) which would fill more of the 19′′ beam of the 1₁₁ − 0₀₀ line than the 38′′ beam of the 1₁₀ −1₀₁ line. Thus the blocking of emission by the outflow may be the cause of the lower intensity rather than a physically lower value of G_cr. As the absorption is both less sensitive to the abundance and takes place in the part of the envelope between the outflow and the observer, the good fit by the models to this part of the line profiles still constrains the velocity profile down to r_mdi.

6. Discussion

Let us now discuss the consistency of our results with those obtained using other tracers. For IRAS4A, Di Francesco et al. (2001) calculated an infall velocity of 0.68 km s^-1 at 2700 AU from interferometric observations of H₂CO using a 2-slab model. Our model gives 0.67 km s^-1 at this radius, fully consistent with the results of Di Francesco et al. (2001) given their assumptions. Our mass infall rate is also consistent with the results of Belloche et al. (2006) using who used radiative transfer of compression-induced collapse models to fit various single-dish observations towards IRAS4A. The value of b required to fit the water observations of IRAS4A is lower than the 0.8 km s^-1 used by Yıldız et al. (2012) for various C¹⁸O lines. While those authors use the same source structure as used here, they only considered a static model which leaves turbulence as the dominant component to the line width. On the other hand, Jørgensen et al. (2004) show that a pure infall (i.e. b = 0) can also be ruled out, so the observed C¹⁸O linewidths must be due to a combination of infall and turbulence. For the sample as a whole we find best-fit values of b in the range 0.2 to 0.9 km s^-1, corresponding to σ_NT values in the range 0.15 to 0.64 km s^-1.

At the minimum detectible infall radius of 1000 AU, the temperature in the model is 20 K which corresponds to a sound speed of 0.26 km s^-1 while the infall velocity is 1.1 km s^-1, so the infall has a Mach number of 4.3. This flow will not result in a shock as long as the mass infall rate does not decrease with decreasing radius. However, the outflow timescale derived by Yıldız et al. (2012) (~6 − 10 × 10⁴ yr) and the free-fall time of the envelope (2 × 10⁵ yr) are approximately an order of magnitude larger than the infall timescale (≳30 000 yr), suggesting that something must slow the infalling material in the inner ~1000 AU. Whether there is a shock formed as material falls onto a flattened disk-like structure will depend on how gradually the radial infall is mediated. Conservation of angular momentum could help to enhance circular motions in the infalling material. In addition, IRAS4A is known to have an hour-glass shaped magnetic field on ~400 AU scales (Girart et al. 2006) which could cause a decrease in the component of motion perpendicular to the axis of the magnetic field. While evidence for a disk shock due to infall onto a disk was claimed due to the existence of hot water in the neighbouring source NGC 1333-IRAS4B by Watson et al. (2007), this was subsequently refuted due to the spatial offset of the emission from the central source by Herczeg et al. (2012). CH₃OH masers have been claimed to be excited in disk shocks in high-mass star forming regions (Torstensson et al. 2011), but have not been detected towards low-mass sources. Evidence for such a shock is therefore limited or otherwise difficult to disentangle from shocks related to the outflow.

The stability of any central condensation depends on how well it can mediate the competing processes which serve to increase or decrease its mass over time. The instantaneous mass infall rate we determine for IRAS4A (2.0 × 10^-4M_⊙ yr^-1) is higher by approximately one order of magnitude than the maximum mass accretion rate in simulations of forming protostars (e.g. Dunham & Vorobyov 2012). From conservation of energy, the maximum radius that a star can have is: $\begin{array}{ccc}{\mathit{R}}_{\mathrm{\ast }}\mathrm{=}\mathrm{61}{\mathit{R}}_{\mathrm{\odot }}\hspace{0.17em}\left(\frac{{\mathit{M}}_{\mathrm{\ast }}}{{\mathit{M}}_{\mathrm{\odot }}}\right)\mathit{,}& & \end{array}$(7)(Palla 1999), where M_∗ is the stellar mass of the protostar. However, using this equation means that the mass drops out of the mass accretion rate calculation. Thus if we assume that all the luminosity of the star is due to accretion, the maximum possible mass accretion rate for IRAS4A is 2 × 10^-5M_⊙ yr^-1, which is an order of magnitude lower than the measured mass infall rate.

The mass infall rate is also 3.5 orders of magnitude higher than the mass-loss rate of 5 × 10^-8M_⊙ yr^-1 in the atomic jet obtained from the [OI] emission by Karska et al. (2013) and at least one order of magnitude higher than the mass outflow rate from low-J CO observations of swept-up gas in the molecular outflow (Yıldız et al. 2012, and in prep.). This suggests that whatever flattened structure exists around IRAS4A, it is gaining mass significantly faster than it is likely to be accreted onto the central protostar or ejected via winds and/or outflows and so is unlikely to be gravitationally stable. How long this has been the case is unclear.

Indeed, IRAS4A is also a known binary with a separation of ~400 AU (Choi et al. 2004), likely formed through turbulent fragmentation of the initial prestellar core (e.g. Offner et al. 2010). The presence of a binary provides additional forces on material in the inner regions, possibly enhancing instability of any circumbinary or circumstellar disk(s). Such instability could lead to episodic accretion, which has been suggested for some time (e.g. Kenyon et al. 1990) and has been predicted to occur due to gravitational instability of a disk caused by mass-loading (e.g. Vorobyov 2009; Vorobyov & Basu 2010) in a similar scenario to that taking place in IRAS4A. It could even lead to the formation of additional companions through rotational fragmentation of the disk (e.g. Burkert & Bodenheimer 1993).

L1527 has similar mass infall and outflow rates, though the maximum mass accretion rate is an order of magnitude smaller. This is primarily due to the small infall velocity resulting in a small gravitational mass though this source shows one of the strongest asymmetries in HCO⁺ (J = 4 − 3) and (J = 3 − 2) observations (Hogerheijde & Sandell 2000). The insensitivity of the water observations to the infall velocity in this source, primarily due to the relatively flat density profile, therefore makes it likely that the infall velocity obtained for L1527 is a lower limit. IRAS 15398 is similar to IRAS4A in that the mass infall rate is an order of magnitude higher than the mass outflow and maximum mass accretion rate. BHR71 and L1157 have similar mass infall, accretion and outflow rates and so may be in more stable configurations. Infall in BHR71, L1157 and IRAS 15398 is supersonic at r_mdi. A larger sample and detailed mapping of the velocity field down to smaller radii is required to tell whether this is related to initial conditions or an evolutionary effect.

The gravitational masses (M_g) given in Table 4 were calculated assuming that the infall speed is solely dependant on the mass of the protostar, and range from 1.17 to 0.08 M_⊙. The lowest value is for L1527, which is half an order of magnitude lower than obtained by Tobin et al. (2012b) based on interferometric observations. As discussed above, the infall velocity for L1527 is particularly uncertain, and the velocity corresponding to the result obtained by Tobin et al. (2012b) does lie within our uncertainties. However, there is also the possibility that infall has only recently started in the outer parts of the envelope as probed by our observations or that the outer layers are not gravitationally unstable.

Finally, it is at this point worth considering why we do not observe infall signatures towards all 29 WISH sources. Four sources show regular P-Cygni profiles, indicative of envelope expansion but what of the remaining 18? Low S/N may play a role in a few observations, but in most cases outflow emission is observed and so red-shifted absorption should be detectible. NGC 1333-4B and DK-Cha both have pole-on outflows which would make detection of the envelope almost impossible. Similarly, the broad absorption towards Serpens SMM1 and SMM3 is undoubtedly caused by the same cloud-scale effects as in SMM4 but because neither is in the region where the two cloud components overlap there is no inverse P-Cygni profile. The more evolved sources may simply not have enough envelope remaining to show inverse or regular P-Cygni profiles. For the remaining younger sources, it could be that the outer layers of the envelope have yet to start infalling, or that the outflow blocks or disrupt enough of the emission and absorption to make identification of infall with water impossible. Alternatively, infall on large-scales may have already stopped or been reversed, in which case this phase is short-lived.

7. Conclusions

We have presented 1–D non-LTE radiative transfer models of water in seven protostars which show infall signatures in Herschel WISH observations. Though these models are spherically-symmetric and include a simplified implementation of the outflow, they nevertheless reproduce the observations well while being constrained by them. Water proves to be an excellent tracer of both the dynamics and chemistry of protostellar envelopes, particularly when both ground-state high-τ lines and more moderately optically thick excited lines are observed together. We find that a step abundance model cannot reproduce the observations while a simple chemical network with reasonable parameters is consistent with the observations. When this is done we find best-fit values for the cosmic-ray induced UV flux in the range (1−5) × 10⁴ times lower than the interstellar UV radiation field, which is in line with estimates for standard conditions in molecular clouds.

Infall of material from the core into the envelope is observed in four sources from ~10 000 to 3000 AU scales and in one source, NGC 1333-IRAS4A, down to at least ~1000 AU. This is only possible for IRAS4A due to the detection in the moderately optically thick 2₀₂ − 1₁₁ line, so future studies should include mid-τ lines where possible in preference to higher-τ lines. In four of the five sources this infall is supersonic and infall in all sources must take place at the outer edge of the envelope, which may be evidence that collapse is global or outside-in rather than inside-out. The mass infall rates derived for IRAS4A and IRAS15398 are so high that any disk-like structure cannot be gravitationally stable, which may result in episodic accretion or fragmentation. However, not all observations showing infall signatures in water can be reproduced with 1–D infalling envelope models. In the case of Serpens-SMM4 and GSS30, the observations are more likely tracing large-scale cloud flows. Using slab-model analysis alone cannot distinguish between these two scenarios. While care must be taken when interpreting inverse P-Cygni or infall profiles, the combination of slab and 1–D envelope models is able to constrain the radial velocity distribution, provided that multiple transitions of a molecule which probes a range of temperatures and optical depths are used.

Online material

Appendix A: Supplementary material

Figures A.1 to A.6 show comparisons between the observations and the best-fit models for the larger sample of sources. Table A.1 provides the observation identification numbers for all Herschel observations included in this paper.

Acknowledgments

The authors would like to thank the anonymous referee, whose comments helped improved the clarity and content of the paper. We also thank Floris van der Tak, Mario Tafalla, Paola Caselli, Neal Evans and Audrey Coutens for their valuable comments on the manuscript, and Daniel Harsono, Catherine Walsh, Pamela Klaassen and Tobias Albertsson for stimulating discussions. JCM is funded by grant 614.001.008 from the Netherlands Organisation for Scientific Research (NWO). Astrochemistry in Leiden is supported by the Netherlands Research School for Astronomy (NOVA), by a Spinoza grant and by the European Community’s Seventh Framework Programme FP7/2007-2013 under grant agreement 238 258 (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

All Tables

All Figures