Free Access
Issue
A&A
Volume 575, March 2015
Article Number A98
Number of page(s) 14
Section Interstellar and circumstellar matter
DOI https://doi.org/10.1051/0004-6361/201425142
Published online 03 March 2015

© ESO, 2015

1. Introduction

Molecular shocks are ubiquitous in the interstellar medium (ISM) of our Galaxy. They can be associated with the formation of “ridges” at the convergence region of molecular clouds (e.g. W43, Nguyen-Luong et al. 2013), to the jets and outflow systems related to the formation of low- to high-mass stars (from low, e.g. Gomez-Ruiz et al. 2013; to high, Leurini et al. 2013; through intermediate, Gómez-Ruiz et al. 2012; for examples, also see Arce et al. 2007; Frank et al. 2014; or Tan et al. 2014 for reviews), or to supernova remnants (SNRs) interacting with molecular clouds (e.g. W44, Anderl et al. 2014). However, an analysis of these shocks is challenging because several physical processes are contributing to the observed emission. The W43 ridges are illuminated by an energetic UV radiation field coming from neighbouring H ii regions of star clusters (Bally et al. 2010). Studying young stellar objects (YSO) by pointing on the central protostar leads in practice to a study of the combination of ejection shocks (that can be multiple, e.g. Kristensen et al. 2013), infall processes, and UV illumination (the latter being all the more important when the mass of the forming object is large, e.g. Visser et al. 2012; San José-García et al. 2013). Outflows from massive star-forming regions are also illuminated by strong radiation in the X-ray regime (e.g. W28 A2, Rowell et al. 2010). Similarly, old SNRs are also often subject to X-rays, and also to γ-ray emission, which are the signature of the acceleration of particles that locally took place shortly after the supernova explosion (e.g. Hewitt et al. 2009).

“Pure” molecular shock regions can be defined as regions where the physics and chemistry are dominantly driven by shocks. Studying these regions, therefore, is of crucial importance to investigating the feedback they exert on their environment, whether this feedback is energetic or chemical. From the energetic point of view, these studies allow us to assess the contribution of shocks to the excitation of e.g. CO on galactic scales, as observed by Herschel (in NGC 1068: Hailey-Dunsheath et al. 2012; M 82: Kamenetzky et al. 2012; NGC 6240 and Mrk 231: Meijerink et al. 2013, or NGC 253: Rosenberg et al. 2014a and Arp299: Rosenberg et al. 2014b) from the inside of the Galaxy. From the chemical point of view, these studies allow us to investigate the formation paths of water (e.g. Leurini et al. 2014b) or more complex molecules (e.g. complex organic ones, Belloche et al. 2013), and to understand their presence in planetary systems, for instance. These “pure” shock regions are not numerous. The most remarkable and well-studied pure shock region is the B1 knot of the L1157 bipolar outflow. Sufficiently distant from the central protostar, this region remains uncontaminated by infall or irradiation processes, and has been the subject of a number of studies dedicated to studying its energetics or chemical composition (Gusdorf et al. 2008b; Codella et al. 2010; Flower & Pineau des Forêts 2012; Benedettini et al. 2013; Busquet et al. 2014; Podio et al. 2014). In particular, L1157 was mapped by the German Receiver for Astronomy at Terahertz Frequencies (GREAT) onboard the Stratospheric Observatory for Infrared Astronomy (SOFIA) in CO (1110), but the low signal-to-noise ratio prevented Eislöffel et al. (2012) from a thorough study of its energetics or chemistry. This article focus on the analysis of a similar shock position in the BHR71 outflow, the “southern twin” of L1157. Because of its southern location, BHR71 is an ideal target to be observed with the Atacama Large Millimeter/submillimeter Array (ALMA), which will never be done for L1157 because of its high northern location. This article is organised as follows: Sect. 2 presents our observations. Section 3 presents the new CO data we made use of in our analysis, while the existing H2 and SiO data is described in Sect. 4. The results of shock modelling and their further use is exposed in Sect. 5, and Sect. 6 summarises our findings.

thumbnail Fig. 1

The BHR71 bipolar outflow in its entirety, in left: CO (65) (colours, associated with the wedge, in main beam temperature units, K km s-1) and (3–2) (white contours, from 30% to 100% of the signal in steps of 10%), both observed by the APEX telescope and integrated between −50 and +50 km s-1, with the resolutions indicated by green circles in the lower right corner; right: the 8 μm emission detected by the Spitzer/IRAC receiver (colours, N09), with the H2 0–0 S(5) emission as observed by the Spitzer/IRS receiver in the inner parts of the outflow (white contours, from G11), and the SiO (5–4) emission (red and black contours, G11) in the upper lobe (the green circles in the lower right corner show the respective resolutions of CO (3–2), (11–10) and SiO (5–4)). On both maps, the grey inset is the field shown in Fig. 2, the SiO peak and knot positions are indicated in blue, the knot 5 and HH320 region are in pink, and the IRS 1&2, knots 6, 8, and HH321 region are indicated by grey dots.

Open with DEXTER

2. The BHR71 bipolar outflow

2.1. Previous work

The double bipolar outflow BHR71 (Bourke et al. 1997; Myers & Mardones 1998; Parise et al. 2006) lies close to the plane of the sky. It is powered by two sources IRS1 and IRS2, separated by ~3400 AU (Bourke et al. 2001), of luminosities 13.5 and 0.5 L (Chen et al. 2008), relatively nearby (~200 pc, Bourke et al. 1995a). The present work builds on previous studies of the BHR71 bipolar outflow in H2 (Neufeld et al. 2009; and Giannini et al. 2011, hereafter N09 and Gia11), and H2 plus SiO (Gusdorf et al. 2011, hereafter G11). In their work, N09 mostly described the Spitzer observations of the outflow. The InfraRed Spectrograph (IRS) was used to map the inner part of the outflow in the pure rotational transitions of H2 as well as in Fe ii and S i transitions. A region corresponding to approximately half the length of the outflow was covered by these observations around the driving protostars. The results were compared to previous observations of the entire region by the InfraRed Array Camera (IRAC) onboard the same telescope, showing that 30% and 100% of the luminosity of bands 3 and 4 could be accounted for by H2 lines emission, similar to L1157. The pure rotational H2 luminosity of the flow was estimated to be 4.4 × 10-2L, less than 1/3 of that of L1157, but measured only from the fraction of the BHR71 outflow that was mapped. The H2 mass above 100 K was constrained to be around 2.5 × 10-3 M, 20 times less than in L1157. However the H2 densities for both outflows were constrained to ~103.8 cm-3. In Gia11, these IRS observations of H2 were detailed line by line. An average column density of H2 of around 1020 cm-2 was extracted from the map. Two temperature components were apparent, at ~300 and ~1500 K. A non-local thermodynamical equilibrium analysis was performed and yielded a high H2 average density of a few 106 cm-3. The total H2 luminosity of the flow was estimated to be twice as high as the pure rotational lines, also based on previous observations of the rovibrational lines emission presented in Giannini et al. (2004). In a few positions where this was possible, the observations of pure rotational H2 lines were combined with those of rovibrational lines to generate a more complete excitation diagram that was also compared to shock models, confirming the high value of H2 density.

In G11, the authors presented a shock-model analysis of various positions in the northern lobe of the BHR71 outflow. They used Spitzer observations of H2, and APEX observations of SiO to constrain shock models. Their analysis was hampered by the rather high beam-size of their SiO observations. However, they identified two positions where a shock analysis was favourable: the “SiO peak” and “SiO knot”. These positions lie at the apex of the inner bipolar structure of the outflow, relatively un-contaminated by envelope infall. They are far away from the protostars, and separated from them by the inner outflow cavity, hence as shielded as possible from their potential UV radiation. For all these reasons they are reminiscent of the L1157-B1 position, and are considered in the present study as “pure shock” positions.

