© ESO, 2013

1. Introduction

Stars form in collapsing dense molecular cores deep inside interstellar clouds (see reviews by di Francesco et al. 2007; Bergin & Tafalla 2007; Lada 1999). Star formation is associated with many physical phenomena that occur simultaneously: infall from the envelope, action of jets and winds resulting in shocks, outflows sweeping up surrounding material, and UV heating of outflow cavity walls (Shu et al. 1987; Spaans et al. 1995; Bachiller & Tafalla 1999). In the earliest phases of star-formation (Class 0 and I objects; André et al. 1993, 2000), the interaction between the jet, wind, and the dense envelope is particularly strong and produces spectacular outflows (Arce et al. 2007).

Atomic and molecular tracers are needed to probe the physical conditions and to evaluate and disentangle the energetic processes that occur in the Class 0/I young stellar objects. Low-J (J ≤ 6, E_u/k_B ≤ 116 K) rotational transitions of carbon monoxide (CO) are among the most widely used tracers (Bontemps et al. 1996), but are only sensitive to the cold gas, T ≤ 100 K, from both the envelope and the entrained outflow material. Nevertheless, spectrally resolved profiles of CO and ¹³CO 6–5 allowed van Kempen et al. (2009) and Yıldız et al. (2012) to attribute the narrow emission lines to the heating of the cavity walls by UV photons (see also Spaans et al. 1995). High-density tracers such as SiO (n_H2 ≥ 10⁵ cm^-3) have been used to study fast J-type shocks produced at bow shocks where the jet plunges into the cloud (Bachiller et al. 2001). At the same time, theoretical studies of line cooling from dense cores predict that most of the released energy is produced in between these two extreme physical regimes and emitted mainly in atomic [O i], high-J CO and H₂O rotational transitions in the far-infrared spectral region in addition to H₂ mid-infrared emission (Goldsmith & Langer 1978; Takahashi et al. 1983; Neufeld & Kaufman 1993; Ceccarelli et al. 1996; Doty & Neufeld 1997). Therefore, to study the energetics of young stellar objects (YSOs) and, in particular, the relative importance of different energetic processes as a function of the evolutionary state of a YSO, line observations in the ~50–200 μm spectral region are necessary.

The Long-Wavelength Spectrometer (LWS) onboard the Infrared Space Observatory (ISO) has for the first time offered spectral access to the complete far-infrared (IR) window (Kessler et al. 1996; Clegg et al. 1996). Many CO rotational transitions from J = 14−13 to J = 29−28 (for NGC 1333-IRAS4) and several H₂O lines up to E_u/k_B ~ 500 K have been detected in Class 0 sources (Giannini et al. 2001; Maret et al. 2002). On the other hand, H₂O remained undetected in Class I sources, the exception being the outflow position of HH46. CO emission was generally found to be weaker than H₂O, whereas the fine structure lines of [O i] and [C ii] dominate the ISO spectra (Nisini et al. 2002b). The gas cooling budget calculations show similar contributions from lines of CO, H₂O, [O i], and to a smaller extent OH in Class 0 sources. Moreover, an evolutionary trend toward a gradual decrease in molecular luminosity and total line luminosity was established as the objects evolve from the Class 0 to Class I phases. This trend was interpreted as the result of weaker shocks and less shielded UV radiation in the later phase of protostellar evolution (Nisini et al. 2002b).

The Photodetector Array Camera and Spectrometer (PACS) (Poglitsch et al. 2010) onboard the Herschel Space Observatory (Pilbratt et al. 2010)¹ with 25 $\mathrm{9}\stackrel{\mathrm{″}}{.}\mathrm{4}\mathrm{\times }\mathrm{9}\stackrel{\mathrm{″}}{.}\mathrm{4}$ spatial pixels provides an 8 × improvement in spatial resolution as compared to ISO/LWS. The PACS field of view of ~ 47″ is smaller than the 80″ ISO beam, but in many cases it still covers the full extent of the emission from nearby YSOs. For a typical distance of 200 pc to our objects (Table 1), regions of ~9400 AU are observed and resolved down to ~1880 AU. The higher sensitivity and better spectral resolution provides an important improvement in the quality of the spectra. PACS is thus well suited for studies of atomic and molecular emission in the Class 0/I objects, as demonstrated by PACS results on individual Class 0/I sources and their outflows (van Kempen et al. 2010b,a; Nisini et al. 2010; Herczeg et al. 2012; Benedettini et al. 2012; Goicoechea et al. 2012). These results have already indicated relative differences in the gas cooling budget from different sources and differences in spatial distributions of emission between different molecules. Visser et al. (2012) have modeled these early data with a combination of shocks and UV heating along the cavity wall. The strong [O i] and OH emission also suggests there are dissociative shocks in the close vicinity of the protostar (van Kempen et al. 2010b).

In our paper, we address the following questions. How does a YSO affect its surrounding cloud and on what spatial scales? What are the dominant gas cooling channels for deeply embedded YSOs? What do they tell us about the physical components and conditions that cause excitation of the observed lines? How do all of these processes change during the evolution from the Class 0 to the Class I stage? To this end, we present Herschel-PACS spectroscopy of 18 Class 0/I YSOs targeting a number of CO, H₂O, OH, and [O i] lines obtained as part of the “Water in star-forming regions with Herschel” (WISH) key program (van Dishoeck et al. 2011). WISH observes about 80 protostars at different evolutionary stages (from prestellar cores to circumstellar disks) and masses (low-, intermediate-, and high-mass) with both the Heterodyne Instrument for the Far-Infrared (HIFI; de Graauw et al. 2010) and PACS. Our paper only focuses on low-mass YSOs and is closely associated to other WISH papers. Specifically, Kristensen et al. (2012) studies the spectrally resolved 557 GHz H₂O line observed towards all our objects with HIFI. Wampfler et al. (2013) analyzes the same sample of sources but focuses on the excitation of OH in the Class 0/I sources, whereas full PACS spectral scans of two sources are published by Herczeg et al. (2012; NGC 1333-IRAS4B) and Goicoechea et al. (2012; Ser SMM1). A synthesis paper discussing the HIFI, PACS, and SPIRE data being obtained in WISH and other programs is planned at the final stage of the program.

The paper is organized as follows. Section 2 introduces the source sample and explains the observations and reduction methods; Sect. 3 presents results that are derived directly from the observations; Sect. 4 focuses on the analysis of the data; Sect. 5 provides the discussion of the results in the context of the available models, and Sect. 6 summarizes the conclusions.

2. Observations

2.1. Sample selection

We used PACS to observe 18 out of 29 Class 0/I objects selected in the low-mass part of the WISH key program. The WISH source list consists of nearby (D ≲ 450 pc), well-known young stellar objects for which ample ground-based single dish and interferometer observations are available (for details concerning the WISH program see van Dishoeck et al. 2011). The remaining 11 sources, that were not targeted with PACS within WISH were observed in the “Dust, Gas and Ice in Time” key program (DIGIT, PI: N. Evans; Green et al. 2013; Dionatos et al. 2013; Jørgensen et al., in prep.; Lee et al., in prep.).

Table 1 presents our sample of objects together with their basic properties. Bolometric luminosities and temperatures were derived using our new PACS data supplemented with observations found in the literature (see Sect. 2.4 for spectral energy distribution discussion). Envelope masses are from Kristensen et al. (2012), which includes a discussion of the impact of new PACS measurements on the derived physical parameters.

2.2. Observing strategy

The far-IR spectra were obtained with PACS, an integral field unit with a 5 × 5 array of spatial pixels (hereafter spaxels). Each spaxel covers $\mathrm{9}\stackrel{\mathrm{″}}{.}\mathrm{4}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}\mathrm{9}\stackrel{\mathrm{″}}{.}\mathrm{4}$, providing a total field of view of ~ 47″ × 47″. The full wavelength coverage consisted of three grating orders (1st: 102–210 μm; 2nd: 71–105 μm; or 3rd: 51–73 μm), two of which were always observed simultanously (one in the blue, λ < 105 μm, and one in the red, λ > 102 μm, parts of the spectrum). The velocity resolution ranges from ~75 to 300 km s^-1, depending on the grating order and the wavelength. The highest spectral resolution is obtained at the shortest wavelengths, below 65 μm. Two nod positions were used for chopping 3′ on each side of the source. Typical pointing accuracy is better than 2′′.

Two observing schemes were used in our program: line spectroscopy mode to cover short spectral regions and range spectroscopy mode to cover the full far-IR spectral energy distribution (SED). Line spectroscopy mode uses small grating steps to provide deep integrations and to fully sample the spectral resolution over short (0.5–2 μm) wavelength intervals. This mode was used to observe selected lines for 16 of 18 objects from our sample (Ser SMM1 and NGC 1333-IRAS2A are the exceptions). We targeted 12 H₂O lines (E_u/k_B ~ 100−1320 K), 12 CO lines (E_u/k_B ~ 580−3700 K), and 4 OH doublets (E_u/k_B ~ 120−291 K), as well as the [O i] and [C ii] lines (full list of available lines is included in Table A.1). BHR71 and Ser SMM4 were only observed in a limited number of scans within the WISH program; range spectroscopy observations of those sources are analyzed in DIGIT (Jørgensen et al., in prep.; and Dionatos et al. 2013, respectively).

