© ESO, 2017

1. Introduction

In the early protostellar stages, fast jets powered by the nascent star, possibly surrounded by a wider angle wind, strongly interact with the parental medium through molecular shocks, producing a slower moving molecular outflow cavity. Several episodes of bow shocks are usually present within the outflow lobes (e.g. Bachiller & Pérez Gutiérrez 1997; Bourke et al. 1997; Lee et al. 2006). In several outflows, fast and well collimated molecular bullets, tracing internal shocks within the jet are also seen (e.g. Gueth & Guilloteau 1999; Nisini et al. 2007; Tafalla et al. 2010). Outflow shocks significantly affect the gas chemical composition, with abundance enhancements up to several orders of magnitude for various molecular species such as SiO, CH₃OH, CS, H₂CO, and H₂O (e.g. Bachiller et al. 2001; Garay et al. 1998, 2002; Benedettini et al. 2007; Arce et al. 2007; Codella et al. 2010). Such enrichment reveals the activation of endothermic reactions by shock heating and the mechanical erosion and processing of grains. Therefore, outflow shocks also offer a unique laboratory for testing astrochemical models. However, up to now enough observations to allow detailed chemical studies that include multiple molecular species have been made in only a handful of outflows.

Fig. 1 Double outflow of BHR71 breaking out from its parental Bok globule. The greyscale image is from the DSS2-red acquired with the UK Schmidt Telescope (c◯1993–1995 by the Anglo-Australian Telescope Board). Blue and red contours are respectively CO (3–2) and H2 0–0 S(5) emission (G15). Black stars indicate the position of the two protostars – IRS1 and IRS2 – that drive the outflows. The white and green triangles indicate the position of the peaks of the Herschel line emission and are labelled from A to E in order of decreasing declination. Yellow rectangle shows the approximate extent of the large-scale Herschel-PACS maps.The two cyan crosses indicate the positions observed with additional pointed Herschel observations.

Open with DEXTER

One interesting target for studying shocks associated with outflows is the BHR71 system. Located at a distance of about 175 pc (Bourke et al. 1995), the Bok globule BHR71 is a well-known example of an isolated star formation with a clear bipolar outflow close to the plane of the sky observable from the southern hemisphere (see Fig. 1). The region comprises two distinct outflows (Bourke et al. 1997, 2001; Parise et al. 2006) powered by two protostellar sources, IRS1 and IRS2 which have luminosities of 13.5 and 0.5 L_⊙, respectively (Chen et al. 2008) and are separated by ~3400 AU (assuming a distance of 200 pc, Bourke et al. 2001). A first view of the global energetics and physical conditions of the outflows have been derived by low J_up  CO and H₂ observations (Neufeld et al. 2009; Giannini et al. 2011). Bright HH objects, HH320 and HH321 (Corporon & Reipurth 1997), and several other knots of shocked gas have been found in the outflows. These knots aroused a great deal of interest, and their shock conditions have been studied in several papers (Giannini et al. 2004; Gusdorf et al. 2011, 2015, hereafter G15). The two HH objects have been imaged also in the S ii transition at 6711 Å, indicating that at least part of the outflows are dissociative (Corporon & Reipurth 1997). In fact, suggestions of the presence of an atomic jet have been recently found by [O i]  observations (Nisini et al. 2015).

In this paper we present the first far infrared images of the BHR71 outflows system observed with Herschel-PACS, mapping the region of 200′′  centred on the powering sources, plus some pointed spectra of CO transitions observed with Herschel-HIFI. The images trace the not-yet-studied warm gas, at intermediate excitation conditions with respect the colder gas traced by low J_up  CO lines (Parise et al. 2006; G15) and the hot gas traced by H₂ rotational lines (Neufeld et al. 2009; Giannini et al. 2011). The main aim of this work is to define the global morphology of the warm gas component and to derive general indication on its physical condition. Moreover, in the two brightest infrared knots we performed a detailed study of the shock conditions by fitting a large set of molecular transitions from different species with shock models. The paper is organised as follows. In Sect. 2 we present the new data acquired with Herschel and the reduction method. The results and the analysis performed on the Herschel data are described in Sects. 3 and 4. In Sect. 5 we present the deeper analysis of the large set of CO transitions available for the brightest knot of the northern lobe. Shock models for this knot an another knot of the southern lobe are presented in Sect. 6. The final conclusions are reported in Sect. 7.

2. Observations and data reduction

We present observations carried out towards the BHR71 outflows system with the PACS (Poglitsch et al. 2010) and HIFI (de Graauw et al. 2010) instruments on-board Herschel (Pilbratt et al. 2010). For the scientific analysis we also make use of already published Herschel and APEX observations. In particular, the map in the [O i] ³P₁–³P₂ line (Nisini et al. 2015) and in the CO (6–5)  line (G15). A summary of the new Herschel observations is presented in Table 1.

2.1. Herschel line maps

The BHR71 outflows system was mapped with PACS in CO (14–13), H₂O (2₂₁–1₁₀), H₂O (2₁₂–1₀₁) and [O i] ³P₀–³P₁ by using the line spectroscopy observing mode. The final extended map was built combining three 47′′ × 47′′  partially overlapping, Nyquist-sampled sub-maps (see Fig. 1). Data were reduced with the HIPE¹ v9.0 package, where they were flat-fielded and flux-calibrated by comparison with observations of Neptune. The flux calibration uncertainty amounts to ~20%, while the Herschel pointing accuracy is ~2′′. Post-pipeline reduction steps were performed to locally fit and remove the continuum emission and to construct integrated line maps (see Fig. 2).

Fig. 2 Integrated line maps from Herschel-PACS of the inner region of the BHR71 outflow system. The bottom right panel shows the H2 0-0 S(5) (Neufeld et al. 2009) and the red (from 0 to 20 km s-1) and blue (from −20 to −9 km s-1) wings of the CO (6–5)  emission (G15). First level is at 3σ, namely: 2.1 × 10-5 erg s-1 cm-2 sr-1  for [O i]  63 μm, 3 × 10-6 erg s-1 cm-2 sr-1  for [O i]  145 μm, 5 × 10-6 erg s-1 cm-2 sr-1  for CO (14–13), 7.5 × 10-6 erg s-1 cm-2 sr-1  for H2O 179 μm, 4.5 × 10-6 erg s-1 cm-2 sr-1  for H2O 108 μm  and 6 K km s-1  for CO (6–5). Steps corresponds to 8% of the maximum of each line up to the 70% of the maximum. Cyan circles indicate the position of the peaks of the Herschel line emission, labelled from A to E in order of decreasing declination plus the protostellar sources IRS1 and IRS2. The two yellow crosses indicate the two positions observed with the additional pointed Herschel observations. The black circles in the bottom right part of each panel show the FWHM of the map.

Open with DEXTER

2.2. Herschel single pointing

Two positions (shown in Fig. 1) were observed in additional spectral lines: CO (5–4), ¹³CO (5–4), CO (9–8), CO (14–13) with HIFI and CO (14–13), (18–17), (20–19), (22–21), OH at 119 μm  and [O i] at 145 μm  with PACS. The two positions correspond to local H₂ peaks, namely the SiO knot (Gusdorf et al. 2011) in the northern lobe (RA(J2000) 12^h01^m36.00^s, Dec(J2000) –65°07′34′′) and knot 8 (Giannini et al. 2004) in the southern lobe (RA(J2000) 12^h01^m38.90^s, Dec(J2000) –65°10′07′′). These observations were performed with a single pointing therefore for each line the HIFI data consist of a single spectrum while the PACS data is composed by 25 spectra combined in a 47′′ × 47′′  map, not Nyquist-sampled with respect to the point spread function (PSF). PACS data were reduced in the same way as the extended maps exposed in the previous subsection by using the chopping-nodding pipeline v10.1.0 within HIPE. HIFI data, both the Wide Band (resolution 1.1 MHz) and High Resolution (resolution 0.25 MHz) Spectrograph backends (WBS and HRS, respectively), were reduced with the standard product generation (SPG) v11.1.0 within HIPE. After consistency check-ups, the horizontal and vertical polarisation spectra are averaged in order to increase the signal-to-noise ratio. The observed lines are well resolved already at the WBS resolution therefore in our scientific analysis we will use those spectra. Calibration of the raw data onto the scale was performed by the in-orbit system, while the spectra were converted to a T_mb scale by adopting a forward efficiency of 0.96 and the main beam efficiencies given in Roelfsema et al. (2012).