Figure 1 shows the entire BHR71 outflow as seen through these various datasets and highlights the particular positions as first introduced by Giannini et al. (2004): the so-called knots 5, 6, and 8 are thus indicated, as well as the two Herbig-Haro objects HH320A/B and HH321A/B. The driving protostars IRS 1 and IRS 2 are also marked. The two extra positions used in G11, “SiO peak” (of emission) and neighbouring “SiO knot” are also shown at the apex of the upper inner lobe of the outflow. In this figure, the previous and most recent observations can be compared (APEX, see Sect. 2.2): intensity in the CO (65) (colours) and (32) (white contours) transitions integrated between 50 and 50 km s-1 in the left panel, and Spitzer IRAC (8 μm, colours) and IRS (H2 00 S(5), white contours) maps overlaid with SiO (54) map (red contours in the upper lobe) in the right panel. The white contours then show how much of the outflow was observed by the IRS onboard Spitzer.

2.2. APEX observations of CO

APEX1 observations of the entire BHR71 outflow were conducted in 2013 and 2014. The analysis of the whole maps is out of the scope of the present work, and will be the subject of a forthcoming publication. In the present study, we used information inferred from the full maps in the 13CO (32), 12CO (32), (4–3), (6–5), and (7–6) transitions, which we briefly describe. The APEX observations towards BHR71 were conducted on several days: June 3 and 4, 2013, and June 28, 2014. We used the heterodyne receivers FLASH345, FLASH460 (First Light APEX Submillimeter Heterodyne receiver, Heyminck et al. 2006), and CHAMP+ (Carbon Heterodyne Array of the MPIfR, Kasemann et al. 2006; Güsten et al. 2008), in combination with the MPIfR fast Fourier transform spectrometer backend (XFFTS, Klein et al. 2012). The central position of all observations was , δ[J2000] = −65°08′530. We checked the focus at the beginning of each observing session, after sunrise and sunset on Saturn. We checked line and continuum pointing locally on IRC+10216, 07454-7112, or η Car. The pointing accuracy was better than ~5′′ rms, regardless of which receiver we used. Table 1 contains the main characteristics of the observed lines and corresponding observing set-ups. The observations were performed in position-switching/on-the-fly mode using the APECS software (Muders et al. 2006). We reduced the data with the CLASS software2. This reduction included baseline subtraction, spatial, and spectral regridding. For all observations, the maximum number of channels available in the backend was used (8192). We obtained maps for all considered transitions, covering the field of the whole outflow shown in Fig. 1.

thumbnail Fig. 2

Overlay of the velocity-integrated CO (65) (colour background) with the CO (32) (white contours) emission observed with the APEX telescope in the inner upper lobe of the BHR71 outflow. For both lines, the intensity was integrated between 50 and 50 km s-1. The wedge unit is K km s-1 in main beam temperature. The CO (32) contours are from 50 to 100% of the maximum, in steps of 10%. The half-maximum contours of the CO (32) and (65) maps are indicated in red and black, respectively. The dark blue circles indicate the positions and beam size of the SOFIA/GREAT observations for the CO (1110) transition. The APEX and SOFIA beam sizes of our CO (65), (1615) and (1110) observations are also provided (upper left corner light green circles, see also Table 1). The pink hexagons mark the position of the so-called knot 5 and HH320AB object.

Open with DEXTER

Figure 2 shows the velocity-integrated CO (65) broad-line emission (colours, resolution 9′′) in the upper inner part of BHR71, overlaid with the CO (32) (white contours, resolution 18′′). The higher angular resolution of the CO (65) line emission shows that it traces the walls of the cavity of the outflow associated with IRS1 (see Fig. 1), with a local maximum of emission also corresponding to the outflow driven by IRS2 (the HH320 AB position in the map). The outflow is seen close to the plane of the sky. The two positions of interest identified by G11 are marked by a blue hexagon and circle corresponding to the beam of our CO (1110) observations with GREAT (which will be the focus of our analysis, see Sects. 2.3 and 3). The position of the SiO peak, detected through 28′′-resolution observations of the SiO (5–4) is localised between two emission peaks in CO (65). The most prominent of these CO peaks is the southern one, which corresponds to the so-called SiO knot; it coincides with an H2 emission peak (also see Fig. 5, which overlays the same CO (65) field with the H2 0–0 S(5) data of N09 and Gia11). The half-maximum emission contours are also given in this map in thick black and red contours for the (6–5) and (3–2) lines, respectively, at their nominal resolutions.

Table 1

Observed lines and corresponding telescope parameters for the APEX and SOFIA observations of BHR71.

2.3. SOFIA-GREAT observations of CO

The observations towards BHR71 were conducted with the GREAT3 spectrometer (Heyminck et al. 2012) during SOFIA’s Cycle 1 “southern deployment” flight on July 23rd, 2013. We observed two positions, towards the northern apex of the inner outflow structure: the SiO knot and peak, as defined in G11 (Figs. 1 and 2). We tuned the receiver to the CO (1110) and (1615) lines frequency 1267.014 GHz LSB and 1841.346 GHz USB. We connected the receiver to a digital FFT spectrometer (Klein et al. 2012) providing a bandwidth of 1.5 GHz with respective spectral resolutions of 0.018 and 0.012 km s-1. The observations were performed in double beam-switching mode, with an amplitude of 80′′ (or a throw of 160′′) at the position angle of 135° (NESW) and a phase time of 0.5 s. The nominal focus position was updated regularly against temperature drifts of the telescope structure. The pointing was established with the optical guide cameras to an accuracy of ~5′′. The line and observation parameters are listed in Table 1. The integration time was 13 min ON source for each line, for respective final rms of ~0.70 and 0.75 K. The data were calibrated with the KOSMA/GREAT calibrator (Guan et al. 2012), removing residual telluric lines, and subsequently processed with the CLASS software.

3. Results: the CO data

3.1. Spectra

thumbnail Fig. 3

CO transitions in the SiO peak (left panels) and knot (right panels) positions indicated in Fig. 2: APEX (3–2), black line; (4–3), red line; (6–5), green line and histograms (upper panels); (7–6), dark blue lines; SOFIA (11–10), pink lines and (16–15), grey lines (lower panels, overlaid on the green histograms of the 6–5 transition). The last two were multiplied by five for comparison purposes. Respective spectral resolutions are 0.33, 0.50, 0.64, 0.54, 0.90 and 0.62 km s-1. The vertical dotted line marks the cloud velocity, −4.5 km s-1 (Bourke et al. 1995b).

Open with DEXTER

Figure 3 presents all 12CO spectra obtained in the two targeted positions in the course of our APEX (Jup = 3, 4, 6, 7) and SOFIA (Jup = 11, 16) observations, after convolution to the same angular resolution, that of the CO (1110) transition (see Table 1 for observing parameters related to each of these lines). The only exception to this convolution is the CO (1615) line, for which only a single-point observations is available. As the line was observed at an angular resolution of roughly 16.̋3, we could in principle use the integrated intensity extracted from this line as an upper limit to that convolved to the 24′′ resolution. In fact, we corrected this value for the beam dilution effect in our analysis, see discussion in Sect. 5.1 on how we made use of this line. On the figure, we split the spectra in two groups for visibility purposes (Jup = 3, 4, 6 in the upper panels, and Jup = 6, 7, 11, 16 in the lower panels).

In both positions, all lines exhibit a profile typical of the presence of a pure shock, with wings extending towards red-shifted velocities, up to 30–40 km s-1; up to Jup = 6, self-absorption is also detected at the velocity of the cloud (around 4.5 km s-1). The CO (3–2) and (4–3) profiles coincide very well, suggesting that these lines are optically thick. This is confirmed for the CO (32) line, for which we observed the 13CO isotopologue on both positions (see Sect. 3.2, Table 1 and Fig. 4 for the spectra obtained on both positions). In both positions, the CO (6–5) profiles start to differentiate from the lower-lying profiles. This departure is confirmed in the higher-lying transitions, and reveals that the lines from Jup = 6 are probably excited in different layers of the shocked region than their lower-lying counterparts. Globally, it can be seen that the excitation conditions vary with the position: the Jup = 11 and 16 profiles differ, and the CO (16–15) is not even detected in the SiO peak position. For both positions, the excitation conditions also vary with the velocity: close to the systemic velocity, the lower-J line emission greatly dominates, an effect previously reported in L1157 B1 by Lefloch et al. (2012). Because of the associated weak or absent detection in the SOFIA lines, and because the two selected positions are not independent, we decided to exclude the SiO peak from the shock analysis presented in Sect. 5 (also see an additional argument in Sect. 4.1).

