The shocked gas of the BHR71 outflow observed by Herschel: indirect evidence for an atomic jet

TL;DR

This study uses Herschel infrared observations to analyze shock-excited gas in the BHR71 outflow, revealing complex shock structures and indirect evidence for an atomic jet component.

Contribution

It provides the first detailed infrared imaging and shock modeling of the BHR71 outflow, identifying multiple gas components and suggesting the presence of an atomic jet.

Findings

Detection of warm, shocked gas knots in the outflow

Identification of two gas components: cold extended and warm compact

Evidence for an additional dissociative shock component, possibly a jet remnant

Abstract

In the BHR71 region, two low-mass protostars drive two distinguishable outflows. They constitute an ideal laboratory to investigate the effects of shock chemistry and the mechanisms that led to their formation. We aim to define the morphology of the warm gas component of the BHR 71 outflow and at modelling its shocked component. We present the first far infrared Herschel images of the BHR71 outflow in the CO(14-13), H$_2$O (2$_{21}$-1$_{10}$), H$_2$O (2$_{12}$-1$_{01}$) and [OI] 145 $\mu$m, lines, revealing the presence of several knots of warm, shocked gas associated with fast outflowing gas. In two of these knots we performed a detailed study of the physical conditions by comparing a large set of transitions from several molecules to a grid of shock models. Herschel lines ratios in the outflow knots are quite similar, showing that the excitation conditions of the fast moving gas do not change significantly within the first $\sim$ 0.068 pc of the outflow, apart at the extremity of the southern blue-shifted lobe that is expanding outside the molecular cloud. Rotational diagram, spectral line profile and LVG analysis of the CO lines in knot A show the presence of two gas components: one extended, cold ($T\sim$80 K) and dense ($n$(H$_2$) = 3$\times$10$^5$-4$\times$10$^6$ cm$^{-3}$) and another compact (18 arcsec), warm ($T$ = 1700-2200 K) with slightly lower density ($n$(H$_2$) = (2-6)$\times$10$^4$ cm$^{-3}$). In the two brightest knots (where we performed shock modelling) we found that H$_2$ and CO are well fitted with non-stationary (young) shocks. These models, however, significantly underestimate the observed fluxes of [OI] and OH lines, but are not too far off those of H$_2$O, calling for an additional, possibly dissociative, J-type shock component. Our modelling indirectly suggests that an additional shock component exists, possibly a remnant of the primary jet

Tables

Table 1. Table 1: Observations towards the BHR71 outflow.

220′′$\times$100′′, spatially Nyquist-sampled maps
Species	transition	$E_{\mathrm{up}}/k_B$	Frequency	Wavelength	Instrument	FWHM	Spectral resolution
		(K)	(GHz)	($\mu$m )		(′′)	(km s-1)
H2O	(221–110)	159.8	2773.977	108.07	PACS	9.4	318
$\left[O\text{i}\right]$	${}_{}{}^{3}P_0^{}$–${}_{}{}^{3}P_1^{}$	326.6	2060.069	145.53	PACS	10.2	260
H2O	(212–101)	80.1	1669.905	179.53	PACS	13.2	210
CO	(14–13)	580.5	1611.794	186.00	PACS	13.2	195
47′′$\times$47′′  at 9$\stackrel{\prime \prime }{.}$4 spatial sampling maps
CO	(22–21)	1397.3	2528.172	118.58	PACS	9.4	300
OH	${}_{}{}^2\Pi _{\frac{3}{2}}^{}{-}^2{\Pi }_{\frac{3}{2}}J={\frac{5}{2}}^{-}-{\frac{3}{2}}^{+}$	120.7	2514.298	119.23	PACS	9.4	300
OH	${}_{}{}^2\Pi _{\frac{3}{2}}^{}{-}^2{\Pi }_{\frac{3}{2}}J={\frac{5}{2}}^{+}-{\frac{3}{2}}^{-}$	120.7	2509.934	119.44	PACS	9.4	300
CO	(20–19)	1109.2	2299.570	130.37	PACS	9.4	275
CO	(18–17)	944.9	2070.616	144.78	PACS	10.2	260
$\left[O\text{i}\right]$	${}_{}{}^{3}P_0^{}$–${}_{}{}^{3}P_1^{}$	326.6	2060.069	145.53	PACS	10.2	260
CO	(14–13)	580.5	1611.794	186.00	PACS	13.2	195
Single spectra
CO	(14–13)	580.5	1611.794	186.00	HIFI	13.2	0.093
CO	(9–8)	248.9	1036.912	289.12	HIFI	21.7	0.14
CO	(5–4)	83.0	576.268	520.23	HIFI	39.0	0.26
13CO	(5–4)	79.3	550.926	540.73	HIFI	40.8	0.27

