© ESO, 2014

1. Introduction

Shocks are ubiquitous phenomena where outflow-envelope interactions take place in young stellar objects (YSO). Large-scale shocks are caused by the bipolar jets and protostellar winds impinging on the envelope along the cavity walls carved by the passage of the jet (Arce et al. 2007; Frank et al. 2014). This important interaction needs to be characterized in order to understand and quantify the feedback from protostars onto their surroundings and, ultimately, to explain the origin of the initial mass function, disk fragmentation, and the binary fraction.

Theoretically, shocks are divided into two main types based on a combination of magnetic field strength, shock velocity, density, and level of ionization (Draine 1980; Draine et al. 1983; Hollenbach et al. 1989; Hollenbach 1997). In continuous (C-type) shocks, in the presence of a magnetic field and low ionization, the weak coupling between the ions and neutrals results in a continuous change in the gas parameters. Peak temperatures of a few 10³ K allow the molecules to survive the passage of the shock, which is therefore referred to as non-dissociative. In jump (J-type) shocks, physical conditions change in a discontinuous way, leading to higher peak temperatures than in C shocks of the same speed and for a given density. Depending on the shock velocity, J shocks are either non-dissociative (velocities below about ~30 km s^-1, peak temperatures of a few 10⁴ K) or dissociative (peak temperatures even exceeding 10⁵ K), but the molecules efficiently reform in the post-shock gas.

Shocks reveal their presence most prominently in the infrared (IR) domain, where the post-shock gas is efficiently cooled by numerous atomic and molecular emission lines. Cooling from H₂ is dominant in outflow shocks (Nisini et al. 2010b; Giannini et al. 2011b), but its mid-IR emission is strongly affected by extinction in the dense envelopes of young protostars (Giannini et al. 2001a; Nisini et al. 2002; Davis et al. 2008; Maret et al. 2009). In the far-IR, rotational transitions of water vapor (H₂O) and carbon monoxide (CO) are predicted to play an important role in the cooling process (Goldsmith & Langer 1978; Neufeld & Dalgarno 1989; Hollenbach 1997) and can serve as a diagnostic of the shock type, its velocity, and the pre-shock density of the medium (Hollenbach et al. 1989; Kaufman & Neufeld 1996; Flower & Pineau des Forêts 2010, 2012).

The first observations of the critical wavelength regime to test these models (λ ~ 45−200μm) were taken using the Long-Wavelength Spectrometer (LWS, Clegg et al. 1996) on board the Infrared Space Observatory (ISO, Kessler et al. 1996). Far-IR atomic and molecular emission lines were detected toward several low-mass deeply-embedded protostars (Nisini et al. 2000; Giannini et al. 2001a; van Dishoeck 2004), but its origin was unclear because of the poor spatial resolution of the telescope (~ 80′′, Nisini et al. 2002; Ceccarelli et al. 2002).

The sensitivity and spectral resolution of the Photodetector Array Camera and Spectrometer (PACS, Poglitsch et al. 2010) on board Herschel allowed a significant increase in the number of detections of far-IR lines in young protostars compared with the early ISO results and has revealed rich molecular and atomic line emission both at the protostellar (e.g., van Kempen et al. 2010; Goicoechea et al. 2012; Herczeg et al. 2012; Visser et al. 2012; Green et al. 2013; Karska et al. 2013; Lindberg et al. 2014; Manoj et al. 2013; Wampfler et al. 2013) and at pure outflow positions (Santangelo et al. 2012, 2013; Codella et al. 2012; Lefloch et al. 2012a; Vasta et al. 2012; Nisini et al. 2013).

The unprecedented spatial resolution of PACS allowed detailed imaging of L1157 providing firm evidence that most far-IR H₂O emission originates in the outflows (Nisini et al. 2010a). A mapping survey of about 20 protostars revealed similarities between the spatial extent of H₂O and high-J CO (Karska et al. 2013). Additional strong flux correlations between these species and similarities in the velocity-resolved profiles (Kristensen et al. 2010, 2012; San José-García et al. 2013; Santangelo et al. 2014) suggest that the emission from the two molecules arises from the same regions. This is further confirmed by finely-spatially sampled PACS maps in CO 16−15 and various H₂O lines in shock positions of L1448 and L1157 (Santangelo et al. 2013; Tafalla et al. 2013). On the other hand, the spatial extent of OH resembles the extent of [O i] and, additionally, a strong flux correlation between the two species is found (Karska et al. 2013; Wampfler et al. 2013). Therefore, at least part of the OH emission most likely originates in a dissociative J-shock, together with [O i] (Wampfler et al. 2010, 2013; Benedettini et al. 2012).

To date, comparisons of the far-IR observations with shock models have been limited to a single source or its outflow positions (e.g., Nisini et al. 1999; Benedettini et al. 2012; Vasta et al. 2012; Santangelo et al. 2012; Dionatos et al. 2013; Lee et al. 2013; Flower & Pineau des Forêts 2013). Even in these studies, separate analysis of each species (or different pairs of species) often led to different sets of shock properties that needed to be reconciled. For example, Dionatos et al. (2013) show that CO and H₂ line emission in Serpens SMM3 originates from a 20 km s^-1J shock at low pre-shock densities (~ 10⁴ cm^-3), but the H₂O and OH emission is better explained by a 30−40 km s^-1C shock. In contrast, Lee et al. (2013) associate emission from both CO and H₂O in L1448-MM with a 40 km s^-1C shock at high pre-shock densities (~ 10⁵ cm^-3), consistent with the analysis of the L1448-R4 outflow position (Santangelo et al. 2013). Analysis restricted to H₂O lines alone often indicates an origin in non-dissociative J shocks (Santangelo et al. 2012; Vasta et al. 2012; Busquet et al. 2014), while separate analysis of CO, OH, and atomic species favors dissociative J shocks (Benedettini et al. 2012; Lefloch et al. 2012b). The question remains how to break degeneracies between these models and how typical the derived shock properties are for young protostars. In addition, surveys of high-J CO lines with PACS have revealed two universal temperature components in the CO ladder toward all deeply-embedded low-mass protostars (Herczeg et al. 2012; Goicoechea et al. 2012; Green et al. 2013; Karska et al. 2013; Kristensen et al. 2013; Manoj et al. 2013; Lee et al. 2013). The question is how the different shock properties relate with the “warm” (T ~ 300 K) and “hot” (T ≳ 700 K) components seen in the CO rotational diagrams.

In this paper, far-IR spectra of 22 low-mass YSOs observed as part of the William Herschel Line Legacy (WILL) survey (PI: E. F. van Dishoeck) are compared to the shock models from Kaufman & Neufeld (1996) and Flower & Pineau des Forêts (2010). All sources are confirmed deeply-embedded YSOs located in the well-studied Perseus molecular cloud spanning the Class 0 and I regime (Knee & Sandell 2000; Enoch et al. 2006, 2009; Jørgensen et al. 2006, 2007; Hatchell et al. 2007a,b; Davis et al. 2008; Arce et al. 2010). H₂O, CO, and OH lines are analyzed together for this uniform sample to answer the following questions: Do far-IR line observations agree with the shock models? How much variation in observational diagnostics of shock conditions is found between different sources? Can one set of shock parameters explain all molecular species and transitions? Are there systematic differences between shock characteristics inferred using the CO lines from the “warm” and “hot” components? How do shock conditions vary with the distance from the powering protostar?

This paper is organized as follows. Section 2 describes our source sample, instrument with adopted observing mode, and reduction methods. Section 3 presents the results of the observations: line and continuum maps, and the extracted spectra. Section 4 shows comparison between the observations and shock models. Section 5 discusses results obtained in Sects. 4 and 6 presents the conclusions.

2. Observations

All observations presented here were obtained as part of the “William Herschel Line Legacy” (WILL) OT2 program on Herschel (Mottram et al., in prep.). The WILL survey is a study of H₂O lines and related species with PACS and the Heterodyne Instrument for the Far-Infrared (HIFI, de Graauw et al. 2010) toward an unbiased flux-limited sample of low-mass protostars newly discovered in the recent Spitzer (c2d, Gutermuth et al. 2009, 2010; Evans et al. 2009) and Herschel (André et al. 2010) Gould Belt imaging surveys. Its main aim is to study the physics and chemistry of star-forming regions in a statistically significant way by extending the sample of low-mass protostars observed in the “Water in star-forming regions with Herschel” (WISH, van Dishoeck et al. 2011) and “Dust, Ice, and Gas in Time” (DIGIT, Green et al. 2013) programs.

