© ESO, 2011

1. Introduction

The propagation of shock waves is an ubiquitous phenomenon in the interstellar medium, where the sound speed is low due to the temperature conditions. Observations over the past few decades (Snell et al. 1980) have shown that, at the early stages of the formation of solar-type stars, the process of mass accretion is almost always associated with mass ejection, in the form of collimated jets, extending from a few astronomical units up to parsecs from the exciting source (Dougados et al. 2000; Bally et al. 1996). The supersonic impact between the jet and the parent cloud generates a shock front, which propagates in the collapsing interstellar gas, and also an inverse shock that propagates along the jet itself. Large cavities, called bipolar outflows, appear, which have been extensively studied through the molecular emission that they generate (Bachiller 1996). At the apices of these cavities, the shock wave heats, compresses and accelerates the ambient interstellar gas. As the temperature rises to at least a few thousands degrees, the energy barriers of numerous reactions involving neutral or ionized molecules can be overcome and the chemistry of certain species can be significantly altered (Bachiller et al. 2001). Similarly, processes specific to the propagation of shock waves affect the interstellar grains (Flower & Pineau des Forêts 2003; Guillet et al. 2007, 2009), leading to the injection of molecular and atomic species into the gas phase, such as silicon-bearing species (Schilke et al. 1997; Gusdorf et al. 2008a, hereafter G08a). The time-scales involved in the heating and in some of the shock chemical processes are short (10² to 10⁴ years; Gusdorf et al. 2008b, hereafter G08b), so the shocked region rapidly acquires a chemical composition distinct from that of the quiescent medium. As the gas cools down radiatively, reactions with high energy barriers no longer operate, and the chemistry is dominated again by low-temperature reactions. The resulting molecular emission can be used as a diagnostic tool to study the physical and chemical characteristics of the shocked region.

Recently, G08b placed constraints on shock parameters using H₂, SiO and CO observations in the B1 knot of the L1157 outflow. At the same time, the HIFI and PACS receivers on Herschel have observed the entire L1157 outflow (Nisini et al. 2010) and provided velocity-resolved data of the L1157-B1 knot (Lefloch et al. 2010; Codella et al. 2010). Their studies indicate that the H₂, SiO and H₂O emission regions are spatially coincident on large scales. They also shed light on the necessity to provide reliable models of water emission in order to better understand its production in shocks and its role in the cooling of the outflow. Recently, Flower & Pineau des Forêts (2010, hereafter FP10) presented an extensive study of water emission for a small grid of stationary shock models, and compared it to Spitzer observations of water emission in NGC 2071, with a certain level of agreement. Our aim is to push the confrontation between models and observations one step further, by using the models, which have already demonstrated their ability to reproduce some observations, to predict water line intensities. Using an approach to comparing models with observations that is similar to that of G08b, as applied to H₂ and SiO, we aim to constrain the physical and chemical parameters of the targeted region. In a second step, we incorporate these preliminary constraints and the latest compilation of molecular data for H₂O in an LVG radiative transfer code to provide predictions of water emission, under the assumption that the emission regions of the various species coincide.

The choice of target is of crucial importance, as our candidate must have been the object of prior studies, and must be a good potential target for contemporary and future instruments such as Herschel and ALMA. Located at a distance of about 175 pc, the Bok globule BHR71 (Bourke et al. 1995) is a well-known example of isolated star formation observable in the southern hemisphere, and hence is a reasonable candidate for future ALMA observations. Its observation by the Herschel telescope is guaranteed through the key-program WISH, and also DIGIT. The associated outflow region contains signatures of two distinct molecular jets (Bourke et al. 1997; Myers & Mardones 1998; Parise et al. 2006) centered on two different protostellar sources, IRS1 and IRS2 (Bourke 2001) separated by  ~3400 AU. Bright HH objects are associated with the blue-shifted lobes of each outflow, HH320 and HH321 (Corporon & Reipurth 1997). They have been imaged in the S ii transition at 6711 Å, indicating that at least part of the outflows are dissociative. The dynamical age of HH321 is estimated to be 400 years (Corporon & Reipurth 1997). It has not yet been possible to determine the dynamical age for HH320. NIR JHK-band spectra covering both HH320 and HH321 were obtained by Giannini et al. (2004). Through detailed shock modelling they derived a pre-shock density of 10⁴ cm^-3 and shock velocity of 41 km s^-1 for HH320A; this was found by fitting a non-steady-state shock model to the observed H₂ lines. The age of the non-steady-state shock is 475 years, which is more or less in agreement with the dynamical age of the HH321 flow. The pre-shock density is lower than is predicted on the basis of CO observations, 10⁵ cm^-3 (Parise et al. 2006). However, the latter is the density of the molecular outflow, which is compressed relative to the ambient pre-shock cloud. These inner parts of the bipolar outflow, covering the HH objects, have also been targeted by the Spitzer Space Telescope, and the corresponding results in terms of rotational H₂ excitation diagrams have been made available by Neufeld et al. (2009). Another article is dedicated to the parameters map analysis based on these data (Giannini et al. 2011). Finally, Garay et al. (1998) have presented maps of SiO (2–1) and (3–2) of the whole bipolar outflow.

In Sect. 2, we present a combined analysis of both rovibrational transitions of H₂ from Giannini et al. (2004) and rotational lines, as compiled by Spitzer (Neufeld et al. 2009). We then briefly present the code that is used to simulate the propagation of shock waves through the interstellar medium, and introduce the grid of models that we computed. Finally we describe the constraints that we obtain in several knots of the BHR71 outflow by comparing more or less complete H₂ observations with models. In Sect. 3, we present the results of comparisons between observations of SiO and the models. We describe the observing setups, the actual observations, and the results in terms of simultaneous fits of the H₂ and SiO molecules. Sect. 4 contains our predictions of water emission line intensities. We briefly summarize our methodology and then discuss the appropriateness of the steady-state and statistical equilibrium assumptions, before presenting a comparison of results obtained using each of these approximations. Section 5 contains our concluding remarks.

2. H₂ observations and modelling

2.1. Observations

There are two distinct sets of H₂ observations. The details of the ground-based observations of the rovibrational transitions are reported in Giannini et al. (2004). The Spitzer observations of the pure rotational transitions, described in Neufeld et al. (2009), are the subject of an separate article (Giannini et al. 2011).

The comparison between observations of H₂ emission and model results is usually based on excitation diagrams that are derived for selected emission regions. Such excitation diagrams display ln(N_vj/g_j) as a function of E_vj/k_B, where N_vj (cm^-2) is the column density of the rovibrational level (v,J), E_vj/k_B is its excitation energy (in K), and g_j = (2j + 1)(2I + 1) its statistical weight (with I = 1 and I = 0 in the respective cases of ortho- and para-H₂). If the gas is thermalized at a single temperature, all the points in the diagram fall on a straight line.