Table 2. Table 2: Knots coordinates.

Knots label	RA (J2000)	Dec (J2000)
A	12${}_{}{}^{h}01_{}^{m}36\stackrel{s}{.}0$	-65°07′30$\stackrel{\prime \prime }{.}$0
B	12${}_{}{}^{h}01_{}^{m}32\stackrel{s}{.}0$	-65°08′09$\stackrel{\prime \prime }{.}$0
C	12${}_{}{}^{h}01_{}^{m}36\stackrel{s}{.}5$	-65°08′09$\stackrel{\prime \prime }{.}$0
IRS1	12${}_{}{}^{h}01_{}^{m}36\stackrel{s}{.}6$	-65°08′48$\stackrel{\prime \prime }{.}$5
IRS2	12${}_{}{}^{h}01_{}^{m}34\stackrel{s}{.}0$	-65°08′44$\stackrel{\prime \prime }{.}$5
D	12${}_{}{}^{h}01_{}^{m}36\stackrel{s}{.}0$	-65°09′29$\stackrel{\prime \prime }{.}$0
E	12${}_{}{}^{h}01_{}^{m}37\stackrel{s}{.}2$	-65°10′27$\stackrel{\prime \prime }{.}$0

Table 3. Table 3: Parameters of the LVG fitting of the CO lines towards knot A at a resolution of 24 ′′ .

Component	$N$(CO)	size	$T$	$n$(H2)
	cm-2	′′	K	cm-3
cold	2$\times$1016 – 6$\times$1016	$>$24	$\sim$ 80	3$\times$105 – 4$\times$106
warm	3$\times$1015 – 9$\times$1015	18	1700 – 2200	2$\times$104 – 6$\times$104

Full text

*Aims. *We aim to define the global morphology of the warm gas component of the BHR 71 outflow and at modelling its shocked component.

Methods. We present the first far infrared Herschel images of the BHR71 outflows system in the CO (14–13), H2O (221–110), H2O (212–101)  and $[{O\textsc{i}}]$  145 $\mu$m transitions, revealing the presence of several knots of warm, shocked gas associated with the fast outflowing gas. In two of these knots we performed a detailed study of the physical conditions by comparing a large set of transitions from several molecules to a grid of shock models.

Results. The Herschel lines ratios in the outflow knots are quite similar, showing that the excitation conditions of the fast moving gas do not change significantly within the first $\sim$ 0.068 pc of the outflow, apart at the extremity of the southern blue-shifted lobe that is expanding outside the parental molecular cloud. Rotational diagram, spectral line profile and LVG analysis of the CO lines in knot A show the presence of two gas components: one extended, cold ($T\sim$80 K) and dense ($n$(H2) = 3$\times$105 – 4$\times$106 cm-3) and another compact (18*′′), warm ($T$ = 1700 – 2200 K) with slightly lower density ($n$(H2) = 2$\times$104 – 6$\times$104 cm-3*). In the two brightest knots (where we performed shock modelling) we found that H2 and CO are well fitted with non-stationary (young) shocks. These models, however, significantly underestimate the observed fluxes of $[{O\textsc{i}}]$  and OH lines, but are not too far off those of H2O, calling for an additional, possibly dissociative, J-type shock component.