3.2. CO (32) line opacities

Figure 4 shows our spectra of the 13CO and 12CO lines in the SiO peak and knot positions at the same spectral and spatial resolutions. These spectra allow us to plot the ratio of the line temperature. For each velocity channel, the line temperature ratio 12CO/13CO is hence also shown (right ordinates) with the following colour-code: red dots show the ratios for the velocity channels where the 13CO is detected at more than 3σ, and orange dots where the 13CO detection was obtained with a confidence level between 2 and 3σ. The grey dots are for all the other channels. Horizontal dashed lines show the reference values of 20 and 40 for these ratios. The orange and red dots all lie below 30. 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), these line ratios yield minimum opacity values of 3 for the 12CO lines. We therefore conclude that the 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(s). As our observations yield a minimum optical thickness of 3 up in the wings, we decided to exclude both the CO (3–2) and (4–3) lines from our analysis.

thumbnail Fig. 4

Comparison of the 12CO (black line) and 13CO (pink line) (32) emission profiles obtained in the SiO peak (upper panel) and knot (lower panel) positions. The ratio of main-beam temperatures is also shown in the form of red dots (for the channels where the 13CO detection is over 3σ), orange dots (for the channels where the 13CO detection is between 2 and 3σ), and grey dots (for the remaining channels). Associated with this distribution are the grey horizontal dotted lines, which show values of 20 and 40 for this ratio.

Open with DEXTER

3.3. Dynamical age of the outflow

An important parameter for the modelling of young outflows is their age. As we have mapped the entire outflow, we are able to give an upper limit of its age, simply obtained by dividing the distance between the furthest point from the driving protostar and the protostar by the associated linewidth. The positions we consider belong to the northern lobe of the outflow powered by the IRS 1 protostar. As the furthest point with significant CO emission belonging to this outflow, we adopted the furthest point on the 10% CO (3–2) emission contour or the furthest point with 3σ emission away from the driving protostar, and found that these two approaches were equivalent. Indeed, the offset from both points relative to the IRS1 position is about (72′′, 280′′), corresponding to a distance of ~289′′. At a distance of 200 pc, this translates in ~5.8 × 104 AU. In this position, the full-width at zero intensity of both the CO (3–2) and (6–5) (unshown) lines is about 1520 km s-1, which converts to a dynamical age of (1.4 − 1.8) × 104 years. Of course, this is an upper limit of the real age of the outflow, as it is likely that the apex of the outflow was associated with greater velocities in the past. A similar measurement for the SiO knot position (i.e. dividing its distance to the protostar by the full-width at zero intensity) yields a dynamical age of ~1800 years.

4. Results: existing data

4.1. H2 observations

thumbnail Fig. 5

Overlay of the map of CO (6–5) emission observed by the APEX telescope (colour background) with the H2 0–0 S(5) emission (white contours), observed with the Spitzer telescope. The wedge unit is K km s-1 (main beam temperature) and refers to the CO observations. The H2 0–0 S(5) contours are from 10% to 100%, in steps of 10%. The light blue contour defines the half-maximum contour of this transition. Like in Fig. 2, the blue circles and dots indicate the positions and beam sizes of the SOFIA/GREAT observations. The beam and pixel sizes of the CO (6–5) (11–10), (16–15), and H2 0–0 S(5) observations are the green circles in the upper left corner. The black contour delineates the half-maximum contour of the H2 0–0 S(2) transition. The field is the same as in Fig. 2, and the knot 5 and HH320 region are in pink. The field covered by Spitzer/IRS to observe the H2 emission is indicated in dashed grey line. It excludes the SiO peak position.

Open with DEXTER

thumbnail Fig. 6

Our model-observations comparisons. Upper panel: CO flux diagram over a beam of 24′′; the observations in the SiO knot position are the black squares, and the model results are the coloured circles (see text). The CO (1615) observational point is corrected for beam-size effect. The uncorrected point is the upper limit (see text), indicated by the grey arrow in this panel. Lower panel: H2 excitation diagram for the SiO knot position, extracted following the different procedures described in the text (empty symbols), global average in black squares, and model results in coloured circles (see text, with the same code as in the upper panel).

Open with DEXTER

We used the Spitzer/IRS observations of the H2 pure rotational transitions, 0–0 S(0) up to S(7), reported and analysed in N09, Gia11, and G11. The reduced data kindly communicated to us by David Neufeld contain rotational transition maps, with 3.̋6 angular resolution, centred on , δ[J2000] = −65°08′5302. Figure 5 shows an overlay of our APEX CO (65) map with the H2 00 S(5) region observed by Spitzer, both at their nominal resolution. The figure shows coinciding maxima between the two datasets in the selected position, and a slightly different emission distribution. This might be the effect of the better spatial resolution of the H2 data, which reveal more peaks than in CO. This overlay also shows the similar morphology of the emission of the S(2) (half-maximum contour in black) and S(5) (half-maximum contour in light blue) transitions on the region we chose to analyse. Unfortunately, the figure also shows the coverage of the H2 emission observations, and that the SiO peak position was not covered by the H2 observations, which is another reason to exclude it from our shock analysis presented in Sect. 5.

As part of our modelling (see Sect. 5), we used the excitation diagram derived for the selected emission region around the SiO knot position. The H2 excitation diagram displays ln(Nνj/gj) as a function of Eνj/kB, where Nνj (cm-2) is the column density of the rovibrational level (v,J), Eνj/kB is its excitation energy (in K), and gj = (2j + 1)(2I + 1) its statistical weight (with I = 1 and I = 0 in the respective cases of ortho- and para-H2). If the gas is thermalised at a single temperature, all points in the diagram thus fall on a straight line.

The initial resolution of the maps is 3.̋6, but we wanted to operate at the same resolution as for CO. Unfortunately, the SiO knot position lies at the edge of the field covered in the H2 map, as can be seen in Fig. 5. We consequently used three different methods for extracting the H2 line fluxes. First, we extracted the flux from the initial map at its nominal 3.̋6 resolution; second, we used the first method, but applied to a map that was convolved to the 24′′ resolution; and third, we associated a filling factor of 0.2 to the fluxes obtained through the second method.

Because of the location of the SiO knot and the rather large beam of our analysis, the first two methods provide lower limits to the H2 line fluxes. On the contrary, with method (3) we voluntarily extracted upper limits to these fluxes. We then used an average value between the values inferred from the three methods, and computed rather large errorbars based on the combination of all three methods. The values were corrected from interstellar extinction following the treatment already applied in G11. The result can be seen in the excitation diagram shown in the lower panel of Fig. 6, where points corresponding to methods (1), (2), (3) are shown in black empty diamonds, upward and downward triangle, and the resulting average points are the black squares, respectively. The overall average H2 values we used to build our figure are given in Table B.1.

thumbnail Fig. 7

The integrated intensity diagram generated from the G11 observations of three SiO lines (Jup = 5, 6, 8), corrected for beam size and for filling factor effect: filling factor of 1, 0.5, and 0.2 in black, dark grey, and light grey squares. The result of our best-fit model for the H2 and CO lines (nH = 104 cm-3, b = 1.5, νs = 22 km s-1, and an age of 3800 years) is shown with 1, 2, or 4% of the pre-shock Si placed in the grain mantles in red empty, dotted, or filled circles.

Open with DEXTER