In their original study of conditions in the BHR71 outflow, Giannini et al. (2004) extracted the excitation diagrams corresponding to nine positions (knots 1 to 9) and four HH objects (HH320A/B and HH321A/B), in order to compare their observations with shock models. In the present study, we attempt to make use of most of the available data, both in the pure rotational and in the rovibrational transitions of molecular hydrogen. Unfortunately, combining the two sets of observations is only possible where both sets of data are available and when a cross calibration can be performed to ensure their consistency. We consequently make use of the following excitation diagrams :

• pure rotational transitions for the “H₂ knots” 5, 6, 8 as denoted by Giannini et al. (2004);

• rotational and rovibrational transitions for the two Herbig-Haro objects that are called HH320AB and HH321AB, as defined by Giannini et al. (2011);

• pure rotational transitions on the “SiO knot” (see Sect. 3.2).

The corresponding list of coordinates is provided in Table 1. All these points are presented in Fig. 1, which shows the internal parts of the BHR71 outflow, as mapped by the Spitzer Space Telescope in the 0–0 S(5) spectral line of H₂ (Neufeld et al. 2009). In each case, the excitation diagrams were extracted for a circular emission region of 5′′ radius, centered on these positions, shown for the case of the so-called “SiO peak” (see Sect. 3.2) in Fig. 1 (smallest circle). In the case of the Herbig-Haro objects, the consistency between the rotational and rovibrational data sets was ensured by means of the following method, as explained by Giannini et al. (2011). First, the image of H₂ 1–0 S(1) line at 2.122 μm (retrieved from the ESO archive at http://archive.eso.org/eso/eso_archive_main.html) is degraded to the Spitzer spatial resolution of 2′′/px. Then, this image is calibrated by using 2MASS fluxes of several bright stars in the field, and the photometry is computed in the same areas of 5′′ radius as the Spitzer images. To find the intercalibration factor between imaging and spectroscopy, the measured 2.12 μm flux was compared with the flux measured in the 1′′ slit where spectroscopic observations have been performed. Finally, the same factor was applied to the fluxes of the other NIR lines, adopting averages of the line ratios (relative to the 2.12 μm transition) in the A and B substructures of HH320 and HH321.

Fig. 1 The internal parts of the BHR71 outflow as observed by Spitzer in the H2 0–0 S(5) transition (Neufeld et al. 2009); the colour scale unit is erg cm-2 s-1 sr-1. The contours correspond to the levels of 3, 9, 15, 21, 27, and 33σ. The pink hexagons give the positions of the IRS1 and IRS2 sources as determined by Bourke et al. (2001). The positions of the Herbig-Haro objects HH320AB and HH321AB are plotted in black as well as those of the SiO peak (referring to the point of maximum SiO (5–4) emission, see Sect. 3), and of the SiO knot, which was chosen to perform our simultaneous fit of H2 and SiO emission, and also to make predictions of water emission. The other emission knots that were analysed are plotted in white (knots 5, 6, and 8, following the labelling of Giannini et al. 2004). Around the SiO-peak, four radii are displayed, corresponding to 5, 9, 12, and 14.4′′. The first represents the size of our region of analysis. The latter three correspond to the beam of the APEX telescope for the SiO (8−7), (6−5) and (5−4) transitions, respectively; see text.

The data were corrected for extinction, adopting the upper limit of the visual extinction from Giannini et al. (2004), A_V = 2, and using the interstellar extinction law of Rieke & Lebofsky (1985). We note that the effect of the interstellar extinction on the observational points is minor. The excitation diagrams resulting from the observations are plotted in Fig. 2. The rotational and rovibrational temperatures derived from these observations are also indicated on the plots in the same colour code, along with the corresponding line. Again, the reddening correction has only minor significance on the temperature values. As anticipated, the rovibrational transitions are found to trace warmer gas than the pure rotational transitions (see, for example, Le Bourlot et al. 2002; Giannini et al. 2006; and G08b).

Fig. 2 The H2 excitation diagrams observed for each position indicated in Fig. 1, in black (pure rotational lines), and grey (rovibrational lines, when available). The data were corrected for extinction, using AV = 2 from Giannini et al. (2004). The inferred rotational temperatures are indicated in black, and the rovibrational temperatures, in grey. The best-fitting model results are also shown in each case, in blue and/or red symbols. The shock parameters are given in each panel in the corresponding colour.

2.2. Grid of models

The shock code that we use to simulate the propagation of shock waves in the interstellar medium is described in G08b. It generates one-dimensional stationary models of both continuous or discontinuous kind (respectively, C- or J-type shocks), and also approximations to time-dependent solutions, in the form of a shock wave possessing the characteristics of both C- and J-type models (see Chièze et al. 1998; Lesaffre et al. 2004a,b). The shock code solves the differential equations that determine the evolution of a set of dynamical and chemical variables: the neutral and ionized fluid temperatures, velocities and densities; the fractional abundances of over 125 species, linked by more than 1000 chemical reactions. The populations of the first 150 levels of molecular hydrogen are calculated within the shock code, based on the treatment of Le Bourlot et al. (2002), with the only difference being that we adopted the most recent H–H₂ collisional rate coefficients (Wrathmall et al. 2007). Finally, grain-related processes are also included in the code, following Flower & Pineau des Forêts (2003).

The parameters that we aim to constrain are the inputs of the models, namely:

• ambient-medium-related quantities: pre-shockdensity n_H and magnetic-field parameter b, defined by B(μG) = $\mathit{b} \mathrm{\times }\sqrt{{\mathit{n}}_{\mathrm{H}} \mathrm{\left(}{\mathrm{cm}}^{-3}\mathrm{\right)}}$;

• shock-related quantities: shock velocity v_s, and age in the case of a non-stationary solution;

• chemistry-related quantities, such as the initial value of the ortho-para ratio for molecular hydrogen, or the initial distribution of silicon-bearing material in the grains.

As a first step, we made use of a grid of stationary C-type shock models with the parameters specified in Table 2. We should mention that the parameter coverage is not complete, in the sense that 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. Following the method presented in G08b (Sect. 4.1), this 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 H₂ observational data. The grid of CJ-type models covers the combination of parameters given in Table 2, for ages ranging from a few tens to a few thousands of years.

For each of the above models, we obtained an excitation diagram that can be compared with observations.

2.3. Comparisons and initial constraints

We describe here the results of the comparisons between observations and our models. A few preliminary remarks about our results are applicable to all the knots and regions analysed. First, the filling factor was considered to be equal to 1 for all transitions; this yields H₂ knot sizes consistent with the thickness of the H₂ emitting layer of our models. Secondly, under this assumption, we find that, in all cases, only non-stationary shock models reproduce satisfactorily both the pure rotational and rovibrational (when available) parts of the H₂ excitation diagram, a result that has already been established in similar studies of other objects (e.g. HH54, Giannini et al. 2006; or L1157-B1, see G08b). Finally, the departure of the observational points from a straight line in the purely rotational parts of the excitation diagram (“saw-tooth” pattern, which is weak in our cases) is indicative of deviations of the ortho-to-para ratio from the high-temperature LTE value of 3.0. Such a pattern can be reproduced in the models by slightly modifying the value of the initial ortho-to-para ratio; but fine-tuning this parameter is not a major consideration in the context of the present study.

We give prior focus to the Herbig-Haro objects, for which we can benefit from H₂ excitation diagrams combining pure rotational and rovibrational transitions, and to the so-called “SiO knot”, for which future observations of SiO rotational transitions are expected to provide further constraints on our shock models. Tables A.1 and A.2 present a summary of the restrictions that comparisons between observational and computed excitation diagrams place on the shock model parameters. In these tables, the values of the pre-shock density, magnetic-field parameter and shock velocity are indicated for our best models, as well as the age range for which the excitation diagrams are in satisfactory agreement with the observations. In the case of the Herbig-Haro objects HH320A/B and HH321A/B, the best-fitting models are non-stationary, and the pre-shock density is constrained to 10⁴ cm^-3, which is in agreement with the average density of the whole globule (9 × 10³ cm^-3), as determined by Bourke et al. (1997). Shock velocities in the range 15 − 35 km s^-1 yield good fits to the observations, depending on the associated magnetic-field parameter and age. We note that the shock ages that we obtain are of the same order of magnitude as the observational dynamical ages. The excitation conditions are remarkably similar in both sources, in spite of their location on a different lobe of the outflow, and the probable mixing of the two outflows from IRS1 and IRS2. In the case of HH320A/B, our findings are in agreement with those of Giannini et al. (2004) in terms of pre-shock density, magnetic field, and shock type and age, but our value of the shock velocity is smaller than theirs. Several factors could account for this difference, which might not be a real discrepancy, given the degeneracy of our fitting models. The first factor is that their value lies slightly outside the range of our parameter search. In addition, we studied HH320A/B, rather than HH320A. Furthermore, we added the pure rotational transitions of H₂ to the analysis. Finally, our shock model incorporates an update of the H–H₂ collisional rate coefficients, following Wrathmall et al. (2007).

The parameters of the SiO knot are not as well constrained, as only the pure rotational observations are available to build the H₂ excitation diagram. Similarly to the HH objects, the best fitting models are of non-stationary shocks with velocities in the range of 15 − 35 km s^-1 for a pre-shock density of 10⁴ cm^-3 (with ages dependent on the value of the magnetic-field parameter). Additionally, a few models with a pre-shock density of 10⁵ cm^-3 are found to fit the observations equally satisfactorily, with even more moderate shock velocities (10–20 km s^-1) and lower shock ages (which again depend on the value of the magnetic-field parameter). Examples of the comparisons between observations and models are shown in the three upper panels of Fig. 2, for the cases of the SiO knot, HH320A/B and HH321A/B, respectively. In each panel, the parameters corresponding to the selected shock model are also indicated.

For the three other points studied (“knot 5”, “knot 6”, and “knot 8”), once again the excitation diagrams do not constrain the parameters as tightly as in the HH regions, since only the pure rotational transitions of H₂ can be used. The number of models that fit these data is large, but we show some examples of the best fits that we could derive in the lower panels of Fig. 2. We note that the modelling work conducted for these positions also suggests similar excitation conditions, regardless of their association with one lobe of the outflow or the other.

3. SiO observations and modelling

3.1. Observations

The observations were carried out mainly in two observational runs, in August and October 2009, with the Atacama Pathfinder EXperiment (APEX¹, Güsten et al. 2006). We used the facility’s own APEX–1 and APEX–2 receivers (see respectively Vassilev et al. 2008, Risacher et al. 2006) to observe respectively the SiO (5–4) and (6–5) transitions, and the SiO (8–7) transition, in combination with the MPIfR Fast Fourier Transform Spectrometer backend (FFTS, Klein et al. 2006). Focus was checked at the beginning of each observing session, and after sunrise on Mars. Line pointing was locally checked on 07454–7112 or IRAS 15194–5115. The pointing accuracy was found to be of the order of 5′′ rms. Table 3 contains the main characteristics of the telescope for the observed transitions:  beam sizes, forward and beam efficiencies, typical system temperatures, and the spectral resolution that was used for each transition. The observations were performed in position-switching mode using the APECS software (Muders et al. 2006), with a reference position at offset ( − 600′′, 600′′) from the central position indicated in Table 1. The data were reduced with the CLASS software (see http://www.iram.fr/IRAMFR/GILDAS).

3.2. Results

The observations resulted in one SiO (5−4) map of the inner, upper (red) lobe, as well as in SiO (6−5) and (8−7) spectra, obtained by pointing the telescope at designated positions inside this lobe. The absolute coordinates of the (0′′, 0′′) position are given in Table 1.

Fig. 3 Upper internal part of the BHR71 outflow. The background map is the same as in Fig. 1. The SiO (5–4) contours, which are overlaid in a dark blue colour, correspond to levels of 6 to 30σ, in steps of 3σ, for an integration between  − 10 and 60 km s-1. The half-maximum, 15σ contour is displayed in red. The associated beam size is indicated in the upper left corner. Our SiO (5 − 4) coverage is actually wider than the field displayed in this figure.

The map of the SiO (5–4) line covers a roughly rectangular region, extending from the lower left offset (35′′, 25′′) to the upper right (−60′′, 170′′), and was performed in raster mode, spending 30 seconds on each point. The beam sampling was 10′′, which is less than the half beam size of 14.35′′, meaning the map is fully Nyquist sampled. The result is plotted in Fig. 3: the SiO emission, integrated between  −10 and 60 km s^-1, is displayed in dark blue contours overlaid on the H₂ 0 − 0 S(5) map already shown in Fig. 1. The half-maximum signal contour is shown in red. The rms of the SiO (5−4) map was found to be 0.26 K km s^-1 at the resolution 2.02 km s^-1, and the contours correspond to levels of 6 to 30σ, in steps of 3σ. The contours show that the SiO (5–4) emission peaks outside of the Spitzer spectral line maps, at a position that will be referred to as the “SiO peak”, with offset (−15′′, 87′′). In order to perform a simultaneous study of the H₂ and the SiO emission, the point to be studied must be located inside the Spitzer map: we chose the point of maximum emission of H₂, at an offset of (0′′, 77′′), to which we refer as the “SiO knot”. We consider the shift between the APEX SiO (5−4) and Spitzer H₂ 0 − 0 S(5) emission peaks to be real, as it is far too large to be attributable to pointing errors. Indeed, the distance between the south-eastern and the north-western tips of the half-maximum SiO emission (roughly of respective coordinates (0′′, 70′′) and ( − 30′′, 110′′), see Fig. 3) is about 50′′, which amounts to 40′′ once deconvolved from the SiO (5 − 4) telescope beam size. Consequently, we consider that, in this direction, which coincides with the line between our SiO knot and peak (see Fig. 3), the emitting region is extended enough to allow us to consider that the signal observed at the position of the SiO knot is free of any residual signal from the SiO peak. Moreover, from the extension of the Spitzer maps, it is impossible to exclude the presence of another, potentially stronger peak of H₂ emission, coinciding with the SiO peak. In order to check, we overlaid our map on the four (way more extended) channel maps of IRAC, especially the fourth, at 8 microns, whose broadband intensity is dominated by the 0–0 S(4) and S(5) transitions of H₂ in this region (Neufeld et al. 2009). We found that these overlays fail to confirm the existence of an H₂ peak coinciding with that of the SiO (5–4) emission.

Fig. 4 The spectra obtained in the three observed SiO transitions at the positions of the SiO peak and knot, of respective offsets ( − 15′′, 87′′) and (0′′, 77′′). Only the SiO knot is considered in our modelling. For each spectrum, we adopted the SiO (5–4) spatial resolution. The ambient cloud velocity, –4.5 km s-1 is indicated by the vertical dashed lines.

The spectra that we obtained of the SiO (5–4) transition at both the SiO peak and knot are presented in Fig. 4, along with the spectra of the other two SiO transitions. In the case of the SiO (5–4) transition, the original beam size was 28.7′′. In order to derive the spectra at the two selected positions as exactly as possible, we resampled the data with a spatial resolution of 5′′ and then averaged the four spectra closest to the two positions to obtain the results that are presented in Fig. 4. In the case of the SiO (6–5) and (8–7) transitions, in addition to this spatial resolution problem, we had also to convolve the data to the SiO (5–4) beam size in order to render the whole set of transitions consistent. After this preliminary convolution, we obtained the spectra displayed in Fig. 4. Our mini-maps around the two positions were extended enough to cover the larger beam of the SiO (5–4) observations.

The shape of the line profiles is typical of bipolar outflows (similar, for instance, to those acquired around L1157: see Nisini et al. 2007), peaking near the cloud velocity and with the main signature being the line wings extending to about 50 km s^-1, most visible in the SiO (5–4) line profile of the peak. Such wings seem to be present, though narrower, at higher frequencies. However, it is then more difficult to estimate their width, as the signal-to-noise ratio decreases as the observed frequency increases. The observed SiO probably arises from a complex mixture of shocks with different velocity projections on the line of sight. The velocity range of our shock models does not exceed the terminal velocity of 50 km s^-1. On the SiO (5–4) and (6 − 5) spectra associated to the SiO peak, a velocity bullet seems to arise between 20 and 50 km s^-1. It could arise from a process that alters the SiO emission within the region of the shock considered, such as an abundance effect, or from the presence of another shock, with a different velocity, within the telescope beam. These possibilities are discussed in Sect. 3.3 but cannot be thoroughly investigated, given the large beam sizes associated with these transitions; they call for a more detailed study, with higher spatial resolution.

The last step was to derive integrated intensities from these spectra. The results for the SiO knot, which is considered in our subsequent analysis, are shown in Table 4: integrated intensities and associated rms are given. The velocity interval over which the integration was made was [−10, 60] km s^-1 for all three transitions. Again, the values inferred are consistent with results for other outflows, such as L1448 or L1157 (Nisini et al. 2007), for which the central source is also a young Class 0 proto-star with a similar bolometric luminosity (respectively, around 13.5, 11, and 1.5−3.5 L_⊙ for BHR71, L1157, and L1448; see Neufeld et al. 2009). In particular, the line profiles we observe for SiO (5−4), (6−5), and SiO(8−7) are remarkably similar to those of L1157-B1, peaking near the cloud velocity. The integrated intensities we inferred are also very close to those compiled for each knot of the L1157 bipolar outflow (Nisini et al. 2007).

3.3. Modelling results

We derive emission-related quantities by means of a radiative transfer code used in combination with the shock models presented above. The radiative transfer code is based on the LVG approximation, and has been described and used in G08a and G08b. At each point of the shock model, it uses the computed parameters of relevance, such as the temperature, density, velocity gradient, and collision-partner density, to solve the equations of statistical equilibrium and calculate emission-related quantities, such as the level populations, optical depths, and local emissivities. The set of equations is solved by means of a simple lambda-iteration, that stops when the relative difference between each of the level populations computed in two successive iterations drops below a convergence value of 10^-4. The set of collisional rate coefficients for SiO with H₂ as a collision partner is simply scaled by the square root of the reduced mass ratio from the SiO–He rate coefficients calculated by Dayou & Balança (2006). The LVG code also provides the integrated intensity over the whole shock model for all transitions considered.

The silicon chemistry is extensively discussed in G08a and G08b. Regarding the initial distribution of silicon-bearing material, G08a made the assumption that Si is initially present exclusively in the grain cores, in the form of silicates, with a fractional abundance of 3.35 × 10^-5 relative to total H. G08b introduced an alternative scenario, in which the Si is partly in the grain mantles, in the form of SiO. Such an assumption was found to be necessary to fit both the SiO and H₂ line intensities by means of a single shock wave model. Furthermore, a simultaneous fit can be achieved only by means of a non-stationary, CJ-type shock. When computing CJ-type models in the present study, we adopt one of three values of the fraction of elemental silicon initially in the form of SiO in the grain mantles, namely 0.00, 0.01, or 0.10.

In the first case, when no SiO is initially present in the grain mantles, the gas-phase silicon is produced by the erosion of the grain cores only, which occurs, in the C- and CJ-shock models, where the predominantly charged grains collide with the neutral particles at the ion-neutral drift speed (ambipolar diffusion). Three reaction paths then lead to the formation of SiO in the gas phase, namely Si⁺(OH, H)SiO⁺(H₂, H)HSiO⁺(e, H)SiO, Si(OH, H)SiO, and Si(O₂, O)SiO; the last two reactions predominate. Following G08a, we adopted the rate coefficient that was experimentally measured for the Si+O₂ reaction by Le Picard et al. (2001) for both Si+O₂ and Si+OH. In the second case, in which SiO is initially present in the grain mantles, their rapid erosion by sputtering in the shock wave generates SiO directly in the gas phase (see also the discussion in G08b).

O₂ is a key reaction partner in the process of formation of SiO in the gas phase, and so its abundance must be modelled to the greatest possible accuracy. Furthermore, the question of the pre-shock abundance of O₂ in the gas phase must be considered carefully. In all cases, as in G08a and G08b, we used an initial gas-phase abundance of n(O₂)/n_H = 1.0 × 10^-7, consistent with the O₂ observations of the Odin satellite (see Pagani et al. 2003; Larsson et al. 2007). On the other hand, the initial fractional abundance of O₂ in chemical equilibrium is 1.0 × 10^-5; we allocated the difference between these two values (i.e. most of the initial O₂) to the grain mantles, consistent with G08a and G08b. Given that the O₂ in the mantles is released rapidly into the gas phase, by sputtering processes induced by the shock wave, this scenario does not differ significantly from that in which all of the O₂ is assumed to be present initially in the gas phase (the chemical equilibrium situation). An alternative hypothesis, in which the same fraction of the corresponding elemental oxygen is assumed to be in the form of water ice, rather than O₂ ice, was considered by G08a and found not to modify significantly the results of the models. Indeed, in this case, the rapid erosion of the grain mantles during the passage of the shock wave releases H₂O into the gas phase, and the subsequent reaction of H₂O with H produces OH, which then reacts with Si to form SiO. A third scenario was put to the test in the present study: placing the elemental oxygen in the grain mantles in the form of CO₂ ice. In this case, CO₂ is released quickly into the gas phase during the passage of the shock wave, but no chemical reactions convert it to O₂ or OH. When SiO is present in the grain mantles, direct sputtering is the dominant source of the gas-phase SiO. However, when SiO is absent from the grain mantles, erosion of the silicon-bearing grain cores is the mechanism that releases silicon into the gas phase, and the deficit in O₂ or OH leads to lower abundances of gas-phase SiO and weaker SiO line emission from low-velocity shocks. In high-velocity shock models, or at higher pre-shock densities, O₂ and OH are destroyed in the gas-phase, through the O₂(H, O)OH and OH(H, O)H₂ sequence, and so the initial composition of the grain mantles is less significant.

Fig. 5 Examples of fits of the integrated intensities of our three transitions of SiO. On the y-axis is the integrated intensity (K km s-1), and on the x-axis is the rotational quantum number of the upper level of the transition. The observations are the black squares linked by the thick black lines, on all three panels. The 3σ errorbars are shown in the figure. The model results are displayed in colour (blue and red), with the indication of the specific shock parameters on each panel. a) C-type shock models: best fit shock parameters nH = 104 cm-3, no SiO initially in the grain mantles; b) CJ-type shock models: best fit shock parameters nH = 104 cm-3 and vs = 28 km s-1, with 1% of elemental silicon initially in the grain mantles, as SiO; c) CJ-type shock models: best fit shock parameters nH = 104 cm-3 and vs = 30 km s-1, with 10% of the silicon initially in the grain mantles as SiO.

