© ESO 2020

1 Introduction

Since the discovery of Herbig-Haro objects, emission lines tracing protostellar jets and outflows have acted as beacons witnessing the formation of stars. Over the course of the last couple of decades, a series of infrared as well as millimeter and submillimeter facilities with increasing sensitivity have enabled the detection of many new lines from increasingly complex molecular species. These detections testify that the surroundings of protostars are not only dynamically active, but also chemically complex. The study of emission lines provides us with the means to understand the physical, kinematical, and dynamical conditions of the emitting gas. When combined with spatial information in spectral line maps, they become essential tools to decipher the complex interplay between different physical components that lead to the formation of stars.

Emission lines provide the efficient cooling of gas that can be excited through various energetic processes, such as the interaction with energetic radiation or the action of supersonic shocks. Lines emitted in the mid- and far-infrared part of the spectrum provide efficient cooling mechanisms of the gas in dense environments around protostars where extinction from the protostellar envelope is very high. Surveys of the environment around protostars with infrared space facilities, such as the Herschel Space Observatoryand the Spitzer Space telescope, have revealed that the major coolants of excited gas are lines from molecules (H₂CO, H₂O), atoms (O, S), and ions (Fe⁺, Si⁺) (e.g., Karska et al. 2018; Green et al. 2013). The contribution of each line in different processes, however, remains unclear and it is hindered by instrumental limitations in terms of sensitivity, angular and spectral resolution, and spatial coverage.

Aiming to observe a large sample of sources during the limited lifetime of cryogenically cooled space missions, spectra are taken with a single observation, pointed on the source. In a few cases, small scale maps covering the protostar along with its outflows have provided valuable information on the great diversity of excitation mechanisms that can be revealed by the spatial association of different molecular and atomic tracers. Such an example is the case of HH 211 that was mapped with Spitzer infrared spectrograph (IRS) and Herschel PACS spectrographs (Tappe et al. 2008; Dionatos et al. 2010a, 2018, respectively), where the morphological diversity of the emission patterns from different molecules have helped us to backtrace the prevailing excitation processes.

Large scale spectral maps of star forming regions are even more scarce. However, such observations can help us study star-formation consistently for a number of individual objects and at the same time the influence of the assembly of forming stars onto their environment. The only case for which such maps exist in the mid- and far-IR are for the NGC 1333 star forming region, that has been observed with Spitzer IRS (Maret et al. 2009) and Herschel PACS (Dionatos & Güdel 2017; hereafter Paper I). In this paper we present Herschel HIFI and SPIRE maps of NGC 1333 complemented with additional Spitzer IRS maps not included in the original work by Maret et al. (2009), aiming to understand kinematic and energetic properties of protostellar jets and outflows. The paper is structured as follows. In Sect. 2 we report on the reduction of the data and in Sect. 3 we discuss the morphology of the reconstructed emission line maps. Analysis of the data is then discussed Sect. 4, where we derive the excitation conditions and the outflow kinematical and dynamical properties, along with estimations of the water abundance. In Sect. 5 we discuss our results, aiming to set a reference for the usage and limitations of different tracers in constraining the jet and outflow properties. Finally, a summary is provided in Sect. 6

2 Observations and data reduction

2.1 Spitzer IRS

Spitzer maps employed in this work belong to two different observing programs with identification numbers (ID) 116 (C. Lawrence, PI) and 20378 (E. Bergin, PI). Observations were obtained on the 6th of March 2006 and between March and September 2006 for ID 116 and 20378, respectively, and data were retrieved from the Spitzer Heritage Archive¹ (SHA). Observations belonging to program 116 employed the long-wavelength, low-resolution (LL, 13.9-39.9 μm) module of the infrared spectrograph (IRS) in slit-scan mode to cover an area of 23.′96 × 17.′61. The low-resolution maps are presented in this work for the first time (see Fig. 1). Observations from program 20378 employed the short-wavelength, low-resolution (SL, 5.1–14.2 μm), the short-wavelength high resolution (SH, 9.9–19.5 μm), and long-wavelength high-resolution (LH, 18.8–37.1 μm) IRS modules. These maps have been presented and discussed in Maret et al. (2009). SL module maps cover an area of 13.′9 × 9.′1, while the SH and LH cover 9.′86 × 5.′91 and 9.′36 × 8.′47, respectively. The high resolution IRS modules can in principle produce higher quality maps compared to the low resolution modules. However the SH maps leave an unobserved gap at the center of the cluster and provide a rather limited coverage. The LH maps on the other hand, with the exception of the central cluster region are sparse and do not provide a complete spatial coverage (see Maret et al. 2009). For these reasons, we based our analysis mostly on the low-resolution modules and used the high resolution maps in regions where they could provide a complete coverage. All data were initially reduced using the CUBISM software (Smith et al. 2007) and rogue pixels were removed using spectra from emission-free areas of the maps. Further processing to produce line maps included Gaussian-line fitting after removing a first or second order polynomial baseline.

2.2 Herschel HIFI

Observations were retrieved from the Herschel Science Archive² (HSA) and were obtained on the 2nd of February 2012 (H. Arce, PI). The Heterodyne Instrument for the Far-Infrared (HIFI; de Graauw et al. 2010) onboard the Herschel Space Observatory (Pilbratt et al. 2010) was used in on-the-fly (OTF) mapping mode to scan the star-forming region NGC 1333 at Nyquist spatial sampling. From the Wide-Band Spectrometer (WBS) back-end providing a ~4 GHz bandwidth the sub-band 1a was used to observe the H₂O line at 557 GHz. The resulting map covers a region of 13.′64 × 6.′48 while the pixel size at the frequency of the H₂O line at 557 GHz is 18. ′′6. The corresponding instrumental half-power beam width is ~37. ′′5. HIFI data were reduced using HIPE version 15.0.3244 using calibration file HIFI_CAL_25_0. HIFI line maps were created using home-grownpipelines. In these, the line intensity and line-centroid velocity shifts for each spaxel were determined through fitting a Gaussian after subtracting a first- or second-order polynomial baseline. The coverage of the HIFI maps is presented in Fig. 1.

2.3 Herschel SPIRE