Table 2

Shock model parameters.

4.2. SiO

As part of our study, we also decided to re-analyse the SiO observations already presented in G11. The emission from three lines was mapped in the northern lobe of the outflow: SiO (5–4), (6–5), and (8–7), and the corresponding spectra were extracted in the SiO peak and knot positions. As maps were obtained, all spectra were convolved to the resolution of the SiO (5–4) line, about 28′′. In the present study, we hence corrected the SiO dataset for the slight difference in resolution between this value and that of the CO (1110) observations, 24′′, by simply multiplying the G11 integrated intensities by (28/24)2. Also, given the rather large beam size, we generated an integrated intensity diagram (displaying against Jup for each transition) for three different filling factor values, 1, 0.5, and 0.2. The result can be seen in the form of respective black, dark grey and light grey squares in Fig. 7. The SiO values we used to build our figure are listed in Table B.1.

5. Discussion

In the following we seek to constrain the physical conditions in the shocked regions. Given the large number of observational constraints, covering several species and several associated transitions, the method-of-choice for obtaining these physical conditions is a comparison to sophisticated shock models. The results of the shock modelling forms the foundation of the discussion: what characterises the chemistry of the outflow, in particular with respect to SiO, and what characterises the energetics.

5.1. Constraining shock models

Shock models are used to constrain the physical conditions in the outflow shocks through comparison with H2 and CO, following the methods already used in Gusdorf et al. (2008a, 2012) and Anderl et al. (2014). We generated CO flux diagrams and H2 excitation diagrams from the observations of the SiO knot position, and compared them with results from a grid of models computed with the “Paris-Durham” (e.g. Flower & Pineau des Forêts 2003) 1D shock model. The H2 diagram has already been introduced in Sect. 4.1 and is shown in Fig. 6. Concerning the observable associated with CO, we chose to display the data in the form of a flux diagram (line flux vs. Jup), because of the increasing number of extragalactic studies that use this form. In the present case, we benefit from the fact that our observations were pointed towards a pure shock position. To generate a CO flux diagram, we consequently integrated the signal obtained over the whole line profile at the angular resolution of 24′′, associated with a filling factor of 1 (see Appendix A for an explanation on filling factors for various CO transitions). Because of their opacity, and their likeliness to be contaminated by ambient gas, we excluded the (3–2) lines from our analysis (see Sect. 3.2). Given the opacity values inferred from this line, we also excluded the (43) line from our fitting procedure. Overall, the two observational tools we used can be seen in Fig. 6: a flux diagram for CO (upper panel) and an excitation diagram for H2 (lower panel) on which the observational points for the SiO knot are always shown in black squares. The CO values we used to build our figures are listed in Table B.1. The observational flux diagram for the SiO peak position can be found in Appendix C (Fig. C.1).

The grid of shock models is that of G11, also used in other publications already mentioned (Gusdorf et al. 2012; Anderl et al. 2014). To summarise, we used a grid covering the following input parameters: pre-shock density nH from 103 to 106 cm-3, magnetic field parameter from b = 0.3 to 2 (defined by B(μ G) = b × [nH (cm-3)]1/2, where B is the intensity of the magnetic field transverse to the shock direction of propagation), and shock velocity νs from 15 to 35 km s-1. Our grid contains both stationary C- and J-type models, and also non-stationary, so-called CJ-type models (Lesaffre et al. 2004a,b). A more complete description of the parameters space coverage (nH, b, and νs) can be found in Table 2. The parameter coverage is not complete: not all velocities are present in our grid for all values of the magnetic-field parameter, b. In fact, the velocity of C-type shocks must remain below a critical value that depends mainly on the pre-shock density and magnetic-field parameter (Le Bourlot et al. 2002; Flower & Pineau des Forêts 2003), which explains the decrease of the maximum shock velocity with the pre-shock density in our grid. Additionally, the time necessary to reach a steadystate depends on the pre-shock density, shock velocity, and magnetic field values. Following the method presented in G08b (Sect. 4.1), the set of C-type shock models enabled us to restrict the range of the search in the parameter space for the CJ-type shock models. We then computed a grid of non-stationary shock models around a first estimate of the shock age, making sure that the range of ages was sufficient to include any model likely to fit the H2 observational data. The final considered ages range from a few tens to around fifteen thousand years. Our grid also includes a few stationary models including the effects of grain-grain interactions, as presented in Guillet et al. (2011); Anderl et al. (2013), and already used in Anderl et al. (2014); Leurini et al. (2014a).

We compared models and observations based on the three higher-J CO lines, from CO (65) to (11–10), and on all the H2 pure rotational lines. The CO (32) and (4–3) lines were excluded from our procedure, as stated above. Our initial CO (1615) observation is in principle an upper limit, since it was associated with a slightly smaller beam than the other lines. Nevertheless we assumed that the CO (1615) emission is as extended as the CO (65) emission, and we estimated the resulting flux over a 24′′ beam by multiplying the value observed over a 16.̋3 beam by a factor (16.3/24)2. We then treated this corrected value as any other CO observation. Similar to what we did in previous studies, we used a χ2 routine to compare these observations to the models. The best results can be seen in Fig. 6. On each of these panels, we show the points of our best fit in red circles, plus three other satisfying models in smaller, purple, green, and blue circles. Our best-fit model (in red points) is non-stationary, with the following characteristics: nH = 104 cm-3, b = 1.5, νs = 22 km s-1, and an age of 3800 years. The χ2 value of our 10 and 20 best-fit models are within a factor 1.4 and 2.1 to the best (lowest) value. Broadly speaking, we found the H2 lines in particular very difficult to fit: the exact curvature of the excitation diagram is only approached by our best model. We tried in vain to improve the quality of the fit by changing the ortho-to-para ratio of H2 in our calculations (switching its value from 3 to 1 and 2). Eventually, we can not exclude that a better fit could be found via a better gridding of our parameters space. However, it turned out that this best-fit model also fitted the CO (32) and (43) lines quite in a very satisfying way. These values can be compared to what was constrained by G11: nH = 104 cm-3, b = 1−2, νs = 25−30 km s-1, and an age of 300800 years. A few comments can be made on these shock parameters.

Shock type. The shock type we constrained is the same as found by G11. Indeed, it is a very general result that in low-mass star forming environments, H2 excitation diagrams can be fitted by these kinds of models only (e.g. Giannini et al. 2006; Gusdorf et al. 2008b; Flower & Pineau des Forêts 2013). All of our ten best fits are CJ-type models. Even more, the J- and C-type models are all at the end of our list of best-fit models: the C-type models generally do not fit the pure rotational H2 excitation diagram, while the J-type models do not generate enough CO emission to match the observations.

Pre-shock density. First, the pre-shock density value is the same as that found in previous analyses of BHR71 (Giannini et al. 2004, for the analysis of the HH320 region, G11 for the SiO knot position), and in the resembling pure shock position L1157-B1 (e.g. Gusdorf et al. 2008b; Flower & Pineau des Forêts 2012). It is also the value associated with our ten best fits to the dataset; it is relatively well constrained, as in particular the general shape of the CO diagram is very sensitive to the density, i.e. to the pre-shock density parameter. In our ranking of models by decreasing χ2 value, the first model with a different pre-shock density has a χ2 value that is more than five times the lowest one. Generally speaking, it is a very standard value when modelling shocks from low-mass YSOs (also see HH54, Giannini et al. 2006), which was also found in shocks associated with SNRs that are interacting with the interstellar medium (Cesarsky et al. 1999; Gusdorf et al. 2012; Anderl et al. 2014). This corresponds to a post-shock density of ~105 cm-3, a value that is one or two orders of magnitude smaller than those found by Gia11 in the HH320 and 321 regions of the same outflow. This could be due to the fact that their study only focussed on the brightest pixel associated with these regions, whereas we consider a rather large beam size. On the other hand, our value is just above the H2 density averaged over the whole inner outflow found by N09 (103.8 cm-3). This can be explained by the fact we are focussing on a bright H2 spot, and also that the post-shock value we are providing here is that of a stationary post-shock. In the (realistic) case where the shock has a finite spatial size, the post-shock gas will expand and come into pressure equilibrium with its surroundings, resulting in a globally lower value than ours.