If one ignores the observations of H₂, it is possible to fit the observations of SiO by means of C-type shock models, without any SiO initially in the grain mantles; cf. G08a. In this case, we find that shock models with a pre-shock density of 10⁴ cm^-3 provide better fits than those with pre-shock densities of 10⁵ or 10⁶ cm^-3. The typical shock velocities are from 25 to 30 km s^-1, with values of the magnetic-field parameter from 0.45 to 2. An example of such a fit can be seen in panel (a) of Fig. 5. We note that, in this case, the assumption that most of the O₂ is initially in the form of CO₂ ice on the grain mantles leads to insufficient levels of SiO emission, as compared with the observations. However, our aim is to model simultaneously the emission of SiO and H₂, assuming their emitting regions coincide, to within the spatial resolution of the observations. Accordingly, we calculate the SiO spectrum for the shock models that already fit the H₂ lines, namely the non-stationary shock models presented in Sect. 2.3. Given the the additional flexibility afforded to these models by the initial fractional abundance of SiO in the grain mantles, we found that most of the CJ-type models that already fit the H₂ data can reproduce the SiO emission satisfactorily. The results are presented in Fig. 5, where panels b) and c) correspond to our best fitting models with, respectively, 1% and 10% SiO initially in the mantles. We find that no model in which SiO is absent from the grain mantles is able to reproduce the intensities of the SiO lines simultaneously with those of H₂, thereby confirming a conclusion of G08b and rendering considerations of the initial chemical form of oxygen irrelevant.