Observations were obtained between August 27th and 28th 2013 and data records were also retrieved from the HSA (H. Arce, PI). The Spectral and Photometric Imaging Receiver (SPIRE; Griffin et al. 2010) onboard the Herschel Space Observatory (Pilbratt et al. 2010)was used in spectral raster mode at the intermediate sampling rate resulting in one beam spatial sampling. The SPIRE spectrograph consisted of two modules, the Spectrograph Short Wavelengths (SSW) module covering from 191 to 318 μm and the Spectrograph Long Wavelengths (SLW) module ranging between 294 and 671 μm. The sampling rate results in spectral maps with a spaxel size of 19′′ and 35′′ for SSW and SLW, respectively, while the total coverage of the combined raster maps extends to an area of ~13′× 13′, shown in Fig. 1. SPIRE data were reduced in HIPE version 15.0.3244 using the standard pipeline for extended sources withthe SPIRE calibration dataset spire_cal_13_1, including apodization. Beyond basic reduction, SPIRE data were further processed with home-grown routines. As with the other data, line maps were constructed by fitting a Gaussian to each line of interest after subtracting a first order polynomial baseline.

3 Line-emission morphology

3.1 IRS maps

In Fig. 1 we present the angular distribution of the H₂ S(1) rotational transition at 17 μm. The H₂ map has a much larger coverage and provides higher contrast when compared to the maps obtained with the high resolution IRS modules of the same line (see Maret et al. 2009). In addition, spectral-line maps are not affected by continuum and scattered light emission as seen in the IRAC bands (Gutermuth et al. 2008) providing a clearer view of the excited gas. Morphologically, the H₂ emission peaks at the position of the B star 03260+3111 to the northeast and has a secondary peak surrounding BD +30459 that resembles an extended shell. At lower levels, molecular hydrogen delineates most known protostellar outflows, with the ones clustered around the SVS13 system being the most prominent. The morphology of the H₂ emission is very similar to that of the [O I] maps presented in Paper I indicating that most of the molecular hydrogen is excited in the same shock interfaces as atomic oxygen. Finally, the H₂ emission extends far to the south, beyond the regions covered with any of the Herschel instruments discussed in the following paragraphs.

Fig. 1 Spectral line map of the H2 S(1) pure rotational transition observed with Spitzer IRS tracing gas excited in outflows and photon-dominated regions in NGC 1333. Positions of Class 0 and I sources in the star-forming region are indicated with filled circles and triangles, respectively.Filled star-symbols mark the position of the B stars 03269+3111 and BD+30459 in the field. Dashed frames delineate the extent of the regions mapped with Herschel HIFI and SPIRE that are discussed in the following sections. The extent of the PACS maps presented in Paper I is also shown for reference.

3.2 HIFI maps

The profile of the H₂O 1₁₀ − 1₀₁ line presents a wide morphological diversity when observed at different locations of NGC 1333. Four representative line profiles are shown in Fig. 2. The lower two spectra at the bottom of Fig. 2 present wide, high-velocity line-wings extending to velocities of ± 50 km s⁻¹. In striking contrast, the emission line-profile at the top of Fig. 2 is very narrow with a width of just 2–3 km s⁻¹ residing at the velocity of local standard of rest of the cloud (v_lsr). The second spectrum below the top appears to be combining both high velocity wings and a very narrow peak along with a very narrow absorption feature close to v_lsr. The wide morphological diversity of the 1₁₀−1₀₁ H₂O line at different locations of NGC 1333 was first reported in Bergin et al. (2003), based on observations with the Submillimeter Wave Astronomy Satellite (SWAS). Similar morphological characteristics based on Herschel pointed observations around the protostellar sources IRAS 2A, IRAS 4A and IRAS 4B have also been reported in Kristensen et al. (2012) and Mottram et al. (2014), where these authors further decompose the high velocity wings into a broad and a medium velocity components.

The integrated intensity map of the H₂O 1₁₀ − 1₀₁ line is presented in the left panel of Fig. 3. Water emission appears to bright at the center and the south of NGC 1333, tracing the general outflow morphology as recovered in the H₂ maps of Fig. 1, but also the CO maps of Plunkett et al. (2013), and the jet morphology revealed in [O I] Paper I. To the north, waterdoes not follow the emission pattern seen in the maps of all other molecules and atoms, which peak around the B stars to the northeast, and instead shows intensity maxima in the surroundings of the protostellar system IRAS 4A–IRAS 4B to the southeast. When mapping the emission corresponding only to the very narrow component (Fig. 2), the morphology changes drastically. Presented in the central panel of Fig. 3, the narrow component peaks in the vicinity of the B-stars to the north and shows secondary maxima around the Herbig-Haro objects HH 7-11 and HH 12, in line with the emission pattern traced in the line maps of other species. It is worth noting that the narrow component around the IRAS 4A–IRAS 4B system, appears to be rather weak despite it being present in the line profile (Fig. 2). The narrow component represents a small fraction of the total water emission and the peak intensity in the narrow-component maps is more than an order of magnitude lower than the maximum intensity of the fast component. While the water lines observed with HIFI are well resolved, their compound profiles make it complex to retrieve information on the overall gas motions. To this end, we fit the water lines with Gaussian curves and calculate the velocity shift of the Gaussian centroids, as this method provides, to a high accuracy, the motions for the bulk of the gas in outflows (see also Paper I, for a detailed description of the method).

In the right hand panel of Fig. 3, we present the velocity-shift map based on the Gaussian-line centroids for the broad-line component. Peak velocities for the bulk of the gas reach up to ±15 km s⁻¹, however maxima do not seem correlate with the peak of the broad-H₂O emission shown in the left-most panel of Fig. 3. When emission is binned in velocity channel maps (Fig. A.1), then water in the line-wings can exceed 30 km s⁻¹. At those velocities emission mainly traces fast material in the outflows of IRAS 4A and IRAS 2A to the south.

In Fig. A.2 we present the map of a serendipitous detection of a line centered at 558.9515 GHz which possibly corresponds to emission from acetone ((CH₃)₂CO) or methyl-formate (CH₃OCHO). The line profile is narrow and integrated emission levels are low and comparable to those reached by the narrow H₂O component. Line emisison peaks at the locations of known protostellar sources (e.g., IRAS 4A–IRAS 4B, IRAS 2A, HH 12) and provides evidence of active chemistry in the extended envelopes of young stellar objects.

Fig. 2 Characteristic examples of H2O 110 − 101 line profiles from different regions across NGC 1333 observed with HIFI. Velocity scale is shifted so that the source velocity is at 0 km s−1 From bottom to top: red-shifted lobe to the north of IRAS 4A - IRAS 4B complex (black), blue-shifted emission to the south of IRAS 2A (blue), on-source emission from IRAS 4A (green), and from the surroundings of the B-star 03260+3111 to the north(red).

3.3 SPIRE maps