Magnetic field strength. The strength of the transverse magnetic field is 150 μG in the pre-shock region, that turns into 1.62 mG in the post-shock, given the compression factor and the fact that the magnetic field is frozen in the neutral fluid. It is comparable to what was constrained in G11 (b = 1−2), and also to what was found in the studies of HH320 in BHR71, L1157-B1, or HH54, or in the aforementioned SNR studies. This value is also consistent with the analysis of Zeeman measurements in molecular clouds where magnetic and kinetic energy densities are in approximate equipartition (Crutcher 1999). Our ten best fits all predict a value of the b parameter between 1.0 and 2.0. Figure 6 shows the best model with a b value outside of this range, in purple points. This model is associated to a χ2 value that is 1.8 times the lowest. This model is also a bit older than our ten best-fit ones (see below). It does not fit the lower-J data as satisfyingly as our best-fit model (red points in Fig. 6).

Shock velocity. The shock velocity values of our ten best-fit models all are between 20 and 25 km s-1. In our ranking of models by decreasing χ2 value, the first model with a velocity out of this range has a χ2 value that is about three times the lowest one. More specifically, the shock velocity that we constrain is slower than, e.g. the full-width at zero intensity of the CO lines. The modelled shock velocity is the difference between the pre-shock velocity and the impact velocity. Consequently, if the pre-shock material is already accelerated, the actual shock velocity is expected to be correspondingly lower. Preliminary acceleration of the pre-shock is possible in the SiO knot position, as it has already been encompassed by the shock that is now propagating at the northern tip of the outflow (a bit more than twice as far from the driving protostars). Alternatively or simultaneously, this velocity discrepancy could be the sign of a limitation to our models. Indeed, we are modelling a multi-dimensional complex shock region, which is seen edge-on, by means of a one-dimension model, which is considered to be face-on. Additionally, it is likely that the rather large beam of our observations contains not one shock structure, but a collection of them. This can be seen in the spatial sub-structure of the H2 emission revealed in Fig. 5, where multiple peaks are detected. This can also be seen in the CO (Figs. 3 and 4) or SiO (G11) line profiles, where a “plateau” or a “bump” can be seen between 20 and 40 km s-1, very much resembling the bullets revealed in the CO (in e.g. Cep E, Gómez-Ruiz et al. 2012) or water line profiles (in e.g. L1448-mm, Kristensen et al. 2011). Another expression of the intrinsic multi-dimensional nature of the observed shocks is the impossibility of fitting the different molecules with the same filling factor values. For instance, H2 is only associated with the hottest regions in our beam, as a 500 K temperature is necessary to populate its levels, whereas CO is more easily excited (for Jup = 6, Eup = 116 K). In particular, unlike H2, CO is probably emitting in the post-shock expansion region mentioned in the “pre-shock density” paragraph above, where the primary shock structure progressively mixes with the ambient medium in all spatial dimensions. This effect cannot be accounted for by a 1D shock model.

Shock age. Finally the shock age is rather large: it is larger than what was previously constrained for this region of BHR71 by G11 (300–800 years), which is due to the high amounts of CO that we detected, that are not compatible with younger CJ-type shocks. However, again the age predicted by our ten best-fit models consistently lies between 3500 and 5500 yr, with the model in purple in the Fig. 6 being slightly older. In our ranking of models by decreasing χ2 value, the first model with an age younger than 3000 years has a χ2 value that is about twice the lowest. The age that we find is between the dynamical age of the SiO knot and that of the global northern lobe (see Sect. 3.3), which hints again at the fact that we are probably catching several shock episodes in our beam.

Limitations. The question of the largest source of uncertainty in these results is difficult to assess. The observational problems (the different beam sizes and the fact that the SiO knot lies close to the edge of the H2 observations) make for an important limitation, but we believe we have taken it into account by using conservative errorbars. More problematic is the complex nature of shocks structures. Lacking high angular resolution, we could only assume that we catch a single shock structure in our beam. Our results seem to indicate this is not the case, given the constrained shock velocity and ages, and the different filling factors we had to adopt for each molecule. A shock similar to those we have in our grid of models is probably propagating in one beam of 24′′, but it is very likely to be associated with less violent shocks corresponding to the mixing of its warm, dense, and accelerated post-shock with the surrounding ISM. To summarise, we could only be aware and accept the fact that we are averaging processes out by working over a beam size typical of single-dish observations. From the modelling point of view, solutions do exist to more realistically account for the observations of complex geometrical shock structures. The first consists of adopting a probability density function (PDF) of shocks, which is to consider that the observed shock characteristics follow a statistical distribution. This distribution can unfortunately only be guessed, and hence generates additional free parameters (see Lesaffre et al. 2013 for applications of this method). The second consists of stitching 1D shock layers (similar to the one we consider here) to either a curve (e.g. Kristensen et al. 2008) or a bow-shock structure (e.g. Gustafsson et al. 2010), hence simulating “pseudo-” 2 or 3D shock propagation. This approach also yields free parameters, such as the inclination of the magnetic field with respect to the shock structure. None of these methods would be ideal in the present case, given the lack of knowledge of the small-scale shock structure within the 24′′ beam (no indication on the collection of shocks is available at the moment), and the relatively limited dataset (if more free parameters are to be considered, more constraints are needed to lift the degeneracy in the results). Their application will become relevant when the SiO knot is observed at higher angular resolution (in CO and H2), and if possible with a maximum number of observations in key species (such as H2, CO, and SiO, but also H2O, OH and O i for instance). We note that at this point, a more tightly sampled shock model will prove necessary.

In the following we will present results based on our best fit only. All of our ten best-fit models were found to be non-stationary (CJ-type), with nH = 104 cm-3, b = 1−2, νs = 20−25 km s-1, and an age of 35004500 years.

5.2. Employing shock models: SiO chemistry

After constraining shock models based on CO and H2 observations, it is now possible to fine-tune the SiO chemistry. SiO has long been interpreted as a tracer of C-type shocks with a velocity greater than 2025 km s-1 (e.g. Caselli et al. 1997; Schilke et al. 1997; Gusdorf et al. 2008a). Indeed, in these types of shocks, the drift velocity between the charged grains and the most abundant, neutral molecular species is sufficient to generate the sputtering of the core of the grains, where all the silicon-bearing material was considered to be locked. However Gusdorf et al. (2008b) have demonstrated: 1) that CJ-type shock models were the only ones to fit the H2 emission in L1157-B1; and 2) that this process was not efficient enough to generate levels of SiO emission comparable to the observations in this kind of (CJ-type) shock models. In order to be able to self-consistently account for both SiO and H2 emission, one solution was considered, consisting of transferring a fraction of SiO from the core to the mantle of the grains in the pre-shock phase. The maximum fraction of SiO thus transferred to the mantles was set to 10%. With mantle sputtering being easier than core sputtering, these authors were then able to fit the H2 and SiO emission by means of a single, non-stationary shock model. The presence of silicon-bearing material in the grain mantles has also been assumed in Coutens et al. (2013) to explain the narrow emission component detected in the SiO (2–1) line profile at the systemic velocity around NGC 1333 IRAS 4A.

Interestingly, the transfer of Si from the core to the mantles in the form of SiO was also considered by Anderl et al. (2013) in the context of denser shock models. Indeed, for pre-shock densities above 105 cm-3, the effect of grain-grain interactions cannot be neglected in shock models, as shown by Guillet et al. (2011). Taking these interactions into account results in the significant production of small, charged grains that couple very well with the neutral fluid, in effect reducing the width of the shock layer, and increasing its maximum temperature. In the narrow shock layers thus produced, Anderl et al. (2013) have shown that the grain core sputtering is not efficient enough to produce levels of SiO emission comparable to the observations, hence the recourse to the transfer of a fraction of Si towards the mantles in the pre-shock phase. Leurini et al. (2014a) have validated this approach by successfully confronting these models with observations.