We wish to emphasize that a single, plane-parallel CJ-type shock wave model is unlikely to provide a completely satisfactory simulation of the molecular line emission from an object as complex as BHR71. In particular, a combination of C-type shock models with appropriate parameters and filling factors might also be able to reproduce the H₂ and SiO line intensities and address the problem of fitting the SiO line profiles. The latter could be strongly affected by geometrical effects, associated with the projection of multiple shock velocities on to the line of sight. However, such a study introduces additional free parameters that are insufficiently constrained by the currently available observations. Furthermore, our model does not take into account the consequences of grain-grain collisions, whose inclusion (following Guillet et al. 2007, 2009, 2011) would have major implications for computation times. In view of these limitations, geometrical and physical, it is not feasible currently to investigate the nature of the bullet that we observe in the SiO(5–4) line profile.

4. Water emission predictions

4.1. Water modelling

The abundance of water at each step of the calculation is the result mainly of two formation pathways. One is the O(H₂, H)OH(H₂, H)H₂O sequence, which takes place in the gas phase and becomes significant when the temperature exceeds a few hundred kelvin (e.g. Charnley 1997; Atkinson et al. 2004). The other is the sputtering of water ice that formed in the mantles from O and H atoms, which were adsorbed on to the grains in the parent dark cloud (see for example Tielens & Hagen 1982). In the latter case, the sputtering process is analogous to that considered for SiO when present in the grain mantles.