In Fig. 4 we present maps for a series of CO transitions detected with the SPIRE spectrograph. The top row includes ¹²CO maps ranging from J_up = 4 to J_up = 8 observed with the SLW module while the middle panel continues with maps for transitions with J_up from 9 to 13 detected with the SSW module. The lower panel displays maps of ¹³CO emission with J_up between 5 and 8 (at the resolution of the SLW module). Brightness scale is kept constant across all maps of Fig. 4 in order to facilitate direct comparisons.

Lower J, ¹²CO maps are rather featureless with the exception of the HH 7-11 outflow which stands out as a bright emission region in all cases. Integrated emission maps at low J CO transitions are dominated by the optically thick component and cannot differentiate the fast moving gas in outflows (e.g., Dionatos et al. 2010b). At J_up = 7 two additional peaks of emission become apparent around the B star 03260+3111 in the north and in the vicinity of HH 12. At J_up = 9 and above more features become increasingly apparent and combined with the higher resolution provided by the SSW module, the shapes of most outflows are outlined (see also Fig. 6).

¹³CO maps have a very different morphology when compared with the corresponding ¹²CO maps of the same upper level. ¹³CO shows an outstanding peak of emission around source 03260+3111 with a secondary peak around SVS13. The difference in the morphology revealed by the two CO isotopologues most likely reflects differences in the excitation and optical depth that appears to vary significantly between, for example, the B-star 03260+3111 and the Herbig-Haro object HH 7–11. A visual inspection between the maps at the limit where no new structures appear with increasing J, suggests that most of the ¹²CO emission becomes optically thin at J_up ~ 8. This value is consistent with the findings of Yang et al. (2018), that were based on detailed radiative transfer models, and applied on a number of different star-forming regions.

In Fig. 5 we present the emission-distribution of [C I]. The lowest energy transition (³P₀–³P₁) appears to have a more homogeneous, distributed morphology resembling that of the lowest J CO transitions. The map of the first excited level (³P₁–³P₂) shows higher contrast between the regions around protostars and the surrounding medium. A comparison between the [C I] maps with [C II], presented in Paper I, reveals that spatial distribution of the neutral and ionized carbon is very different. [C II] is clearly tracing the photon-dominated region around the B-star 03260+3111 whereas [C I] highlights regions of higher density, possibly outlining extended envelopes surrounding embedded protostars.

Fig. 3 HIFI maps of the broad and narrow 110−101 H2O line components (left and central panels, respectively). The right-hand panel presents a velocity map corresponding to the Gaussian centroid-shift for the broad component that reflects the motion for the bulk of the outflowing gas.

Fig. 4 SPIRE spectral-line maps of the 12CO ladder observed with the long- and short-wavelength modules (upper and central panels, respectively). 13CO maps observed with the long-wavelength module are presented at the lower panel.

Fig. 5 SPIRE maps of the lowest energy transitions of [C I] observed with the long-wavelength module. Upper panel: map corresponding to the fundamental transition [C I] (3P0–3P1) at 610 μm; lower panel: map corresponding to the first transition (3P1–3P2) centered at 370.4 μm (~492 and 809 GHz respectively). To enable comparisons, the [C II] emission from Paper I is superimposed in the lower panel (orange contours).

Fig. 6 SPIRE map with superposition of the elliptical regions adopted from Plunkett et al. (2013), which are used to calculate the properties of outflows in a consistent manner.

4 Analysis

In the following sections we estimate the prevalent excitation conditions responsible for the molecular line emission, and derive the kinematical,and dynamical properties for the gas carried out, as revealed by the different outflow tracers. In order to enable comparisonsbetween this work and the results from Paper I and Plunkett et al. (2013), we assumed the same driving sources and adopted the elliptical outflow loci defined in Plunkett et al. (2013). These are shown for reference in Fig. 6, superimposed on the CO (9–8) line map.

4.1 Excitation conditions of H₂ and CO

The physical conditions driving the excitation of excited gas can be directly constrained by means of excitation diagrams. These diagrams display the logarithm of the column density over the statistical weight as a function of the upper level energy, for a range of transitions of a species. Assuming that the gas is thermalized and the medium is optically thin, a straight line can be fit to excitation diagram points, thus determining the rotational temperature (T_rot) of the gas from the inverse slope and the total column density from the intercept (e.g., Dionatos et al. 2013; Yang et al. 2018). In the case of molecular hydrogen, we can also estimate the ortho-to-para ratio (OPR) and the optical depth using deviations from the assumedlinear distribution of data points in the diagrams. In particular, deviations of the ortho-to-para ratio from the equilibrium value of 3 would appear as vertical displacements between the ortho and para transitions forming a “saw-tooth” pattern. The ortho-to-para ratio can be best determined by examining the alignment of the S(5) with the neighboring S(4) and S(6) rotational H₂ transitions (Dionatos et al. 2010a), however due to the limited sensitivity of the IRS maps at shorter wavelengthswe compare the alignment of the S(1) to the S(0) and S(2) transitions. In addition, the H₂ S(3) transition at 9.7 μm is sensitive to the amount of dust along the line of sight, being located within a wide-band silicate absorption feature at the same wavelength; consequently, a visual extinction value can be estimated by comparing the displacement of this data point in relation to the S(2) and S(4) transitions (e.g., Dionatos et al. 2009). Using these extinction estimations we dereddened the H₂ lines using the extinction law of Chapman et al. (2009) before deriving any measurements from the excitation diagrams.

An example of an H₂ excitation diagram is provided in Fig. 7; diagrams for all outflow lobes are shown in Appendix A. Excitation temperatures, column densities, OPR, and visual extinction values derived from the excitation analysis are summarized in Table 1. As a general trend, the H₂ emissions can be well approximated with a single temperature component (see also Fig. A.3). The H₂ S(0) line intensity appears to be overestimated in most of the cases, most likely due to residual emission from a persisting rogue pixel in the IRS scans. This resulted in underestimating the H₂ ortho-to-para ratio (see Table 1) and therefore derived ortho-to-para ratios should be considered only for comparisons between different regions. Overestimating the S(0) intensity levels however did not compromise the derivation of excitation temperature or column density as did not impact the quality of linear fits. Few exceptions apply when a low number of H₂ transitions are detected; for those cases the excitation temperature was slightly underestimated and column density overestimated. The derived H₂ excitation temperatures lie between ~400 and 500 K and column densities between 10¹⁹ and 10²⁰ cm⁻². Visual extinction ranges between ~10 and 20 mag, in good agreement with extinction estimations from near-IR and X-ray observations. (Walawender et al. 2008).