In the northern pointing, the CO (14–13)  line has been detected by both PACS and HIFI. There is a difference between the two measurements of ~25% that can be considered as the estimate of the intercalibration error between the two instruments. Therefore, in the following analysis we assume a conservative error of the absolute fluxes of 30% that takes into account both the instrumental calibration error and the error of cross calibration between the different instruments.

3. Results

PACS line maps are displayed in Fig. 2. It is worth noting that the PACS lines are spectroscopically un-resolved, therefore it is not possible to separate the contribution of the red and blue outflow lobes. To this aim, in the same figure we also show the emission of the CO (6–5)  line observed with APEX (G15), integrated over the red and blue wings that shows the relative position of the two outflows. North of the two exciting sources there are the red lobe of the IRS1 outflow and the blue lobe of the IRS2 outflow, and the two lobes are spatially well separated. The situation of the southern lobes is more complicated, with the fainter red lob of the IRS2 outflow overlapping the more extended and brighter blue lobe of the IRS1 outflow.

The emission of the PACS lines appears to be clumpy with several individual peaks surrounded by a low level extended emission. The strongest peak is observed towards the protostellar source IRS1 (L_bol = 13.5L_⊙), the exciting source of the more extended outflow. The second protostellar source IRS2 (L_bol = 0.5L_⊙) is weaker but still visible as a local maximum in the maps. We named the five knots located in the outflow lobes from A to E in decreasing value of declination. A first series of knots, namely A, C, D and E are well aligned along a north-south straight line crossing IRS1 and have a separation of ~40′′  in the northern lobe and ~60′′  in the southern lobe. They likely might indicate the direction of the driving jet emitted by IRS1. Knot B is the only knot clearly associated with the blue lobe of the IRS2 outflow. No clear, strong peak along the straight line connecting B to IRS2 is present in the southern lobe suggesting that in this lobe the emission in the PACS lines should be dominated by the blue lobe of the IRS1 outflow, as also suggested by the CO (6–5)  map. The emission peaks are spatially coincident in all the observed PACS lines and are roughly located at the positions of shocked spots previously identified through H₂ observations (Neufeld et al. 2009). We measured the position of the knots from the peak of CO (14–13)  emission, apart from knot E whose coordinates are measured from the [O i] 63 μm  map. In Table 2 the coordinates of the five outflow knots and the two protostars are listed.

The additional single pointing observations with Herschel PACS and HIFI were acquired towards knot A in the northern lobe and about 15′′  north-east of knot E in the southern lobe, indicated by two crosses in Figs. 1 and 2. The spectra of the four lines observed with HIFI are shown in Fig. 3. Towards knot A, all the four lines are detected with good signal-to-noise ratio, with wings extending up to ~60 km s^-1, a velocity higher than previously reported (~45 km s^-1, G15). The CO (5–4)  line presents a self-absorption feature at the velocity of the environmental cloud (−4.5 km s^-1, Bourke et al. 1995), observed in all CO lines up to the (6–5) (G15). From the ¹²CO (5–4)/¹³CO (5–4)  ratio we measured the optical depth of the ¹²CO (5–4) line that results to be optically thin only in the red wing at km s^-1. The intensity peak of the CO lines moves towards higher velocity for lines with higher J_up, indicating that the low and high J_up  lines are emitted from different layers of the shocked region. In the southern lobe pointing the lines profiles are different, despite the blue wings extending up to a velocity (ν~−60 km s^-1) similar to that observed in the red lobe. In particular, no self-absorption and no peak shifting is observed in the southern tip of the blue lobe. Moreover, the CO (14–13)  was not detected and the line ratio CO (9–8)/CO (5–4)  is lower, indicating that the excitation condition at this position is lower than in knot A. This is also confirmed by the fact that the northern pointing corresponds to a peak of the high excitation transitions observed with PACS while the southern pointing is not centred on a peak. The ¹²CO (5–4)/ ¹³CO (5–4) ratio shows that also in the southern blue lobe the ¹²CO (5–4) line is optically thin only in the wing, at v< −10 km s^-1. The CO (9–8)  line in this point of the southern lobe is at velocity bluer than systemic (ν<−4.5 km s^-1) while in the CO (5–4)  line a small red wing is present. There are several possible explanations for this. The first possibility is the larger beam of the CO (5–4)  line with respect to the CO (9–8), secondly, the excitation condition of the gas in the red lobe of the IRS2 outflow could be insufficiency to populate the CO J_up = 9 level, and a third possiblilty is that the emission at red velocity of the CO (9–8)  line is under the sensitivity level of HIFI spectrograph. In all cases we can conclude that the red lobe of the IRS2 outflow does not contribute to the measured flux of the transitions with high excitation levels.

	Fig. 3 Plot of the lines observed with HIFI towards knot A in the northern lobe (upper panel) and about 15′′  north-east of knot E in the southern lobe (lower panel). Spectra were rebinned at a common spectral resolution of 0.5 km s-1. In the insets a zoom in the lines wings is shown.
Open with DEXTER

Towards the two outflow positions (indicated in Figs. 1 and 2) we also acquired PACS maps of high J_up  CO transitions, namely CO (14–13), CO (18–17), CO (20–19), CO (22–21), and of the OH 119 μm  and [O i]  145 μm  lines. Towards knot A, the CO lines are bright, while the OH and [O i]  are less intense (Fig. 4). All the lines have a similar, partially resolved, roundish form with the peak at the centre of the map and a deconvolved circular size of 18′′, indicating that they are emitted from the same extended gas component. Towards the southern pointing only a ~3σ level CO (14–13)  emission is detected while the other lines are not detected, confirming once again the lower excitation condition of this position.

	Fig. 4 Maps of the four CO lines, OH 119 μm  and [O i]  145 μm  observed with PACS towards knot A in the red lobe of the IRS1-driven outflow. First level and level steps are 9 × 10-18 W m-2, corresponding to 3σ. The crosses mark the central position of the 25 spatial pixels of the PACS FOV. The arrow in the top left panel indicates the direction of the IRS1 protostar.
Open with DEXTER

4. Analysis of the extended outflow

4.1. The distribution of the warm gas

The far-infrared PACS lines with excitation temperatures between 80 and 580 K allow for the first time to have a complete view of the distribution of the warm gas component in the BHR71 outflow system, tracing the intermediate excitation condition between the cold gas traced by low J_up  CO lines (Parise et al. 2006; G15) and the hot gas traced by H₂ (Neufeld et al. 2009; Giannini et al. 2011). The emission in the PACS lines shows the presence of seven bright knots surrounded by a low level emission. Two knots correspond to the two protostars that generate the outflows system and the other five knots are close to clumps of high temperature gas ( K) already observed in H₂ rotational transitions (see Fig. 2). They are also associated with the gas moving at the highest velocity as shown by the comparison of the [O i]  63 μm  map with the map of the reddest (from 10.5 km s^-1  to 40 km s^-1) and bluest (from −30 km s^-1  to −19.5 km s^-1) velocities of the CO (6–5)  line wings (Fig. 5). They most likely represent shock spots where the fast jet that drives the outflows impact against the lower expanding gas.

Since the PACS lines have coarse spectral resolution, we cannot spectroscopically separate the contribution of the two outflows to the total emission. This is not a problem for the northern lobes since the two outflows are spatially well separated but it could be a problem for the southern lobes that partially overlap. However, the spectral profile of the CO lines observed with HIFI at the tip of the southern lobe close to knot E, shows that the percentage of the emission at velocity redder than the systemic velocity (−4.5 km s^-1) with respect to the total emission decreases for increasing J_up, being only the 3% for the CO (9–8)  transition. Therefore we can expect an even lower percentage for the CO (14–13)  line and for the other PACS lines from [O i] and water, that trace the same warm gas component, as shown from Fig. 2 for BHR71 and as also observed in other outflows (Santangelo et al. 2013; Busquet et al. 2014). In conclusion, at the tip of the southern lobe the contribution of the red lobe of the IRS2 outflow to the PACS lines is negligible. The situation in the southern positions closer to the two protostars could be different since in the uppest J_up  available CO line, the CO (6–5)  observed with APEX, the emission in the red wing is still significant. However, towards the bright knot D the similar morphology of the blue wing of the CO (6–5)  and the PACS lines (Fig. 2) shows that even in this knot the contribution of the red lobe of the IRS2 outflow to the PACS lines should be negligible.