The modelling of water emission was performed in a similar manner to that of SiO. Again, we used an external radiative transfer module based on the LVG approximation; the molecular data were taken from the LAMDA database (http://www.strw.leidenuniv.nl/~moldata/). Regarding the rate coefficients for excitation of H₂O by He, we used the data of Green et al. (1993). For collisions of H₂O with H₂, we adopted the same approach as FP10 and used the CTMC results of Faure et al. (2007), rather than the more precise but less complete (in terms of their temperature range) quantal calculations of Dubernet & Grosjean (2002), Grosjean et al. (2003), and Dubernet et al. (2009). As no data exist for the excitation of H₂O in collisions with H, we used the same rate coefficients as for excitation by ortho-H₂, following FP10.

The more complex structure of H₂O, as compared to SiO, results in some changes being required to the computational method in the radiative transfer module. First, the code must be modified to take account of the more complicated radiative selection rules. Second, in the more complex case of water as compared to SiO, a value of 10^-10 is required for the convergence parameter. Finally, we found it necessary to adapt the expression that is used for the escape probability. For a plane-parallel geometry, such as is assumed in our shock models, the escape probability is ${\mathit{\beta }}_{\mathrm{\perp }}\mathrm{=}\frac{\mathrm{1}\mathrm{-}{\mathrm{e}}^{\mathrm{-}{\mathit{\tau }}_{\mathrm{\perp }}}}{{\mathit{\tau }}_{\mathrm{\perp }}}\mathit{,}$(1)where τ _⊥  is the LVG optical depth in the z-direction, normal to the shock front. On the other hand, Neufeld & Kaufman (1993) proposed the following approximation to the angle-averaged escape probability, $\mathit{\beta ̅}\mathrm{av}\mathrm{=}\frac{\mathrm{1}}{\mathrm{1}\mathrm{+}\mathrm{3}{\mathit{\tau }}_{\mathrm{\perp }}}\mathit{,}$(2)in the case of water, we found that the use of the second formula gave rise to numerical instabilities and longer computation times. Consequently, we adopted the first formula, both when solving the equations of statistical equilibrium and when evaluating the local line emissivity at each point in the shock wave. This approach is consistent with that of FP10 and enables us to make direct comparisons with results obtained using their method of treating radiative line transfer in shock waves. The implications of the different forms of the escape probability for the intensities of the rotational transitions of SiO are discussed in Appendix A of G08a. We recall that the geometry of the shock wave is not well known.

Before comparing with the method of FP10, we benchmarked our LVG code against the RADEX code (van der Tak et al. 2007), at four representative points, covering different conditions of temperature, density, and velocity gradient, namely (respectively in K, cm^-3, and s^-1): (25, 10⁴, 2 × 10^-13), (50, 10⁵, 10^-13), (300, 10⁴, 10^-12), and (467, 3.71 × 10⁵, 9 × 10^-11). For the purpose of this comparison, we adopted the same escape probability formula as in the RADEX code. For ortho-H₂O, the computed excitation temperatures showed no mean relative discrepancy above 6.3% for the ten levels emitting the strongest lines at each of the four points. In the case of para-H₂O, the relative discrepancy was occasionally larger, reaching 43.9%. Although it remains unclear why the discrepancies were larger for para- than for ortho-H₂O, the overall level of agreement with the RADEX code was reassuring.

4.2. Steady-state versus statistical-equilibrium

Fig. 6 Upper panels: comparison between profiles obtained when the computation of the water line emission is internal to the shock code (steady-state approximation: blue curves) or external to the shock code (statistical equilibrium approximation: red curves); the temperature of the neutral fluid (continuous line) and the fractional abundance of water (dashed lines) are shown for our two reference, C-type (left-hand panel), and CJ-type (right-hand panel) shock models. In both cases, the shock parameters are displayed in the corresponding panels. Lower panels: integrated line intensities for all para-water transitions; steady-state approximation to the computation of the water emission (dark blue open circles); same shock model, but with the statistical-equilibrium approximation to the computation of the water emission (large light-blue open circles); present shock model, with the statistical-equilibrium assumption (small red filled circles).

Fig. 7 Evolution of the computed local optical depths of two transitions as a function of the flow time of the neutral fluid in our reference models, the stationary C-type (left panel) and non-stationary CJ-type (right panel). The o-H2O 110–101 (at 557 GHz) and p-H2O 111–000 (at 1113 GHz) lines are considered, and their optical depths are shown in blue and red, respectively. The corresponding total optical depths are indicated in the same colour code. The neutral temperature profile shown in black. The shock parameters are specified in each panel.

In this section, we compare the results of the present calculations with those obtained using the approach of FP10. We present the intensities of all ortho- and para-transitions, integrated over the whole extent of the shock wave. We compared one stationary C-shock model, with input parameters n_H    =    2    ×    10⁴ cm^-3, v_s    =    20 km s^-1, and b = 1, and one non-stationary CJ-shock model, with input parameters n_H    =  10⁴ cm^-3, v_s    =    28 km s^-1, b = 1.5, and age 500 yrs, which fits both the H₂ and SiO observations (see Fig. 5, panel b). The principal difference between the two versions of the shock code is that FP10 included the population densities of the levels of water as variables of the model, and the differential equations describing the evolution of the level populations were solved simultaneously with those for the other dynamical and chemical variables, in the steady-state approximation. Thus, the water line emissivities were computed within the shock code, at each point of the shock profile, and could be integrated internally. In the present approach, the level population densities (in statistical-equilibrium) and the line intensities are calculated in a separate LVG module, to which the shock profiles are provided as input data.

In addition to the different approximations used when calculating the population densities of the water levels (steady-state² and statistical-equilibrium³), there are two further, potentially significant differences between these two sets of calculations. In the steady-state case, the radiative cooling function that is calculated at each point includes the contribution of the water emission, determined from the level population densities and the line escape probabilities. This contribution has an effect on the temperature profile in regions where H₂O is a significant coolant. In the statistical-equilibirum assumption, the contribution of water emission to the cooling of the medium is approximated by means of the tables provided by Neufeld & Kaufman (1993). In addition, the effect of adsorption of gas-phase species on the radius of the grains is treated differently in the two versions of the shock code; this affects the gas–grain interactions in the tail of the cooling flow, where adsorption takes place. As a consequence, the coupling of the gas to the grain temperature, as the post-shock region is approached, occurs differently in the two versions of the shock code.

In both the C-type and CJ-type cases, we have three sets of LVG calculations. We denote by “FP10” the results of the steady-state (s-s) shock code. Additionally, the dynamical and chemical profiles generated by either version of the shock code may be combined with a separate, statistical-equilibrium (s-e) treatment of the water level populations. The upper panels of Fig. 6 compares the neutral temperature and water fractional abundance profiles from the two versions of the shock code. In the case of the C-type model, differences are noticeable only in the cooling flow, where water is a significant coolant (see top right panel of Fig. 2). Nonetheless, the differences in the thermal profiles are modest, and the water fractional abundances in the two versions remain indistinguishable. In the case of the CJ-type model, even a small displacement of the location of the J-discontinuity can modify the maximum neutral temperature, which in turn has an effect on the fractional abundance profiles: see the r.h.s. upper panel of Fig. 6.

The lower panels of Fig. 6 show the integrated intensities of the total of 156 transitions of para-H₂O that result from the three different calculations. In the case of the C-shock, the similarity of the results is striking and extends to very weak lines. The time-scale for the level populations to approach equilibrium is given by  [n(H₂)q(J + 1 → J)]^-1 (where q(J + 1 → J) is the de-excitation rate coefficient for collisions with molecular H₂) in the limit of high densities; it is even shorter in the limit of low densities, when radiative decay dominates collisional de-excitation. To validate the assumption of s-e, one must compare this population transfer time-scale with a dynamical time-scale, such as the flow time of the neutral particles through the shock region, that corresponds to the width of the neutral temperature profile shown in the upper panels. Adopting the representative values of 10⁴ cm^-3 for the density and 10^-10 cm^-3 for the collisional de-excitation rate coefficient results in a population transfer time-scale of 0.03 yr, much less than the flow time of the neutral particles through the C-type shock considered. In the case of the reference CJ-type shock model, this condition is verified for most of the transitions, as the overall agreement between light and dark blue open circles in Fig. 6 confirms. For the brightest transitions, our independent LVG calculation is consistent with that denoted “FP10”. In this case, the discrepancies with the red filled circles are larger, owing to the differences in the water abundance profiles, which are attributable in turn to the temperature differences arising from: the small displacement of J-discontinuity; the treatment of the cooling due to water; and the strength of the grain–gas interaction as the post-shock region is approached. We note that the pattern created by each set of points is similar for the brightest, and hence observable, transitions, and that the pattern is shifted vertically from one set of results to the other. Such a shift is indistinguishable from that arising from the uncertainty in the filling factor, when comparisons are made with observations.

We conclude from the above analysis that the use of the s-e approximation is justified in the present context of the study of water emission.

4.3. Water emission predictions for BHR71

As a preamble to this section – aimed at providing predictions of water line intensities for use in the interpretation of observations – we study the evolution of the optical depth of the fundamental transitions of both o-H₂O and p-H₂O in each of our reference models. This evolution, together with the total optical depths, integrated through the shock wave, are shown in Fig. 7. Owing to the one-dimensional nature of our approach, we probably overestimate the optical depths and find that the transitions are optically thick in most parts of our shock models; this is especially the case of the C-type shock, as the water accumulates in the large post-shock region. The example of the CJ-type shock shows the spread of values of the optical depths in the emitting region – low when water is dissociated in the high-temperature regime, followed by higher values when water reforms. This figure underlines the requirement for an adequate treatment of the radiative transfer.

As already mentioned, H₂, SiO and H₂O are believed to have coincident regions of emission, at least to the resolution of currently available observations (e.g. Nisini et al. 2010). For completeness, we provide predictions of water line intensities for all the models that were found to fit the H₂ pure rotational excitation diagram acceptably at position of the SiO knot. Given the relative uniformity of the shock conditions in the inner parts of the outflow, revealed by the similarity of the H₂ excitation diagrams at the different positions considered, our models should be applicable to water observations at any position. Stationary, C-type shock models are excluded, for the reasons given in Sect. 3.3.

The results are given in Tables B.1 and B.2 for o-H₂O and p-H₂O, respectively. The transitions that appear in the Tables are the ones that will be the target of Herschel observations. The first three transitions in each table will be observed by HIFI, and the remainder (nine transitions of o-H₂O and three of p-H₂O) will be observed by PACS, in the framework of the WISH key-program (van Dishoeck et al. 2011). The calculations comprise the entire set of transitions that are listed on the LAMDA website, and the complete results are available upon request.

The first general comment concerns the dependence of the water emission line intensities on the pre-shock density. We find that the integrated intensities calculated for a pre-shock density of 10⁵ cm^-3 are systematically around one order of magnitude larger than those computed in models with pre-shock densities of 10⁴ cm^-3, regardless of the values of the other parameters, such as the shock velocity or magnetic-field strength. In addition, and with all other parameters remaining equal, we find that the computed integrated intensities vary by factors of 2 − 3, and sometimes more, over the range of shock ages considered; this variation may be sufficient, in certain circumstances, to enable the observations to constrain the age of the shock wave.

Fig. 8 Integrated intensities of the strongest transitions of o- and p-H2O to be targeted by HIFI and PACS, for the eight shock models presented in panels b) and c) of Fig. 5. The o-H2O and p-H2O transitions are indicated by circles and squares, respectively. The shock model parameters are abbreviated such that “n4b1.50v28t500”, for example, corresponds to a pre-shock density of 104 cm-3, a magnetic-field parameter of 1.5, a shock velocity of 28 km s-1, and an age of 500 years.