A typical excitation diagram for carbon monoxide is shown in Fig. 7; more diagrams for all other lobes are included in Appendix A. The distribution of CO lines can be well approximated with two linear segments, corresponding to two temperature components. These segments correspond to lines detected with the SSW and the SLW modules, and in some cases, appear to have a small offset. While this offset may be due to a miscalibration between the two modules, the offset does not appear to be consistent for all lobes. CO lines from J_up = 9 upwards can be optically thin as it is shown for a large number of embedded sources in Yıldız et al. (2012) and Yang et al. (2018). Indeed, optical depth affects mostly SLW lines (as discussed in more detail in the previous section) and as a result, the derived number of molecules in the upper state level of each transition for J_up < 9 is underestimated, creating the apparent mismatch in the excitation diagrams. Therefore column density effects may result in slightly underestimating the total column density of the cold component, however thederived temperature should remain unaffected. Column densities range between ~10¹⁶ and 10¹⁷ cm⁻² and differ by a factor of ~2–3 between the cold and the warm component, while temperatures range between ~50 and ~100 K for the two components. The results from the excitation analysis are detailed for each molecule, source and outflow lobe, in Table 1.

Fig. 7 Left: excitation diagram for H2; observations(open circles) are fit to calculate the temperature and column density after dereddening (open squares) The distribution of data points can be well approximated with a single temperature component. Right: excitation diagram for CO, where the distribution of data points is best fit with a two distinct temperature components. Excitation diagrams for all outflow regions are presented in Fig. A.3.

Fig. 8 Left: [CI] line ratio vs. temperature diagnostic calculated with RADEX for two hydrogen atom densities. Right: temperaturemap reconstructed from the [CI] maps presented in Fig. 4 and the use of the diagnostic.

4.2 Excitation of [C I]

Given the quite different morphology revealed in the maps of the two [C I] lines detected with SPIRE, we attempt to determine the possible underlying excitation mechanisms. We employ the ratio of the ³P₂–³P₁ transition at370 μm over the ³P₁–³P₀ transition at610 μm as a temperature probe, since these lines arise from levels that are separated well in energies and therefore their population depends mainly on the efficiency of the excitation process. We analyzed the [C I] observations employing RADEX, which is a radiative transfer code based on the escape probability approximation (van der Tak et al. 2007). Collisional rates for [C I] (Schroder et al. 1991) were taken from the LAMDA³ database. We ran a grid of RADEX models for a range of densities n between 10² and 10⁶ cm⁻³, kinetic temperatures between 10 and 200 K and total carbon column densities ranging from 10¹⁵ to 10¹⁷ cm⁻². The [C I] line ratio appears insensitive to variations of the total density as shown in the left panel of Fig. 8 and depends mostly on the temperature. Using this diagnostic, we attribute a single approximate temperature to the line ratio corresponding in each SPIRE spaxel to create the temperature map shown in the right hand panel of Fig. 8. [C I] appears warmer at the base of outflows and close to known protostars, reaching temperatures of ~100 K; further out it traces more diffuse, colder gas at ~30–50 K.

4.3 H₂, CO and, H₂O mass flux carried by outflows

Being the most abundant molecule, the observations of H₂ provide a unique opportunity to directly derive the total mass flux carried out by the outflows acting in NGC 1333. We employ the method first described in Dionatos et al. (2009) where the mass flux is derived from the H₂ column density through the following relation: (1)

where μ is the mean atomic weight, m_H is the mass of the hydrogen atom, N(H₂) is the molecular hydrogen column density, A corresponds to the sampled area, du_t is the tangential velocity and dl_t corresponds to the assumed length of the outflow. The factor 2 in Eq. (1) accounts for the two hydrogen atoms comprising a hydrogen molecule. The area A and length dl_t of the outflow correspond to the ellipses circumscribing the outflows in NGC 1333 as described in Plunkett et al. (2013). The column density for these areas is calculated in the excitation diagrams of the previous section. Outflow tangential velocities are taken from Raga et al. (2013), based on proper motion studies of mid-infrared knots traced in multi-epoch Spitzer IRAC images. Knot-emission in the IRAC bands is dominated by H₂ lines tracing shocks, and therefore provide a reliable measurement of the same gas recorded in the IRS spectro-images presented here.

Equation (1) can be modified to provide the mass flux carried out by any molecule X from its column density N(X) subtended in an area A of length dl_t along the path of motion, for gas moving with a velocity as measured on the plane of the sky d u_t : (2)

In the case of CO, column densities were calculated by means of the excitation diagrams in the previous section. From the two temperature components we employed the measurements corresponding to the warmer one (~100 K) as it suffers less from self-absorption and can therefore record more reliably the bulk of the outflow gas. In order to account for the total amount of gas we assumed a “canonical” CO abundance value of 10⁻⁴.

The physical conditions dominating H₂O excitation cannot be determined from a single transition, as is the case with the data available in this study. However an estimation of the water column density can be derived using the relation derived in Tafalla et al. (2013), which is based on detailed local velocity gradient (LVG) model solutions applied to a number of different protostellar outflows: (3)

where I is the line intensity of the water transition at 557 GHz. In addition, we employed an abundance ratio for water-ice of 5 × 10⁻⁵ (van Dishoeck et al. 2014), assuming that all water is efficiently released into the gas phase even at moderate outflow velocities (Suutarinen et al. 2014). Lacking any proper motion observations for either CO or H₂O, we employ the tangential velocities from Raga et al. (2013) that were used also in mass flux estimations of H₂. The same tangential velocities were also used for the mass flux calculations of the [O I] jet in Paper I. The assumption of a common velocity for all species is a simplification that introduces uncertainties in our calculations, albeit it enables direct comparisons. Results from the mass flux calculations for H₂, CO, and H₂O are reported in Table 2 separately for each outflow lobe, along with the length for each outflow, the tangential velocity employed and the corresponding dynamical timescale that is calculated from the relation: (4)

In Table 3 we report the total mass flux corresponding to each protostar for all molecules discussed here, along with the source properties such as the bolometric temperature and luminosity from Plunkett et al. (2013). Dynamical timescales in Table 4 are the average for the two lobes reported in Table 3. To enable direct comparisons, in the last two columns of Table 4 we report the mass flux corresponding to the [O I] jet estimated using two different methods in Paper I.

As a general trend, the CO mass flux is a factor of ~10 higher when compared to H₂, which in its turn is another factor of ~10 higher when compared tothe mass flux of H₂O. The water mass flux is of the same order of magnitude as the mass flux derived for the [O I] emission.

4.4 H₂O momentum and energy

Water lines observed with HIFI are resolved and therefore the momentum carried out by the water excited in outflows can be calculated as discussed in Paper I using the following relation: (5)