The range spectroscopy mode uses large grating steps to quickly scan the full 50–210 μm wavelength range with Nyquist sampling of the spectral resolution. This mode achieves a spectral resolution of R = λ/Δλ ≈ 1000–1500 over the full spectral range, which includes 37 high-J CO transitions, as well as  140 H₂O transitions (J < 10, E_u/k_B < 2031 K) and 11 OH doublets. NGC 1333-IRAS2A, 4A, 4B, and Ser SMM1 were observed with full range spectroscopy within WISH. The NGC 1333-IRAS2A data were taken during the science demonstration phase, when the optimal PACS settings were not yet known. The data are therefore of poorer quality than the other full range spectroscopy observations.

2.3. Reduction methods

Both line spectroscopy and range spectroscopy basic data reduction were performed with the Herschel Interactive Processing Environment v.8 (HIPE, Ott 2010). The flux was normalized to the telescopic background and calibrated using Neptune observations. Spectral flatfielding within HIPE was used to increase the signal-to-noise ratio (for details, see Herczeg et al. 2012; Green et al. 2013). The overall flux calibration is accurate to ~30%, based on the flux repeatability for multiple observations of the same target in different programs, cross-calibrations with HIFI and ISO, and continuum photometry. The 5 × 5 datacubes were further processed with IDL.

Since the spaxel size stays fixed, whereas the Herschel beam size increases with wavelength, the wavelength-dependent loss of radiation in a spaxel for a well-centered point source is observed to be ~30% in the blue to ~60% in red parts of the spectra (see PACS Observers Manual). Most of the radiation that leaks outside a given spaxel is captured by the adjacent ones. However, the far-IR emission from many Class 0/I objects is spatially extended on scales of >10′′, which are resolvable by PACS. For these sources the central spaxel fluxes corrected for the point spread function (PSF) using the standard wavelength-dependent values provided by the Herschel Science Center largely underestimate the total emission from the source. Thus, in this paper, either a sum of 25 spaxels (for lines at λ ≥ 100 μm) or a sum of the spaxels with detected emission (for weak lines at λ ≤ 100 μm) are taken to calculate line fluxes used for most of the analysis. The only exceptions are in Sects. 3.1 and 4.3, where central spaxel fluxes corrected for the PSF using the standard factors are calculated in order to study the emission in the direct vicinity of the YSOs. All line fluxes are listes in Tables A.2 and A.3.

The approach to use the sum of the fluxes of all spaxels results in a lower signal-to-noise ratio of the detected lines; some of the weak lines become undetected. Therefore, we developed the “extended source correction” method, which provides wavelength-dependent correction factors for the brightest spaxel(s). This method is well suited for the extended, Class 0/I sources. The details of the method are given in Appendix B; it is primarily applied to sources for which full line scans are available.

The PACS maps show that Herschel was mispointed for some of our objects. Continuum emission of BHR71, IRAS 15398 and TMR1 peaks in between a few spaxels. In the case of IRAS 15398, our observations were centered on the 2MASS position, which is offset by ~10″ from the far-infrared source position as determined from SCUBA maps by Shirley et al. (2000). Continuum emission from TMC1 and TMC1A peaks off-center in the PACS array, but is well confined to a single spaxel.

2.4. Spectral energy distributions

Radiation from the inner regions of a YSO is absorbed by dust in the envelope and re-emitted in the far-IR. As the evolution proceeds, the SED due to cold dust of a young Class 0 source evolves to a warmer Class I SED with emission also observed at shorter wavelengths. Therefore, SEDs are a useful starting guide for estimating the evolutionary stage of a YSO (Lada 1999; André et al. 2000). PACS continuum observations cover the SED peak of these embedded sources and thus provide a more accurate determination of bolometric luminosity (L_bol) and bolometric temperature (T_bol) than previously available.

Based on our PACS continuum measurements and the literature measurements, new values of L_bol and T_bol have been calculated. They are included in Table 1 and were also presented by Kristensen et al. (2012). The details of the calculations, the continuum values at different PACS wavelengths, and the actual SEDs are presented in Appendix C.

3. Results

3.1. Emission spectra

PACS spectroscopy of our sources reveals rich emission line spectra superposed on the dust continuum emission. Several transitions of the CO, H₂O, and OH molecules, as well as atomic emission from [O i] are detected. Emission in the [C ii] line is only rarely detected and associated with the young stellar object.

Figure 1 presents a line inventory at the on-source position for Class 0 sources (central spaxel²). The Class 0 spectra show detections of at least one line of H₂O, CO, OH, and [O i] each for every object (all except NGC 1333-IRAS2A, which has a high upper limit). The H₂O 2₁₂–1₀₁ line at 179.5 μm is the strongest observed water line and often the strongest far-IR line in general, only comparable with CO 14–13 and [O i] 63.2 μm lines. CO transitions from J = 14 − 13 to J = 48 − 47 are detected in the richest spectra; typically CO emission from transitions higher than J = 31 − 30 is either weak or undetected. The OH doublet at 84 μm is detected for all sources, except NGC 1333-IRAS2A. The discussion of other OH transitions can be found in Wampfler et al. (2013). The [O i] ${}^{\mathrm{3}}\mathrm{P}_{\mathrm{1}}{\mathrm{-}}^{\mathrm{3}}{\mathrm{P}}_{\mathrm{2}}$ and ${}^{\mathrm{3}}\mathrm{P}_{\mathrm{0}}{\mathrm{-}}^{\mathrm{3}}{\mathrm{P}}_{\mathrm{1}}$ lines at 63.2 μm and 145.5 μm are detected for all sources except NGC 1333-IRAS2A (both lines undetected) and NGC 1333-IRAS4B (the 145.5 μm line undetected).

For Class I objects, on-source spectra are presented in Fig. 2. At least one water line is detected in all Class I sources except TMC1A; H₂O in RNO91 is detected when a few lines are co-added. Unlike the case of the Class 0 sources in our sample, the H₂O 2₁₂–1₀₁ (E_u/k_B = 114 K) line at 179.5 μm is no longer the strongest water or molecular line. For all sources except HH46, the H₂O 2₂₁–1₁₀ line (E_u/k_B = 194 K) at 108.07 μm or the H₂O 3₀₃–2₁₂ line (E_u/k_B = 196 K) at 174.63 μm is the strongest water line, whereas CO 16–15 or CO 18–17 and OH 84.6 μm lines are the strongest molecular lines. The CO lines are typically weaker from Class I than from Class 0 objects, up to a factor of 10 compared with the brightest Class 0 sources, and the CO 24–23 is even undetected for one Class I object, RNO91. On the other hand, the OH 84.6 μm line and both fine-structure [O i] lines are seen in all sources. The [O i] line at 63.2 μm is always the strongest emission line in the far-IR spectrum of the Class I sources. The profiles of the [O i] line at 63.2 μm are discussed in Sect. 3.3.

3.2. Spatial extent of line emission

PACS maps of the line emission in the detected species show a variety of patterns and thus allow us to spatially resolve the emission from different components of a young stellar object.

The Class 0 source NGC 1333-IRAS4A and the Class I source L1489 are used here to demonstrate the differences in spatial distributions of the emission from the objects in our sample. Figures 3 and 4 show PACS 5 × 5 maps for the two sources in the [O i] 63.2 μm, H₂O 2₁₂–1₀₁, H₂O 2₂₁–1₁₀, CO 14–13, CO 24–23, and OH lines. In each map the CO 6–5 blue and red outflow directions are overplotted for comparison (Yıldız et al., in prep.). The same figures for the rest of our objects are included in the Online Material.

The NGC 1333-IRAS4A emission in [O i], CO, and H₂O cover the outflow direction over the entire map, corresponding to a radius of 25″ or 5900 AU from the protostar. The [O i] emission peaks at the red outflow position. The CO and H₂O maps also show a pattern of extended emission but are less concentrated than [O i], although some of this apparent extent can be attributed to the larger PSF at longer wavelengths. The CO emission, however, is rather symmetric and peaks in the center/red outflow position, whereas H₂O, contrary to [O i], is more pronounced in the blue outflow lobe, including the peak of the emission. The OH 84.6 μm line is detected both on-source and off-source, but with a pattern that is difficult to compare with other lines because of low signal-to-noise (as well as at the neigboring IRAS4B position in the SE corner of the map). OH follows the [O i] emission by peaking at the center/red outflow position. On the other hand, the maps of L1489 show that the emission from all species peaks strongly on-source, i.e. within a 5″ radius corresponding to 700 AU distance from the protostar. Weaker molecular and atomic emission is detected along the outflow direction and is more pronounced in the blue outflow position.

These differences are shown further in Fig. 5, which illustrates the extent of line emission from various species and transitions in NGC 1333-IRAS4A, HH46, and L1489, including higher excited H₂O lines. The distributions are normalized to the emission in the central spaxel. For HH46, [O i] and H₂O are strong in the red-outflow position R_I, whereas OH and CO 24−23 are observed only on-source. L1489 shows some extended emission in the blue outflow (in particular in [O i] and OH), but clearly most of the emission originates in the central spaxel. Since L1489 is much closer to us than the other two sources (see Table 1), the extended emission in NGC 1333-IRAS4A and HH46 indeed covers a much larger area on the sky.

NGC 1333-IRAS4A and L1489 are thus the prototypes for the two morphologically different groups of objects: sources with extended emission and sources with compact emission. Figures D.1 and D.2 show the spectra in the four discussed species in the blue outflow, on-source and red outflow positions for those two groups. The adopted selection rule is based on the ratio of the on-source and the outflow [O i]: the sources where the outflow [O i] emission (in a selected position) accounts for more than the half of the on-source emission form the extended group, whereas the sources for which the off-source emission is ≤50% compared with the on-source emission form the compact group.