In Sect. 3 we already noticed that towards knots A the peak velocity of the CO lines moves to higher velocity for lines with higher J_up. In Fig. 6 we compare the spatial morphology of CO (3–2), CO (6–5) and CO (14–13) emission in knot A and we find that the peak of the CO emission moves towards north for CO lines with higher J_up, with a separation of about 20′′  between the CO (3–2)  and the CO (14–13)  peaks. Moreover, the peak of the highest excitation CO line is the closer to the peak of H₂ 0-0 S(5) (E_up = 4586 K) emission tracing hot gas. The peak shifting is along the straight-line connecting knots A, C and D and crossing IRS1, suggesting that the shift is linked to the bow shock generated by the impact of the high-velocity jet against the slower moving ambient shell. The situation in the bright knot D of the southern lobe is different (see Fig. 6), with the low J_up  CO (3–2) line having the classical morphology of the cavity walls and the higher J_up  CO (6–5) and (14–13) showing a roundish structure with spatially coincident peaks.

A displacement in the emission peak of the CO lines with increasing J_up  was also observed towards the flow containing the Herbig-Haro object HH54 (Bjerkeli et al. 2014). In HH54, however, the gradient goes in the direction of the exciting source, that is, in the opposite direction with respect what we observe in knot A of BHR71. Bjerkeli et al. (2014) interpreted the observed displacement of the CO lines peak as the consequence of the presence of an under-dense, hot clump upstream of the HH54 shock wave, traced by a well defined distinct gas component in the spectra of the low J_up  CO lines. Although no distinct velocity components are present in the CO spectra of the BHR71 outflow, the observed displacement of the location of the emission maximum of CO with increasing J_up  in our source could be attributed to different layers of shocked gas impacting an ambient medium with a significant stratification of particle density. The difference in the morphology of the CO lines emission between knots A and D is probably due to the different environmental conditions in which the two shocks have been originated and/or to geometrical effects that could change the angle of view under which the two shocks are observed. In fact, as clearly visible in Fig. 1 that compares the outflow traced by CO (3–2)  with an optical image of the region at 0.54 μm, the outflow is breaking out from its parental Bok globule with more than half of the blue southern lobe projected outside the obscured region while in the red northern lobe only the very final part is projected outside the globule. Although it is more difficult to understand what is really happening on the rear side of the globule, the measured much larger mass and kinetic energy of the red lobe with respect the blue lobe (Bourke et al. 1997) confirms what Fig. 1 suggests. Our Herschel-PACS map covers the part of the northern lobe still fully embedded in the obscured and denser medium inside the molecular dark cloud while the final part of the mapped portion of the southern lobe is propagating outside the cloud. In this respect knots A and D are in a different environment: the northern knot A is located completely inside the cloud while the southern knot D is located at the border of cloud, likely in a less dense medium. Moreover, while knot A is associated to only the red lobe of the IRS1 outflow in the position of knot D two different lobes coexist.

Remarkably, towards knot A the peak of the high J_up  CO corresponds to the peak of atomic lines such as [O i]  63 μm  and [Fe ii] 26 μm  and of the H₂ rotational lines (Fig. 7) as also observed in other bright, outflow shocks such as L1157-B1 and HH54 (Benedettini et al. 2012; Bjerkeli et al. 2014). This implies that the presence of a dissociative and ionising shock but also that the jet material is already partly molecular, or molecules reformation is efficient after the passage of the shock front. A model of the shock towards knot A and D is presented in Sect. 6.

	Fig. 5 Map of the high velocity wings of CO (6–5)  overlaid to [O i] 63 μm. The red wing (red line) is integrated from 10.5 km s-1 to 40 km s-1 and the blue wing (blue line) from –30 km s-1 to –19.5 km s-1. CO levels are at 50%, 70% and 90% of the maximum. The position of the two protostars are indicated with a triangle.
Open with DEXTER

	Fig. 6 Map of the H2 0-0 S(5) line (grey-scale map) compared with CO (3–2) (red line), CO (6–5) (black line) CO (14–13)  (green line). Upper panel is centred on knot A, lower panel on knot D. Contours start at 90% of the peak value in order to highlight the position of the peak of the emission. The circles show the HPBW of the different lines. The arrows indicate the direction of IRS1.
Open with DEXTER

	Fig. 7 Map of the H2 0-0 S(5) line compared with CO (14–13)  (red line), [O i]  63 μm  (green line), [Fe ii] 26 μm  (blue line) towards knot A. Contour levels is 70%, 80% and 90% of the local maxima of each line. The circles show the HPBW of the different lines.
Open with DEXTER

4.2. Line ratios

With only two lines per species it is impossible to constrain the physical parameters of the emitting gas. However, we can derive some general indications on if and how the physical and excitation condition of the gas change along the outflow, comparing the ratios of lines of the same species in the seven infrared knots. In order to avoid different beam filling effects we smoothed all the PACS maps at the common spatial resolution of 13.̋2 and extracted the brightness in the position of the seven Herschel knots. In Fig. 8 we show the H₂O 179 μm/108 μm, the CO (6–5)/(14–13) and the [O i]  63 μm/145 μm  ratios in the seven knots along the outflow. The ratio between the two water lines is quite constant in all the knots but IRS1 where it is higher of about a factor of two. A higher H₂O 179 μm/108 μm  is an indication of a lower particle density or/and lower water column density. IRS1 is also the knot with the lower CO (6–5)/(14–13) ratio which indicates a much higher temperature (T> 1200 K) for this source. This is not unexpected since towards very young protostars a hot CO component was usually observed, traced by CO lines with excitation temperature higher than 1500 K (e.g. Dionatos et al. 2013; van Kempen et al. 2010; Karska et al. 2013). In particular, towards IRS1 Karska et al. (2013) detected the CO (24–23) line, indicating high temperature for the gas surrounding the protostar. The CO (6–5)/(14–13) ratio shows a decreasing trend from the northern lobe knots to the southern lobe knots and this is indicative of a lower temperature at the southern tip of the outflow that confirms our results on the lower excitation condition presented in the previous sections.

The ratio of the two [O i]  lines gives an estimate of the hydrogen particle density. We compared our observed ratios with theoretical line intensity ratios calculated from the non-LTE code also used in Nisini et al. (2015) that considers the first five levels of oxygen. The three main peaks of the IRS1 outflow have ratio higher than 22, indicative of density between 10⁴–10⁵ cm^-3  while on sources and at the southern peak of the outflow (knot E) we measured a slightly lower ratio (15), indicating that in the latter positions the density could be lower. It is worth noting that the lowest ratio of 6.5 observed towards IRS2 is not compatible with model predictions. Unpredicted low [O i]  63 μm/145 μm  ratios have been already observed in several sources (Liseau et al. 2006) and interpreted as due to self-absorption of the 63 μm  line. This interpretation was corroborated by the recent observation of a spectroscopically resolved [O i]  63 μm  line towards a high-mass source which shows a very prominent absorption (Leurini et al. 2015). The self-absorption of the [O i]  63 μm  line could be the reason of the extremely low 63 μm/145 μm  ratio observed towards IRS2. Indeed IRS2 has been classified as a Class 0 protostar and it is likely at an earlier evolutionary stage with respect to IRS1 (Chen et al. 2008) therefore it should be deeply embedded in its dusty envelope whose extinction capability is still efficient at 63 μm  but it does not affect lines at higher wavelengths. We cannot exclude the possibility that the [O i]  63 μm  line towards IRS1 is also partially absorbed therefore the density inferred form the [O i]  lines ratio should be considered a lower limit.

5. Modelling of the CO lines in knot A

For the brightest point of the red lobe of the IRS1 outflow (knot A) we have measurements of CO lines from J_up = 3 to J_up = 22 that allows us to make a deep analysis of the physical conditions present in this pure shock position. The region has already been modelled in several papers (Gusdorf et al. 2011; G15), here we add the new contribution of the high J_up  CO lines observed by PACS that trace a gas component not traced by the transitions modelled by previous works.