where the sum acts over the mass per velocity bin, v is the mean velocity for each velocity bin and v_r is the rest velocity for the water molecule. The mass per velocity bin () was derived from the column density (Eq. (3)) but in this case calculations were applied per velocity bin. (6)

The outflow energy traced by the H₂O line can be derived according to the following relation: (7)

where the mass per velocity bin is defined in Eq. (6) and velocities follow the same notation as for the momentum. The uncertainties in the mass derivation discussed above also affect the estimation of the momentum and energy deposited by the outflows. Derived water momentum and energy is reported in in Table 4, along with the embedded-source properties. In the same table we provide the corresponding momentum and energy values for the [O I] jet and the CO outflow (from Paper I and Plunkett et al. (2013), respectively), for comparisons. The momentum traced by water is in good agreement with the values derived for [O I] but much lower than the momentum derived for CO. In addition water appears to trace a less energetic outflow component compared to CO or [O I].

4.5 H₂O abundance

Water can be formed and excited following very different pathways. A direct determination of the water abundance requires comparison to another molecule of known abundance which is excited through the same mechanism and tracing similar conditions to those traced by water. Following these lines, we calculated the water abundance based on comparisons between the H₂O and H₂ column densities measured in individual outflows, based on the observation that the two molecules are co-spatial and appear to have similarline strengths (Tafalla et al. 2013). We also derived water abundances comparing the column density of H₂O to that of the warm component traced by CO, based on the observation that intermediate J CO and water line profiles appear to closely correlate and therefore trace the same volume of gas (Kristensen et al. 2017). To convert CO-based water abundances to H₂ we assumed a nominal CO abundance of 10⁻⁴, however direct estimations in outflows indicate that this can be as low as 10⁻⁵ (Dionatos et al. 2013, 2018). We report the results of these comparisons for the outflows driven by each source in Table 5. H₂O abundances based on comparisonsto H₂ average to ~ 5.3 × 10⁻⁶ while those based on comparisons to warm COaverage to ~7 × 10⁻⁷.

4.6 Molecular abundances around 03260+3111

In addition to outflow regions, we measured column densities for all molecules examined in this work within a circular region extending to a 30′′ radius around the B-star 03260+3111 to the NE. Excitation of gas in that region is dominated by energetic photons as also evidenced from the detection of very bright [O I] and [C II] emission presented in Paper I (Hollenbach & Tielens 1997), and can therefore be used as a reference to compare with the excitation of gas occurring along outflows. Using the same line of arguments as in the previous paragraph, we derived a water abundance of 4.6 × 10⁻⁸ based on comparisonswith CO (Table 5). The abundance based on comparisons with molecular hydrogen is another factor of ~ 6 lower than the CO-derived ones and 10³ times lower when compared to H₂-based abundances in outflows. This large discrepancy can be explained by a number of physical mechanisms acting in parallel. Photodissociation of molecular hydrogen may occur with absorption of a photon in the Lyman and Werner bands while at moderate column densities far-ultraviolet (FUV) lines become optically thick and H₂ begins to self-shield against radiation at A_V ≃ 2 (Hollenbach & Tielens 1997). Water in contrast absorbs UV radiation less selectively than H₂ and shields less efficiently, so it can be dissociated at higher column densities even for moderate radiation field strengths. Furthermore, when H₂ is dissociated, then the reaction H₂O + H → OH + H₂ leads in further water dissociation through collisions, while at the same time, replenishing the H₂ reservoir (Kristensen et al. 2017). Similar to water, CO shields less efficiently and at higher densities compared to H₂, which can explain the smaller discrepancy of the CO-based water abundance measured around 03260+3111 when compared to outflows.

Interestingly, when comparing CO and H₂ column densities we derive a CO abundance of 1.5 × 10⁻⁵. This measurement provides additional evidence that lower CO abundances than the canonical value of 10⁻⁴ may also apply in conditions very different than the ones prevailing in outflows, and hints that a lower value may be more appropriate in dense regions independently of the local excitation conditions.

5 Discussion

In the following discussion we attempt to shed light on physical processes responsible for the excitation the various molecules, and on the different kinematical and dynamical components they may trace. Our aim is to set a reference for the use and limitations of different tracers in constraining jet and outflow properties. Water emission in particular is discussed in the context of possible formation and excitation pathways, which can explain the observed abundance variations.

5.1 H₂ and CO

Molecular hydrogen appears brightest around the intermediate-mass sources to the NE of the Spitzer maps indicating heating under the influence of energetic ultraviolet (UV) radiation from those sources. In particular, the H₂ emission morphology around BD+30459 traces a shell-like structure with the source positioned at the center. In the central regions of the map, H₂ is excited around HH6, HH12, and SVS13; It appears brightest at the location of the HH 7-11 outflow and also at a symmetrical outflow structure to the northwest of the driving source SVS13. To the south, H₂ emission delineates a web of interlaced linear and bow-shaped features, tracing extended emission from molecular jets and outflows. Despite the limited sensitivity of the Spitzer maps, it is clear that H₂ emission does not trace the bulk of the outflowing material but is rather confined in certain regions. This is consistent with the result from Paper I, where [O I] maps closely follow the H₂ morphology,indicating that both tracers are excited in shocks. Inline with our analysis, H₂ traces warm gas with temperatures of ~450 K while in more sensitive maps a hot component at T ~ 1000 K is also commonly observed (e.g., Maret et al. 2009; Dionatos et al. 2009, 2010a). This interpretation is also consistent with the intrinsic properties of H₂, that, as a very light molecule, has its first excited rotational levels at energies above 300 K.

The mass flux traced by H₂ is at least an order of magnitude higher compared to that derived from [O I], under the assumption that the two elements move at the same velocity. When such shock interfaces are examined in detail, H₂ appears to be excited at the flanks of bow-shocks while atomic lines become bright at the shock-heads (Dionatos et al. 2014) and therefore the assumption of a common velocity between the two species is likely an oversimplification. If the mass flux traced by the two components is preserved, then [O I] should be moving a factor between 10 and 20 faster than H₂. The mass flux traced by H₂ therefore corresponds to the mass of the gas excited in energetic shocks (internal working surfaces) along the jet-propagation path.

Depending on the excitation conditions ¹²CO may trace diverse excitation processes. While low-J transitions are in general optically thick, emission from the line wings can trace high velocity gas in outflows. At the limit the emission becomes optically thin (for J_up ~ 8–9 Yang et al. 2018), ¹²CO appears to trace gas excited in shocks or gas heated under the influence of ultraviolet radiation. This transition is readily evidenced in the CO maps of Fig. 4, especially for the transitions falling in the short-wavelength SPIRE module.