This paper presents the Herschel/PACS spectra of 22 low-mass deeply-embedded YSOs located exclusively in the Perseus molecular cloud (see Table 1) to ensure the homogeneity of the sample (similar ages, environment, and distance). The sources were selected from the combined SCUBA and Spitzer/IRAC and MIPS catalog of Jørgensen et al. (2007) and Enoch et al. (2009), and all contain a confirmed embedded YSO (Stage 0 or I, Robitaille et al. 2006, 2007) in the center.

The WILL sources were observed using the line spectroscopy mode on PACS which offers deep integrations and finely sampled spectral resolution elements (minimum 3 samples per FWHM depending on the grating order, PACS Observer’s Manual¹) over short wavelength ranges (0.5−2 μm). The line selection was based on the prior experience with the PACS spectra obtained over the full far-infrared spectral range in the WISH and DIGIT programs and is summarized in Table A.1. Details of the observations of the Perseus sources within the WILL survey are shown in Table A.2.

PACS is an integral field unit with a 5×5 array of spatial pixels (hereafter spaxels) covering a field of view of ~ 47′′ × 47′′. Each spaxel measured ~ 9.4′′ × 9.4′′, or 2×10^-9 sr, and at the distance to Perseus (d = 235 pc, Hirota et al. 2008) resolves emission down to ~2300 AU. The total field of view is about 5.25 × 10^-8 sr and ~11 000 AU. The properly flux-calibrated wavelength ranges include: ~55−70 μm, ~72−94 μm, and ~105−187 μm, corresponding to the second (<100 μm) and first spectral orders (>100 μm). Their respective spectral resolving power are R ~ 2500−4500 (velocity resolution of Δν ~ 70−120 km s^-1), 1500−2500 (Δν ~ 120−200 km s^-1), and 1000−1500 (Δν ~ 200−300 km s^-1). The standard chopping-nodding mode was used with a medium (3′) chopper throw. The telescope pointing accuracy is typically better than 2′′ and can be evaluated to first order using the continuum maps.

The basic data reduction presented here was performed using the Herschel Interactive Processing Environment v.10 (hipe, Ott 2010). The flux was normalized to the telescopic background and calibrated using observations of Neptune. Spectral flatfielding within HIPE was used to increase the signal-to-noise (for details, see Herczeg et al. 2012; Green et al. 2013). The overall flux calibration is accurate to ~ 20%, based on the flux repeatability for multiple observations of the same target in different programs, cross-calibrations with HIFI and ISO, and continuum photometry.

Custom idl routines were used to extract fluxes using Gaussian fits with fixed line width (for details, see Herczeg et al. 2012). The total, 5×5 line fluxes were calculated by co-adding spaxels with detected line emission, after excluding contamination from other nearby sources except Per 2, 3, 10, and 18, where spatial separation between different components is too small (see Sect. 3.1). For sources showing extended emission, the set of spaxels providing the maximum flux was chosen for each line separately. For point-like sources, the flux is calculated at the central position and then corrected for the PSF using wavelength-dependent correction factors (see PACS Observer’s Manual). Table A.3 shows the line detections toward each source, while the actual fluxes will be tabulated in the forthcoming paper for all the WILL sources (Karska et al., in prep.). Fits cubes containing the spectra in each spaxel will be available for download in early 2015 at http://www.strw.leidenuniv.nl/WISH/.

3. Results

In the following sections, PACS lines and maps of the Perseus YSOs are presented. Most sources in this sample show emission in just the central spaxel. Only a few sources show extended emission and those maps are compared to maps at other wavelengths to check for possible contamination by other sources and their outflows. In this way, the spaxels of the maps with emission originating from our objects are established and line fluxes determined over those spaxels. The emergent line spectra are then discussed.

3.1. Spatial extent of line emission