Consequently, we ran our best-fit model for the H2 and CO data (nH = 104 cm-3, b = 1.5, νs = 22 km s-1, and an age of 3800 years), including 0 to 10% of the pre-shock silicon-bearing material in the form of SiO in the grain mantles. We then generated the corresponding integrated intensity diagrams, which we compared with the observations. The result can be seen in Fig. 7. For a filling factor value conservatively varied between 0.2 and 1, we show that no more than 4% of SiO needs to be placed in the grain mantles to reproduce the observations. The presence of silicon-bearing material on the grain mantles could be explained by the fact that the SiO knot position has already been processed in the past by the shocks that are now located at the apex of the northern outflow lobe, further north, as can be seen in Fig. 1. Under this assumption, we have demonstrated that it is possible to self-consistently fit H2, and velocity-resolved CO and SiO observations in the “SiO knot”, pure shock position of BHR71.

5.3. Employing shock models: energetics

In this section, we discuss the energetics associated with the shocks along two axes: the CO excitation generated from our best-fit model, and the derivation of the outflow parameters.

CO excitation. Indeed, the excitation of CO has been intensively used in the literature to demonstrate the presence of various processes at work in the regions of star formation. For example, van Kempen et al. (2010b) and Visser et al. (2012) obtained CO flux diagrams (similar to our Fig. 6) with PACS, centred on low-mass protostars at the origin of various outflows (HH46, NGC 1333 IRAS 2A, and DK Cha). Based on these diagrams, they evidenced the existence of three physical components corresponding to three distinct processes contributing to the excitation of CO: a passively heated envelope, UV-heated outflow cavity walls, and small-scale shocks along the cavity walls. The SiO knot position is distant from the central protostars IRS 1 and IRS 2 (see e.g. Fig. 1), which most likely prevents any UV heating from the central star to be operating. Furthermore, it lies outside the envelope, as revealed by IRAC (Tobin et al. 2010), NH3 (1,1) (Bourke et al. 1995a, 1997) or N2H+ (Chen et al. 2008) emission. In other words, The SiO knot position offers the opportunity to study pure small-scale shocks along cavity walls. However, the comparison of our shock models with those presented in van Kempen et al. (2010b) or Visser et al. (2012) is not really meaningful. These authors indeed used a multiple-shock layers model, with a pre-shock density of 106.5 cm-3 close to the protostar, and 104.5 cm-3 further away along the envelope (see Fig. 4 of Visser et al. 2012). Furthermore the 1D shocked layers they used were outputs of the Paris-Durham model we are also using in the present study. However they did not include the tip of the envelope that would correspond to the SiO knot position, as they solely focus on the outflow cavity in the vicinity of the protostar. The conclusion from our analysis is that our constrained pre-shock density of 104 cm-3 is consistent with the range of values they used (as they considered a decreasing density further away from the protostar).

Another, more classical approach to discuss energetics and physical conditions based on CO excitation consists of using excitation diagrams, as reviewed in Visser (2014). Figure 8 presents the excitation diagram built from our best-fit model for CO. As extensively described in Goldsmith & Langer (1999), in local thermodynamical equilibrium conditions and under optically thin regime, the points in this type of diagram are expected to fall on a straight line, whose inverse of the slope is the excitation temperature of the transitions (also see the description of H2 excitation diagrams in Sect. 4.1). As pointed out in Visser (2014), numerous studies can be found in the literature describing the building of this kind of excitation diagrams based on PACS observations centred on outflow-driving protostars of low- to high-masses, and their subsequent decomposition in at least two gas components, one warm (Tex = 320 ± 50 K), and one hot (Tex = 820 ± 150 K). Usually the breakpoint between these two components is around 1800 K, between the (2524) and (2625) transitions. We have produced a similar diagram from our best-fit model, as can be seen in Fig. 8. For comparison purposes, we have also applied a two-component linear fit of the lines accessible to PACS, with the warm component (Tex = 216 K) fitting the level populations for Jup= 12 to 25, and a hot component (Tex = 422 K) fitting the level populations for Jup = 26 to 40. The values we constrained for the excitation temperatures within these two components are systematically below those inferred from PACS observations centred on outflow-driving protostars of all possible mass (van Kempen et al. 2010a; Herczeg et al. 2012; Goicoechea et al. 2012; Manoj et al. 2013; Karska et al. 2013, 2014; Green et al. 2013a,b; Dionatos et al. 2013; Lee et al. 2013; Lindberg et al. 2014). This shows that pure shocks contribute to both the so-called “warm” and “hot” components. Although the excitation temperature associated with the warm component is close to the observed values, the figure also shows that pure shocks fail to entirely account for any of these warm or hot components. Indeed the excitation temperatures that we find here are less than those derived from the observations in those papers. This means that close to the protostar, additional mechanisms, or different kinds of shocks should account for the CO observations. It also shows to which extent the collection of continuous temperature components that is generated by our shock models can be interpreted as a two-excitation component, through the analysis of CO excitation diagrams, for transitions with Jup = 12 to 40.

thumbnail Fig. 8

The CO excitation diagram produced from our best-fit model (red squares). Two temperature components can be fitted to our shock model over the PACS range of observations: a warm component (Tex = 216 K, dark blue line) fitting the level populations for Jup = 12 to 25 (light blue squares), and a hot component (Tex = 422 K, light green line) fitting the level populations for Jup = 26 to 40 (dark green squares).

Open with DEXTER

Outflow parameters. Finally, we studied the energetics associated with this pure shock position. Our first method is to follow the procedure presented in Anderl et al. (2014) for the W44 SNR shock study. Indeed, we can infer the mass, momentum, and energy associated with our best-fit model, under the assumption of a filling factor equal to one. The results are presented in the second column of Table 3. First, the total mass contained in the beam is 1.8 × 10-2 M. This value is far greater than that determined for the mapped area (half of the inner part of the outflow) by N09 (2.5 × 10-3 M), or Gia11 (0.6 × 10-2 M). This might be partly explained by the fact that our beam size is rather large and that we decided to focus on the brightest H2 spot in the map, contrary to those studies that operated on values averaged over larger areas. More convincingly, this is probably due to the fact that our observations–models approach gives us access to the population of the (v = 0, J = 0, 1) H2 levels, contrary to the observations. The contribution of these two levels to the total column density is significant: a linear fit to the Spitzer observations presented in Fig. 6 yields a column density N(H2) ~ 6.9 × 1019 cm-2, whereas linearly fitting the (v = 0, J = 0, 1) part of the modelled rotational diagram yields N(H2) ~ 1.8 × 1021 cm-2, which is 26 times larger than that measured based only on the observations. Moreover, our value is consistent with the total mass of 1 M for the northern lobe determined by Bourke et al. (1997) based on CO (10) and (21) line observations (combined with 13CO and C18O (1–0)). From a different perspective, the value that we found is simultaneously ~50 times smaller than that found in the bright CO positions of the W44 SNR studied in the same way by Anderl et al. (2014), where molecular shocks also propagate. The corresponding momentum is 0.4 M km s-1, also a good factor ~100 below the values inferred using similar methods in W44. This is also a factor 10 below the value found based on more observational methods (described in the next paragraph) over a 12.̋5 beam in the massive star-forming region W28 A2 that encompasses three different outflows at least (Gusdorf et al. 2015). Finally, the dissipated energy is 4.2 × 1043 erg, typically two orders of magnitudes below the equivalent quantity in W44. This is also two orders of magnitudes below the energy dissipated in a 12.̋5 beam in W28 A2 but comparable to what was found by Gomez-Ruiz et al. (2013) for the entire outflows (associated with low-mass YSOs) L1448-IRS3 and HH211-mm (based on the use of more observational methods). From this method, it seems that BHR71 can be defined as an energetic low-mass outflow, but the relatively high numbers we found could be the effect of the method we used, based on a sophisticated shock model. From the energetic point of view, the impact of the whole BHR71 outflow on the Galactic ISM is much less than that, e.g. of the whole W44 SNR, because the dissipated energy is smaller, but also because of its much smaller size: the BHR71 outflow is roughly covered by a rectangle of ~0.1 × 0.5 pc2, whereas a SNR such as W44 resembles a circle of ~26 pc radius.