The H₂ mass flux is found an order of magnitude lower when compared to the mass flux derived from the CO lines mapped by SPIRE. When the abundances of the two molecules are considered, this result appears counterintuitive and can only be explained in the context of the H₂ excitation in shocks. In this respect, low-excitation CO lines pick up the bulk of the outflowing material in a turbulent layer surrounding the protostellar jet. H₂ in contrast, despite being the most abundant molecule, does not capture the total mass of the gas entrained in outflows. It is only CO lines of higher excitation energy (above J_up ~ 14) that start tracing gas excited in shocks rather than the turbulent mixing layer and therefore it is only those, higher-J lines, that can trace the same excitation conditions as H₂ (e.g., Yıldız et al. 2013; Kristensen et al. 2017; Yang et al. 2018).

5.2 H₂O

Continuing the discussion from Sect. 3.2, the water-line profiles can be decomposed into two distinct components: a narrow, weaker component moving close to the systemic velocity of the cloud with a width of just ~2–3 km s⁻¹, and a faster component forming line-wings which extend to ~25 km s⁻¹, associated to outflows. The distribution of the narrow component in Fig. 3, shows that it strongly peaks around the intermediate mass sources to the north and especially around the B-star 03260+3111. Strong emission from [O I], [C II], H₂, [C I], and CO at the same location provide evidence that UV radiation from the star creates a photon dominated region (PDR) in its surroundings (Hollenbach & Tielens 1997). The narrow emission component is therefore associated with the influence of UV radiation and when traced, it can provide evidence of the underlying excitation processes.

Other regions in the narrow component map (Fig. 3), where water emission can attest the influence of UV heating, are found to the west of HH 12, along HH 7–11, and in the vicinity of the IRAS 4 system. Maps of ¹²CO transitions with J_up ≥ 9 display emission maxima at exactly the same locations, providing evidence that intermediate excitation CO lines can have an important component that is driven by UV heating. These findings are also in agreement with the analysis of Yıldız et al. (2012) where the narrow component of the CO 6–5 line along the outflow cavities from IRAS 4A–IRAS 4B is interpreted to be resulting from the UV heating of the gas. When comparing the excitation conditions and outflow kinematical properties, no apparent difference is found between outflows with and without signatures of UV heating, indicating that the influence of UV radiation is rather localized and smoothed out at the outflow scales used for our calculations. This is also consistent with the distribution of [C II] which appears prominent in dense photodissociation regions but is not detected along outflows.

The morphology of the faster H₂O component is very similar to the emission maps of jets and outflows traced by [O I], CO, and H₂. The bulk of the H₂O outflow velocity correlates well with the velocity distribution for the bulk of [O I], with a major difference being that H₂O emission traces velocities upto 15 km s⁻¹ as opposed to 50 km s⁻¹ reached in the [O I] maps. At those velocities we expect that most of the water should be released into the gas phase through sputtering (Suutarinen et al. 2014; Kristensen et al. 2017). The distribution of the integrated H₂O intensity maps does not appear to follow the bulk H₂O outflow velocity pattern (see Fig. 3), suggesting that water emission is deficient in fast shocks and predominantly originates from material coasting in outflows after having been accelerated in shocks. Line-wings of water trace velocities up to ~30 km s⁻¹, which is in agreement with the predictions for the dissociation of the molecule at higher velocities in J-shocks (Suutarinen et al. 2014). Water reformation through warm gas-phase reactions is possible in the post-shock cooling zone, however estimated timescales are very long for typical outflow densities (Kristensen et al. 2017).

Based on pointed HIFI observations of water lines around a sample of embedded protostars, Mottram et al. (2014) proposed a decomposition of the line profile to three components: (i) a very narrow component close to the cloud systemic velocity corresponding to emission from the envelope, (ii) a intermediate component arising from shocks along the cavity walls and, (iii) fast components corresponding to molecular bullets or all other locations where the low-velocity wind directly impacts the envelope. The narrow component map in Fig. 3 shows that energetic radiation can escape inside outflow cavities but it can also be generated in-situ from the action of fast shocks (e.g., HH 7–11). The maximum velocity reached for the bulk of the H₂O gas in the line-centroid maps is ± 15 km s⁻¹, which reflects exactly the same limits set for the “cavity shocks” in Mottram et al. (2014) and possibly traces gas that is coasting in outflows after being accelerated. These boundaries also reflect a possible a physical limit where most of water is dissociated in shocks and not replenished at the same rate in the post-shock cooling zone. Water emission at higher velocities is rather localized, as demonstrated in the velocity-channel maps of Fig. A.1, which is consistent with the formation and excitation in molecular bullets.

Water emission peaks in the surroundings of the IRAS 4 system as seen in the integrated emission maps of Fig. 3, while the maximum is found at the location of IRAS 4A. Higher excitation water emission lines in the outflows of IRAS 4B have been previously reported by Herczeg et al. (2012) where a ~50% of the cooling was measured to take place through water lines (Karska et al. 2018). The water line profile around IRAS 4A (Fig. 2) displays signatures of both radiative and outflow components. High-velocity outflows are also prominent in the water velocity channel maps of Fig. A.1. Signatures of UV-heated gas extending in the outflow cavities of IRAS 4A have been isolated in the narrow component of ¹²CO and ¹³CO 6–5 lines (Yıldız et al. 2012). Such radiation can alternatively be produced in accretion shocks very close to the forming star and then scatter to larger distances in the outflow cavities before being absorbed, heating up the material in the cavity walls (Spaans et al. 1995).

The water abundance with respect to CO for a circular region of a 10′′ radius centered on IRAS 4A measures 1.7 × 10⁻⁶, which is a factor of 2 higher than the abundance measured in the outflows of the same source and at least 50% higher than any other outflow. The angular resolution of the CO and H₂O maps does not allow us to isolate the source from the outflow positions, however the observed abundance enhancement must be resulting from the area very close to IRAS 4A which is not included in the abundance calculations of the outflows. These findings are consistent with the best-fitting water abundance profile used to model multiple water transitions around IRAS 4A (Mottram et al. 2013), where the abundance is found to sharply increase in the inner 100 au from the source. Energetic radiation from the central source cannot account for the observed water enhancement, as the example of 03260+3111demostrates that the net effect of UV-heating would be to decrease rather than to increase the water abundance. Gas phase chemistry in outflows could increase the water abundance up to a factor of 2 (Suutarinen et al. 2014), however it is hard to explain why such an enhancement would selectively take place at the base of outflows and not further out. Yet the outflows of IRAS 4A appear to be an exceptional case; following the Gaussian decomposition of Mottram et al. (2014), these outflows possess a uniquely fast “spot-shock” component which, however is not shifted from the source velocity and therefore appears to arise from the envelope cavities rather than a fast jet.