Table 2 provides a summary of the patterns of H₂O 2₁₂−1₀₁ line and continuum emission at 179 μm for all WILL sources in Perseus (maps are shown in Figs. A.1 and A.2). The mid-infrared continuum and H₂ line maps from Jørgensen et al. (2006) and Davis et al. (2008, including CO 3−2 observations from Hatchell et al. 2007a are used to obtain complementary information on the sources and their outflows. For a few well-known outflow sources, large-scale CO 6−5 maps from Yıldız et al. (2014) are also considered.

As shown in Table 2, the majority of the PACS maps toward Perseus YSOs do not show any extended line emission. The well-centered continuum and line emission originates from a single object and an associated bipolar outflow for 12 out of 22 sources. Among the sources with spatially-resolved extended emission on the maps, various reasons are identified for their origin as illustrated in Fig. 1. In the map of Per 2, contribution from three nearby protostars and a strong outflow from the more distant L1448-MM source cause the extended line and continuum pattern. Emission in the Per 12 map is detected away from the continuum peak, but the emission originates from a large-scale outflow from NGC 1333-IRAS4A, not the targeted source. The H₂O emission in the Per 20 map is detected in the direction of the strong outflow and its extent is only slightly affected by the small mispointing revealed by the asymmetric continuum emission.

Extended emission beyond the well-centered continuum, as in the case of Per 20, is seen clearly only in Per 9 and Per 21-22. Additionally, the continuum peaks for Per 3, 8, and 11 are off-center, whereas the line emission peaks on-source, suggesting that some extended line emission is associated with the source itself and not only due to the mispointing². Similar continuum patterns are seen in Per 6 and Per 10, but here the line emission peaks a few spaxels away from the map center. In both cases, contribution from additional outflows/sources is the cause of the dominant off-source line emission.

To summarize, when the contamination of other sources and their outflows can be excluded, Perseus YSOs show that the H₂O 2₁₂−1₀₁ line emission is either well-confined to the central position on the map or shows at best weak extended line emission (those are marked with “e” in Table 2). In total, 7 out of 22 sources show extended emission in the H₂O 2₁₂−1₀₁ line associated with the targeted sources. Emission in CO, OH, and other H₂O lines follows the same pattern (see Fig. A.3). Similarly compact emission was seen in a sample of 30 protostars surveyed in the DIGIT program (Green et al. 2013). In contrast, the WISH PACS survey (Karska et al. 2013) revealed strong extended emission in about half of the 20 low-mass protostars. There, the analysis of patterns of molecular and atomic emission showed that H₂O and CO spatially co-exist within the PACS field of view, while OH and [O i] lines are typically less extended, but also follow each other spatially and not H₂O and CO.

3.2. Line detections

In the majority of our sources, all targeted rotational transitions of CO, H₂O, and OH are detected, see Fig. 2 and Tables A.1 and A.3. The [O i] line at 63 μm and the [C ii] line at 158 μm will be discussed separately in a forthcoming paper including all WILL sources (Karska et al., and in prep.) and are not included in the figure and further analysis.

The CO 16−15 line with an upper level energy (E_u/k_B) of about 750 K is seen in 17 (~80%), the CO 24−23 (E_u/k_B ~ 1700 K) in 16 (~70%), and the CO 32−31 (E_u/k_B ~ 3000 K) in 8 (~30%) sources. The most commonly detected ortho-H₂O lines are: the 2₁₂−1₀₁ line at 179 μm (E_u/k_B ~ 110 K) and the 4₀₄−3₁₃ line at 125 μm (E_u/k_B ~ 320 K), seen in 15 sources (~70%), whereas the 6₁₆−5₀₅ line at 82 μm (E_u/k_B ~ 640 K) is detected in 13 sources (~ 60%). The para-H₂O line 3₂₂−2₁₁ at 90 μm (E_u/k_B ~ 300 K) is seen toward 9 sources (~50%).

Three OH doublets targeted as part of the WILL survey, the OH ${\mathrm{2}}^{}{\mathrm{\Pi }}_{{}^{\mathrm{1}}\mathit{/}\mathrm{2}}\mathit{,J}\mathrm{=}{}^{\mathrm{1}}\mathit{/}\mathrm{2}\mathrm{-}{\mathrm{2}}^{}{\mathrm{\Pi }}_{{}^{\mathrm{3}}\mathit{/}\mathrm{2}}\mathit{,J}\mathrm{=}{}^{\mathrm{3}}\mathit{/}\mathrm{2}$ at 79 μm, ${}^{\mathrm{2}}\mathrm{\Pi }_{{\mathrm{3}}^{}\mathit{/}\mathrm{2}}$J = ⁷/ 2−⁵/ 2 doublet at 84 μm (E_u/k_B ~ 290 K), and ${}^{\mathrm{2}}\mathrm{\Pi }_{{\mathrm{1}}^{}\mathit{/}\mathrm{2}}$J = ³/ 2−¹/ 2 doublet at 163 μm (E_u/k_B ~ 270 K), are detected in 14, 15, and 13 objects, respectively (~60−70%).

Sources without any detections of molecular lines associated with the targeted protostars are Per 4, 7, 12, and 17 (see Table A.3). In Per 13 only two weak H₂O lines at 108 and 125 μm are seen, whereas in Per 16 only a few of the lowest-J CO lines are detected. A common characteristic of this weak-line group of objects, is a low bolometric temperature (all except Per 4) and a low bolometric luminosity (see Table 1), always below 1.3 L_⊙. However, our sample also includes a few objects with similarly low values of L_bol, that show many more molecular lines (in particular Per 3, but also Per 15 and 21), so low luminosity by itself is not a criterion for weak lines. On the other hand, low bolometric temperature and high luminosity is typically connected with strong line emission (Kristensen et al. 2012; Karska et al. 2013).

3.3. Observed line ratios

Observed line ratios are calculated using the fluxes obtained from the entire 5×5 PACS maps in cases of extended emission where contamination by a nearby source and/or outflows is excluded. For point sources, a wavelength-dependent PSF correction factor is applied to fluxes obtained from the central spaxel (for details, see Sect. 2). In principle, ratios using lines that are close in wavelength could be calculated using smaller flux extraction regions, and no PSF correction would be required. However, the transition wavelengths of H₂O are not proportional to the upper energy levels, as is the case for CO, and comparisons of lines tracing similar gas have to rely on lines that lie far apart in wavelength. We explore to what extent the size of the extraction region affects the inferred line ratios. Since our aim is to understand the influence of the extended emission associated with the source(s), the ratios of lines close in wavelength are studied to avoid the confusion caused by the PSF variations.

Fluxes of the nearby H₂O 4₀₄−3₁₃ and CO 21−20 lines located at 124−125 μm are calculated first using only the central spaxel (compact region) and then using all the spaxels with detected line emission (extended region). Figure 3 illustrates that the line ratios calculated in these two regions are fully consistent, both for sources with compact and extended emission. Therefore, for the subsequent analysis the extended region is used for comparisons with models. Table A.4 shows the minimum and maximum values of the observed line ratios, their mean values, and the standard deviations for all sources with detections.

Even taking into account the uncertainties in flux extraction, the line ratios span remarkably narrow ranges of values, see Table 3 for a selection of H₂O, CO, and OH line ratios. The largest range is seen in the H₂O 2₁₂−1₀₁/CO 16−15 line ratio, which spans an order of magnitude. In all the other cases the line ratios are similar up to a factor of a few. The most similar are the OH line ratios which differ only by a factor of two, consistent with previous studies based on a large sample of low-mass YSOs in Wampfler et al. (2013). The observed similarities also imply that the line ratios do not depend on protostellar luminosity, bolometric temperature, or envelope mass (see Fig. B.1).

Our line ratios for Perseus sources are consistent with the previously reported values for other deeply-embedded protostars observed in the same way (“on source”) as tabulated in Table 3. Some differences are found for PACS observations of shock positions away from the protostar (“off source”). Most notably, the ratios using the low excitation H₂O 2₁₂−1₀₁ line at 179 μm are up to a factor of two larger than those observed in the protostellar vicinity. Such differences are not seen when more highly-excited H₂O lines are compared with each other, for example H₂O 2₂₁−1₁₀ and 4₀₄−3₁₃ lines, or with the high-J CO lines, for example CO 24−23. The ratios of two CO lines observed away from the protostar, e.g., the CO 16−15 and CO 24−23 ratios, are at the low end of the range observed toward the protostellar position.

Spectrally-resolved profiles of the H₂O 2₁₂−1₀₁ line observed with HIFI toward the protostar position reveal absorptions at source velocity removing about 10% of total line flux (e.g., Kristensen et al. 2010; Mottram et al. 2014). Our unresolved PACS observations therefore provide a lower limit to the H₂O emission in the 2₁₂−1₀₁ line. This effect, however, is too small to explain the differences in the line ratios at the “on source” and “off source” positions.

4. Analysis

Multiple molecular transitions over a wide range of excitation energies are detected, which points to the presence of hot, dense gas and can be used to constrain the signatures of shocks created as a result of outflow-envelope interaction. In particular, the line ratios of H₂O, CO, and OH are useful probes of various shock types and parameters that do not suffer from distance uncertainties.

Similarities between the spatial extent of different molecules (Sect. 3.1) coupled with similarities in the velocity-resolved line profiles among these species (Kristensen et al. 2010; Yıldız et al. 2013; San José-García et al. 2013; Mottram et al. 2014) strongly suggest that all highly-excited lines of CO and H₂O arise from the same gas. Some differences may occur for OH, which is also associated with dissociative shocks and can be affected by radiative excitation (see Sect. 5). Modeling of absolute line fluxes requires sophisticated two-dimensional (2D) physical source models for the proper treatment of the beam filling factor (Visser et al. 2012). Those models also show that UV heating alone is not sufficient to account for the high excitation lines. Hence, the focus in this analysis is on shocks. Since the absolute flux depends sensitively on the assumed emitting area, in the subsequent analysis only the line ratios are compared.

In the following sections, properties of shock models and the predicted line emission in various species are discussed (Sect. 4.1) and observations are compared with the models, using line ratios of the same species (Sect. 4.2), and different species (Sect. 4.3). Special focus will be given to C-type shocks where grids of model results are available in the literature. Observations suggest that most of the mass of hot gas is in C-type shocks toward the central protostellar positions, at least for H₂O and CO with J< 30 (Kristensen et al. 2013, and in prep.); higher-J CO, OH and [O i] transitions, on the other hand, will primarily trace J-type shocks (e.g., Wampfler et al. 2013; Kristensen et al. 2013). The excitation of OH and [O i] will be analyzed in a forthcoming paper; the J> 30 CO emission is only detected toward ~30% of all sources and so is likely unimportant for the analysis and interpretation of the data presented here. Only a limited discussion of J-type shocks is therefore presented below.

4.1. Model line emission

Models of shocks occurring in a medium with physical conditions typical for the envelopes of deeply-embedded young stellar objects provide a valuable tool for investigating shock characteristics: shock type, velocity, and the pre-shock density of (envelope) material.

Model grids have been published using a simple 1D geometry either for steady-state C and J type shocks (Hollenbach et al. 1989; Kaufman & Neufeld 1996; Flower & Pineau des Forêts 2010) or time-dependent C−J type shocks (Gusdorf et al. 2008, 2011; Flower & Pineau des Forêts 2012). The latter are non-stationary shocks, where a J-type front is embedded in a C-type shock (Chieze et al. 1998; Lesaffre et al. 2004a,b). These shocks are intermediate between pure C- and J-type shocks and have temperatures and physical extents in between the two extremes.

C-J shocks may be required for the youngest outflows with ages less than 10³ yrs, (Flower & Pineau des Forêts 2012, 2013, for the case of IRAS4B, Per 11). Here the dynamical age of the outflow is taken as an upper limit of that of the shock itself, which may be caused by a more recent impact of the wind on the envelope. The age of our sources is on the order of 10⁵ yrs (Sadavoy et al. 2014, for Class 0 sources in Perseus) and they should have been driving winds and jets for the bulk of this period, so this timescale is long enough for any shocks close to the source position to have reached steady state. While we cannot exclude that a few individual shocks have been truncated, our primary goal is to examine trends across the sample. Invoking C-J type shocks with a single truncation age as an additional free parameter is therefore not a proper approach for this study. The focus is therefore placed on comparing C-type shock results from Kaufman & Neufeld (1996, KN96 from now on) and Flower & Pineau des Forêts (2010, F+PdF10 from now on) with the observations.

All models assume the same initial atomic abundances and similarly low degrees of ionization, x_i ~ 10^-7 for C shocks. The pre-shock transverse magnetic field strength is parametrized as ${\mathit{B}}_{\mathrm{0}}\mathrm{=}\mathit{b}\mathrm{\times }\sqrt{{\mathit{n}}_{\mathrm{H}}\mathrm{\left(}{\mathrm{cm}}^{-3}\mathrm{\right)}}$μgauss, where n_H is the pre-shock number density of atomic hydrogen and b is the magnetic scaling factor, which is typically 0.1−3 in the ISM (Draine 1980). The value of b in KN96 and F+PdF10 is fixed at a value of 1.

The main difference between the two shock models is the inclusion of grains in the F+PdF10 models (Flower & Pineau des Forêts 2003). The latter models assume a standard MRN distribution of grain sizes (Mathis et al. 1977) for grain radii between 0.01 μm and 0.3 μm and a fractional abundance of the PAH in the gas phase of 10^-6 (the role of PAHs and the sizes of grains are discussed in Flower & Pineau des Forêts 2003, 2012). As the electrons are accelerated in the magnetic precursor, they attach themselves to grains thereby charging the grains and thus increasing the density of the ionized fluid significantly. This increase in density has the effect of enhancing the ion-neutral coupling (Draine 1980), thereby effectively lowering the value of b compared to the KN96 models. As a consequence, the maximum kinetic temperature is higher in the F+PdF10 models for a given shock velocity, ν. The stronger coupling between the ions and neutrals results in narrower shocks (Flower & Pineau des Forêts 2010), with shock widths scaling as ∝ b²(x_in_Hν)^-1 (Draine 1980). This proportionality does not capture the ion-neutral coupling exactly, as, for example, the grain size distribution influences the coupling (Guillet et al. 2007, 2011). The compression in C shocks also changes with the coupling since the post-shock density depends on the magnetic field, n_post ~ 0.8νn_Hb^-1 (e.g., Karska et al. 2013). The column density of emitting molecules is a function of both shock width and compression factor, and as a zeroth-order approximation the column density is $\mathit{N}\mathrm{~}{\mathit{n}}_{\mathrm{post}}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}\mathit{L}\mathrm{~}\mathit{b}{\mathit{x}}_{\mathrm{i}}^{-1}$. The ionization degree is not significantly different between the KN96 and F+PdF10 models because it is primarily set by the cosmic ray ionization rate (ζ = 5 × 10^-17 s^-1, F+PdF10), and thus the F+PdF10 models predict lower column densities than the KN96 models for a given velocity and density.