5.1. Rotational diagram

We calculated the rotational temperature from the rotational diagram of the CO lines, which gives a lower limit to the kinetic temperature if the gas is not in local thermodynamic equilibrium (LTE). In order to enlarge our data set as much as possible we also used the low J_up  CO line observed with APEX and SOFIA (G15). We smoothed all the CO maps at the common resolution of 24′′. Four CO lines, namely CO (5–4), CO (9–8), CO (11–10)  and CO (16–15), however, have been observed only with a single pointing and with a different beam (except the CO (11–10)  that was observed with a beam of 24′′) and are not used in the fitting. In building the rotational diagram we took into consideration the measured size of the CO emission, assuming that lines with J_up  from 3 to 11 are extended over the beam of 24′′ while the lines with J_up  from 14 to 22 have a size of 18′′  (see Sect. 3). For the spectrally resolved lines observed with APEX and SOFIA the line flux was calculated integrating the emission in the velocity range between –4.5 km s^-1  and 50 km s^-1  as in G15, while for the unresolved Herschel-PACS lines we used a Gaussian fit to the line profile. We attribute to each line flux an uncertainty of 30% in order to take into consideration the intercalibration error between different instruments. The rotational diagram is shown in Fig. 9.

Two components are clearly present: a low excitation component fitting the low J_up  lines and a higher excitation component fitting the high J_up  lines. To derive the column density and the rotational temperature of the two components we fitted two sets of lines separately, grouping the lines measured with the same instrument, i.e. the APEX (J_up = 3, 4, 6, 7) and the PACS lines (J_up = 14, 18, 20, 22), respectively. Indeed the point at which the slope of the fitting stray-line changes is between the two sets of transitions. The best fit parameters are: T_rot = 68 K and N(CO)= 6 × 10¹⁶cm^-2  for the low J_up  APEX lines and T_rot = 267 K and N(CO)= 6 × 10¹⁵cm^-2  for the high J_up  PACS lines. Note that the HIFI and the SOFIA measurements, not used in the fitting itself, are in agreement with the fitting results. It is worth noting that the low J_up  lines, at least the ones with J_up ≤ 5, are not optically thin at all velocities so that column density derived from the rotational diagram represents a lower limit and the rotational temperature an upper limit. The two rotational temperatures are similar to those found in other shock positions along protostellar outflows as L1157 and Cep E (Benedettini et al. 2012; Lefloch et al. 2012, 2015).

	Fig. 8 PACS lines ratios of H2O 179 μm/108 μm  (left panel), CO (6–5)/(14–13) (middle panel) and [O i]  63 μm/145 μm  (right panel) in the seven infrared knots. Arrows represent lower limits.
Open with DEXTER

	Fig. 9 Rotational diagram of CO lines in knot A. For the low Jup  lines the derived rotational temperature and column density are values averaged in the 24′′  beam, assuming a filling factor 1. For the high Jup  lines a size of 18′′  is assumed.
Open with DEXTER

These two different gas components show a different line profile. In fact, the spectral profile of the CO (14–13)  and CO (16–15)  lines is well described by a single exponential law I_CO(ν) = I_CO(0) exp(–ν/ ν₀) with the same slope ν₀ = 10 km s^-1  for the two lines (Fig. 10), while the lower J_up lines have a more complex profile indicating that different gas components contribute to the total emission. A similar slope (ν₀ = 12 km s^-1) has been found for the profile of the high J_up CO lines in the B1 shock of the L1157 outflow (Lefloch et al. 2012) and also these authors attribute the emission of the CO line with J_up  13 to shock excited gas. Another indication that the high J_up  CO lines are dominated by the high excitation gas related to the shock is that also the line profile of SiO (Gusdorf et al. 2011), a species introduced in the gas phase by the shock through sputtering from the dust grains or vaporisation in grain-grain collisions (e.g. Guillet et al. 2009) is well described by the same exponential law (Fig. 10).

5.2. LVG analysis

	Fig. 10 Exponential fit of the line profile for CO and SiO lines towards knot A. The line represent the fit to the line profile with the exponential low I(ν) = I(0) exp(–ν/10 km s-1).
Open with DEXTER

We compared the CO lines emission with a grid of Large Velocity Gradient (LVG) models calculated with the Ceccarelli et al. (2003) code. The molecular data were taken from the BASECOL² database (Dubernet et al. 2013) and the line width (FWHM of the line profile) was set to a fixed value of 9 km s^-1, as measured in the HIFI spectra. We consider the first 41 levels corresponding to energy levels up to 2000 K. The explored range in CO column density is from 4 × 10¹⁵ to 1 × 10¹⁸ cm^-2, the range in particle density is from 5 × 10³ to 3 × 10⁸ cm^-3, the range in temperature is from 25 to 2000 K.

It is impossible to fit all the observed CO lines from J_up = 3 to J_up = 22 with a single gas component, while a two-component model gives a good fit to the data. In order to reduce the free parameters we fixed all the parameters that are well constrained by the observations. In particular, the size of the two components can be measured from the maps: the low J_up component is extended, that is, we assume a filling factor of 1, and the size of the high J_up transitions is measured from the PACS maps to be 18′′. The CO column density of the low J_up lines is well constrained from the fitting of the APEX lines at N(CO)~ 6 × 10¹⁶ cm^-2, a value consistent with that derived from the rotational diagram. If we fix these three parameters the fitting with the LVG grid of models constrains the other model parameters, as given in Table 3. The first gas component fitting the low J_up lines is extended (size > 24′′), cold (T ~ 80 K) and dense (n(H₂) = 3 × 10⁵–4 × 10⁶ cm^-3); the second component fitting the high J_up lines is compact (size = 18′′), warm (T = 1700–2200 K) and with a slightly lower density (n(H₂) = 2 × 10⁴–6 × 10⁴ cm^-3). The comparison between the observed CO line fluxes and the best fit model is shown in Fig. 11. The temperature of the warm gas derived from CO is in agreement with the temperature derived from the H₂ lines (Giannini et al. 2011) confirming that high J_up CO lines and rotational H₂ lines trace the same gas component as also suggested from the spatial coincidence of their emission (Fig. 7). The column density and temperature of the first colder component is in agreement with the estimates derived in the previous section with the rotational diagram, indicating that it is not far from LTE. The temperature is also compatible with the lower limits of the kinetic temperature of the lower excitation gas derived from low J_up  CO lines in other positions along the outflow (Parise et al. 2006).

6. Comparison with shock models

In this section, we present comparisons of our observations with outputs of the Paris-Durham shock code (Flower & Pineau des Forêts 2015) aimed at constraining the physical and chemical conditions of the two brightest shock positions: knots A and D. To this purpose we used a large dataset of emission lines: we started from H₂Spitzer-IRS observations (Neufeld et al. 2009; Giannini et al. 2011) and CO lines from APEX (G15) and our new Herschel-PACS and HIFI observations, gradually including the new Herschel-PACS lines from other species: [O i], OH, and H₂O. The comparison method requires the combination of a shock model with a radiative transfer code based on the LVG approximation and was originally introduced in Gusdorf et al. (2008). Because the datasets available differ in knot A and D, and also because of the different shock structures, we used a slightly different comparison method for each position.

Fig. 11 Best fit of the LVG model for CO lines in knot A considering two components. The parameters of the first component (dotted line) are: N(CO)= 6 × 1016 cm-2, size = 24′′, T = 70 K, n(H2) = 4 × 106 cm-3: they are the physical parameters of the low excitation gas component averaged inside the 24′′  beam centred towards the knot A position; this component is, however, more extended than the beam. The parameters of the second component (dashed line) are: N(CO)= 3 × 1015 cm-2, size = 18′′, T = 2200 K, n(H2) = 4 × 104 cm-3. The solid line represents the sum of the two components.

Open with DEXTER