In this respect, a water reservoir may preexist frozen onto the dust grains in the protostellar envelope, and is evaporated at the locations where temperatures reach above ~100 K. In such casesimilar water reservoirs should exist in the envelope of most, if not all other YSOs, however no similar behavior is traced anywhere else in NGC 1333. This suggests that IRAS 4A is a particular case unlike any other YSO in our maps. Indeed, when observed with interferometers IRAS 4A is resolved into a binary system, where the weaker (A2) component is displaying complex chemical signatures of a hot corino (Bottinelli et al. 2004; Sahu et al. 2019). Interferometric observations of HO (Persson et al. 2012) confirm that water emission is coincident with IRAS 4A2 and a series of shocks in the NS direction while no water emission is seen to arise from IRAS 4A1. In this scenario, hot corinos display significant emission from molecules that were inherited from the prestellar phase and survive frozen on the dust grains within their envelopes until they are heated and released in the gas phase for the first time. Deuterium fractionation and a very low ortho-to-para ratio of water in particular favors grain-surface formation at very low temperatures during the prestellar phase (Ceccarelli et al. 2014; Dionatos 2020). Maintaining primordial material in their envelopes may indicate that such sources are younger compared to other YSOs seen in the field, which appear to have processed and exhausted the reservoir of primordial ices in their envelopes. This scenario is also consistent with the abundance increase in the inner envelope implied by models (Mottram et al. 2013), and can also explain the exceptionally high-velocity water which is not observed in other sources (e.g., Fig. A.1, Mottram et al. 2014).

Assuming all volatile oxygen which is not locked up in CO, is bound in water, and that all water is in the gas phase, then the water abundance with respect to H₂ in thermodynamic equilibrium is estimated to be 5−6 × 10⁻⁴ (van Dishoeck et al. 2014). Abundance estimations in this work are at best a factor of ten lower compared to the expected values, however they are in agreement with determinations based on a variety of other observations and methods (e.g., Bjerkeli et al. 2012; Tafalla et al. 2013; Kristensen et al. 2017). Water abundance estimations based on CO are a factor of ~10 lower compared to the ones based on H₂ with the exception of the PDR extending around 03260+3111 where this relation is inverted mainly due to the very different responses of water and H₂ to UV heating (photodissociation vs self-shielding). In some cases water has been reported to closely follow the emission pattern of H₂ (Tafalla et al. 2013), however, the spatial correlation between the two molecules is not clear in the maps presented here. In addition, when UV-heating is present, then water abundances derived from comparisons with H₂ can be severely underestimated. In contrast, CO appears to be a better proxy for estimating water abundances, as both molecules are affected in a similar way under the influence of UV radiation, while their resolved velocity profiles indicate that they trace the same volume of gas (Kristensen et al. 2017).

The main mechanism leading to lower water abundances than expected is not yet clear. As demonstrated in the PDR region, UV-heating can efficiently reduce water abundance and in fact, even modest radiation-field strengths are sufficient to account for the observed abundance values (Kristensen et al. 2017). The narrow-component water map, however, confirms that UV heating may have an influence along HH 7–11 and the outflow of IRAS 4A but not in other outflows. In these two particular cases the water abundance does not appear to vary significantly when compared to other outflows, and therefore, it is not clear what is the net effect of energetic radiation along the outflow cavities. On the other hand, UV heating may be ubiquitous in outflows but not traced in our maps due to the available sensitivity. Fast moving material in outflows can produce dissociative J-shocks, which seem to provide a convincing alternative that can explain the observed water abundance deficiency. The mass flux and momentum carried by water in outflows is comparable with the values measured for [O I] that traces fast shocks, and this link provides additional support to the shock-dissociation scenario to explain the lower abundances of water.

6 Conclusions

We presented Spitzer and Herschel spectral mapping observations of the NGC 1333 star-forming region. Spitzer maps comprise of slit-scan observations using the long-slit IRS modules and record emission of molecular hydrogen. Herschel observations provide a velocity resolved map of water at 557 GHz obtained with HIFI and a series of CO and [C I] line maps observed with SPIRE. H₂ and CO emission was analyzed by means of excitation diagrams and mass flux was calculated for each molecule and outflow. The water mass flux was in addition calculated using a model-derived description relating the water line intensity and the column density. Thevelocity-resolved lines for water were employed to calculate the momentum and energy of the outflows and compare it with previous estimations of other outflow tracers. Finally we calculated the abundance of water based on comparisons with CO and H₂ for all outflows and a few more positions of particular interest. The main results of our analysis can be summarized as follows:

• Morphologically, molecular hydrogen traces emission related to sources of energetic radiation and, at lower levels, to outflows.

• Water lines are resolved to a narrow (slow) and a wide (fast) component. When present, the narrow component contributes up to ~5% of the total line luminosity. Its spatial distribution indicates that it is excited in regions where UV heating is active. The fast component appears to trace gas coasting in outflows after being accelerated in fast shocks as suggested from comparisons between the velocity centroid and integrated emission maps.

• Integrated line emission of water peaks on IRAS 4A showing at the same time an abundance increase of at least a factor of two compared to all other regions in NGC 1333. While outflow and UV-heating signatures are present in the line profile, these processes alone are not sufficient to explain the observed water enhancement. The fact that IRAS 4A2 is a hot corino can provide reasonable grounds to interpret the difference in the water emission between IRAS 4A and other protostars in NGC 1333. The nature of IRAS 4A2 is also confirmed with the serendipitous detection of acetone or methyl-formate in the HIFI maps.

• SPIRE maps of carbon monoxide trace progressively deeper layers into NGC 1333 with increasing upper level energy of the transitions. Emission becomes optically thin at J_up ~ 9. Higher excitation CO transitions have similar morphological characteristics to the H₂ maps, being highly excited around the B-stars in the north but also tracing the complex web of outflows in the south. Intermediate J CO lines appear to peak at the locations where the narrow H₂O component becomes brighter, suggesting that a component of the intermediate J CO emission is sensitive to UV radiation. Radiative heating must be acting on small scales as it produces no observable differences in the excitation and kinematics between outflows with and without UV signatures.

• A similar optical depth effect to CO is observed also for the two [C I] emission lines detected with SPIRE, as the higher excitation transition appears to be excited around known protostellar sources. Temperature maps based on RADEX models show that [C I] traces warm (~100 K) material around protostellar envelopes.