Table 2 summarizes the results of using the same criterion for the [O i] 63.2 μm, CO 14–13, H₂O 2₁₂–1₀₁, and OH lines in all objects³. The general trends are: (1) in the compact group, [O i] and OH emission dominate the central spaxel, whereas CO and H₂O either follow the same pattern or are off-source-dominated; (2) in the extended group OH is often strong off-source (except L1527 and HH46 where it dominates on-source), similar to CO and H₂O; (3) Most objects in the extended group are Class 0 objects, with the exception of TMC1A and HH46; (4) Class 0 and I sources are almost equally represented in the compact group.

In a few cases both H₂O and CO are extended but in a different manner. For example, L1527 and NGC 1333-IRAS4A show a brighter CO line and a weaker H₂O line in the red outflow position and the opposite in the blue outflow position (Fig. 5). L483 shows similar differences, but with the brighter CO and weaker H₂O line in the blue outflow position. In those three cases the [O i] line is stronger at the position of weak H₂O; the same holds for the OH in case of NGC 1333-IRAS4A (OH is not detected off-source in L483 and L1527). These differences are further discussed in Sect. 5.2.

For all objects, the [O i] emission is seen from the young stellar object and associated outflows rather than extended cloud emission. In the NGC 1333-IRAS4A, 4B, Ser SMM3 and SMM4 regions, spaxels where the emission originates from the nearby sources are omitted. When detected, the [C ii] emission is usually spread across the entire detector and seen in different strengths in the two nods, which both indicate that the emission is primarily produced by the parent cloud. TMC1 is the only source with [C ii] detected from the central source (maps in both nods are shown in the Online Material). In Ser SMM1, [C ii] emission follows the pattern of other species along the outflow direction (Goicoechea et al. 2012). The [C ii] emission is not discussed further in this paper.

3.3. Velocity shifts in O i and OH lines

Figure 6 shows the [O i] line at 63.2 μm towards the Class 0 sources NGC 1333-IRAS4A, L1527, Ser SMM1, and SMM4 as well as the Class I sources TMC1A and HH46 (for comparison between PACS and ISO fluxes of the [O i] lines see Appendix E). The blue and red outflow profiles show significant line velocity shifts and, in particular in the case of HH46 and Ser SMM4, high-velocity line wings⁴.

Early results by van Kempen et al. (2010b) showed that the bulk of the [O i] emission in HH46 comes from low-velocity gas. On top of this “quiescent” profile, high-velocity gas was detected in the blueshifted jet with a centroid velocity of about –170 km s^-1 and in the redshifted jet with a centroid velocity of +100 km s^-1. Such velocity shifts, indicative of an optically invisible “hidden” atomic jet, are seen towards at least a third of our objects (shown in Fig. 6). High-velocity tails are detectable in a few [O i]-bright sources, but have a minor contribution to the total line emission.

The [O i] 63.2 μm velocity shifts and profile wings may be associated with similar features of the OH line at 84.6 μm. For NGC 1333-IRAS4A, the OH 84.6 μm line from the source spaxel is redshifted by 90 km s^-1, compared with 50 km s^-1 for the [O i] profile shift (see also Fig. D.3). HH46 also shows a tentative detection of blue-shifted high-velocity OH material that resembles the [O i] pattern. Within our sample, no other molecular lines have significant centroid velocity shifts, with typical limits of ~40 km s^-1 at <100 μm and ~100 km s^-1 at >100 μm. Such velocity shifts are at the velocity calibration limit and may be introduced by emission that is spatially offset within the slit(s). Inclination affects the projected velocity of the jet but is not likely to be the explanation for why a majority of sources do not show a velocity shift.

4. Analysis

4.1. Rotational diagrams

Boltzmann (or rotational) diagrams are used to determine the rotational temperatures T_rot from level populations for the Class 0/I objects from our sample (see Goldsmith & Langer 1978, for Boltzmann diagrams). For optically thin thermalized lines, the natural logarithm of the column density of the upper level N_u of a given transition over its degeneracy g_u is related linearly to the energy E_u of that level: $\ln \frac{{\mathit{N}}_{\mathrm{u}}}{{\mathit{g}}_{\mathrm{u}}}\mathrm{=}\ln \frac{{\mathit{N}}_{\mathrm{T}}}{\mathit{Q}\mathrm{\left(}{\mathit{T}}_{\mathrm{rot}}\mathrm{\right)}}\mathrm{-}\frac{{\mathit{E}}_{\mathrm{u}}}{{\mathit{k}}_{\mathrm{B}}{\mathit{T}}_{\mathrm{rot}}}$(1)where Q(T_rot) denotes the rotational partition function at a temperature T_rot for a given molecule, N_T is the total column density and k_B is the Boltzmann constant.

The emitting region is unresolved in the PACS data due to the low spatial resolution, thus the number of emitting molecules, , is calculated for each transition, defined as: ${\mathrm{𝒩}}_{\mathrm{u}}\mathrm{=}\frac{\mathrm{4}\mathit{\pi }{\mathit{d}}^{\mathrm{2}}{\mathit{F}}_{\mathit{\lambda }}\mathit{\lambda }}{\mathit{hcA}}$(2)F_λ denotes the flux of the line at wavelength λ, d is the distance to the source, A is the Einstein coefficient, c is the speed of light and h is Planck’s constant.

Figure 7 shows CO and H₂O rotational diagrams calculated using the fluxes measured over the entire 5 × 5 PACS array for the Class 0 source Ser SMM3 and the Class I source L1489. Diagrams for all objects are included in the appendix (Figs. F.1 and F.2).

Full range scan observations cover many more CO transitions than our targeted line scans (van Kempen et al. 2010a; Herczeg et al. 2012; Goicoechea et al. 2012; Manoj et al. 2013; Green et al. 2013; Dionatos et al. 2013). In those observations, two excitation temperature components are clearly present. The lower-T_rot component of ~250−300 K dominates mid-J transitions with E_u/k_B below ~1000−2000 K. A higher-T_rot component of ~500−1000 K dominates high-J transitions with E_u/k_B above ~2000 K. We call these components warm and hot, respectively, in order to distinguish them from the cool component, T_rot ~ 100 K, observed in the J < 14 lines (van Kempen et al. 2009; Yıldız et al. 2012; Goicoechea et al. 2012). Motivated by these observations of complete CO ladders, we fit two linear components to our more limited set of CO data. The exceptions are L723, L1489, TMR1, and TMC1A, where there is no indication of the hot component in our dataset. The physical interpretation of these two components is discussed in Sect. 5.

As an example, the CO diagram of the Class 0 object Ser SMM3 in Fig. 7 show a break around E_u/k_B ~ 1200−2000 K, with a rotational temperature for the warm component, T_rot(warm), of 292 ± 14 K and a rotational temperature for the hot component, T_rot(hot), of 670 ± 50 K. The error bars reflect the uncertainties given by the fit and include the uncertainties in individual line fluxes as given in Table A.2. The temperature fits include only relative flux uncertainties between lines and not the absolute flux uncertainty, which would shift the total luminosity up and down but would not change the temperature. The corresponding values for the Class I object L1489 are 405 ± 20 K and 480 ± 55 K for the warm and hot components, respectively. However, since a single temperature fit to these latter data is valid within the errors, one component fit is used to derive the temperature of 471 ± 20 K. As a result, the ratio of the number of emitting molecules of the hot CO over the warm CO is ~0.1 for Ser SMM3, whereas no value for hot CO can be given for L1489.

The rotational temperatures of warm and hot CO and H₂O, the numbers of emitting CO molecules and the ratio of hot over warm CO for all sources are tabulated in Table 3. The uncertainties due to the limited number of observed lines, associated with the fits and the selection of the break point are discussed in Appendix H.

Median rotational temperatures of the CO warm component are 325 and 420 K for Class 0 and Class I sources, respectively (calculated using the unbracketed values from Table 3). The hot CO average temperature is ~700 K for the Class 0 sources; for the Class I sources the temperatures seem to be 100–200 K lower, but in general they are poorly constrained due to the limited detections. The median logarithm of the number of emitting CO molecules for the warm and hot components of Class 0 sources is 49.4 and 48.7, respectively. Therefore, about 16% of the CO molecules observed by PACS in Class 0 sources are hot.

The H₂O rotational diagrams of Class 0/I sources in Fig. 7 show scatter in the single-temperature fits that significantly exceeds the measurement errors and is due to subthermal excitation and optical depth effects (Herczeg et al. 2012; see also Johnstone et al. 2003 for the case of CH₃OH). We refrain from calculating of the number of emitting H₂O molecules because the high optical depths require orders of magnitude correction factors. Highly excited H₂O emission from at least some Class 0 and Class I sources is seen in the H₂O 8₁₈–7₀₇ line at 63.3 μm (E_u/k_B = 1071 K). The single rotational temperatures obtained from the fit to the H₂O diagrams in Fig. 7 are ~120 K for the Class 0 source and ~170 K for the Class I source. Similarly low values of H₂O rotational temperatures are also obtained from full spectroscopy observations (Herczeg et al. 2012; Goicoechea et al. 2012). The fact that the H₂O rotational temperature is 100 K or higher already indicates that H₂O cannot be in the entrained outflow gas seen in ¹²CO low-J lines.

4.2. Far-infrared line cooling

The CO rotational temperatures of Class 0/I sources presented in Sect. 4.1 are used to estimate the flux in non-observed lines and to calculate the total far-infrared CO cooling. The extrapolation of the fluxes is limited to the transitions in the PACS range, from J = 14–13 to J = 49–48. This accounts for ~80% of the CO cooling calculated for the first 60 rotational transitions of CO (J ≤ 60, E_u/k_B ≤ 10006 K), used in CO cooling calculations by Nisini et al. (2002b) (see Appendix I). Additional line emission arising from the T_rot ~ 100 K component seen in the CO J ≤ 13 with SPIRE (Goicoechea et al. 2012) and ground-based data (Yıldız et al. 2012) is not included in this estimation, because these lines probe a different physical component, the entrained outflow gas.