Finally, we consider the computed data pertaining to the eight shock models that provide the best simultaneous fits for H₂ and SiO; these models are shown in Fig. 5, in panels (b) and (c). The results of the calculations of the integrated intensities of the brightest o- and p-H₂O transitions, likely to be targetted by HIFI and PACS, are shown in Fig. 8. The filling factor is assumed to be the same and equal to 1 for all the transitions and too small to discriminate between these models. Nevertheless, forthcoming Herschel observations should enable us to enhance our understanding of the chemistry of water in interstellar shock waves.

Water is believed to be the main interstellar reservoir of oxygen, which is the most abundant element in the universe after H and He and plays a major role in interstellar chemistry. Given its rotational energy level structure, water is an important gas coolant and a convenient tracer of shocked regions. Forthcoming observations, through their comparison with the predictions of shock models, should enable us to investigate both the chemical and the physical processes at work in bipolar outflows.

5. Concluding remarks

We have presented a detailed analysis of observations at a few positions in the the bipolar outflow BHR71, assuming that the emission from these regions arises in shock waves. In two cases, namely the Herbig-Haro objects HH320 and HH321, we are able to make use of both the recent Spitzer observations of pure rotational transitions and previous, ground-based observations of rovibrational transitions of molecular hydrogen. In order to compare with the predictions of our shock models, we used excitation diagrams. The physical conditions were found to be similar in both cases, with a pre-shock density of 10⁴ cm^-3 and moderate shock velocities, in the range 15 − 35 km s^-1. The only shock models that fit these data satisfactorily are non-stationary, with ages that depend on the magnetic-field strength and do not exceed 2000 years. Three other positions in the upper and lower lobes of the outflow were studied, using only the pure rotational transitions of H₂ observed by Spitzer. Our analysis suggests that the excitation conditions at these positions are similar to those in the HH objects and can only be reproduced by non-stationary models.