Table 3

Outflow parameters over the beam of our observations in H2, CO, and SiO.

We also used a second method, following Bontemps et al. (1996); Beuther et al. (2002), to calculate the outflow parameters related to the considered species (CO, H2 and SiO). These parameters (dynamical age td, mass M, mass entrainment rate , momentum P, mechanical force Fm, kinetic energy Ek, and mechanical luminosity Lmech), are calculated using the equations: where N is the column density of the considered species, R the radius of our analysis region, and νmax the approximate zero-base linewidth of the considered species. We used CO, H2, and SiO column densities extracted from our best-fit model associated with a filling factor of 1. The results are indicated in Table 3. Lacking the necessary spectral resolution for the H2 data, we assumed that the δνmax and td values for H2 are identical identical as for CO. This leads to overestimating the mass, momentum, and energy calculated based on the H2 data with respect to the results purely obtained from the model (first method). This is because the global shock velocity of our best-fit model, 22 km s-1, is less than the zero base linewidth that we attributed to H2 in this second method. Regardless, the CO mechanical luminosity we constrain is again consistent with that measured by Bourke et al. (1997) for the whole outflow (0.5 L), whereas the corresponding value we calculated for H2 exceeds that given by N09 for the fraction of the outflow that they mapped (0.05 L).

We can now compare the outflow parameters as inferred from our observations of BHR71 with similar studies found in the literature for outflows associated with protostars of various masses. An impressive number of such studies have been aimed at isolated targets, making it difficult to give some perspective on our results based on a one-to-one comparison (e.g. Leurini et al. 2006; Lebrón et al. 2006; Fontani et al. 2009; Liu et al. 2010; Guzmán et al. 2011; Cyganowski et al. 2011). We consequently focus on a comparison with “survey studies”, i.e. studies that are aimed at extracting outflow parameters from a sample of sources. From this respect, we found that the outflow parameters derived from CO observations of BHR71 are obviously systematically less than the values inferred in more massive environments (e.g. Duarte-Cabral et al. 2013; Klaassen et al. 2011), let alone in the outflows associated with O-type stars (López-Sepulcre et al. 2009). This is also true for SiO-related energetic parameters: BHR71 is less energetic from this point of view than its massive counterparts (Sánchez-Monge et al. 2013). Indeed these authors computed the mass, momentum, and energy associated with the SiO (21) and (54) lines. Their study was focussed on the whole outflows, resulting in values considerably larger than those inferred here (4-110 and 235  M; 262130 and 11440  M km s-1; (0.275) and (0.06 − 16) × 1046 ergs). On the contrary, the energetic parameters inferred from our study of one bright beam are similar to those inferred over the extent of the whole outflows associated with similar low-mass stars as IRS1, for instance. This is the case for the outflows in the Perseus cloud (NGC 1333, L1448, HH211, L1455, e.g. Curtis et al. 2010; Gomez-Ruiz et al. 2013). The conclusion is that either BHR71 is indeed more energetic than its low-mass counterparts, or that our method partly based on shock modelling naturally yields higher numbers than in all these studies, that often make use of fewer CO lines, for instance, than in the present work.

6. Conclusions

We have presented new observations of the BHR71 bipolar outflows, obtained with SOFIA in 12CO (1110), and (16–15), and with APEX in 12CO (32), (4–3), (6–5), (7–6), and 13CO (32).

We combined these data with existing datasets: in H2 (Spitzer/IRS) and SiO (APEX), and we have compared the observations in the form of a CO flux diagram, an H2 excitation diagram, and an SiO integrated intensity diagram to a grid of sophisticated shock models.

Our best fit is non-stationary (CJ-type) and has the following parameters: nH = 104 cm-3, b = 1.5, νs = 22 km s-1, and an age of 3800 years. The age and velocity of the shock model hint at the presence of more than one shock structure within the rather large beam of our observations, 24′′. This was also suggested by the fact that we had to assume different filling factor values for the different considered species. From the analysis of our ten best-fit models, we consider that the constrained values are quite robust, as these models all had nH = 104 cm-3, b = 1−2, νs = 20−25 km s-1, and an age of 35004500 years.

However this modelling can still be used to discuss the feedbacks of the shocks encompassed in our observations. We studied its chemical feedback in terms of SiO chemistry, placing an upper limit of 4% of the total silicon-bearing material in the form of SiO in the grain mantles in the pre-shock region.

We quantified the contribution of shocks to the excitation of CO around low-mass protostars surrounded by outflows, and shown that the CO excitation diagram from a shocked layer where the gas temperature is calculated point by point can be interpreted as a two-component one for levels from Jup = 12 to 40.

We inferred global outflow parameters from our shock model: a mass of 1.8 × 10-2 M was measured in our beam, in which a momentum of 0.4 M km s-1 was dissipated, corresponding to an energy of 4.2 × 1043 erg. We also analysed the energetics of the outflow species by species. Both methods suggest that BHR71 is a rather energetic outflow.

Three perspectives lie ahead of the present study. The first is to generalise our analysis to the whole outflow, and to observationally constrain the outflow energetic parameters over the whole mapped area. This will yield meaningful comparisons with observational studies of various similar objects. The second is to observe the SiO knot position with ALMA to resolve the multiple shock structures caught in our present single-dish beam. This should allow us to isolate a single shock structure, which would help us to get closer to the geometrical configuration that the Paris-Durham model was tailored to fit, and lead the shock analysis to a new level of understanding, both in terms of chemistry and energetics. Finally, at this point, more emission lines will also have to be included in our analysis, especially H2O lines.


1

This publication is partly based on data acquired with the Atacama Pathfinder EXperiment (APEX). APEX is a collaboration between the Max-Planck-Institut für Radioastronomie (MPIfR), the European Southern Observatory, and the Onsala Space Observatory.

3

GREAT is a development by the MPI für Radioastronomie and the KOSMA/Universität zu Köln (Kölner Observatorium für SubMillimeter Astronomie), in cooperation with the Max-Planck-Institut für Sonnensystemforschung and the Deutsches Zentrum für Luft- und Raumfahrt Institut für Planetenforschung.

Acknowledgments

We thank the anonymous referee for comments that helped to improve the quality of this work. We thank the SOFIA operations and the GREAT instrument teams, whose support has been essential for the GREAT accomplishments, and the DSI telescope engineering team. Based [in part] on observations made with the NASA/DLR Stratospheric Observatory for Infrared Astronomy. SOFIA Science Mission Operations are conducted jointly by the Universities Space Research Association, Inc., under NASA contract NAS2-97001, and the Deutsches SOFIA Institut, under DLR contract 50 OK 0901. This work was partly funded by grant ANR-09- BLAN-0231-01 from the French Agence Nationale de la Recherche as part of the SCHISM project. It was also partly supported by the CNRS programme “Physique et Chimie du Milieu Interstellaire”.

References

Appendix A: Filling factors

The half-maximum emission contours of the CO (65) and (32) lines were already shown in Fig. 2, in thick black and red contours, at their nominal resolutions of 9′′ and 18.̋1. For both the SiO peak and knot, the filling factor of the emission with respect to a beam of 24′′ (that of the CO (1110) line observations), inferred from the red, CO (32) contours is of the order of 1. As we performed our analysis over a 24′′ beam size, and not at the resolution of the CO (32) line, we decided to verify that this filling factor assumption was correct for all CO transitions at a 24′′ resolution. We hence convolved all the APEX data to this resolution. The result is shown in Fig. A.1. The dataset is remarkably consistent in terms of the size of the emitting region: the half-maximum emission contour for all lines is broadly the same, except for the (7–6) transition, which could be due to insufficient signal-to-noise values. Based on this figure, we chose to operate with a filling factor of 1 for all CO transitions observed in the SiO knot position over a beam of 24′′.