The far-IR cooling in H₂O lines is calculated by scaling the total H₂O flux observed over the full PACS range towards NGC 1333-IRAS4B and Serpens SMM1 (Herczeg et al. 2012; Goicoechea et al. 2012) to the limited number of lines observed here in the line spectroscopy mode. These two sources, even though both classified as Class 0, have very different water spectra, with IRAS4B showing numerous high-excitation water lines that are absent in Ser SMM1.

In the water rich source NGC 1333-IRAS4B, the total luminosity of the water lines equals 2.6 × 10^-4L_⊙, whereas the luminosity calculated from the selected lines equals 1.0 × 10^-4L_⊙. For Serpens SMM1, where H₂O lines are much weaker as compared to CO, the total water luminosity observed in the PACS range equals 2.4 × 10^-4L_⊙, of which 1.1 × 10^-4L_⊙ is detected in the small set of lines observed in the line scan observations. Thus, the scaling factor from the line scan observations to the total far-IR water cooling, based on these two sources, is ~2.4⁵.

Despite the obvious limitations of the method, which assumes similar gas properties for all the sources, it provides more reasonable values of the cooling than the extrapolation using the H₂O rotational temperature (see Appendix H and I). Indeed, as argued in Appendix I, the adopted scaling of the H₂O luminosity should be robust for a broad range of objects within the quoted uncertainties of ~30%. As a further validation of our approach, the values derived for Class I sources in Taurus agree within 30% with the full range spectroscopy observations obtained in the DIGIT program (Lee et al., in prep.). H₂O cooling in the PACS range accounts for ~ 86% of the total cooling in this molecule (from non-LTE large velocity gradient model of Serpens SMM1, Goicoechea et al. 2012).

For OH, a scaling factor of ~1.5 is derived based on the full scan observations of Ser SMM1 and IRAS4B, calculated in the same way as for H₂O. Fluxes of OH for all our sources are from (Wampfler et al. 2013).

The cooling in the [O i] lines is calculated as the sum of the observed fluxes in the 63 μm and 145 μm lines. Cooling in [C ii] is omitted because the line emission is a minor contributor to the cooling budget (Goicoechea et al. 2012) and usually originates from extended cloud emission.

We use the definition of the total far-IR line cooling, L_FIR, adopted by Nisini et al. (2002b): L_FIR = L_O   i + L_CO + L_H2O + L_OH. The results of the calculations are summarized in Table 4 and illustrated in Fig. 8, which shows the relative contributions to the total far-IR gas cooling, L_FIR, by different atomic and molecular species, individually for each object from our sample. Note that this definition does not include the contribution of H₂ near- and mid-infrared lines to the cooling. At off-source shock positions, H₂ can be a significant, and even dominant, gas coolant (Neufeld et al. 2009; Nisini et al. 2002a). However, at the much more obscured central source position, data on H₂ are generally not available. In case where they are, the H₂ mid-infrared emission (after extinction correction) is found to have a negligible contribution to the total cooling for Class 0 sources (Herczeg et al. 2012; Goicoechea et al. 2012).

H₂O, CO and [O i] are the most important gas cooling channels among the considered species. For the sources with the largest L_FIR in our sample (most of them are young Class 0 sources, except HH46), H₂O and CO clearly dominate the far-IR emission. [O i] is the most important coolant for those objects with relatively low L_FIR (and more evolved, except Ced110-IRS4). The relative contribution of H₂O is the highest for NGC 1333-IRAS4A and IRAS4B, where it accounts for ~50% of L_FIR. The relative contribution of CO is the largest for Ser SMM1 and SMM3, where ~40% of L_FIR is emitted in this channel. [O i] dominates the far-IR emission in HH46 and RNO91 (both are Class I sources) with a relative contribution of ~60%. The OH contribution to the total cooling is the highest for the Class I Taurus sources, in particular L1489, TMR1 and TMC1, where ~20% of the total cooling is done by OH. For the rest of the sources the typical OH contribution is ~5−10%.

The Class I source HH46 was studied by van Kempen et al. (2010b), who used the central spaxel emission to determine relative contributions of 28, 22, 40, and 10% for CO, H₂O, [O i], and OH. Our results for this object, but using the full PACS array and correcting for the missing lines, yield: 16, 10, 63, and 11% for those species. This example illustrates that depending on the extent of emission in different molecular and atomic species, the derived values of the relative cooling can vary significantly, both between species and across the protostellar envelope. In HH46, a significant amount of the [O i] emission is found off-source and therefore in our calculations it is an even more significant cooling channel as compared to van Kempen et al. results (see Table 2 and Fig. D.1 in Appendix).

4.3. Flux correlations

To gain further insight into the origin of the emission of the various species, the correlation of various line fluxes is investigated. In addition, correlations between the far-infrared line cooling and other source parameters are studied. The strength of the correlations is quantified using the Pearson coefficient, r. The 3σ correlation corresponds to the Pearson coefficient of $\mathit{r}\mathrm{=}\mathrm{3}\mathit{/}\sqrt{\mathit{N}\mathrm{-}\mathrm{1}}$, where N is the number of sources used for the calculation. For our sample of N = 18 sources, | r |  < 0.7 denotes a lack of correlation, 0.7 <  | r |  < 0.8 a weak (3σ) correlation and | r |  > 0.8 a strong correlation.

Figure 9 compares the luminosities in the H₂O 2₁₂–1₀₁ line as a function of the CO 14–13, CO 24–23, [O i] 63.2 μm, and OH 84.6 μm line fluxes. Both the integrated luminosities over the 5 × 5 array and the central spaxels corrected for the missing flux according to standard PSF corrections are presented. Data for BHR71 and Ser SMM4 are not included since some lines are missing for these two sources.

A very tight correlation is found between the H₂O 2₁₂–1₀₁ and the CO lines, both on large and on small scales (r ~ 0.95 for 5 × 5 fluxes and r = 0.90 for central spaxel luminosities). Thus, H₂O and CO – both in CO 14–13 and CO 24–23 – most likely arise in the same gas. OH also strongly correlates with H₂O, in particular on the larger scales (r = 0.87), where it very likely probes the same gas as CO and H₂O. However, on the central spaxel, the correlation between H₂O and OH weakens. The same results are found when the H₂O 3₀₃–2₁₂ line is used (E_u/k_B ~ 200 K).

On the other hand, the [O i] 63 μm line luminosities do not correlate with the H₂O line luminosity, either on large or on small scales (r = 0.42 and r = 0.21, respectively). Significant variations in the [O i] luminosities are found in particular for the H₂O-bright Class 0 sources. Similar results are obtained for the [O i] 145 μm line.

A strong correlation is also found between the two [O i] lines, at 63.2 and 145.5 μm (r = 0.96, see Fig. 11). The observed ratios of the [O i] 63 μm/145 μm lines vary from ~6 to ~26 with a median value of ~13 (see Table 5). Such a median line ratio corresponds to n_H > 10³ cm^-3 for optically thin lines excited by collisions with H₂ (see Fig. 4 in Liseau et al. 2006). Line ratios of ~20, as observed for example in IRAS 15398, L1489 and RNO91, are indicative of J-type shocks (Hollenbach & McKee 1989).

The median [O i] 63.2 to 145.5 flux ratio of ~10 for the Class 0 sources and ~16 for the Class I sources could imply stronger shocks for the more evolved sources. However, the [O i] 63.2 flux may be suppressed for Class 0 sources because they have higher envelope masses and therefore higher extinctions than Class I sources (see also Herczeg et al. 2012; Goicoechea et al. 2012).

Mass loss rates can be derived directly from the luminosity of [O i] at 63 μm using the Hollenbach (1985) equation, dM/dt = 10 L(O i 63), in units of 10^-5M_⊙ yr^-1. This formula assumes that the [O i] line traces the primary wind shock and that the shock is a major coolant, which is valid for pre-shock densities of ≤10⁵ cm^-3. Mass loss rates calculated from this simple formula for our sources vary from 10^-8 to 10^-6M_⊙ yr^-1 and are the highest for Serpens SMM1 and HH46. Median values of the mass loss rate are 1 × 10^-7M_⊙ yr^-1 and 0.5 × 10^-7M_⊙ yr^-1 for Class 0 and Class I sources, respectively, which reflects the lower [O i] 63 μm line luminosities for the more evolved sources in Fig. 11. In general, these values are in the same range as those of the entrained outflow material derived from low-J CO observations (Yıldız et al. 2013). However, for some deeply-embedded and heavily extincted sources, e.g. NGC 1333-IRAS4A, the low-J CO mass-loss rates are up to 2 orders of magnitude higher than those determined from the [O i] 63.2 line (see Yıldız et al. 2012 for low-J CO spectroscopy of NGC 1333-IRAS4A).

Figure 10 shows the line luminosities of the CO 14–13, H₂O 2₁₂–1₀₁, and [O i] lines and the total far-IR line cooling L_FIR as functions of physical parameters: bolometric luminosity, L_bol, and envelope mass, M_env (for correlations with bolometric temperature and density at 1000 AU see Appendix G). Values for these parameters are taken from the physical models of Kristensen et al. (2012) and the results of this paper. The values of L_bol and M_env are not fully independent in our sample (a weak correlation exists, see Kristensen et al. 2012).