• Excitation analysis shows that H₂ traces warm (~500 K) gas excited in outflows. The derived A_v value of ~15 mag along the outflows is consistent with previous estimations in the same region. Excitation of CO lines detected with SPIRE is best approximated with two temperatures of ~50 K and 100 K while column densities range from ~10¹⁹ cm⁻² for H₂ to ~10¹⁶ cm⁻² for CO.

• CO mass flux is measured to be an order of magnitude higher compared to H₂ and two orders of magnitude relative to H₂O. At the limit that the emission becomes optically thin, CO appears to trace the bulk of the outflow mass while H₂ and H₂O probe the mass of the gas excited in shocks. The mass flux discrepancy between H₂ and H₂O can be explained by dissociation of water in shocks, which is also consistent with the low abundance measurements. This scenario is further supported by the observations that the H₂O momentum and energy correlate strongly with the [O I] measurements from Paper I.

• Water abundance in outflows averages between 7 × 10⁻⁷ and 5 × 10⁻⁶, based on comparisons with CO and H₂, respectively. CO- and H₂-based water abundances are a factor of 10 and 10³ lower in the PDR around the B-star 03260+3111, primarily due to the different responses of these molecules to energetic radiation. Therefore UV heating can explain the lower water abundances, however the narrow H₂O component indicates that radiation may be active in only a couple of cases, and for those outflows significant abundance variations are observed. Dissociation of water in shocks on the other hand is consistent with the results of the morphological, kinematical and excitation analysis and provides the only alternative pathway to interpret the observed water abundance deficiency.

• CO atthe limit that it becomes optically thin is a safer proxy to estimate water abundances than H₂, as it hasa closer response to UV fields and traces the same velocity components in outflows as water.

The observations presented here but also in Paper I demonstrate the power of large-scale maps in interpreting emission from complex systems such as low-mass protostars in a consistent way. They also underline the importance of such maps as opposed to pointed observations in performing consistent comparisons between sources and identifying the ones that appear to have particular properties. While the role of the various molecules in tracing diverse processes starts taking shape, observations of mid-infrared and far infrared lines at higher spectral and angular resolution in the future will be important to understand the detailed mechanisms governing the excitation and cooling of protostellar jets. To this end, observations with SOFIA and soon JWST along with possible future missions such as SPICA will play an important role in deciphering the role of protostellar ejecta as a feedback mechanism in star-formation.

Acknowledgements

This research was supported by the Austrian Research Promotion Agency (FFG) under the framework of the Austrian Space Applications Program (ASAP) project PROTEUS (FFG-866005). The research of LEK is supported by a research grant (19127) from VILLUM FONDEN. Research at the Centre for Star and Planet Formation is funded by the Danish National Research Foundation.

Appendix A Additional maps and figures

Fig. A.1 Channel maps of the water line displaying the integrated emission at different velocity intervals which are indicated onthe top right corner of each panel.

Fig. A.2 Map of a narrow line, possibly from either acetone ((CH3)2CO) or methyl formate (CH3OCHO). Being excited in the vicinity of embedded sources, it may trace hot-corino activity.

Fig. A.3 Excitation diagrams for the H2 fluxes extracted from Spitzer within the elliptical regions.

Fig. A.3 continued.

Fig. A.4 Excitation diagrams for CO fluxes extracted from the SPIRE maps.

Fig. A.4 continued.

References

All Tables

All Figures

Fig. 1 Spectral line map of the H2 S(1) pure rotational transition observed with Spitzer IRS tracing gas excited in outflows and photon-dominated regions in NGC 1333. Positions of Class 0 and I sources in the star-forming region are indicated with filled circles and triangles, respectively.Filled star-symbols mark the position of the B stars 03269+3111 and BD+30459 in the field. Dashed frames delineate the extent of the regions mapped with Herschel HIFI and SPIRE that are discussed in the following sections. The extent of the PACS maps presented in Paper I is also shown for reference.

In the text

Fig. 2 Characteristic examples of H2O 110 − 101 line profiles from different regions across NGC 1333 observed with HIFI. Velocity scale is shifted so that the source velocity is at 0 km s−1 From bottom to top: red-shifted lobe to the north of IRAS 4A - IRAS 4B complex (black), blue-shifted emission to the south of IRAS 2A (blue), on-source emission from IRAS 4A (green), and from the surroundings of the B-star 03260+3111 to the north(red).

In the text

	Fig. 3 HIFI maps of the broad and narrow 110−101 H2O line components (left and central panels, respectively). The right-hand panel presents a velocity map corresponding to the Gaussian centroid-shift for the broad component that reflects the motion for the bulk of the outflowing gas.
In the text

	Fig. 4 SPIRE spectral-line maps of the 12CO ladder observed with the long- and short-wavelength modules (upper and central panels, respectively). 13CO maps observed with the long-wavelength module are presented at the lower panel.
In the text

Fig. 5 SPIRE maps of the lowest energy transitions of [C I] observed with the long-wavelength module. Upper panel: map corresponding to the fundamental transition [C I] (3P0–3P1) at 610 μm; lower panel: map corresponding to the first transition (3P1–3P2) centered at 370.4 μm (~492 and 809 GHz respectively). To enable comparisons, the [C II] emission from Paper I is superimposed in the lower panel (orange contours).

In the text

	Fig. 6 SPIRE map with superposition of the elliptical regions adopted from Plunkett et al. (2013), which are used to calculate the properties of outflows in a consistent manner.
In the text

Fig. 7 Left: excitation diagram for H2; observations(open circles) are fit to calculate the temperature and column density after dereddening (open squares) The distribution of data points can be well approximated with a single temperature component. Right: excitation diagram for CO, where the distribution of data points is best fit with a two distinct temperature components. Excitation diagrams for all outflow regions are presented in Fig. A.3.

In the text

	Fig. 8 Left: [CI] line ratio vs. temperature diagnostic calculated with RADEX for two hydrogen atom densities. Right: temperaturemap reconstructed from the [CI] maps presented in Fig. 4 and the use of the diagnostic.
In the text

	Fig. A.1 Channel maps of the water line displaying the integrated emission at different velocity intervals which are indicated onthe top right corner of each panel.
In the text

	Fig. A.2 Map of a narrow line, possibly from either acetone ((CH3)2CO) or methyl formate (CH3OCHO). Being excited in the vicinity of embedded sources, it may trace hot-corino activity.
In the text

	Fig. A.3 Excitation diagrams for the H2 fluxes extracted from Spitzer within the elliptical regions.
In the text

	Fig. A.3 continued.
In the text

	Fig. A.4 Excitation diagrams for CO fluxes extracted from the SPIRE maps.
In the text

	Fig. A.4 continued.
In the text