*Conclusions. *Our modelling indirectly suggests that an additional shock component exists, possibly a remnant of the primary jet. Direct, observational evidence for such a jet must be searched for.

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, CH3OH, CS, H2CO, and H2O (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.

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 $\sim$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_{\rm up}$  CO and H2 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\textsc{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_{\rm up}$  CO lines (Parise et al., 2006; G15, ) and the hot gas traced by H2 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] ${}^{3}P_{1}$–${}^{3}P_{2}$ 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), H2O (221–110), H2O (212–101) and [O i] ${}^{3}P_{0}$–${}^{3}P_{1}$ by using the line spectroscopy observing mode. The final extended map was built combining three 47*′′$\times$47′′*  partially overlapping, Nyquist-sampled sub-maps (see Fig. 1). Data were reduced with the HIPE111HIPE is a joint development by the Herschel Science Ground Segment Consortium, consisting of ESA, The NASA Herschel Science Center, and the HIFI, PACS and SPIRE consortia. v9.0 package, where they were flat-fielded and flux-calibrated by comparison with observations of Neptune. The flux calibration uncertainty amounts to $\sim$ 20$\%$, while the Herschel pointing accuracy is $\sim$2*′′*. Post-pipeline reduction steps were performed to locally fit and remove the continuum emission and to construct integrated line maps (see Fig 2).

2.2 Herschel single pointing