Another important difference is that the F+PdF10 models take into account that molecules frozen out onto grain mantles can be released through sputtering when the shock velocity exceeds ~15 km s^-1 (Flower & Pineau des Forêts 2010, 2012; Van Loo et al. 2013). Therefore, the gas-phase column densities of molecules locked up in ices increase above this threshold shock velocity with respect to the KN96 models, an effect which applies to both CO and H₂O. Furthermore, H₂O forms more abundantly in the post-shock gas of F+PdF10 models, because H₂ reformation is included, unlike in the KN96 models (Flower & Pineau des Forêts 2010).

Molecular emission is tabulated by KN96 for a wide range of shock velocities, from ν= 5 to 45 km s^-1 in steps of 5 km s^-1, and a wide range of pre-shock densities n_H, from 10⁴ to 10^6.5 cm^-3 in steps of 10^0.5 cm^-3. The F+PdF10 grid is more limited in size, providing line intensities for only two values of pre-shock densities, namely 10⁴ and 10⁵ cm^-3, and a comparable range of shock velocities, but calculated in steps of 10 km s^-1. Calculations are provided for CO transitions from J = 1−0 to J = 60−59 in KN96 and only up to J = 20−19 in F+PdF10. The two sets of models use different collisional rate coefficients to calculate the CO excitation. F+PdF10 show line intensities for many more H₂O transitions (in total ~120 lines in the PACS range, see Sect. 2) than in the older KN96 grid (18 lines in the same range), which was intended for comparisons with the Submillimeter Wave Astronomy Satellite (SWAS, Melnick et al. 2000) and ISO data. KN96 use collisional excitation rates for H₂O from Green et al. (1993) and F+PdF10 from Faure et al. (2007). Line intensities for OH are only computed by KN96, assuming only collisional excitation and using the oxygen chemical network of Wagner & Graff (1987). The reaction rate coefficients in that network are within a factor of 2 of the newer values by Baulch et al. (1992) and tabulated in the UMIST database (www.udfa.net, McElroy et al. 2013), see also a discussion in van Dishoeck et al. (2013).

We also use CO fluxes extracted from the grid of models presented by Kristensen et al. (2007), since high-J CO lines are missing in F+PdF10. This grid is based on the shock model presented in Flower & Pineau des Forêts (2003) and covers densities from 10⁴−10⁷ cm^-3 and velocities from 10−50 km s^-1 (denoted as F+PdF* from now on). The main difference compared to the results from F+PdF10 is that CO and H₂O level populations are not calculated explicitly through the shock; instead analytical cooling functions are used to estimate the relevant line cooling and afterwards are line fluxes extracted (Flower & Gusdorf 2009). Models with b = 1 are used. The CO line fluxes presented here are computed using the 3D non-LTE radiative transfer code LIME (Brinch & Hogerheijde 2010), for levels up to J = 80−79. The CO collisional rate coefficients from Yang et al. (2010) extended by Neufeld (2012) are used.

In the following sections, the model fluxes of selected CO, H₂O, and OH lines are discussed for a range of shock velocities and three values of pre-shock densities: 10⁴, 10⁵, and 10⁶ cm^-3. We note that the post-shock densities traced by observations are related to the pre-shock densities via the compression factor dependent on the shock velocity and magnetic field. In C shocks, the compression factor is about 10 (e.g., Karska et al. 2013).

Figure 4 compares model fluxes of various CO, OH, and H₂O lines from the KN96 models (panels a and d) and, for a few selected CO and H₂O lines, compares the results with the F+PdF10 or F+PdF* models (panels b, c, e, and f).

4.1.1. CO

The KN96 model line fluxes for CO 16−15, CO 21−20, and CO 29−28 are shown in panel a of Fig. 4. The upper energy levels of these transitions lie at 750 K, 1280 K, and 2900 K, respectively, while with increasing shock velocity, the peak C-shock temperature increases from about 400 K to 3200 K (for 10 to 40 km s^-1) and is only weakly dependent on the assumed density (see Fig. 3 of KN96). Therefore, the CO 16−15 line is already excited at relatively low shock velocities (ν~ 10 km s^-1, for n_H = 10⁴ cm^-3), whereas the higher-J levels become populated at higher velocities. At a given shock velocity, emission from the CO 16−15 line is the strongest because of its lower critical density, n_cr ~ 9 × 10⁵ cm^-3 at T = 1000 K (Neufeld 2012). This situation only changes for the highest pre-shock densities, when the line becomes thermalized and the cooling in other lines dominates.

The CO 16−15 flux from the F+PdF10 C-type shock models is comparable to the KN96 flux for the ν~ 10 km s^-1 shock, but increases less rapidly with shock velocity despite the sputtering from grain mantles (panel b of Fig. 4). Because of the lower magnetic scaling factor b in the models with grains (see Sect. 4.1), it is expected that the column density and the corresponding line fluxes are lower in the F+PdF10 models. For slow shocks the higher temperatures in the latter models compensate for the smaller column density resulting in a similar CO 16−15 flux.

In J-type shocks, the peak temperatures of the post-shock gas are 1400, 5500, 12 000, and 22 000 K for the shock velocities of 10, 20, 30, and 40 km s^-1 respectively (Neufeld & Dalgarno 1989; Kaufman & Neufeld 1996). For 10−20 km s^-1 shocks, these high temperatures more easily excite the CO 16−15 line with respect to C-type shock emission. For shock velocities above 20 km s^-1, such high temperatures can lead to the collisional dissociation of H₂ and subsequent destruction of CO and H₂O molecules, resulting in the decrease in CO fluxes. This effect requires high densities and therefore the CO flux decrease is particularly strong for the pre-shock densities 10⁵ cm^-3.

4.1.2. H₂O