Fig. 12 Comparison between the best-fit model (red points) obtained in G15 from comparisons of a grid of models from the Paris-Durham shock code with observations (black squares) of H2 (from the Spitzer-IRS telescope) and CO (from the APEX telescope). The target position is the northern shock called SiO-knot in G15 and roughly corresponding to knot A of the current study. Upper panel: H2 excitation diagram. Middle panel: CO integrated intensity diagram. Lower panel: [O i]  line flux. The newly presented CO (5–4, 9–8, 14–13, 18–17, 20–19, and 22–21) transitions observed with Herschel HIFI and PACS have been added to the previously fitted APEX CO lines in the intensity diagram and are satisfactorily fitted by the model.

Open with DEXTER

6.1. Knot A

The shock structure in knot A was already analysed in Gusdorf et al. (2011) and G15. The observational dataset consisted of lines from H₂ 0–0 S(i) (with i = 0 to 7), CO (J_up = 3, 4, 6, 7, 11, 16) and SiO (J_up = 5, 6, 8). More specifically, the observable associated with each molecule was an excitation diagram for H₂, and an integrated intensity diagram for CO. These were both used to constrain the input parameters of the shock models. The integrated intensity diagram of SiO was used in a second step to fine-tune our constraints on the pre-shock relative distribution of silicon in the mantle and core of interstellar grains. The CO and SiO spectra pointed towards the existence of red-shifted gas, consistent with a single shock structure. In this region, one non-stationary shock model was consequently found to reproduce observations of H₂, CO and SiO, with the following input parameters: pre-shock density of n(H)= 10⁴ cm^-3, shock velocity of ν_s = 22 km s^-1, magnetic field parameter b = 1.5 (defined by B(μG)= b × [ n(H)(cm^-3)]^1/2, where B is the intensity of the magnetic field transverse to the shock direction of propagation), and age of ~3800 yr.

The first step of our new analysis in knot A was to add the integrated intensities of newly presented CO (5–4, 9–8, 14–13, 18–17, 20–19, and 22–21) transitions observed with Herschel HIFI and PACS to the CO modelled in G15. We then repeated the analysis done in G15 by comparing the CO and H₂ observables to outputs from the same grid of models and calculating a chi-square. We found that the best-fit model from G15 still provides the best fit to this more complete dataset. The corresponding results for CO and H₂ are shown in Fig. 12.

As we also dispose of [O i]  observations by Herschel-PACS (³P₁–³P₂ at ~63 μm and ³P₀–³P₁ at ~145 μm) and we have compared their emissivities to the predictions of the best-fit model found fitting the CO and H₂ emission. We did not include the [O i]  lines in the best fitting process because contrary to CO lines they are not velocity-resolved and contrary to H₂ they might be subject to self-absorption (especially the lower-lying one, see e.g. Leurini et al. 2015, for a discussion based on velocity-resolved observations in a different region). In particular, the fact that self-absorption might affect the lower lying transition justifies that the corresponding flux should be considered as a lower limit. An additional reason to justify why we did only a posteriori comparison between models and observations of [O i]  lines is the relatively simple assumptions under which their emissivities are modelled: optically thin and in LTE conditions. While it is difficult to diagnose the first optical thickness condition, an indication for the LTE assumption is given by the critical density of the [O i]  63 μm, of the order of 10⁶ cm^-3, likely difficult to fulfil in the analysed shock region. The result of the comparison of observed [O i]  lines with the best fit model can be seen in the lower panel of Fig. 12. For both lines, we found a huge discrepancy between observations and models. Indeed the modelled fluxes are 1% and 2% of the observed values for the 63 μm  and 145 μm  line, respectively. Two explanations can be advanced for this discrepancy: either (i) there is a single shock structure in the observed beam, and the physical ingredients of the Paris-Durham model are not appropriate to describe the physics and chemistry at work in this region; or (ii) a second shock structure exists, strongly emitting in [O i]  and possibly negligibly in CO or H₂, such as an atomic jet component. In the first case (i), for example, an enhancement of [O i]  and OH can be the consequence of the photo-dissociation of H₂O. The photo-dissociation of molecules can be due to the presence of an UV radiation field that is not taken into account in the version of the Paris-Durham model that we used. However, we would expect this effect to have also a strong impact on the predicted CO and H₂ line emission (namely significantly reducing this emission, as these molecules would also undergo dissociation to some extent). Moreover knot A is about 1.3 × 10⁴ AU away from IRS1 making the possible proto-stellar radiation field practically non influential at this position. This is why the second assumption (ii) is a more convincing scenario. Indeed, it is possible that an atomic jet component exists, that would simultaneously be responsible for the [O i]  emission and would not modify the H₂ and CO comparisons that have been presented in G15 and here. For instance, a dissociative J-type shock propagating in a less dense medium could generate the observed amounts of [O i]  emission, bring a negligible addition in terms of CO line emission, and a weak one in terms of H₂ emission. This would be the case, for example, of the stationary model with n(H)= 10³ cm^-3, b = 1.5, and ν_s = 30–40 km s^-1  presented in Hollenbach & McKee (1989). In such models, the dissociation of molecules is generated by the shock itself. Alternatively, a purely atomic component could be considered, for which models are lacking because of the difficulty to include the presence of self-generated or externally-generated radiation field. A situation similar to what we observe in knot A of the BHR71 outflow has been observed also in other shocked clumps of the L1157 and L1148 outflows (Benedettini et al. 2012; Santangelo et al. 2013) where it was found that the bright emission of the [O i]  lines cannot be reproduced by only a C-type shock or a low-velocity, non-dissociative J-type shock but requires the presence of an additional dissociative shock that could be associated to the primary jet. Additional evidence of the existence of a fast driving jet were found by Herschel near the protostar where high velocity gas components were detected, likely produced by a collimated jet very close to the ejecting protostar (Kristensen et al. 2013; Nisini et al. 2015). More observations, with a better spatial and/or spectral resolution, in particular of atomic lines such as [O i], are necessary to unveil the primary jet in outflows and to confirm our current interpretation of [O i]  and OH emission in the shocked spots along the lobes.

Finally, we also compared the prediction of the best-fit model to our Herschel-PACS observations of OH (lines at 119.2 μm  and 119.4 μm) and H₂O (lines at 108 μm  and 179 μm). We note that the OH lines and the o-H₂O 179 μm  line might suffer from self-absorption, hence the observed values should be treated as lower limits (e.g. Leurini et al. 2015). We found that the best-fit model significantly underestimates the integrated intensity of the OH 119.2 μm  and 119.4 μm  lines (models: 1.0 × 10^-3 and 1.8 × 10^-3 K km s^-1, observations: 0.25 and 0.36 K km s^-1, respectively), similarly as what we obtained for [O i]. Conversely, the 108 μm  and 179 μm  water lines comparisons provided a much smaller difference: 2.5 K km s^-1  and 41.0 K km s^-1  respectively for the models, and 1.2 K km s^-1  and 8.4 K km s^-1  for the observations. We might be able to further decrease the discrepancy between observations and models for water comparisons by simultaneously fine-tuning the physical and chemical networks in our models, and by using more, preferably velocity-resolved observations of H₂O lines. This will be the subject of a dedicated work in a forthcoming publication.

6.2. Knot D

	Fig. 13 CO ladder obtained in knot D with the APEX telescope by G15, convolved to the 24′′  resolution. The nominal spectral resolution was given in G15.
Open with DEXTER

	Fig. 14 Channel-by-channel optical depth of the CO (3–2) line (right-hand y-axis), obtained under two assumption for the [12CO]/[13CO] relative abundance: 50 and 60 (empty and full red circles), for all channels with at least a 3σ detection in 13CO (3–2). Also shown the 12CO and 13CO (3–2) spectra (in black and pink respectively, left-hand y-axis), obtained in knot D with the APEX telescope.
Open with DEXTER