The CO and H₂O luminosities, as well as L_FIR, correlate with both L_bol and M_env. The correlations are particularly strong when the envelope mass is used. However, no correlation is found for the [O i] fluxes and the [O i] line ratios with physical parameters.

5. Origin of far-IR line emission

The protostellar environment contains many physical components that can give rise to far-infrared line emission. As illustrated in Fig. 5 of van Dishoeck et al. (2011), these include (i) the warm quiescent inner part of the envelope passively heated by the luminosity of the source (the “hot core”); (ii) the entrained outflow gas; (iii) UV-heated gas along the cavity walls; (iv) shocks along the outflow cavity walls where the wind from the young star directly hits the envelope; (v) bow shocks at the tip of the jet where it impacts the surrounding cloud; (vi) internal working surfaces within the jet; and (vii) a disk embedded in the envelope. In the case of shocks, both C- and J-type shocks are possible. Spatially disentangling all of these components is not possible with the resolution of Herschel, but our data combined with velocity information from HIFI and physical-chemical models of the molecular excitation provide some insight into which components most likely dominate the emission (Visser et al. 2012; Herczeg et al. 2012).

Velocity-resolved spectra of ¹²CO up to J_u = 10 (Yıldız et al. 2010, 2012; Bjerkeli et al. 2011; Benedettini et al. 2012) and of several H₂O (Kristensen et al. 2010; Codella et al. 2010; Kristensen et al. 2011, 2012; Santangelo et al. 2012; Vasta et al. 2012) and OH (Wampfler et al. 2011) lines indicate that the bulk of the far-infrared emission has broad line wings with Δν > 15 km s^-1. The quiescent warm inner envelope is primarily seen in the narrow (Δν ≲ 5 km s^-1) ¹³CO and C¹⁸O isotopolog mid-J spectra centered at the source position but this component does not contribute to H₂O emission. The broad line wings combined with the spatial extent of both CO and H₂O beyond the central spaxel (i.e., beyond 1000 AU) argue against the hot core (option (i)) and the disk (option (vii)) being the main contributors, at least for the extended sources. Internal jet shocks (option (vi)) have high-velocity shifts of 50 km s^-1 or more. While such fast-moving gas is seen in some [O i] (van Kempen et al. 2010b, and Sect. 3.3) and HIFI H₂O spectra (Kristensen et al. 2011), it is only a minor fraction of the total emission for our sources (based on a lack of extremely-high velocity (EHV)/ bullet emission in HIFI spectra) and is ignored here. Our maps are generally not large enough to cover the bow shocks at the tip of the outflow, ruling out option (v).

For these reasons, the far-infrared line emission seen here most likely originates in the currently shocked and UV-heated gas along the cavity walls in the protostellar environment (options (iii) and (iv)). This warm and hot gas should be distinguished from the cooler entrained outflow gas with rotational temperatures of ~100 K traced by ¹²CO line wings up to J_u ≈ 10 (Yıldız et al. 2012, option (ii)). Visser et al. (2012) developed a physical model that includes a passively heated envelope, UV-heated outflow cavity walls and heating by small-scale C-shocks. The model successfully reproduces emission from the central spaxel in three sources (emission from other spaxels was not included in the models). The cavity walls, carved by the outflows, are illuminated by the UV from the protostellar accretion shocks at the star-disk boundary although UV produced by shocks inside the cavity itself can also contribute (Spaans et al. 1995; van Kempen et al. 2009). The small-scale C-type shocks are produced when the protostellar wind hits the outflow cavity walls and locally increases the gas temperature to more than 1000 K, driving much of the oxygen into water (Kaufman & Neufeld 1996). The UV-heated gas can provide ¹²CO emission up to J_u ~ 20 depending on source parameters. Evidence for UV-heated gas also comes from strong extended ¹³CO mid-J emission along the cavity walls (van Kempen et al. 2009; Yıldız et al. 2012). In contrast, shocks are predicted to produce the bulk of high-J CO (J ≥ 20) and all of the H₂O in the PACS range (Visser et al. 2012). The importance of the UV-heating was suggested to increase with evolutionary stage in the three sources studied by Visser et al. (2012), whereas the shock emission weakens for the more evolved sources. Qualitative analysis of the [O i] and OH distribution in one case, HH46, argued for a separate dissociative J-type shock centered on the source to explain the detected emission (van Kempen et al. 2010b).

In the following subsections, we seek to test these model predictions qualitatively against the new observations and for a larger sample. As we will see, the UV-heated component contributes less to the observed PACS lines than thought before. Also, evolutionary trends will be explored and compared to previous ISO observations (Giannini et al. 2001; Nisini et al. 2002b).

5.1. Spatial extent of line emission and correlations

The molecular and atomic far-IR emission is commonly observed to be extended along the outflow direction which is consistent with the above picture of the emission being associated with the cavity walls. The physical scales of the extent are ~10⁴ and ~10³ AU for the extended and compact groups, respectively.

Most of the extended sources are well-known Class 0 objects primarily known for their large-scale outflows. The extended emission seen in high-J CO and H₂O transitions is therefore most likely associated directly with shocks that are currently taking place in the outflows. A few Class 0 sources famous for their large-scale outflows, e.g., BHR71 and NGC 1333-IRAS2A, are found in the compact group. Those outflows extend well beyond the PACS field of view mapped here and their active outflow hot spots are not covered in our data. The excitation conditions for high-J CO and H₂O emission in the outflow shocks have been determined recently in several studies based on PACS and HIFI observations (e.g., Santangelo et al. 2012; Vasta et al. 2012; Benedettini et al. 2012; Bjerkeli et al. 2012; Lefloch et al. 2012; Codella et al. 2012; Tafalla et al. 2013). All studies conclude that only a combination of high density and temperature can account for the observed intensities. In that respect, the dense, extended envelopes surrounding the Class 0 objects provide the ideal excitation conditions for extended sources (Kristensen et al. 2012).

The compact group contains most Class I objects in our sample and objects with outflows seen nearly face-on (e.g. Ced110-IRS4). The envelopes surrounding Class I objects are less dense than around Class 0 objects (Jørgensen et al. 2002; Kristensen et al. 2012), and thus less conducive to high-J CO and H₂O excitation except close to the source itself. Their outflows are generally compact.

The fact that both the CO 14–13 and 24–23 fluxes as well as the H₂O 2₁₂–1₀₁ and CO 14–13 fluxes correlate very strongly with each other, both in the 5 × 5 and central spaxels (Fig. 9), and also show a similar spatial distribution, strengthens the conclusion that the H₂O and mid/high-J CO emission originate in the same physical component. Such strong H₂O emission can only originate in non-dissociative shocks (Neufeld & Dalgarno 1989). Thus, the observed correlation suggests that most of the H₂O and CO emission arises in shocks and that only a minor fraction of the CO 14–13 emission originates in the UV-heated gas. Velocity-resolved line profiles of the CO 14–13 (or neighboring CO lines such as CO 16–15) are needed to quantify the relative contributions of UV heating and shocks to the mid-J CO line profiles. For the single case for which such HIFI data are available in the WISH program, Ser SMM1, the line profile decomposition indicates that ~25 ± 4% of the flux comes from UV heated gas, but this number may vary from source to source (Kristensen et al., in prep.).

The fact that H₂O and CO emission is sometimes extended in a different manner (L1527, L483 and NGC 1333-IRAS4A) (see Sect. 3.2) could be due to a variety of reasons, including asymmetric densities and UV fields across the core. Different densities on the two sides could be caused by irregularities or gradients of the surrounding cloud material distribution (Liu et al. 2012; Tobin et al. 2011) and can lead to different excitation of the molecules. The density also affects the penetration of the UV field, which in turn results in abundance variations. For example, the UV field from the star-disk boundary can readily photodissociate the H₂O molecule but is not hard enough to dissociate CO. Higher extinction towards the red outflow lobe, as found for NGC 1333-IRAS4B (Herczeg et al. 2012), can also result in different fluxes towards the two outflow lobes, but this should affect all species similarly as long as lines at similar wavelengths are compared.

The [O i] emission is usually extended along the outflow direction and resembles the pattern of the OH extent. In the two exceptions, L1527 and HH46, the [O i] is extended but the OH emission is seen only at the source position (see Table 2). For the compact sources, both [O i] and OH consistently peak together at the central position. Together with a weaker correlation of OH with H₂O and no correlation of [O i] with H₂O (Fig. 9) as well as a strong correlation between [O i] and OH (Wampfler et al. 2013), this confirms the early conclusion that [O i] and at least part of the OH have a common origin. The [O i] and OH emission is most likely produced by dissociative shocks, in contrast to the non-dissociative shocks which explain the CO and H₂O emission (see Sect. 5.3).

5.2. Excitation

A remarkable observational finding is the rather constant CO rotational temperature of ~300 K for the warm component, as well as the common presence of a hotter component (see also Manoj et al. 2013; Green et al. 2013). In order to assess the physical conditions that could reproduce the observed T_rot of the warm component, non-LTE molecular excitation models were run for hydrogen densities of n(H₂) = 10⁴ − 10⁷ cm^-3 and kinetic temperatures of T_kin = 10² − 10³ K using RADEX (van der Tak et al. 2007) (see Fig. 12). The ratio of column density and the line width N(CO)/ΔV, which enters the escape probability calculation, is chosen such that the emission is optically thin. CO transitions from J = 15 − 14 to J = 25 − 24 are modeled using the high-T collisional rate coefficients from Yang et al. (2010), extended by Neufeld (2012). The transitions are plotted as rotational diagrams and fitted with one rotational temperature, corresponding to the warm component.