The H₂O fluxes show a strong increase with shock velocities above ν~ 10−15 km s^-1 in both models, especially at low pre-shock densities (panels d-f of Fig. 4). At this velocity, the gas temperature exceeds ~400 K, so the high-temperature route of H₂O formation becomes efficient (Elitzur & Watson 1978; Elitzur & de Jong 1978; Bergin et al. 1998; KN96), which quickly transfers all gas-phase oxygen into H₂O via reactions with H₂ (KN96).

In contrast to CO, the upper level energies of the observed H₂O lines are low and cover a narrow range of values, E_up ~ 200−600 K. As a result, the effect of peak gas temperature on the H₂O excitation is less pronounced (Fig. 3 of KN96) and after the initial increase with shock velocity, the H₂O fluxes in the KN96 models stay constant for all lines. For high pre-shock densities, the higher lying levels are more easily excited and, as a consequence, the fluxes of the H₂O 6₁₆−5₀₅ line become larger than those of the H₂O 2₁₂−1₀₁ line. The critical densities of these transitions are about two orders of magnitude higher than for the CO 16−15 line and the levels are still sub-thermally excited at densities of 10⁶–10⁷ cm^-3 (and effectively optically thin, Mottram et al. 2014).

The H₂O 2₁₂−1₀₁ fluxes in the C-type F+PdF10 models are remarkably similar to those found by KN96 (panel e of Fig. 4), while the H₂O 4₀₄−3₁₃ fluxes are lower by a factor of a few over the full range of shock velocities and pre-shock densities in the F+PdF10 models (panel f of Fig. 4). Lower H₂O fluxes are expected because of the smaller column of H₂O in the models with grains (Sect. 4.1). For the lower-J lines, the various factors (lower column density through a shock but inclusion of ice sputtering and H₂ reformation) apparently conspire to give similar fluxes as for KN96.

In J-type shocks, the fluxes of H₂O 2₁₂−1₀₁ and 4₀₄−3₁₃ lines increase sharply for shock velocities 10−20 km s^-1 (panels e and f of Fig. 4). The increase is not as steep at 30 km s^-1 shocks for n_H = 10⁵ cm^-3, when the collisional dissociation of H₂ and subsequent destruction of H₂O molecules occurs (fluxes for larger shock velocities are not computed in F+PdF10 and hence not shown). Below 30 km s^-1, line fluxes from J-type shocks are comparable to those from C-type shock predictions except for the H₂O 2₁₂−1₀₁ fluxes at high pre-shock densities, which are an order of magnitude higher with respect to the C-type shock predictions. The difference could be due to smaller opacities for the low-excitation H₂O line in the J shocks.

4.1.3. OH

The fluxes of the ${}^{\mathrm{2}}\mathrm{\Pi }_{{\mathrm{3}}^{}\mathit{/}\mathrm{2}}$J = ⁷/ 2−⁵/ 2 doublet at 84 μm (E_u/k_B ~ 290 K) calculated with the KN96 models are shown with the H₂O lines in panel d of Fig. 4. Not much variation is seen as a function of shock velocity, in particular beyond the initial increase from 10 to 15 km s^-1, needed to drive oxygen to OH by the reaction with H₂. At about 15 km s^-1, the temperature is high enough to start further reactions with H₂ leading to H₂O production. The trend with increasing pre-shock density is more apparent, with OH fluxes increasing by two orders of magnitude between the 10⁴ cm^-3 to 10⁶ cm^-3, as the density becomes closer to the critical density of the transition.

4.2. Models versus observations – line ratios of the same species

Comparison of observed and modeled line ratios of different pairs of CO, H₂O, and OH transitions is shown in Fig. 5. The line ratios are a useful probe of molecular excitation and therefore can be used to test whether the excitation in the models is reproduced correctly, which in turn depends on density and temperature, and thus shock velocity.

4.2.1. CO line ratios

Given the universal shape of the CO ladders observed toward deeply-embedded protostars (see Sect. 1 and the discussion on the origin of CO ladders in Sect. 5), three pairs of CO lines are compared with models: (i) CO 16−15 and 21−20 line ratio, corresponding to the “warm”, 300 K component (panel a); (ii) CO 16−15 and 29−28 line ratio, combining transitions located in the “warm” and “hot” (>700 K) components (panel b); (iii) CO 24−23 and 29−28, both tracing the “hot” component (panel c). All ratios are consistent with the C-type shock models from both KN96 and F+PdF10 for pre-shock densities above n_H = 10⁴ cm^-3. For the CO 16−15/21−20 ratio, a pre-shock density of n_H = 10⁵ cm^-3 and shock velocities of 20−30 km s^-1 best fit the observations. Shock velocities above ~ 25 km s^-1 are needed to reproduce the observations of the other two ratios at the same pre-shock density. Alternatively, higher pre-shock densities with velocities below 30 km s^-1 are also possible.

The KN96 C-shock CO line ratios for lower-to-higher-J transitions (panel b of Fig. 5) decrease with velocity, because of the increase in peak temperature that allows excitation of the higher-J CO transitions. The effect is strongest at low pre-shock densities (see Sect. 4.1.1) and for the sets of transitions with the largest span in J numbers. The CO 16−15/29−28 line ratio (ΔJ_up = 13) decreases by almost three orders of magnitude between shock velocities of 10 and 40 km s^-1 over the range of pre-shock densities. In contrast, the CO 16−15/21−20 and CO 24−23/29−28 line ratios show drops of about one order of magnitude with increasing velocity (panel a and c of Fig. 5). These model trends explain why the observed CO line ratios are good diagnostics of shock velocity.

In absolute terms, the line ratios calculated for a given velocity are inversely proportional to the pre-shock density. The largest ratios obtained for n_H = 10⁴ cm^-3 result from the fact that the higher-J levels are not yet populated at low shock-velocities, while the lower-J transitions reach LTE at high shock-velocities and do not show an increase of flux with velocity. This effect is less prominent at higher pre-shock densities, where the higher-J lines are more easily excited at low shock velocities.

The F+PdF* CO line ratios, extending the Flower & Pineau des Forêts (2003) grid to higher-J CO lines, are almost identical to the KN96 predictions for pre-shock densities n_H = 10⁴ cm^-3. For higher densities, the low-velocity C-shock models from F+PdF* are systematically lower than the KN96 models, up to almost an order of magnitude for 10−15 km s^-1 shocks at n_H = 10⁶ cm^-3. Therefore, the pre-shock density is not as well constrained solely by CO lines.

For densities of 10⁵ cm^-3, shock velocities of 20−30 km s^-1 best reproduce the ratios only using transitions from the “warm” component, while shock velocities above 25 km s^-1 match the ratios using the transitions from the “hot” component. Velocities of that order are observed in CO J = 16−15 HIFI line profiles (Kristensen et al. 2013, and in prep.), but higher-J CO lines dominated by the hot component have not been obtained with sufficient velocity resolution.

4.2.2. H₂O line ratios

Two ratios of observed H₂O lines are compared with the C- and J-type shock models: (i) the ratio of the low excitation H₂O 2₁₂−1₀₁ and moderate excitation 4₀₄−3₁₃ lines (panel d of Fig. 5) and (ii) the ratio of the highly-excited H₂O 6₁₆−5₀₅ and 4₀₄−3₁₃ lines (panel e). Similar to the CO ratios, C-type shocks with pre-shock densities of 10⁵ cm^-3 reproduce the observations well. Based on the observations of ratio (i), C shocks with a somewhat larger (F+PdF10) or smaller (KN96) pre-shock density are also possible for a broad range of shock velocities. On the other hand, no agreement with the J-type shocks is found for this low-excitation line ratio. Observations of ratio (ii) indicate a similar density range as ratio (i) for the KN96 models, but extend to 10⁴ cm^-3 for the F+PdF10 models, with agreement found for both C- and J-type.

The model trends can be understood as follows. For 10−20 km s^-1 shocks, increasing temperature in the C shock models from KN96 allows excitation of high-lying H₂O lines and causes the H₂O 2₁₂−1₀₁/4₀₄−3₁₃ line ratio to decrease and the H₂O 6₁₆−5₀₅/4₀₄−3₁₃ line ratio to increase. At higher shock velocities, the former ratio shows almost no dependence on shock velocity, while a gradual increase is seen in the ratio using two highly-excited lines in the KN96 models. At high pre-shock densities (n_H = 10⁶ cm^-3), the upper level transitions are more easily excited and so the changes are even smaller.