The analysis of knot D is slightly different. Indeed, in this region, the southern, blue-shifted lobe of the outflow driven by IRS1 coexists with the southern, red-shifted lobe of the outflow driven by IRS2 (see Fig. 2). This can also be seen in Fig. 13, that shows all spectrally resolved spectra obtained with the APEX telescope (G15) towards knot D, convolved to a spatial resolution of 24′′  (for consistency reasons with the analysis in knot A). All spectra indicate the presence of blue-shifted gas up to about –35 km s^-1  and of red-shifted gas to about 15 km s^-1. Also shown in Fig. 13 is the APEX ¹³CO (3–2) line G15. We used this line in combination with ¹²CO (3–2) to estimate the optical thickness of the latter line. Figure 14 shows spectra of the ¹³CO and ¹²CO (3–2) lines at the knot D position at the same spectral and spatial resolutions. These spectra allow us to calculate the opacity for each velocity channel where the ¹³CO is detected at more than 3σ. Assuming an identical excitation temperature for the ¹²CO and ¹³CO lines, and a typical interstellar abundance ratio of 50–60 (e.g. Langer & Penzias 1993), we conclude that the ¹²CO (3–2) emission is optically thick at least in the low-velocity regime of the spectral wings. The large optical thickness is consistent with the constancy of the line integrated intensities (in the non-absorbed components) from CO (3–2) and (4–3), when convolved to the same resolution.

Fig. 15 Comparison between the best-fit model (coloured points) from comparisons of a grid of models from the Paris-Durham shock code with observations (black squares) of H2 (from the Spitzer-IRS telescope) and CO (from the APEX telescope). The target position is the southern knot D shock where the blue lobe of the outflow driven by IRS1 overlaps with that driven by IRS2. We consequently fitted two shock layers in this position: one on the red-shifted, one on the blue-shifted gas component. Upper panel, left: CO integrated intensity diagram for the blue component. Upper panel, right: CO integrated intensity diagram for the red component. Lower panel, left: H2 excitation diagram, the three empty symbols correspond to fluxes extracted with different methods (see text for details), the filled squares indicates the average values. Lower panel, right: [O i] line flux.

Open with DEXTER

The presence of blue- and red-shifted gas in knot D compelled us to vary the modelling method with respect to the analysis performed in knot A. Indeed, instead of trying to fit only one shock layer onto the observations, we attempted to model each blue and red lobe independently, using an approach that was already used in Gusdorf et al. (2012) and Anderl et al. (2014). We hence generated two integrated intensity diagrams for CO by using velocity-resolved line spectra (available for J_up = 3, 4, 6, 7): one for the blue-shifted component (from –35 km s^-1  to –4.5 km s^-1), and one for the red-shifted component (from –4.5 km s^-1  to 20 km s^-1). We then compared each of these excitation diagrams to our grid of models to select the best-fit models for knot D. Before doing this, we first extended our grid of models to a larger parameters space.

Indeed, as can be seen in Fig. 1, the position of knot D lies at the edge of the parental Bok globule of the outflow. For this reason, we ran additional shock models (stationary and non-stationary ones) for a lower pre-shock density value: 10³cm^-3. Additionally, we noticed that the linewidth of the red-shifted gas component was smaller than that of the blue-shifted one. We consequently extended our grid towards lower values of the shock velocity (towards 15 km s^-1 included).

Comparing the velocity-resolved CO lines emission with the extended grids of models, we found that no model with a lower pre-shock density is be able to fit either the blue- or red-shifted gas observations. This is because non-stationary, young models at this pre-shock density are J-type shocks, which are too narrow to generate significant levels of CO emission. Physically, the fact that shocks at the edge of the cloud still propagate in a relatively dense medium can be explained by their propagation in a medium that has already been shock-processed, hence densified by, for example, the primary ejection episodes of the protostar(s).

The best-fit models for each spectral range are shown in the top panels of Fig. 15. The input parameters of the best-fit models are: n(H)= 10⁴ cm^-3, b = 1.0, and ν_s = 20 km s^-1, and an age of about 3900 yr for the blue-shifted CO emission, and n(H)= 10⁴ cm^-3, b = 1.5, and ν_s = 20 km s^-1, and an age of about 1000 yr for the red-shifted CO emission. For both spectral ranges, the best-fit models are non stationary, with pre-shock density and the b value similar to the values found in knot A. For the red-shifted gas component, we chose to display in Fig. 15 two solutions: the absolute best-fit, with a velocity of 20 km s^-1  (red symbols and lines), and another good solution with a velocity closer to the relatively small observed linewidth (15 km s^-1, orange symbols and lines). The difference between the two solutions is not significant. We remind the reader that when comparing 1D shock models with CO observations, one should not necessarily expect that the observed linewidth should be equal to the shock velocity. Indeed it is possible that shocks propagate in an already moving medium, or are orientated towards the plane of the sky resulting respectively in a lower or higher shock velocity.

Finally, we found an evolutionary trend in the best-fit solution: the best-fit for the blue-shifted emission has an age comparable to what was found for the red-shifted gas in knot A, around 3900 yr. Looking at the relative position of the two knots A and D with respect to the exciting source IRS1, one would expect an older age for the more distant knot A than knot D. However, both lines profiles and the best fit shock model suggest that the shock velocity in the southern, blue-shifted lobe should be lower than that in the northern red-shifted lobe. Moreover, the age estimate, as well as all the other input parameters of our shock model, suffer of an uncertainty given by the imperfect coverage of the input parameters space. Therefore, despite the different distance of A and D form the protostar, their age similarity is not too surprising even because both structures belong to the same outflow.

On the other hand, the age associated with the red-shifted gas (from the outflow driven by IRS2) is younger, 1000–1500 yr old. This value is consistent with the dynamical age of this outflow that we estimated from the CO (3–2) map. Indeed, for the outflow driven by IRS2, we found CO (3–2) emission with 10 km s^-1  linewidth up to ~56′′  from the driving source in the northern lobe and ~113′′  in the southern lobe, translating in a dynamical age comprised between 5000 and 11 000 yr (also slightly younger than the dynamical age inferred by G15 for the outflow driven by IRS1). The uncertainty associated to these values is important, however both the observational dynamical ages and results from models seem to suggest a younger evolutionary stage for IRS2. This is consistent also with the findings of Chen et al. (2008) based on the analysis of the spectral energy distribution of both IRS1 and IRS2 and with our results shown in Sect. 4.2.

In order to validate these results, we also plotted the integrated intensity of the CO (14–13) line inferred from our Herschel-PACS measurements. As this observation was not velocity resolved, this value should not be compared to the blue- or red-shifted best-fit, but to their combination. We found that this value is indeed consistent (within the errorbars) with the sum of the integrated intensity of the two best-fit. Similarly, we used the H₂ observations by Spitzer-IRS to further check our model results. We extracted the level populations in various methods already discussed in G15, resulting in the set of symbols shown in the lower left-hand panel of Fig. 15. Since these observations are not velocity resolved, we summed the contribution from the two best-fits and compared them with the observed excitation diagram. The result can be seen in colours in the lower left-hand panel of Fig. 15. Two sets of models are plotted, since we chose to display the sum of the two best-fit contributions (purple symbols and line) and the sum of the contributions from the blue-shifted best-fit with that from the red-shifted best-fit with a low velocity (pink symbols and line). In both cases, the results are fitted with a similar quality as in knot A. We note that the fit is not perfect, and we are working on solutions to improve our models, specially from the point of view of the geometrical description of the observed region. Similarly, we summed either the [O i]  fluxes from our two best-fit models or the blue-shifted gas best-fit with the red-shifted gas best fit with the lower velocity and compared these values to the observations. The results can be seen in the right-hand lower panel of Fig. 15. Similarly to what found for knot A, also for knot D the best fit models drastically underestimate the emission of both [O i]  transitions. We refer to the previous subsection for an explanation of this result.

Finally, in this knot, we have observed H₂O emission from the lines at 108 μm  and 179 μm  with PACS, that is, without spectrally resolving the line profile. We could then compare the values predicted by the sums of our best models: 39.8−45.2 K km s^-1 and 2.5–2.9 K km s^-1, respectively with the observed integrated intensities: 11.7 K km s^-1 and 1.0 K km s^-1. Given the possible self-absorption (see discussion in the above subsection), the observed values are lower limits. We note that a discrepancy between models and observations does exist, but less impressive than for the [O i]  lines, similarly as the conclusions reached in the previous paragraph for knot A, with the same conclusions.

7. Summary and conclusions