We present a new map of the SiO(5 − 4) transition in the inner part of the bipolar outflow BHR71, together with new observations of the SiO(6–5) and SiO(8–7) lines, at one position that is also covered by the Spitzer observations and at another position that lies outside of the Spitzer maps. At the first position, referred to as “the SiO knot”, we make use of the combination of data (pure rotational transitions of both H₂ and SiO) in order to constrain the parameters of our shock models more tightly than what is possible with the rotational transitions of H₂ alone. The modelling based on H₂ leads to the conclusion as for the other positions in the outflow, namely non-stationary shocks provide the best fits. In particular, shocks with a pre-shock density of 10⁵ cm^-3, low velocities (10–20 km s^-1) and ages (less than 1000 years) were found to be good candidates. Including the SiO emission places constraints on the initial repartition of silicon-bearing material but is of little help in discriminating between the models which fit the H₂ rotational excitation diagram.

We have validated our approach to calculating the water line intensities, which uses the results of shock model as input to a separate radiative transfer code, based on the LVG approximation. To this end, we provide a full comparison between the statistical-equilibrium and the steady-state approximations, based on reference shock models of the stationary, C-type, and non-stationary, CJ-type. We make predictions of water emission line intensities for the SiO knot in the BHR71 outflow. As we expect emission by water to coincide with that by SiO and H₂, we provide the results of the LVG calculations of the water line intensities for models that fit satisfactorily both the SiO and H₂ emission⁴. Given that the excitation conditions are relatively uniform in the internal parts of the outflow, we consider that our calculations can be used to compare with observations of water throughout this region.