The H₂O 2₁₂−1₀₁/4₀₄−3₁₃ line ratios calculated using the C-type shock models from F+PdF10 are a factor of a few larger than the corresponding ratios from the KN96 models (see the discussion of absolute line fluxes in Sect. 4.1.2). As a result, when compared to observations, the F+PdF10 models require pre-shock densities of at least n_H = 10⁵ cm^-3, while the KN96 models suggest densities that are lower by a factor of a few. Overall, the best fit to both the CO and H₂O line ratios is for pre-shock densities around 10⁵ cm^-3.

4.2.3. OH line ratios

Comparison of the observed OH 84 and 79 μm line ratio with the KN96 C-type models (panel f of Fig. 5) indicates an order of magnitude higher pre-shock densities, n_H = 10⁶ cm^-3, with respect to those found using the CO and H₂O ratios. However, the KN96 models do not include any far-infrared radiation, which affects the excitation of the OH lines, in particular the 79 μm (Wampfler et al. 2010, 2013). Additionally, part of OH most likely originates in a J-type shock, influencing our comparison (Wampfler et al. 2010; Benedettini et al. 2012; Karska et al. 2013; Kristensen et al. 2013).

Similar to the absolute fluxes of the 84 μm doublet discussed in Sect. 4.1.3, not much variation in the ratio is seen with shock velocity. The ratio increases by a factor of about two between the lowest and highest pre-shock densities.

4.3. Models and observations – line ratios of different species

Figure 6 compares observed line ratios of various H₂O and CO transitions with the C and J-type shock models. The line ratios of different species are sensitive both to the molecular excitation and their relative abundances.

4.3.1. Ratios of H₂O and CO

Comparison of observations to the KN96 and F+PdF10 models shows that the C-type shocks at pre-shock density n_H = 10⁵ cm^-3, which best reproduces the line ratios of same species, fail to reproduce the observed line ratios of different species (Fig. 6). There are only a few cases where the observations seem to agree with the models at all. For n_H = 10⁵ cm^-3, a few H₂O/CO line ratios fit at low velocities (<20 km s^-1; panels a, c and e) but this does not hold for all ratios. Moreover, such low shock velocities have been excluded in the previous section. Higher densities, n_H = 10⁶ cm^-3, are needed to reconcile the observations of the H₂O 2₁₂−1₀₁/CO 29−28 line ratio. Observations of all the other ratios, using more highly-excited H₂O lines, are well below the model predictions.

The patterns seen in the panels in Fig. 6 can be understood as follows. The KN96 C type shock models show an initial rise in the H₂O 2₁₂−1₀₁ and CO 16−15 line ratios from 10 to 15 km s^-1 shocks (panel a), as the temperature reaches the 400 K and enables efficient H₂O formation. Beyond this velocity, the line ratios show no variations with velocity. The decrease in this line ratio for higher densities, from 10 at 10⁴ cm^-3 to 0.1 at n_H = 10⁶ cm^-3 is due to the larger increase of the column of the population in the J_u = 16 level with density compared to the increase in the H₂O 2₁₂ level (see Fig. 4 above).

Line ratios of H₂O 2₁₂−1₀₁ and higher-J CO lines (e.g., 29−28, panel b) show more variation with velocity. A strong decrease by about an order of magnitude and up to two orders of magnitude are seen for the ratios with CO 24−23 and CO 29−28, respectively (the ratio with CO 24−23 is not shown here). These lines, as discussed in Sects. 4.1.1 and 4.2.1, are more sensitive than H₂O to the increase in the maximum temperature attained in the shock that scales with shock velocities and therefore their flux is quickly rising for higher velocities (Fig. 4). The decrease is steeper for models with low pre-shock densities, since n_H ~ 10⁶ cm^-3 allows excitation of high-J CO lines at lower temperatures. At this density, the H₂O/CO line ratios are the lowest and equal about unity.

Because of the lower CO 16−15 fluxes in the C-type shock models from F+PdF10 and similar H₂O 2₁₂−1₀₁ fluxes (Fig. 4), the H₂O-to-CO ratios are generally larger than in the KN96 models. The exceptions are the ratios with higher-J CO which are more easily excited, especially at low shock velocities, in the hotter C-type shocks from F+PdF*.

For the same reason, the increasing ratios seen in the J shock models are caused by the sharp decrease in CO 16−15 flux for shock velocity ν = 30 km s^-1, rather than the change in the H₂O lines. At such high-velocities for J shocks, a significant amount of CO can be destroyed by reactions with hydrogen atoms (Flower & Pineau des Forêts 2010; Suutarinen et al. 2014). Since the activation barrier for the reaction of H₂O with H is about 10⁴ K, the destruction of H₂O does not occur until higher velocities.

Similar trends to the line ratios with H₂O 2₁₂−1₀₁ are seen when more highly-excited H₂O lines are used (panels c-f of Fig. 6), supporting the interpretations that variations are due to differences in CO rather than H₂O lines.

4.3.2. Ratios of CO and H₂O with OH

Figure 7 shows line ratios of CO or H₂O with the most commonly detected OH doublet at 84 μm. The ratios are calculated for three values of pre-shock densities (10⁴, 10⁵, and 10⁶ cm^-3) using exclusively the KN96 models, because the F+PdF10 grid does not present OH fluxes.

In general, the observed CO/OH, and H₂O/OH ratios are similar for all sources but much lower than those predicted by the models assuming that a significant fraction of the OH comes from the same shock as CO and H₂O (see Sect. 5.1). The only exception is the CO 24−23/OH 84 μm ratio where models and observations agree for densities 10⁴−10⁵ cm^-3 and shock velocities below 20 km s^-1. For any other set of lines discussed here, the observations do not agree with these or any other models.

As discussed in previous sections, the trends with shock velocity are determined mostly by the changes in the CO lines, rather than the OH itself, as seen in Fig. 4 (panel d). At shock velocities below 20 km s^-1 the OH model flux exceeds that of CO owing to the abundance effect: not all OH has been transferred to H₂O yet at low temperatures. Because of the lower critical densities of the CO lines (n_cr ~ 10⁶−10⁷ cm^-3) compared with the OH line (n_cr ~ 10⁹ cm^-3), the lines for various pre-shock densities often cross and change the order in the upper panels of Fig. 7. The corresponding trends in the H₂O/OH line ratios are similar to those of CO/OH, except that the variations with shock velocity are smaller and the critical densities are more similar.

5. Discussion

5.1. Shock parameters and physical conditions

Spectrally resolved HIFI observations of the CO 10−9 and 16−15 line profiles (Kristensen et al. 2013 and in prep.; Yildiz et al. 2013) as well as various H₂O transitions (Kristensen et al. 2012; Mottram et al. 2014) reveal at least two different kinematic shock components: non-dissociative C-type shocks in a thin layer along the cavity walls (so-called “cavity shocks”) and J-shocks at the base of the outflow (also called “spot shocks”), both caused by interaction of the wind with the envelope. Both shocks are different from the much cooler entrained outflow gas that is observed in the low-J CO line profiles (Yıldız et al. 2013).

One possible physical explanation for our observed lack of variation is that although the outflow structure depends on the mass entrainment efficiency and the amount of mass available to entrain (the envelope mass), the wind causing the shocks does not depend on these parameters. Instead the cavity shock caused by the wind impinging on the inner envelope depends on the shock velocity and the density of the inner envelope (Kristensen et al. 2013; Mottram et al. 2014). Thus, the lack of significant variation in the line ratios suggests that the shock velocities by the oblique impact of the wind are always around 20−30 km s^-1.

In Sect. 4 the observed emission was compared primarily to models of C-type shock emission. Although J-type shocks play a role on small spatial scales in low-mass protostars (Kristensen et al. 2012, 2013; Mottram et al. 2014) their contribution to CO emission originating in levels with J_up ≲ 30 is typically less than ~50%. Since higher-J CO emission is only detected toward 30% of the sources, the J-type shock component is ignored for CO. For the case of H₂O, spectrally resolved line profiles observed with HIFI reveal that the profiles do not change significantly with excitation up to E_up = 250 K (Mottram et al. 2014); J-type shock components typically contribute <10% of the emission. It is unclear if the trend of line profiles not changing with excitation continues to higher upper-level energies, in particular all the way up to E_up = 1070 K (J = 8₁₈−7₀₇ at 63.32 μm). OH and [O i], on the other hand, almost certainly trace dissociative J-type shocks (e.g., van Kempen et al. 2010; Wampfler et al. 2013) but a full analysis of their emission will be presented in a forthcoming paper. Thus, in the following the focus remains on comparing emission to models of C-type shocks.