thumbnail Fig. A.1

Same field as in Fig. 2, shown in CO (65) (colours) convolved at the 24′′ resolution. The half-maximum contours of the emission of the CO (32) (black line), (4–3) (red), (6–5) (green), and (7–6) (blue) lines is also shown at the same 24′′ resolution. The SiO peak and knot, HH320 AB and knot 5 positions, as well as the beam sizes of the (6–5), (16–15) and (11–10) transitions are also indicated as in Fig. 2.

Open with DEXTER

Appendix B: CO, SiO, and H2 observables

Table B.1

CO and SiO integrated intensity values, and H2 level populations measured over a beam of 24′′ centred on the SiO knot position.

Appendix C: The CO flux diagram in the SiO peak position

thumbnail Fig. C.1

Observational CO flux diagram over a beam of 24′′ obtained in the SiO peak position. The CO (1615) line is not detected: an upper limit estimate is displayed in the form of a black arrow.

Open with DEXTER

All Tables

Table 1

Observed lines and corresponding telescope parameters for the APEX and SOFIA observations of BHR71.

Table 2

Shock model parameters.

Table 3

Outflow parameters over the beam of our observations in H2, CO, and SiO.

Table B.1

CO and SiO integrated intensity values, and H2 level populations measured over a beam of 24′′ centred on the SiO knot position.

All Figures

thumbnail Fig. 1

The BHR71 bipolar outflow in its entirety, in left: CO (65) (colours, associated with the wedge, in main beam temperature units, K km s-1) and (3–2) (white contours, from 30% to 100% of the signal in steps of 10%), both observed by the APEX telescope and integrated between −50 and +50 km s-1, with the resolutions indicated by green circles in the lower right corner; right: the 8 μm emission detected by the Spitzer/IRAC receiver (colours, N09), with the H2 0–0 S(5) emission as observed by the Spitzer/IRS receiver in the inner parts of the outflow (white contours, from G11), and the SiO (5–4) emission (red and black contours, G11) in the upper lobe (the green circles in the lower right corner show the respective resolutions of CO (3–2), (11–10) and SiO (5–4)). On both maps, the grey inset is the field shown in Fig. 2, the SiO peak and knot positions are indicated in blue, the knot 5 and HH320 region are in pink, and the IRS 1&2, knots 6, 8, and HH321 region are indicated by grey dots.

Open with DEXTER
In the text
thumbnail Fig. 2

Overlay of the velocity-integrated CO (65) (colour background) with the CO (32) (white contours) emission observed with the APEX telescope in the inner upper lobe of the BHR71 outflow. For both lines, the intensity was integrated between 50 and 50 km s-1. The wedge unit is K km s-1 in main beam temperature. The CO (32) contours are from 50 to 100% of the maximum, in steps of 10%. The half-maximum contours of the CO (32) and (65) maps are indicated in red and black, respectively. The dark blue circles indicate the positions and beam size of the SOFIA/GREAT observations for the CO (1110) transition. The APEX and SOFIA beam sizes of our CO (65), (1615) and (1110) observations are also provided (upper left corner light green circles, see also Table 1). The pink hexagons mark the position of the so-called knot 5 and HH320AB object.

Open with DEXTER
In the text
thumbnail Fig. 3

CO transitions in the SiO peak (left panels) and knot (right panels) positions indicated in Fig. 2: APEX (3–2), black line; (4–3), red line; (6–5), green line and histograms (upper panels); (7–6), dark blue lines; SOFIA (11–10), pink lines and (16–15), grey lines (lower panels, overlaid on the green histograms of the 6–5 transition). The last two were multiplied by five for comparison purposes. Respective spectral resolutions are 0.33, 0.50, 0.64, 0.54, 0.90 and 0.62 km s-1. The vertical dotted line marks the cloud velocity, −4.5 km s-1 (Bourke et al. 1995b).

Open with DEXTER
In the text
thumbnail Fig. 4

Comparison of the 12CO (black line) and 13CO (pink line) (32) emission profiles obtained in the SiO peak (upper panel) and knot (lower panel) positions. The ratio of main-beam temperatures is also shown in the form of red dots (for the channels where the 13CO detection is over 3σ), orange dots (for the channels where the 13CO detection is between 2 and 3σ), and grey dots (for the remaining channels). Associated with this distribution are the grey horizontal dotted lines, which show values of 20 and 40 for this ratio.

Open with DEXTER
In the text
thumbnail Fig. 5

Overlay of the map of CO (6–5) emission observed by the APEX telescope (colour background) with the H2 0–0 S(5) emission (white contours), observed with the Spitzer telescope. The wedge unit is K km s-1 (main beam temperature) and refers to the CO observations. The H2 0–0 S(5) contours are from 10% to 100%, in steps of 10%. The light blue contour defines the half-maximum contour of this transition. Like in Fig. 2, the blue circles and dots indicate the positions and beam sizes of the SOFIA/GREAT observations. The beam and pixel sizes of the CO (6–5) (11–10), (16–15), and H2 0–0 S(5) observations are the green circles in the upper left corner. The black contour delineates the half-maximum contour of the H2 0–0 S(2) transition. The field is the same as in Fig. 2, and the knot 5 and HH320 region are in pink. The field covered by Spitzer/IRS to observe the H2 emission is indicated in dashed grey line. It excludes the SiO peak position.

Open with DEXTER
In the text
thumbnail Fig. 6

Our model-observations comparisons. Upper panel: CO flux diagram over a beam of 24′′; the observations in the SiO knot position are the black squares, and the model results are the coloured circles (see text). The CO (1615) observational point is corrected for beam-size effect. The uncorrected point is the upper limit (see text), indicated by the grey arrow in this panel. Lower panel: H2 excitation diagram for the SiO knot position, extracted following the different procedures described in the text (empty symbols), global average in black squares, and model results in coloured circles (see text, with the same code as in the upper panel).

Open with DEXTER
In the text
thumbnail Fig. 7

The integrated intensity diagram generated from the G11 observations of three SiO lines (Jup = 5, 6, 8), corrected for beam size and for filling factor effect: filling factor of 1, 0.5, and 0.2 in black, dark grey, and light grey squares. The result of our best-fit model for the H2 and CO lines (nH = 104 cm-3, b = 1.5, νs = 22 km s-1, and an age of 3800 years) is shown with 1, 2, or 4% of the pre-shock Si placed in the grain mantles in red empty, dotted, or filled circles.

Open with DEXTER
In the text
thumbnail Fig. 8

The CO excitation diagram produced from our best-fit model (red squares). Two temperature components can be fitted to our shock model over the PACS range of observations: a warm component (Tex = 216 K, dark blue line) fitting the level populations for Jup = 12 to 25 (light blue squares), and a hot component (Tex = 422 K, light green line) fitting the level populations for Jup = 26 to 40 (dark green squares).

Open with DEXTER
In the text
thumbnail Fig. A.1

Same field as in Fig. 2, shown in CO (65) (colours) convolved at the 24′′ resolution. The half-maximum contours of the emission of the CO (32) (black line), (4–3) (red), (6–5) (green), and (7–6) (blue) lines is also shown at the same 24′′ resolution. The SiO peak and knot, HH320 AB and knot 5 positions, as well as the beam sizes of the (6–5), (16–15) and (11–10) transitions are also indicated as in Fig. 2.

Open with DEXTER
In the text
thumbnail Fig. C.1

Observational CO flux diagram over a beam of 24′′ obtained in the SiO peak position. The CO (1615) line is not detected: an upper limit estimate is displayed in the form of a black arrow.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.