The inferred physical conditions cover a broad range of T_kin and n(H₂) (Fig. 12). The range of possible scenarios can be described by two limiting solutions: (i) CO is subthermally excited in hot (T_kin ≥ 10³ K), low-density (n(H₂) ≤ 10⁵ cm^-3) gas (Neufeld 2012; Manoj et al. 2013); or (ii) CO is close to LTE in warm (T_kin ~ T_rot) and dense (n(H₂) > n_crit ~ 10⁶ cm^-3) gas. We favor the latter interpretation based on the strong association with H₂O emission and the high densities needed to excite H₂O (see Sect. 5.1 and Herczeg et al. 2012).

If H₂O and CO indeed co-exist, the low rotational temperatures of H₂O (T_rot ~ 100−200 K) as compared to those of CO (T_rot ~ 300−750 K) indicate sub-thermal excitation of H₂O, because the critical density of H₂O is higher than that of high-J CO by 1−2 orders of magnitude (depending on J level; e.g. n_crit ~ 10⁶ for CO 14–13 and ~ 10⁸ cm^-3 for H₂O 2₁₂–1₀₁). RADEX simulations of optically thick H₂O (see Herczeg et al. 2012 and Fig. B.1 in Kristensen et al. 2011) indicate densities n_H > 10⁶ cm^-3 to reproduce our observed H₂O rotational temperatures. This clearly favors the high-density, moderately high-temperature solution (ii) for the CO excitation. The actual gas density very likely lies between the CO and H₂O critical densities if CO is assumed to be close to thermally excited, i.e., between a few 10⁶ cm^-3 to 10⁸−10⁹ cm^-3. Such high densities are only found within a few hundred AU of the protostar or at shock positions where the gas is strongly compressed (Santangelo et al. 2012; Vasta et al. 2012; Tafalla et al. 2013). The hot CO component is not modeled in this work because it is less well defined than in other studies, which have larger samples of sources with complete 55–200 μm PACS spectra (see Manoj et al. 2013). However, as is clear from Fig. 12, T_rot ~ 450−800 K requires even warmer and higher density gas in our scenario (ii).

In the alternative solution (i) in which the entire CO far-infrared ladder comes from low-density, very high-temperature, high column density material (Neufeld 2012), the question arises where this gas is located in the protostellar environment. The only region close to the protostar (within the PACS central spaxel) where gas densities are lower than ~10⁵ cm^-3 is inside the outflow cavities. If the CO and H₂O emission were to originate from the jet itself, velocity shifts of the high-J CO and H₂O lines would be expected in the short-wavelength parts of PACS spectra, where the spectral resolution is the highest. As discussed in Sect. 3.3, we do not see clear evidence for velocity shifts of 100 km s^-1 or more in short-wavelength H₂O and CO lines. Low velocity emission from inside the cavity can originate from dusty disk winds (Panoglou et al. 2012), but more modeling is required to determine whether the line fluxes and profiles can be reproduced in this scenario. Moreover, this explanation would raise the question why there is no contribution from shocks along the higher density cavity walls, which surely must be present as well.

5.3. Comparison with shock models

The association of the CO and H₂O emission with warm dense shocked gas can be further strengthened by comparison with shock models. H₂O and CO are efficiently excited in non-dissociative C-type shocks (Kaufman & Neufeld 1996; Bergin et al. 1998; Flower & Pineau des Forêts 2010). The models by Visser et al. (2012) used the output of the Kaufman & Neufeld (1996) models to compute line fluxes along the cavity walls. To test our full data set against these models and further explore parameter space, rotational temperatures were computed from the grid of Kaufman & Neufeld (1996) model fluxes and presented in Fig. 13. The results show that CO excitation temperatures T_rot ≈ 300 K for J = 15–25 are readily found for pre-shock densities of ~3 × 10⁴−10⁶ cm^-3 and a wide range of shock velocities. Similar results are obtained for the H₂O excitation temperatures: the observed values of 100–200 K correspond to the pre-shock densities of ~10⁵ − 10⁶ cm^-3 and do not constrain the shock velocities.

The density probed by the CO and H₂O lines is not the pre-shock density; by the time CO and H₂O cooling becomes important, the gas has been compressed by the shock front and a more relevant density is the post-shock density (Flower & Pineau des Forêts 2010). The compression factor, n_post/n_pre, is $\sqrt{\mathrm{2}}\mathrm{\times }{\mathit{M}}_{\mathrm{A}}$, where M_A is the Alfvén Mach number, ν_shock/ν_A, and ν_A is the Alfvén velocity (Draine & McKee 1993). The Alfvén velocity is ${\nu }_{\mathrm{A}}\mathrm{=}\mathit{B}\mathit{/}\sqrt{\mathrm{4}\mathit{\pi \rho }}\mathrm{=}\mathrm{2.18}$ km s^-1$\mathit{b}\mathit{/}\sqrt{{\mathit{\mu }}_{\mathrm{H}}}$ where ρ is the pre-shock density, ρ = μ_Hm_Hn_pre; the latter relation comes from the assumption that the magnetic field is frozen into the pre-shock gas and is B = b$\sqrt{{\mathit{n}}_{\mathrm{pre}}\mathrm{\left[}{\mathrm{cm}}^{-3}\mathrm{\right]}}$μGauss. For a mean atomic weight, μ_H, of 1.4, the compression factor is $\begin{array}{ccc}\frac{{\mathit{n}}_{\mathrm{post}}}{{\mathit{n}}_{\mathrm{pre}}}& \mathrm{=}& \sqrt{\mathrm{2}}{\mathit{M}}_{\mathrm{A}}\mathrm{=}\sqrt{\mathrm{2}}\frac{{\nu }_{\mathrm{shock}}}{{\nu }_{\mathrm{A}}}\\ & \mathrm{=}& \sqrt{\mathrm{2}}\frac{{\nu }_{\mathrm{shock}}}{\mathrm{2.18}\hspace{0.17em}\mathit{b}\mathit{/}\sqrt{{\mathit{\mu }}_{\mathrm{H}}}}\mathrm{=}\mathrm{0.78}\mathrm{\times }{\nu }_{\mathrm{shock}}\mathit{/}\mathit{b}\mathit{.}\end{array}$For a standard value of b = 1, the compression factor is thus ≳10 for all the shock velocities considered here, implying that the relevant post-shock densities are all ≳10⁶ cm^-3 as also shown by the RADEX results (Fig. 13). At high velocities, ~50 km s^-1, the post-shock densities are greater than 10⁷ cm^-3 and the emission is in LTE.

Density is clearly the critical parameter: higher T_rot require higher pre-shock densities which in turn lead to higher post-shock densities. Thus, this shock model analysis is consistent with that presented in Sect. 5.2 and with the typical pre-shock densities found in the inner envelope. The presence of UV radiation in the outflow cavity may affect the shock structure, however, and new irradiated models are required to fully test this scenario (Kaufman et al., in prep.). With the addition of UV, some OH emission is very likely also produced by photodissociation of H₂O.

A large fraction of the [O i] and some OH emission must originate in a different physical component. Because the emission from these species does not correlate with that of CO and H₂O, especially in the central spaxel, their origin is most likely in dissociative shocks rather than non-dissociative C- or J-type shocks. In particular, the absolute flux of [O i] is generally too high to be produced in non-dissociative, C-type shocks (Flower & Pineau des Forêts 2010). Dissociative shocks are typically found at the terminal bow shock where the >100 km s^-1 jet rams into the surrounding cloud, but they can also be located closer to the protostar where the protostellar wind impacts directly on the dense inner envelope near the base of the wind (van Kempen et al. 2010b). In dissociative shocks, [O i] and OH emission greatly dominate over that of H₂O as molecules are gradually reformed in the dense post-shock gas (Neufeld & Dalgarno 1989). Some emission may also arise in the jet itself, as evidenced by the velocity shifts of ~100 km s^-1 of the [O i] line in a few sources (Fig. 6), but this is usually a minor component according to our data.

5.4. Evolution from Class 0 to Class I

Class 0 and Class I objects are relatively short phases of protostellar evolution (0.16 and 0.54 Myr, respectively, see Evans et al. 2009), yet our far-IR observations reveal a significant change in their spectral line characteristics. In the following, evolutionary trends of line fluxes, spatial extent, excitation and far-infrared gas cooling are discussed.

Figure 14 shows the median fluxes of the CO 14–13, H₂O 2₁₂–1₀₁, and [O i] ${}^{\mathrm{3}}\mathrm{P}_{\mathrm{1}}{\mathrm{-}}^{\mathrm{3}}{\mathrm{P}}_{\mathrm{2}}$ lines for the 12 Class 0 and 6 Class I sources. Although the Class I sample is small, a factor of 5 decrease is observed in the CO and H₂O fluxes for the more evolved sources. These trends indicate that the molecular emission decreases quickly with evolution, whereas atomic emission remains similar. One explanation could be that molecular abundances are lower due to enhanced photodissociation as the envelope becomes more dilute. When the average density of the envelope drops (see Kristensen et al. 2012 for our sources), UV photons penetrate more readily and on larger spatial scales, thus dissociating more H₂O and OH to form O. This scenario was originally put forward by Nisini et al. (2002b).