Figure 8 summarizes the different line ratios as a function of pre-shock density discussed in the previous sections. General agreement is found between the observations and models when line ratios of different transitions of the same species are used (see top row for H₂O, CO, and OH examples), indicating that the excitation of individual species is reproduced well by the models. The H₂O line ratios are a sensitive tracer of the pre-shock gas density since they vary less with shock velocity than those of CO. The C shock models from KN96 with pre-shock gas densities in the range of 10⁴−10⁵ cm^-3 are a best match to the observed ratios, consistent with values of 10⁵ cm^-3 from the C shock models of F+PdF10. For the considered range of shock velocities, the compression factor in those shocks, defined as the ratio of the post-shock and pre-shock gas densities, varies from about 10 to 30 (Neufeld & Dalgarno 1989; Draine & McKee 1993; Karska et al. 2013). The resulting values of post-shock densities, traced by the observed molecules, are therefore expected to be ≥ 10⁵−10⁶ cm^-3.

The CO line ratios, on the other hand, are not only sensitive to density, but also to the shock velocities, because of their connection to the peak temperature attained in the shock. In the pre-shock density range of ≥ 10⁴−10⁵ cm^-3, indicated by the H₂O line ratios, shocks with velocities above 20 km s^-1 best agree with the CO observations. Within this range of densities, the predictions from both the KN96 and F+PdF10 C shock models show a very good agreement with each other.

The ratio of two OH lines from the KN96 models compared with the observations suggest higher pre-shock densities above 10⁵ cm^-3, but this ratio may be affected by infrared pumping (Wampfler et al. 2013). In addition, some OH emission traces (dissociative) J-shocks, based on its spatial connection and flux correlations to [O i] emission (Wampfler et al. 2010, 2013; Karska et al. 2013). The single spectrally-resolved OH spectrum towards Ser SMM1 (Fig. 3, Kristensen et al. 2013) suggests that the contribution of the dissociative and non-dissociative shocks is comparable. Thus, observed CO/OH and H₂O/OH line ratios are only affected at the factor ~2 level and the discrepancy in the bottom row of Fig. 8 remains.

Overall, the observed CO and H₂O line ratios are best fit with C-shock models with pre-shock densities of ~10⁵ cm^-3 and velocities ≳20 km s^-1, with higher velocities needed for the excitation of the highest-J CO lines.

The shock conditions inferred here can be compared to the temperatures and densities found from single-point non-LTE excitation and radiative-transfer models, e.g., radex (van der Tak et al. 2007) and from non-LTE radiative transfer analysis of line intensity ratios (Kristensen et al. 2013; Mottram et al., and in prep.) toward various sources. Typically, densities ≳ 10⁶ cm^-3 and temperatures of ~ 300 K and ≳ 700 K are required to account for the line emission (Herczeg et al. 2012; Goicoechea et al. 2012; Santangelo et al. 2012, 2013; Vasta et al. 2012; Karska et al. 2013). Within this range of densities, the predictions from both the KN96 and F+PdF10 C-shock models reproduce CO observations. However, the disconnect between predicted pre-shock conditions required to reproduce H₂O and CO is puzzling (see below).

A small number of individual sources have been compared directly to shock models (Lee et al. 2013; Dionatos et al. 2013) and the conclusions are similar to what is reported here: pre-shock conditions of typically 10⁴−10⁵ cm^-3, and emission originating in C-type shocks. None of the sources analyzed previously were therefore special or atypical, rather these shock conditions appear to exist toward every embedded protostar.

At shock positions away from the protostar, dissociative or non-dissociative J-type shocks at the same pre-shock densities are typically invoked to explain the FIR line emission (Benedettini et al. 2012; Santangelo et al. 2012; Busquet et al. 2014). Differences between the protostar position and the distant shock positions are revealed primarily by our line ratios using the low-excitation H₂O 2₁₂−1₀₁ (Table 3) and can be ascribed to the differences in the filling factors and column densities between the immediate surrounding of the protostar and the more distant shock positions (Mottram et al. 2014).

5.2. Abundances and need for UV radiation

In contrast with the ratios of two H₂O or CO lines, the ratios calculated using different species do not agree with the shock models (Fig. 8, bottom row). The ratios of H₂O-to-CO lines are overproduced by the C shock models from both the KN96 and F+PdF10 grids by at least an order of magnitude, irrespective of the assumed shock velocity. Although there are a few exceptions (e.g., the ratio of H₂O 2₁₂−1₀₁ and CO 16−15), the majority of the investigated sets of H₂O and CO lines follow the same trend. Observations agree only with slow, <20 km s^-1, J shock models, but as shown above, those models do not seem to reproduce the excitation properly (Fig. 5).

The discrepancy between the models and observations is even larger in the case of the H₂O-to-OH line ratios, as illustrated in Fig. 8. The two orders of magnitude disagreement with the C shock models cannot be accounted by any excitation effects for any realistic shock parameters. Additional comparison to J shock models is not possible because of a lack of OH predictions for J shocks in the F+PdF10 models.

The CO-to-OH ratios are overproduced by about an order of magnitude in the C shock models, similar to the H₂O-to-CO ratios. The agreement improves for fast (ν= 40 km s^-1) shocks in high density pre-shock medium (~ 10^6.5 cm^-3), but those parameters are not consistent with the line ratios from the same species.

An additional test of the disagreement between models and observations is provided by calculating the fraction of each species with regard to the sum of CO, H₂O and OH emission. For that purpose, only the strongest lines observed in our program are used. As seen in Table 4, the observed percentage (median) of H₂O is about 30% and OH is about 25%. In contrast, KN96 models predict typically 70−90% of flux in the chosen H₂O lines and only up to 2% in the OH lines.

The only possible way to reconcile the models with the observations, after concluding that the excitation is treated properly in the models, is to reconsider the assumed abundances. The fact that the H₂O-to-CO and CO-to-OH ratios are simultaneously overestimated suggests a problem with the abundances of H₂O and OH, rather than that of CO. The scenario with the overestimated H₂O abundances and underestimated OH abundances would translate into a too large H₂O-to-OH abundance ratio in the models. A possible and likely solution is photodissocation of H₂O to OH and subsequently to atomic oxygen. As noted above, some OH also comes from the dissociative shock seen in [O i].

A significant shortcoming of all these shock models lies in their inability to account for grain-grain interactions, which has been shown to significantly alter the structure of the shocks propagating in dense media (n_H> 10⁵ cm^-3, Guillet et al. 2007, 2009, 2011). These grain-grain interactions mostly consist of coagulation, vaporization, and shattering effects affecting the grains. Their inclusion in shock models necessitates a sophisticated treatment of the grains, especially following their charge and size distribution (Guillet et al. 2007). Most remarkably, such interactions eventually result in the creation of small grain fragments in large numbers, which increases the total dust grain surface area and thereby changes the coupling between the neutral and the charged fluids within the shock layer. The net effect is that the shock layer is significantly hotter and thinner (Guillet et al. 2011), which in turn affects the chemistry and emission of molecules (Guillet et al. 2009).

Unfortunately, at the moment these models are computationally expensive and are not well-suited for a grid analysis; moreover, the solutions found by Guillet et al. do not converge for preshock densities of 10⁶ cm^-3 or higher. A recent study by Anderl et al. (2013) shows that it is possible to approximate these effects in a computationally efficient way, and subsequently evaluated the line intensities of CO, OH, and H₂O on a small grid of models. When including grain-grain interactions, CO lines were found to emit significantly less than in “simpler” models, while a smaller decrease was found for H₂O lines, and a small increase for OH. These trends probably still need to be systematically investigated on larger grids of models before they can be applied to our present comparison efforts.