We shall discuss in a future publication ground-based observations of CO transitions, undertaken with the APEX facility, that could bring further insights into the physical conditions in the BHR71 outflow. Complementary observations of sulphur-bearing molecules and other shock tracers, such as HNCO, are planned. One of the aims of these studies is to investigate the suitability of “chemically active outflows”, such as BHR71, for observations with the new generation of instruments, such as ALMA.

Online material

Appendix A: Model constraints for the HH objects and the SiO knot

Appendix B: Water emission predictions for the SiO knot

Acknowledgments

B.P. acknowledges the support of the Deutsche Forschungsgemeinschaft through the Emmy Noether grant PA1692/1–1. We also wish to thank Malcolm Walmsley and an anonymous referee for useful comments that contributed to improving the completeness and clarity of our paper.

References

All Tables

All Figures

Fig. 1 The internal parts of the BHR71 outflow as observed by Spitzer in the H2 0–0 S(5) transition (Neufeld et al. 2009); the colour scale unit is erg cm-2 s-1 sr-1. The contours correspond to the levels of 3, 9, 15, 21, 27, and 33σ. The pink hexagons give the positions of the IRS1 and IRS2 sources as determined by Bourke et al. (2001). The positions of the Herbig-Haro objects HH320AB and HH321AB are plotted in black as well as those of the SiO peak (referring to the point of maximum SiO (5–4) emission, see Sect. 3), and of the SiO knot, which was chosen to perform our simultaneous fit of H2 and SiO emission, and also to make predictions of water emission. The other emission knots that were analysed are plotted in white (knots 5, 6, and 8, following the labelling of Giannini et al. 2004). Around the SiO-peak, four radii are displayed, corresponding to 5, 9, 12, and 14.4′′. The first represents the size of our region of analysis. The latter three correspond to the beam of the APEX telescope for the SiO (8−7), (6−5) and (5−4) transitions, respectively; see text.

In the text

Fig. 2 The H2 excitation diagrams observed for each position indicated in Fig. 1, in black (pure rotational lines), and grey (rovibrational lines, when available). The data were corrected for extinction, using AV = 2 from Giannini et al. (2004). The inferred rotational temperatures are indicated in black, and the rovibrational temperatures, in grey. The best-fitting model results are also shown in each case, in blue and/or red symbols. The shock parameters are given in each panel in the corresponding colour.

In the text

Fig. 3 Upper internal part of the BHR71 outflow. The background map is the same as in Fig. 1. The SiO (5–4) contours, which are overlaid in a dark blue colour, correspond to levels of 6 to 30σ, in steps of 3σ, for an integration between  − 10 and 60 km s-1. The half-maximum, 15σ contour is displayed in red. The associated beam size is indicated in the upper left corner. Our SiO (5 − 4) coverage is actually wider than the field displayed in this figure.

In the text

	Fig. 4 The spectra obtained in the three observed SiO transitions at the positions of the SiO peak and knot, of respective offsets ( − 15′′, 87′′) and (0′′, 77′′). Only the SiO knot is considered in our modelling. For each spectrum, we adopted the SiO (5–4) spatial resolution. The ambient cloud velocity, –4.5 km s-1 is indicated by the vertical dashed lines.
In the text

Fig. 5 Examples of fits of the integrated intensities of our three transitions of SiO. On the y-axis is the integrated intensity (K km s-1), and on the x-axis is the rotational quantum number of the upper level of the transition. The observations are the black squares linked by the thick black lines, on all three panels. The 3σ errorbars are shown in the figure. The model results are displayed in colour (blue and red), with the indication of the specific shock parameters on each panel. a) C-type shock models: best fit shock parameters nH = 104 cm-3, no SiO initially in the grain mantles; b) CJ-type shock models: best fit shock parameters nH = 104 cm-3 and vs = 28 km s-1, with 1% of elemental silicon initially in the grain mantles, as SiO; c) CJ-type shock models: best fit shock parameters nH = 104 cm-3 and vs = 30 km s-1, with 10% of the silicon initially in the grain mantles as SiO.

In the text

Fig. 6 Upper panels: comparison between profiles obtained when the computation of the water line emission is internal to the shock code (steady-state approximation: blue curves) or external to the shock code (statistical equilibrium approximation: red curves); the temperature of the neutral fluid (continuous line) and the fractional abundance of water (dashed lines) are shown for our two reference, C-type (left-hand panel), and CJ-type (right-hand panel) shock models. In both cases, the shock parameters are displayed in the corresponding panels. Lower panels: integrated line intensities for all para-water transitions; steady-state approximation to the computation of the water emission (dark blue open circles); same shock model, but with the statistical-equilibrium approximation to the computation of the water emission (large light-blue open circles); present shock model, with the statistical-equilibrium assumption (small red filled circles).

In the text

Fig. 7 Evolution of the computed local optical depths of two transitions as a function of the flow time of the neutral fluid in our reference models, the stationary C-type (left panel) and non-stationary CJ-type (right panel). The o-H2O 110–101 (at 557 GHz) and p-H2O 111–000 (at 1113 GHz) lines are considered, and their optical depths are shown in blue and red, respectively. The corresponding total optical depths are indicated in the same colour code. The neutral temperature profile shown in black. The shock parameters are specified in each panel.

In the text

Fig. 8 Integrated intensities of the strongest transitions of o- and p-H2O to be targeted by HIFI and PACS, for the eight shock models presented in panels b) and c) of Fig. 5. The o-H2O and p-H2O transitions are indicated by circles and squares, respectively. The shock model parameters are abbreviated such that “n4b1.50v28t500”, for example, corresponds to a pre-shock density of 104 cm-3, a magnetic-field parameter of 1.5, a shock velocity of 28 km s-1, and an age of 500 years.

In the text