An alternative explanation could be that the molecules are less excited because of the lower densities. The strong correlation of the line fluxes with M_env would be consistent with this option. Also, comparison of CO rotational temperatures for Class 0 and I sources shows that the former are characterized by higher temperatures of the hot component, possibly due to the higher densities. The warm component rotational temperatures of CO are remarkably constant for all sources, however. The rotational temperature of H₂O tends to increase for the Class I sources and is mainly due to the hot water detections in the Taurus objects. Since there are only 6 Class I sources in this sample it is difficult to examine how robust the high H₂O rotational temperatures are for the Class I sources, although these lines are now also seen in some Class II sources (Riviere-Marichalar et al. 2012). Dust extinction could hide the highly-excited H₂O lines at shorter wavelengths in the more embedded Class 0 sources (Herczeg et al. 2012). Disks may also contribute to the far-IR emission seen from some Class I sources but are most likely entirely obscured at short wavelengths by the more massive envelopes around Class 0 sources.

The evolution of the far-infrared gas cooling in CO, H₂O, OH, and [O i] is presented in Fig. 15. Our calculations for Class 0 and I sources are supplemented by results from Podio et al. (2012) for the Class II sources (average cooling of DG Tau A and RW Aur) that drive prominent jets. Similar to the selected line flux trends, a factor of 10 decrease in the CO and H₂O cooling is detected for the Class I sources as compared to the Class 0 objects. Further decrease, although much less, is seen for the Class II objects. The absolute values of the OH cooling are similar in the Class 0 and I objects, but they decrease by a factor of 10 for the Class II objects. Cooling in [O i] does not change that dramatically in the three classes: a factor of 2 lower for Class I and II sources with respect to the Class 0 sources.

As a result, the relative contribution to the total cooling by different species changes significantly in the course of the embedded phase of evolution. In the earliest phase, CO and H₂O dominate the total gas cooling with only a minor contribution from OH and [O i]. In the more evolved Class I phase, all species considered here contribute almost equally to the total gas cooling. The evolution of the total molecular cooling, L_mol, where L_mol = L_CO + L_H2O + L_OH with respect to atomic cooling, L_OI, shows a significant decrease by a factor of ~5 between Class 0 and I sources. In the Class II sources, [O i] is the main coolant among the considered species, followed by CO and H₂O. These conclusions are robust within the uncertainties of our calculations of the total luminosities. As discussed in Sect. 4.2, this discussion excludes any potential contribution from H₂ cooling.

Similarly, the total far-infrared gas cooling and the total cooling over the bolometric luminosity decrease with evolutionary stage, most drastically between the Class 0 and Class I phases. These results agree with the ISO observations of 17 Class 0 and 11 Class I sources (Giannini et al. 2001; Nisini et al. 2002b). The decrease in molecular emission was attributed to the less powerful jet impacting the lower-density envelope of the more evolved sources (Nisini et al. 2002b). Our results show that density may be more important than shock velocity in controlling the emission (see Fig. 13).

6. Conclusions

We have characterized the Herschel/PACS spectra of 18 deeply embedded low-mass protostars. The spatially resolved information allows us to link the emission in different species with the physical components of a young stellar object and study the spatial scales over which the object interacts with its surroundings. The conclusions are as follows:

• 1.

  Emission from the H₂O 2₁₂ − 1₀₁ line at 179.5 μm (E_u/k_B = 114 K) is detected towards all Class 0 and I YSOs from our sample. The highly excited H₂O 8₁₈–7₀₇ line (E_u/k_B = 1071 K) at 63.3 μm is detected in 7 out of 18 sources, in both Class 0 and I objects. CO lines are detected from J = 14 − 13 up to J = 48 − 47.

• 2.

  Emission in the H₂O, CO, OH and [O i] lines is extended along the outflow directions mapped in lower-J CO emission. The extent of emission covers the PACS array in the outflow direction for half of our sources (extended sample), which accounts for a region of ~ 10⁴ AU. The rest of the sources show compact emission limited mostly to the central spaxel (a region of ~ 10³ AU) (compact sample). The extended sample is dominated by Class 0 sources with evidence for active shocks (“hot spots”) currently taking place along the outflow.

• 3.

  Fluxes of the H₂O 2₁₂ − 1₀₁ line and the CO 14–13 and 24−23 lines are strongly correlated with each other as well as with the physical source parameters L_bol and M_env suggesting they arise in the same physical component. In contrast, H₂O and CO fluxes correlate less strongly with the OH, in particular on small scales, and not at all with the [O i] line fluxes.

• 4.

  Rotational diagrams of CO show two distinct components – the warm component with T_rot1 ~ 350 K and the hot component with T_rot2 ~ 700 K. The hot component is weaker for Class I sources. Rotational diagrams of H₂O show scatter due to subthermal excitation and optical depth effects and are characterized by T_rot1 ~ 150 K for all sources, with a tendency for higher temperatures for Class I sources.

• 5.

  CO and H₂O are argued to arise in non-dissociative shocks along the outflow walls with a range of pre-shock densities. UV heating appears to play a minor role in the excitation of these PACS lines but affects the shock structure and chemistry by dissociating H₂O to OH. These shocks are most likely responsible for the warm CO component. The origin of the hot CO component requires further data and analysis.

• 6.

  [O i] and part of the OH originate largely in dissociative shocks at the point of direct impact of the wind on the dense envelope. Only a few sources show high-velocity [O i] emission tracing a hidden atomic jet.

• 7.

  H₂O is a major coolant among the observed far-IR lines and accounts for 25 to 50% of the total far-IR gas cooling, L_FIR. CO is the other important coolant with a contribution of 5 to 50% to L_FIR. The [O i] contribution to the total cooling accounts for 5-30% and increases for the more evolved sources, whereas 1–15% is radiated in the OH lines.

• 8.

  Weaker emission in the molecular lines for the more evolved sources can result either from lower envelope densities resulting in less excitation or from lower abundances of the molecules due to enhanced photodissociation. Both aspects very likely play a role, but density is plausibly the critical factor in controlling the line emission.

Online material

Appendix A: Targeted lines and measurements

Table A.1 lists the species and transitions observed in the PACS range with the line spectroscopy mode in the decreasing wavelength order. Information about the upper level energies, E_u/k_B, the Einstein coefficients, A, the weights, g_up, and frequencies, ν, are obtained from the Cologne Database for Molecular Spectroscopy (Müller et al. 2001, 2005), the Leiden Atomic and Molecular Database (Schöier et al. 2005) and the JPL Catalog (Pickett et al. 1998).

Tables A.2 and A.3 list line fluxes for all our sources in units of 10^-20 W cm^-2. The uncertainties are 1σ measured in the continuum on both sides of each line; calibration uncertainty of 30% of the flux should be included for comparisons with other modes of observations or instruments.

Appendix B: Extended source correction method

To account for the combination of real spatial extent in the emission and the wavelength-dependent PSF, we developed an “extended source correction” method. We first inspected the 5 × 5 spectral (or contour) maps. Contributions from NGC 1333-IRAS4A, Ser SMM6, and an unlabeled object were subtracted from the observations of NGC 1333-IRAS4B, Ser SMM3, and SMM4 (e.g. Yıldız et al. 2012; Dionatos et al. 2013). In the next step, we used two long-wavelength lines of CO and H₂O (CO 14–13 at 185.999 μm and H₂O 2₁₂ − 1₀₁ at 179.527 μm) to visualize the spatial extent of the line emission attributed to each object and summed all the spaxels that contained emission. Since all of our lines, except the [O i] line at 63.2 μm, are spectrally unresolved by PACS, single or double (for OH doublets, closeby and blended lines) Gaussian fits to the resulting spectra are used to calculate the line fluxes of the detected lines.

Summing the spectra from all 25 spaxels increases the noise and often prevents detecting weak lines. Using only those spaxels that contain most of the emission results in a much higher signal-to-noise and ultimately a higher detection rate for lines but fails to include emission that leaks out of those spaxels because of real spatial extent in the line and the instrumental PSF. Therefore, our “extended source correction” method provides a wavelength-dependent correction factor to account for the missing flux. The main idea of the method is to use the brightest spaxels, which contain most of the emission, to measure the line fluxes and then correct the value for the missing flux, contained in the omitted spaxels. This method assumes that weak lines are similarly distributed to the strong ones, with observed differences in spatial distributions caused only by the wavelength dependence in the PSF.

The correction factors are derived using the strongest lines, i.e. those that can be measured in the brightest spaxels as well as in all spaxels that contain emission from the object (usually 25 of them). The ratio of flux in the small, bright extraction region and the large extraction region yields a wavelength-dependent correction curve (see Fig. B.1). A 0th-order fit (horizontal line) is used for the short-λ part of the spectrum, whereas a 1st- or 2nd-order polynomial is used for the long-λ part of the spectrum (≥100–110 μm). All line fluxes are then measured in only the brightest spaxels and multiplied by this correction factor. This method was used primarily for the full spectral scans.

Appendix C: Spectral energy distributions

Figure C.1 shows the SEDs for all of our sources obtained from our PACS spectroscopy and literature measurements from Spitzer-IRAC and MIPS (Evans et al. 2009), 2MASS (Skrutskie et al. 2006), SCUBA (Shirley et al. 2000; Di Francesco et al. 2008), as well as APEX/LABOCA, Bolocam, SEST, ISO and IRAS telescopes. The PACS measurements cover the peak of dust emission and are in good agreement with the previous observations (for more details, see Kristensen et al. 2012).

For PACS, the overlap regions of different orders cause the regions of 70–73 μm, 98–105 μm and 190–220 μm to be less reliable in terms of continuum shapes and flux densities. These regions were thus excluded from our SED analysis. New values of L_bol and T_bol are calculated following the standard definitions of the two physical quantities (e.g., Dunham et al. 2010). Several methods of interpolation were tested for consistency of the results including linear interpolation, midpoint, prismodial method and trapezoidal summation. Among them the trapezoidal summation offered the most stable values and is used in this study.