Regardless of the effect of grain-grain interactions, Snell et al. (2005) invoked several scenarios to reconcile high absolute H₂O fluxes with the shock models for the case of supernova remnants. These include (i) the high ratio of atomic to molecular hydrogen, which drives H₂O back to OH and O; (ii) freeze-out of H₂O in the post-shock gas; (iii) freeze-out of H₂O in the pre-shock gas, and (iv) photodissociation of H₂O in the pre- and post-shock gas. Because of the high activation barrier of the H₂O + H → OH + H₂ reaction (~ 10⁴ K), the first scenario is not viable. The freeze-out in the post-shock gas (ii) is not effective in the low density regions considered in Snell et al. (2005), but can play a role in the vicinity of protostars, where densities above ~ 10⁶ cm^-3 are found (e.g., Kristensen et al. 2012, this work). However, this mechanism alone would not explain the bright OH and high-J H₂O lines seen toward many deeply-embedded sources (e.g., Karska et al. 2013; Wampfler et al. 2013). A similar problem is related to the freeze-out in the pre-shock gas (iii), which decreases the amount of e.g., O, OH, and H₂O in the gas phase for shock velocities below 15 km s^-1.

Therefore, the most likely reason for the overproduction of H₂O in the current generation of shock models, at the expense of OH, is the omission of the effects of ultraviolet irradiation (scenario iv) of the shocked material. The presence of UV radiation is directly seen in Ly-α emission both in the outflow-envelope shocks (Curiel et al. 1995; Walter et al. 2003) and at the protostar position (Valenti et al. 2000; Yang et al. 2012). Additionally, UV radiation on scales of a few 1000 AU has been inferred from the narrow profiles of ¹³CO 6−5 observed from the ground toward a few low-mass protostars (Spaans et al. 1995; van Kempen et al. 2009; Yıldız et al. 2012). H₂O can be photodissociated into OH over a broad range of far-UV wavelengths, including by Ly-α, and this would provide an explanation for the disagreement between our observations and the models. Photodissociation of CO is less likely, given the fact that it cannot be dissociated by Ly-α and only by very hard UV photons with wavelengths <1000 Å. The lack of CO photodissociation is consistent with weak [C i] and [C ii] emission observed toward low-mass YSOs (Yıldız et al. 2012; Goicoechea et al. 2012; Karska et al. 2013). At the positions away from the protostars, on the other hand, the bow-shocks at the tip of the protostellar jets can produce significant emission in the [C i] (van Kempen et al. 2009). Therefore, it is unlikely that lower line excitation at those positions is due to the weaker UV. The differences seen in the resolved line profiles (e.g., Santangelo et al. 2012; Vasta et al. 2012; Mottram et al., in prep.) indicate that the lower column densities involved are the more likely reason for differences in the excitation.

Visser et al. (2012) proposed a scenario in which the lower-lying CO transitions observed with PACS (14 <J_u< 23) originate in UV-heated gas and higher-J transitions (J_u> 24) in shocked material. All water emission would be associated with the same shocks as responsible for the higher-J CO emission with less than 1% of the water emission coming from the PDR layer. Although not modeled explicitly, the UV irradiation from the star-disk boundary impinging on the shocks naturally accounts for the lower H₂O abundance and exceeds by at least two orders of magnitude the H₂O destruction rate by He⁺ and H$\begin{array}{c}\mathrm{+}\\ \mathrm{3}\end{array}$, assuming a normal interstellar radiation field (G₀ = 1). The authors predict that while the dynamics of the hot layers where both shocks and UV irradiation play a role will be dominated by the shocks, only the UV photons penetrate farther into the envelope, where the dynamics would resemble the quiescent envelope. The lower-temperature UV-heated gas has indeed been observed to be quiescent on the spatial scales of the outflow cavity through observations of medium-J¹³CO lines (Yıldız et al. 2012, 2014).

Flower & Pineau des Forêts (2013) proposed a model where all emission originates in a non-stationary shock wave, where a J-type shock is embedded in a C-type shock. Without a detailed modeling of individual sources based on different source parameters it is not possible to rule out any of these solutions. However, the trends reported here suggest that it is possible to find a pure shock solution, in agreement with Flower & Pineau des Forêts (2013), as long as UV photons are incorporated to provide dissociation of H₂O. Complementary observations, preferably at higher angular resolution, are required to break the solution degeneracy and determine the relative role the shocks and UV photons play on the spatial scales of the thickness of the cavity wall. Models whose results depend sensitively on a single parameter such as time, are ruled out because the observed line ratios are so similar across sources.

6. Conclusions

We have compared the line ratios of the main molecular cooling lines detected in 22 low-mass protostars using Herschel/PACS with publicly available one-dimensional shock models. Our conclusions are the following:

• Line ratios of various species and transitions are remarkably similar for all observed sources. No correlation is found with source physical parameters.

• Line ratios observed toward the protostellar position are consistent with the values reported for the positions away from the protostar, except for some ratios involving the low-excitation H₂O 2₁₂−1₀₁ line. Coupled with the larger absolute fluxes of highly-excited H₂O and CO lines at the protostellar positions, this indicates that lines at distant off-source shock positions are less excited.

• General agreement is found between the observed line ratios of the same species (H₂O, CO, and OH) and the C shock models from Kaufman & Neufeld (1996) and Flower & Pineau des Forêts (2010). Ratios of H₂O are particularly good tracers of the density of the ambient material and indicate pre-shock densities on the order of ≥ 10⁵ cm^-3 and thus post-shock densities on the order of 10⁶ cm^-3. Ratios of CO lines are more sensitive to the shock velocities and, for the derived range of pre-shock densities, indicate shock velocities above 20 km s^-1.

• Ratios of CO lines located in the “warm” component of CO ladders (with T_rot ~ 300 K) are reproduced with shock velocities of 20−30 km s^-1 and pre-shock densities of 10⁵ cm^-3. The CO ratios using the lines from the “hot” component (T_rot ≳ 700 K) are better reproduced by models with shock velocities above 25 km s^-1.

• A lack of agreement is found between models and the observed line ratios of different species. The H₂O-to-CO, H₂O-to-OH, and CO-to-OH line ratios are all overproduced by the models by 1−2 orders of magnitude for the majority of the considered sets of transitions.

• Since the observed molecular excitation is properly reproduced in the C shock models, the most likely reason for disagreement with observations is the abundances in the shock models, which are too high in case of H₂O and too low in case of OH. Invoking UV irradiation of the shocked material, together with a dissociative J shock contribution to OH and [O i], would lower the H₂O abundance and reconcile the models and observations.

New UV-irradiated shock models will allow us to constrain the UV field needed to reconcile the shock models with observations (M. Kaufman, priv. comm.) Those models should also account for the grain-grain processing, which affects significantly the shocks structure at densities ~ 10⁶ cm^-3 (Guillet et al. 2011). The effects of shock irradiation as a function of the distance from a protostar will help to understand the differences in the observed spectrally-resolved lines from HIFI at “on source” and distant shock-spot positions.

Online material

Appendix A: Supplementary material

Table A.1 provides molecular/atomic information about the lines observed in the WILL program.

Table A.2 shows the observing log of PACS observations including observations identifications (OBSID), observation day (OD), date of observation, total integration time, and pointed coordinates (RA, Dec).

Table A.3 informs about which lines are detected toward the Perseus sources. The full list of line fluxes for all WILL sources including Perseus will be tabulated in the forthcoming paper (Karska et al., in prep.).

Table A.4 summarizes observed and modeled line ratios used in the Analysis section.

Figures A.1 and A.2 show line and continuum maps around 179.5 μm for all the Perseus sources in the WILL program.

Figure A.3 show maps in the H₂O 4₂₃−3₁₂ line at 78.74 μm, OH 84.6 μm, and CO 29−28 at 90.16 μm for Per1, Per5, Per9, and Per20, all of which show bright line emission and centrally peaked continuum. The lines are chosen to be located close in the wavelength so that the variations in the PSF does not introduce significant changes in the emission extent.

Appendix B: Correlations with source parameters

Figure B.1 shows selected H₂O-to-CO line ratios as a function of source physical parameters (bolometric luminosity, temperature, and envelope mass). Lack of correlation is seen in all cases.

---

Acknowledgments

Herschel is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA. A.K. acknowledges support from the Polish National Science Center grant 2013/11/N/ST9/00400. 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). N.J.E. was supported by NASA through an award issued by the Jet Propulsion Laboratory, California Institute of Technology.

References

All Tables

All Figures