We presented Herschel maps in CO (14–13), H₂O (2₂₁–1₁₀), H₂O (2₁₂–1₀₁)  and [O i]  145 μm  of the BHR71 outflows system plus additional CO lines measurements towards two outflow positions. These far infrared maps show seven bright knots surrounded by a low level extended emission. Two knots correspond to the two protostars IRS1 and IRS2 exciting the two outflows composing the BHR71 system while the other knots are located in the outflow lobes and represent shock episodes along the outflow jet. Indeed, the far infrared emission peaks, spatially coincident in all the observed PACS lines, are roughly coincident with the positions of the local maxima of H₂ lines and of the fastest outflowing gas traced by the CO (6–5)  line. In the southern lobe the blue-shifted gas of the IRS1 outflow and the red-shifted gas of the IRS2 outflow partially overlap, however several evidences indicate that the emission in the observed PACS lines is dominated by the flux emitted from gas associated with the blue lobe of the IRS1 outflow.

The ratio of the Herschel-PACS lines intensity, between lines of the same molecular species, have been calculated in the seven knots. The ratios are quite similar within the errors, showing that the excitation conditions of the fast moving gas do not change significantly along the outflow, apart a lower density and temperature at the extremity of the southern blue lobe where the outflow is flowing outside the parental molecular cloud. More peculiar ratios are found towards the two protostars from which we deduced that IRS1 have a much higher temperature (T> 1200 K) and that a significant extinction is present towards IRS2 able to absorb part of the [O i]  63 μm  emission.

In the brightest position of the red lobe of the IRS1 outflow (knot A) we observed additional CO lines with J_up up to 22. Combining our Herschel observations with measurements of the low J_up CO, H₂ and SiO lines (Gusdorf et al. 2011; G15; Neufeld et al. 2009) we were able to make a deep analysis of the physical and shock conditions present in this pure shock position. Rotational diagram, spectral profile shape and LVG analysis of the CO lines showed the presence of two gas components: one extended (size > 24′′), cold (T ~ 80 K) and dense (n(H₂) = 3 × 10⁵–4 × 10⁶ cm^-3) and another compact (18′′), warm (T = 1700–2200 K) and slightly less dense (n(H₂) = 2 × 10⁴–6 × 10⁴ cm^-3). The CO column density of the cold component is about one order of magnitude higher than that of the warm component.

In the two infrared brightest knots, knot A in the northern lobe and knot D in the southern lobe, we compared observations with the Paris-Durham shock code. We found that non-stationary shock models well reproduce observations of CO and H₂ in both knots. In knot A, the best-fit model is non stationary with n(H)= 10⁴ cm^-3, b = 1.5, and ν_s = 22 km s^-1, and an age of about 3800 yr. In knot D, two shocks coexist the blue-shifted lobe of the IRS1 outflow and the red-shifted lobe of the IRS2 outflow. The blue-shifted CO emission can be described well with a non-stationary shock with n(H)= 10⁴ cm^-3, b = 1.0, and ν_s = 20 km s^-1, and an age of about 3900 yr, conditions very similar to that found in the red-shifted counterpart. On the other hand, the red-shifted CO emission from the IRS2 outflow can be reproduced by a non-stationary shock with similar parameter, n(H)= 10⁴ cm^-3, b = 1.5, and ν_s = 20 km s^-1, but a younger age of about 1000 yr. Because of the smaller linewidth associated with this component, we also investigated the possibility to fit the observations with a lower velocity shock with n(H)= 10⁴ cm^-3, b = 1.0, and ν_s = 15 km s^-1, and an age of about 1450 yr. Interestingly, the younger age of the outflow powered by IRS2 is in agreement with the evolutionary estimate derived from the spectral energy distribution fitting of the two driving protostars (Chen et al. 2008) and with the higher obscuration of the IRS2 protostar deduced from the [O i]  lines ratio.

The non-stationary shock models that well reproduce the CO and H₂ observations underestimate by about two orders of magnitude the observed [O i]  lines at 63 μm  and 145 μm, and the OH lines at 119.2 μm  and 119.4 μm, when available. We interpret this discrepancy as the evidence for the existence of a primary, mostly atomic, jet component in the outflow. In fact a shock with a radiative precursor, with n(H)= 10³ cm^-3, b = 1.5, and ν_s = 30–40 km s^-1  as presented in Hollenbach & McKee (1989) can account for the observed [O i]  emission knot A without significantly altering our conclusions with respect to CO and H₂. Our Herschel-PACS observations of [O i]  and OH, however, are both spatially and spectroscopically unresolved, therefore they do not directly probe the jet. Moreover, our shock model does not account for the effects of a possible UV radiation field affecting the shock (whether this UV field is created by the shock itself or emanates from the protostar). It is worth noting that several other Herschel observations both close to protostars and in the outflow lobes have been interpreted as evidence of dissociative shocks produced by the primary jet (Benedettini et al. 2012; Kristensen et al. 2013; Nisini et al. 2015) but only additional observations of species such as OH, [C ii] and [O i], possibly spectroscopically resolved, for example with SOFIA, could directly probe the primary collimated jet in outflows.

Finally, we compared the predictions of the models that best fit the CO and H₂ emission in knots A and D to water observations. We found small discrepancies that might be resolved by further investigations on the shock modelling of water, and their applications to a larger number of velocity-resolved observations of water lines. This will be the object of a forthcoming publication.

Our results in the BHR71 outflow are consistent with other recent studies of outflows that have carried out global analysis of emission lines with different excitation temperatures, for example the CO lines with low J_up observable from the ground together with the intermediate and high J_up lines observed with Herschel (e.g. Lefloch et al. 2012; Benedettini et al. 2012) or large samples of H₂O lines observed with Herschel (e.g. Santangelo et al. 2013; Busquet et al. 2014). These studies have shown the presence of multiple components in the molecular outflowing gas at different temperatures and densities. So far, three components have been detected in most outflows and have been associated by the authors to the following physical origin: a cold component with temperature of about 60–80 K associated to the walls of the cavity formed by the entrained gas, a warm component with temperature of about 200–500 K associated to a non-dissociative shock and a hot component with temperature of about 800–2000 K, also emitting in H₂ rotational lines and atomic lines, associated with a dissociative J-type shock in the compact region where the protostellar jet impacts the less fast gas. In this respect we found that the chemically active BHR71 outflow, as well as L1157, the other well studied chemically active outflow (Benedettini et al. 2012; Lefloch et al. 2012; Busquet et al. 2014), does not show evident peculiarity with respect to less chemically active outflows. This is also true when comparing the luminosities of the main far-infrared coolants and the mass loss rate that instead are correlated to the bolometeric luminosity and mass of the envelope of the exciting protostar (Neufeld et al. 2009; Karska et al. 2013; Tafalla et al. 2013; Nisini et al. 2015).

The results presented in this paper show that the BHR71 outflow is an ideal astrophysical laboratory for studying shocks in protostellar outflows. However, the major limitation of the analysis performed on the large dataset already acquired for this source is the exiguous spatial resolution of the single dish observations. Improving the spatial resolution with interferometric data will allow us to make a step forward in modelling the physics and the chemistry of protostellar shocks. This makes the BHR71 outflow one of the best candidates for future ALMA observations.

Acknowledgments

We thank Prof. D. Neufeld for providing us with the reduced Spitzer data. G.B. acknowledges the support of the Spanish Ministerio de Economia y Competitividad (MINECO) under the grant FPDI-2012-18204. G.B. is supported by the Spanish MICINN grant AYA2014-57369-C3-1-P. PACS has been developed by a consortium of institutes led by MPE (Germany) and including UVIE (Austria); KU Leuven, CSL, IMEC (Belgium); CEA, LAM (France); MPIA (Germany); INAF-IFSI/OAA/OAP/OAT, LENS, SISSA (Italy); IAC (Spain). This development has been supported by the funding agencies BMVIT (Austria), ESA-PRODEX (Belgium), CEA/CNES (France), DLR (Germany), ASI/INAF (Italy), and CICYT/MCYT (Spain). HIFI has been designed and built by a consortium of institutes and university departments from across Europe, Canada and the United States under the leadership of SRON Netherlands Institute for Space Research, Groningen, The Netherlands and with major contributions from Germany, France and the US. Consortium members are: Canada: CSA, U.Waterloo; France: CESR, LAB, LERMA, IRAM; Germany: KOSMA, MPIfR, MPS; Ireland, NUI Maynooth; Italy: ASI, IFSI-INAF, Osservatorio Astrofisico di Arcetri-INAF; Netherlands: SRON, TUD; Poland: CAMK, CBK; Spain: Observatorio Astronómico Nacional (IGN), Centro de Astrobiología (CSIC-INTA). Sweden: Chalmers University of Technology – MC2, RSS & GARD; Onsala Space Observatory; Swedish National Space Board, Stockholm University – Stockholm Observatory; Switzerland: ETH Zurich, FHNW; USA: Caltech, JPL, NHSC. The Digitized Sky Surveys were produced at the Space Telescope Science Institute under US Government grant NAG W-2166. The images of these surveys are based on photographic data obtained using the Oschin Schmidt Telescope on Palomar Mountain and the UK Schmidt Telescope. The plates were processed into the present compressed digital form with the permission of these institutions. The digitized images of the southern sky are copyrightc ◯1993-5 by the Anglo-Australian Observatory Board, and are distributed herein by agreement.