Appendix D: Spatial extent of line emission

Figures D.1 and D.2 show the spectra in the on-source and outflow positions for objects with extended emission and objects with compact emission, respectively (see Sect. 3.2).

Figure D.3 presents the spectra of the [O i] 63.2 μm line, the OH 84.6 μm line, the H₂O 7₁₆–6₀₇ line, and the CO 30–29 line for the central target, two spaxels in the red-shifted outflow, and two spaxels in the blue-shifted outflow of NGC 1333-IRAS4A.

Appendix E: Comparing PACS and ISO far-IR spectra

Higher sensitivity of Herschel/PACS compared to ISO-LWS allows us to improve the detection rate of H₂O (15 out of

16 Class 0/I sources) and higher-J CO transitions (14 out of 16 sources detected in CO 24–23). In particular, water detections in Class I sources are now possible for the majority of the sources. Detection of the more highly excited CO transitions allows to distinguish the hot component on the rotational diagram, which was not possible with ISO (Nisini et al. 2010).

The chopping capabilities and the spatial resolution of Herschel at the distances of our objects allow us to distinguish the YSO–related atomic emission from the emission from the nearby objects or the surrounding cloud. Herschel observations show that only two objects (Ser SMM1 and TMC1) show [C ii] emission associated with the YSO. The [O i] emission, on the other hand, is clearly extended in the outflow direction and most probably traces the hidden jet. For some of the outflow-dominated emission sources, the Herschel beam does not cover the full extent of the [O i] emission, whereas the ISO beam can suffer from the cloud or nearby sources emission. Table E.1 shows the comparison between the [O i] emission for the sources observed with both instruments.

Appendix F: Rotational diagrams

Figures F.1 and F.2 show CO and H₂O rotational diagrams for all sources from our sample. Two-component fits are used for the CO diagrams and one-component fits for the H₂O diagrams. These fits are used to determine the total cooling budget in the two molecules for each objects, as discussed in Sect. 4. The errors in the temperatures reflect the statistical error of the fit taking the uncertainties in individual line fluxes as listed in Table A.2 into account. They do not include the absolute flux uncertainties since the relative fluxes between lines within a single spectrum have much lower uncertainties.

Appendix G: Correlations

Figure G.1 shows correlations between selected line luminosities and bolometric temperature (T_bol) and density at 1000 AU (n_H2). Strong correlations are found with the latter quantity.

Appendix H: Rotational temperature uncertainties

Figure H.1 shows the CO and H₂O rotational diagrams for the NGC 1333-IRAS4B and Serpens SMM1, using the data published in Herczeg et al. (2012) and Goicoechea et al. (2012). The full spectroscopy data is shown, with the full line coverage in the PACS range, as well as the selected lines only, typically observed in our line spectroscopy mode for 16 sources in our sample.

Rotational temperatures calculated from the rotational diagrams constructed using the full and limited line configurations are in good agreement for the CO. For the assumed ranges of the two components, the warm component T_rot error of the fit is ±15 K and the hot component T_rot error is ±50 − 100 K.

The change of the energy break in a wide range of transitions results in ±20 K error for the T_rot(warm) and ±40 K for the T_rot(hot) for the full spectroscopy data. Those ranges are not well determined by line spectroscopy data only and thus in this work we always use 1700 K for this kind of observations.

H₂O rotational temperatures, on the other hand, are less accurately determined for the line spectroscopy observations than the formal error of the fit would imply. The scatter due to the subthermal excitation and very likely high opacities of the water lines result in significant differences in T_rot calculation, depending on the choice of observed lines. The fit to the water lines chosen in our program underestimates the resulting temperature by about 50–80 K for the NGC 1333-IRAS4B and Serpens SMM1.

Appendix I: Cooling budget calculations

Since different methods of cooling budget calculation exist in the literature and will appear due to the availability of new Herschel observations, we perform here a comparison between the methods and estimate the differences between the resulting budgets.

Appendix I.1: Carbon monoxide

In order to calculate the CO cooling, Nisini et al. (2002b) calculated LVG models that reproduced the detected transitions and used them to determine the fluxes for the first 60 transitions of CO. Due to the limited ISO sensitivity, CO transitions from J = 14–13 to J = 22–21 (Class I) and J = 29–28 (Class 0) were available for the brightest sources only which did not allow them to distinguish the hot component. Their method corresponds to a single-component fit to the excitation diagrams.

We use the entire PACS array line fluxes of NGC 1333-IRAS4B and Ser SMM1 from Herczeg et al. (2012) and Goicoechea et al. (2012) to compare the observed CO total luminosities with those calculated using fits to the excitation diagrams (see Table I.1). In particular, we show the results of the two-components fits for PACS data and three-components fits for PACS and SPIRE data, available for Ser SMM1. Three wavelength ranges are included: (i) the PACS range from 54–60 to 190 μm; (ii) the Nisini et al. 2002 range, namely 44–2601 μm; (iii) the PACS + SPIRE range, from 60 to 650 μm (for SMM1 only). Additionally, we include the fitting results to the observed lines only, in order to check how good our two/three-component linear fits reproduce the observations.

The calculations in Table I.1 show that multi-component linear fits to the excitation diagrams agree well with the observed values of the total CO luminosity from the detected lines (rows 1+2 and 5+6). The uncertainties correspond to different choice of the break energy for the two components. The fits are then used to extrapolate the fluxes of the lines which are either blends or fall in the region of ~100 μm, where the measured line fluxes are less reliable (row 3 for PACS range and row 7 for PACS+SPIRE range). The resultant total CO luminosities are a good measure of the far-IR CO cooling (for the PACS range) and total CO cooling (for PACS+SPIRE range).

The example of SMM1 shows that the additional CO emission from the entrained outflow gas increases the CO luminosity by a factor of 1.3 with respect to the extrapolated values from the warm CO component (5.23 × 10^-2L_⊙ versus 6.93 × 10^-2L_⊙).

The relative CO luminosity in different spectral regions (A: J = 4–3 to J = 13–12 B: J = 14–13 to J = 24–23 and C: J = 25–24 to J = 44–43) for SMM1 is ~2:2:1 (A:B:C) and for IRAS4B is ~2:1 (B:C), when the additional cold CO component is included. Thus, approximately 80% of the CO luminosity comes from transitions lower than CO 24–23, roughly equally in both SPIRE and PACS ranges. The poorer determination of the hot component rotational temperature does not affect significantly the total cooling determination.

Appendix I.2: Water

Nisini et al. (2002) have calculated H₂O luminosities based on the assumption that H₂O emission arises from the same gas as CO, for which large velocity gradient (LVG) models determined the gas temperature and density. The model prediction were used to extrapolate H₂O line fluxes for rotational transitions with J < 10 and E_u/k_B < 2031 K.

Total H₂O cooling in this work (see Sect. 4.2) is calculated based on the full spectroscopy data for the NGC 1333-IRAS4A and Serpens SMM1 (Herczeg et al. 2012 and Goicoechea et al. 2012). The average scaling factor than transfers the luminosity observed in the selected lines in the line spectroscopy mode to the total water luminosities (as observed in the range spectroscopy) is 2.4 ± 0.3.

An alternative method considered for the water cooling calculation is the extrapolation of the non-observed line fluxes based on the H₂O rotational temperature. Table I.2 compares the results of this method with the values obtained when the scaling factor was used. The extrapolation is done for (i) the 328 lowest rotational transitions of water (J < 10, E_u/k_B < 2031 K, so called “Nisini et al. range”); (ii) the same transitions but only for PACS range. Molecular information is obtained from the JPL and CDMS catalogs (Pickett et al. 1998; Müller et al. 2001, 2005).

Calculations for Serpens SMM3 show that the extrapolation of the fluxes based on the fitted rotational temperature results in a factor of ~ 2 higher total water luminosities than the value calculated using the scaling factor of 2.4. Even higher values are obtained when we extend the range used in Nisini et al. (2002). Additionally, the rotational temperature of H₂O derived from the line scan data is very likely underestimated by a factor of ~ 1.3 − 1.6 (see Appendix H). This uncertainty has an effect on the derived, extrapolated, H₂O total cooling.

In summary, from the comparisons it is concluded that the total luminosity of both CO and H₂O in the PACS range is accurate to 30% or better.

Appendix J: PACS maps for all sources

Figures J.1–J.16 show PACS 5 × 5 maps for all the sources, except NGC 1333-IRAS4A and L1489 presented in the main text, in the [O i] ${}^{\mathrm{3}}\mathrm{P}_{\mathrm{1}}{\mathrm{-}}^{\mathrm{3}}{\mathrm{P}}_{\mathrm{2}}$, H₂O 2₁₂–1₀₁, CO 14–13 and OH lines (unless stated otherwise in the captions). In each map the CO 6–5 blue and red outflow direction are overplotted for comparison (Yıldız et al., in prep.). Color frames show the blue and red outflow positions used to create Figs. D.1 and D.2.

---

Acknowledgments

The authors thank Linda Podio for making her results concerning cooling in Class II sources available for us prior to the publication and the referee for numerous suggestions which improved the paper. A.K. acknowledges support from the Christiane Nüsslein-Volhard-Foundation, the L’Oréal Deutschland and the German Commission for UNESCO via the “Women in Science” prize. J.R.G. is supported by a Ramón y Cajal research contract and thanks the Spanish MICINN for funding support through grants AYA2009-07304 and CSD2009-00038. D.J. acknowledges support from an NSERC Discovery Grant. Herschel is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA. Astrochemistry in Leiden is supported by the Netherlands Research School for Astronomy (NOVA), 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).

References

All Tables

All Figures