Two positions (shown in Fig. 1) were observed in additional spectral lines: CO (5–4), 13CO (5–4), CO (9–8), CO (14–13) with HIFI and CO (14–13), (18–17), (20–19), (22–21), OH at 119 $\mu$m  and [O i] at 145 $\mu$m  with PACS. The two positions correspond to local H2 peaks, namely the SiO knot (Gusdorf et al., 2011) in the northern lobe (R.A.(J2000) $12^{\rm h}01^{\rm m}36.00^{\rm s}$, Dec.(J2000) -65°07*′34′′) and knot 8 (Giannini et al., 2004) in the southern lobe (R.A.(J2000) $12^{\rm h}01^{\rm m}38.90^{\rm 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′′$\times$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 T${}_{\rm A}^{*}$ scale was performed by the in-orbit system, while the spectra were converted to a Tmb 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 $\sim$ 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_{\rm bol}$=13.5 L*⊙), the exciting source of the more extended outflow. The second protostellar source IRS2 ($L_{\rm bol}$=0.5 L⊙) 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 $\sim$ 40′′*  in the northern lobe and $\sim$ 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 H2 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 $\mu$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 $\sim$60 km s*-1*, a velocity higher than previously reported ($\sim$ 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 12CO (5–4)/13CO (5–4)  ratio we measured the optical depth of the 12CO (5–4) line that results to be optically thin only in the red wing at $\varv\geqslant$ 2 km s*-1*. The intensity peak of the CO lines moves towards higher velocity for lines with higher $J_{\rm up}$, indicating that the low and high $J_{\rm 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 ($\varv\sim-$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 12CO (5–4)/ 13CO (5–4) ratio shows that also in the southern blue lobe the 12CO (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 ($\varv<-$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_{\rm 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.

Towards the two outflow positions (indicated in Figs. 1 and 2) we also acquired PACS maps of high $J_{\rm up}$  CO transitions, namely CO (14–13), CO (18–17), CO (20–19), CO (22–21), and of the OH 119 $\mu$m  and $[{O\textsc{i}}]$  145 $\mu$m  lines. Towards knot A, the CO lines are bright, while the OH and $[{O\textsc{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 $\sim 3\sigma$ level CO (14–13)  emission is detected while the other lines are not detected, confirming once again the lower excitation condition of this position.

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_{\rm up}$  CO lines (Parise et al., 2006; G15, ) and the hot gas traced by H2 (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 ($T\geqslant$ 1000 K) already observed in H2 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\textsc{i}}]$  63 $\mu$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 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_{\rm 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_{\rm 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_{\rm 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_{\rm 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 H2 0-0 S(5) ($E_{\mathrm{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_{\rm up}$  CO (3–2) line having the classical morphology of the cavity walls and the higher $J_{\rm 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_{\rm 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_{\rm up}$  CO lines. Although no distinct velocity components are present in the CO spectra of the BHR 71 outflow, the observed displacement of the location of the emission maximum of CO with increasing $J_{\rm 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 $\mu$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_{\rm up}$  CO corresponds to the peak of atomic lines such as $[{O\textsc{i}}]$  63 $\mu$m  and $[{Fe\textsc{ii}}]$ 26 $\mu$m  and of the H2 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.

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$\aas@@fstack{\prime\prime}$2 and extracted the brightness in the position of the seven Herschel knots. In Fig. 8 we show the H2O 179$\mu$m/108$\mu$m, the CO (6–5)/(14–13) and the $[{O\textsc{i}}]$  63$\mu$m/145$\mu$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 H2O 179$\mu$m/108$\mu$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\textsc{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 104 – 105 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\textsc{i}}]$  63$\mu$m/145$\mu$m  ratios have been already observed in several sources (Liseau et al., 2006) and interpreted as due to self-absorption of the 63 $\mu$m  line. This interpretation was corroborated by the recent observation of a spectroscopically resolved $[{O\textsc{i}}]$  63 $\mu$m  line towards a high-mass source which shows a very prominent absorption (Leurini et al., 2015). The self-absorption of the $[{O\textsc{i}}]$  63 $\mu$m  line could be the reason of the extremely low 63$\mu$m/145$\mu$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 $\mu$m  but it does not affect lines at higher wavelengths. We cannot exclude the possibility that the $[{O\textsc{i}}]$  63 $\mu$m  line towards IRS1 is also partially absorbed therefore the density inferred form the $[{O\textsc{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_{\rm up}$= 3 to $J_{\rm 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_{\rm 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_{\rm 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_{\rm up}$  from 3 to 11 are extended over the beam of 24′′while the lines with $J_{\rm 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_{\rm up}$  lines and a higher excitation component fitting the high $J_{\rm 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_{\rm up}$=3, 4, 6, 7) and the PACS lines ($J_{\rm 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_{\mathrm{rot}}$=68 K and $N$(CO)=6$\times$1016cm*-2*  for the low $J_{\rm up}$  APEX lines and $T_{\mathrm{rot}}$=267 K and $N$(CO)=6$\times$1015cm*-2*  for the high $J_{\rm 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_{\rm up}$  lines, at least the ones with J_{\rm up}$$\leq5, 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).

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}(\varv)=I_{CO}$(0) exp(-$\varv/\varv_{0}$) with the same slope $\varv_{0}$= 10 km s*-1*  for the two lines (Fig. 10), while the lower $J_{\rm up}$  lines have a more complex profile indicating that different gas components contribute to the total emission. A similar slope ($\varv_{0}$= 12 km s*-1*) has been found for the profile of the high $J_{\rm 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_{\rm up}$$\geqslant13 to shock excited gas. Another indication that the high $J_{\rm 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

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 BASECOL222http://basecol.obspm.fr 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$\times$1015 to 1$\times$1018 cm*-2*, the range in particle density is from 5$\times$103 to 3$\times$108 cm*-3*, the range in temperature is from 25 to 2000 K.

It is impossible to fit all the observed CO lines from $J_{\rm up}$= 3 to $J_{\rm 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_{\rm up}$  component is extended, that is, we assume a filling factor of 1, and the size of the high $J_{\rm up}$  transitions is measured from the PACS maps to be 18*′′. The CO column density of the low $J_{\rm up}$  lines is well constrained from the fitting of the APEX lines at $N$(CO) $\sim$ 6$\times$1016* 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_{\rm up}$  lines is extended (size$>$ 24*′′), cold ($T\sim$80 K) and dense ($n$(H2) = 3$\times$105 – 4$\times$ 106 cm-3*); the second component fitting the high $J_{\rm up}$  lines is compact (size=18*′′), warm ($T$ = 1700 – 2200 K) and with a slightly lower density ($n$(H2) = 2$\times$104 – 6$\times$104 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 H2 lines (Giannini et al., 2011) confirming that high $J_{\rm up}$  CO lines and rotational H2 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_{\rm 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 H2 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\textsc{i}}]$, OH, and H2O. 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.

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 H2 0–0 S($i$) (with $i$ = 0 to 7), CO ($J_{\rm up}$= 3, 4, 6, 7, 11, 16) and SiO ($J_{\rm up}$= 5, 6, 8). More specifically, the observable associated with each molecule was an excitation diagram for H2, 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 H2, CO and SiO, with the following input parameters: pre-shock density of $n$(H) = 104cm*-3*, shock velocity of $\varv_{\rm s}$ = 22 km s*-1*, magnetic field parameter $b$ = 1.5 (defined by $B$($\mu$G) = $b\times$ [$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 $\sim$ 3800 years.

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 H2 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 H2 are shown in Fig. 12.

As we also dispose of $[{O\textsc{i}}]$  observations by Herschel-PACS (${}^{3}P_{\rm 1}$–${}^{3}P_{\rm 2}$ at $\sim$63 $\mu$m and ${}^{3}P_{\rm 0}$–${}^{3}P_{\rm 1}$ at $\sim$145 $\mu$m) and we have compared their emissivities to the predictions of the best-fit model found fitting the CO and H2 emission. We did not include the $[{O\textsc{i}}]$  lines in the best fitting process because contrary to CO lines they are not velocity-resolved and contrary to H2 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\textsc{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\textsc{i}}]$  63 $\mu$m, of the order of 106 cm*-3*, likely difficult to fulfil in the analysed shock region. The result of the comparison of observed $[{O\textsc{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 $\mu$m  and 145 $\mu$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\textsc{i}}]$  and possibly negligibly in CO or H2, such as an atomic jet component. In the first case (i), for example, an enhancement of $[{O\textsc{i}}]$  and OH can be the consequence of the photo-dissociation of H2O. 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 H2 line emission (namely significantly reducing this emission, as these molecules would also undergo dissociation to some extent). Moreover knot A is about 1.3$\times$104 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\textsc{i}}]$  emission and would not modify the H2 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\textsc{i}}]$  emission, bring a negligible addition in terms of CO line emission, and a weak one in terms of H2 emission. This would be the case, for example, of the stationary model with $n$(H) = 103 cm*-3*, $b=1.5$, and $\varv_{\rm 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\textsc{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\textsc{i}}]$, are necessary to unveil the primary jet in outflows and to confirm our current interpretation of $[{O\textsc{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 $\mu$m  and 119.4 $\mu$m) and H2O (lines at 108 $\mu$m  and 179 $\mu$m). We note that the OH lines and the o-H2O 179 $\mu$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 $\mu$m  and 119.4 $\mu$m  lines (models: 1.0$\times 10^{-3}$ and 1.8$\times 10^{-3}$ K km s*-1*, observations: 0.25 and 0.36 K km s*-1*, respectively), similarly as what we obtained for $[{O\textsc{i}}]$. Conversely, the 108 $\mu$m  and 179 $\mu$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 H2O lines. This will be the subject of a dedicated work in a forthcoming publication.

6.2 Knot D

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 13CO (3–2) line G15 . We used this line in combination with 12CO (3–2) to estimate the optical thickness of the latter line. Fig. 14 shows spectra of the 13CO and 12CO (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 13CO is detected at more than 3$\sigma$. Assuming an identical excitation temperature for the 12CO and 13CO lines, and a typical interstellar abundance ratio of 50–60 (e.g. Langer & Penzias, 1993), we conclude that the 12CO (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.

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_{\rm 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: 103cm*-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) = 104 cm*-3*, $b$ = 1.0, and $\varv_{\rm s}$ = 20 km s*-1*, and an age of about 3900 years for the blue-shifted CO emission, and $n$(H) = 104 cm*-3*, $b$ = 1.5, and $\varv_{\rm s}$ = 20 km s*-1*, and an age of about 1000 years 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 years. 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 years 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 $\sim$56*′′*  from the driving source in the northern lobe and $\sim$113*′′*  in the southern lobe, translating in a dynamical age comprised between 5000 and 11000 years (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 H2 observations by Spitzer-IRS to further check out 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\textsc{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\textsc{i}}]$  transitions. We refer to the previous subsection for an explanation of this result.

Finally, in this knot, we have observed H2O emission from the lines at 108 $\mu$m  and 179 $\mu$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\textsc{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), H2O (221–110), H2O (212–101)  and $[{O\textsc{i}}]$  145 $\mu$m  of the BHR 71 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 BHR 71 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 H2 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\textsc{i}}]$  63$\mu$m  emission.

In the brightest position of the red lobe of the IRS1 outflow (knot A) we observed additional CO lines with $J_{\rm up}$  up to 22. Combining our Herschel observations with measurements of the low $J_{\rm up}$  CO, H2 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\sim$80 K) and dense ($n$(H2) = 3$\times$105 – 4$\times$106 cm-3*) and another compact (18*′′), warm ($T$ = 1700 – 2200 K) and slightly less dense ($n$(H2) = 2$\times$104 – 6$\times$104 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 H2 in both knots. In knot A, the best-fit model is non stationary with $n$(H) = 104 cm*-3*, $b$ = 1.5, and $\varv_{\rm s}$ = 22 km s*-1*, and an age of about 3800 years. 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) = 104 cm*-3*, $b$ = 1.0, and $\varv_{\rm s}$ = 20 km s*-1*, and an age of about 3900 years, 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) = 104 cm*-3*, $b$ = 1.5, and $\varv_{\rm s}$ = 20 km s*-1*, but a younger age of about 1000 years. 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) = 104 cm*-3*, $b$ = 1.0, and $\varv_{\rm s}$ = 15 km s*-1*, and an age of about 1450 years. 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\textsc{i}}]$  lines ratio.

The non-stationary shock models that well reproduce the CO and H2 observations underestimate by about two orders of magnitude the observed $[{O\textsc{i}}]$  lines at 63 $\mu$m  and 145 $\mu$m, and the OH lines at 119.2 $\mu$m  and 119.4 $\mu$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) = 103 cm*-3*, $b$ = 1.5, and $\varv_{\rm s}$ = 30–40 km s*-1*  as presented in Hollenbach & McKee (1989) can account for the observed $[{O\textsc{i}}]$  emission without significantly altering our conclusions with respect to CO and H2. Our Herschel-PACS observations of $[{O\textsc{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\textsc{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 H2 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_{\rm up}$  observable from the ground together with the intermediate and high $J_{\rm up}$  lines observed with Herschel (e.g. Lefloch et al. 2012; Benedettini et al. 2012) or large samples of H2O 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 H2 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 BHR 71 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 BHR 71 outflow one of the best candidates for future ALMA observations.

Acknowledgements.

We thank Prof. D. Neufeld for providing us with the reduced Spitzer data. GB acknowledges the support of the Spanish Ministerio de Economia y Competitividad (MINECO) under the grant FPDI-2012-18204. GB 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 U.S. 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 copyright ©1993-5 by the Anglo-Australian Observatory Board, and are distributed herein by agreement.