References

All Tables

All Figures

	Fig. 1 Double outflow of BHR71 breaking out from its parental Bok globule. The greyscale image is from the DSS2-red acquired with the UK Schmidt Telescope (c◯1993–1995 by the Anglo-Australian Telescope Board). Blue and red contours are respectively CO (3–2) and H2 0–0 S(5) emission (G15). Black stars indicate the position of the two protostars – IRS1 and IRS2 – that drive the outflows. The white and green triangles indicate the position of the peaks of the Herschel line emission and are labelled from A to E in order of decreasing declination. Yellow rectangle shows the approximate extent of the large-scale Herschel-PACS maps.The two cyan crosses indicate the positions observed with additional pointed Herschel observations.
Open with DEXTER
In the text

	Fig. 2 Integrated line maps from Herschel-PACS of the inner region of the BHR71 outflow system. The bottom right panel shows the H2 0-0 S(5) (Neufeld et al. 2009) and the red (from 0 to 20 km s-1) and blue (from −20 to −9 km s-1) wings of the CO (6–5)  emission (G15). First level is at 3σ, namely: 2.1 × 10-5 erg s-1 cm-2 sr-1  for [O i]  63 μm, 3 × 10-6 erg s-1 cm-2 sr-1  for [O i]  145 μm, 5 × 10-6 erg s-1 cm-2 sr-1  for CO (14–13), 7.5 × 10-6 erg s-1 cm-2 sr-1  for H2O 179 μm, 4.5 × 10-6 erg s-1 cm-2 sr-1  for H2O 108 μm  and 6 K km s-1  for CO (6–5). Steps corresponds to 8% of the maximum of each line up to the 70% of the maximum. Cyan circles indicate the position of the peaks of the Herschel line emission, labelled from A to E in order of decreasing declination plus the protostellar sources IRS1 and IRS2. The two yellow crosses indicate the two positions observed with the additional pointed Herschel observations. The black circles in the bottom right part of each panel show the FWHM of the map.
Open with DEXTER
In the text

	Fig. 3 Plot of the lines observed with HIFI towards knot A in the northern lobe (upper panel) and about 15′′  north-east of knot E in the southern lobe (lower panel). Spectra were rebinned at a common spectral resolution of 0.5 km s-1. In the insets a zoom in the lines wings is shown.
Open with DEXTER
In the text

	Fig. 4 Maps of the four CO lines, OH 119 μm  and [O i]  145 μm  observed with PACS towards knot A in the red lobe of the IRS1-driven outflow. First level and level steps are 9 × 10-18 W m-2, corresponding to 3σ. The crosses mark the central position of the 25 spatial pixels of the PACS FOV. The arrow in the top left panel indicates the direction of the IRS1 protostar.
Open with DEXTER
In the text

	Fig. 5 Map of the high velocity wings of CO (6–5)  overlaid to [O i] 63 μm. The red wing (red line) is integrated from 10.5 km s-1 to 40 km s-1 and the blue wing (blue line) from –30 km s-1 to –19.5 km s-1. CO levels are at 50%, 70% and 90% of the maximum. The position of the two protostars are indicated with a triangle.
Open with DEXTER
In the text

	Fig. 6 Map of the H2 0-0 S(5) line (grey-scale map) compared with CO (3–2) (red line), CO (6–5) (black line) CO (14–13)  (green line). Upper panel is centred on knot A, lower panel on knot D. Contours start at 90% of the peak value in order to highlight the position of the peak of the emission. The circles show the HPBW of the different lines. The arrows indicate the direction of IRS1.
Open with DEXTER
In the text

	Fig. 7 Map of the H2 0-0 S(5) line compared with CO (14–13)  (red line), [O i]  63 μm  (green line), [Fe ii] 26 μm  (blue line) towards knot A. Contour levels is 70%, 80% and 90% of the local maxima of each line. The circles show the HPBW of the different lines.
Open with DEXTER
In the text

	Fig. 8 PACS lines ratios of H2O 179 μm/108 μm  (left panel), CO (6–5)/(14–13) (middle panel) and [O i]  63 μm/145 μm  (right panel) in the seven infrared knots. Arrows represent lower limits.
Open with DEXTER
In the text

	Fig. 9 Rotational diagram of CO lines in knot A. For the low Jup  lines the derived rotational temperature and column density are values averaged in the 24′′  beam, assuming a filling factor 1. For the high Jup  lines a size of 18′′  is assumed.
Open with DEXTER
In the text

	Fig. 10 Exponential fit of the line profile for CO and SiO lines towards knot A. The line represent the fit to the line profile with the exponential low I(ν) = I(0) exp(–ν/10 km s-1).
Open with DEXTER
In the text

	Fig. 11 Best fit of the LVG model for CO lines in knot A considering two components. The parameters of the first component (dotted line) are: N(CO)= 6 × 1016 cm-2, size = 24′′, T = 70 K, n(H2) = 4 × 106 cm-3: they are the physical parameters of the low excitation gas component averaged inside the 24′′  beam centred towards the knot A position; this component is, however, more extended than the beam. The parameters of the second component (dashed line) are: N(CO)= 3 × 1015 cm-2, size = 18′′, T = 2200 K, n(H2) = 4 × 104 cm-3. The solid line represents the sum of the two components.
Open with DEXTER
In the text

	Fig. 12 Comparison between the best-fit model (red points) obtained in G15 from comparisons of a grid of models from the Paris-Durham shock code with observations (black squares) of H2 (from the Spitzer-IRS telescope) and CO (from the APEX telescope). The target position is the northern shock called SiO-knot in G15 and roughly corresponding to knot A of the current study. Upper panel: H2 excitation diagram. Middle panel: CO integrated intensity diagram. Lower panel: [O i]  line flux. The newly presented CO (5–4, 9–8, 14–13, 18–17, 20–19, and 22–21) transitions observed with Herschel HIFI and PACS have been added to the previously fitted APEX CO lines in the intensity diagram and are satisfactorily fitted by the model.
Open with DEXTER
In the text

	Fig. 13 CO ladder obtained in knot D with the APEX telescope by G15, convolved to the 24′′  resolution. The nominal spectral resolution was given in G15.
Open with DEXTER
In the text

	Fig. 14 Channel-by-channel optical depth of the CO (3–2) line (right-hand y-axis), obtained under two assumption for the [12CO]/[13CO] relative abundance: 50 and 60 (empty and full red circles), for all channels with at least a 3σ detection in 13CO (3–2). Also shown the 12CO and 13CO (3–2) spectra (in black and pink respectively, left-hand y-axis), obtained in knot D with the APEX telescope.
Open with DEXTER
In the text

	Fig. 15 Comparison between the best-fit model (coloured points) from comparisons of a grid of models from the Paris-Durham shock code with observations (black squares) of H2 (from the Spitzer-IRS telescope) and CO (from the APEX telescope). The target position is the southern knot D shock where the blue lobe of the outflow driven by IRS1 overlaps with that driven by IRS2. We consequently fitted two shock layers in this position: one on the red-shifted, one on the blue-shifted gas component. Upper panel, left: CO integrated intensity diagram for the blue component. Upper panel, right: CO integrated intensity diagram for the red component. Lower panel, left: H2 excitation diagram, the three empty symbols correspond to fluxes extracted with different methods (see text for details), the filled squares indicates the average values. Lower panel, right: [O i] line flux.
Open with DEXTER
In the text