A&A 482, 809-829 (2008)
DOI: 10.1051/0004-6361:20078900

SiO line emission from C-type shock waves: interstellar jets and outflows

A. Gusdorf1,2 - S. Cabrit3 - D. R. Flower1 - G. Pineau des Forêts2,3,4


1 - Physics Department, The University, Durham DH1 3LE, UK
2 - Institut d'Astrophysique Spatiale (IAS), Bâtiment 121, 91405 Orsay, France
3 - LERMA (UMR 8112 du CNRS), Observatoire de Paris, 61 avenue de l'Observatoire, 75014 Paris, France
4 - Université Paris-Sud 11 and CNRS (UMR 8617), France

Received 23 October 2007 / Accepted 12 February 2008

Abstract
We study the production of SiO in the gas phase of molecular outflows, through the sputtering of Si-bearing material in refractory grain cores, which are taken to be olivine. We calculate also the rotational line spectrum of the SiO. The sputtering is driven by neutral particle impact on charged grains, in steady-state C-type shock waves, at the speed of ambipolar diffusion. The emission of the SiO molecule is calculated by means of an LVG code. A grid of models, with shock speeds in the range $20 < v_{\rm s} < 50$ km s-1 and preshock gas densities $10^4 < n_{\rm H} < 10^6$ cm-3, has been generated. We compare our results with those of an earlier study (Schilke et al. 1997). Improvements in the treatment of the coupling between the charged grains and the neutral fluid lead to narrower shock waves and lower fractions of Si ($\la$10%) being released into the gas phase. Erosion of grain cores is significant ($\ga$1%) only for C-type shock speeds $v_{\rm s} > 25$ km s-1, given the adopted properties of olivine. More realistic assumptions concerning the initial fractional abundance of O2 lead to SiO formation being delayed, so that it occurs in the cool, dense postshock flow. Good agreement is obtained with recent observations of SiO line intensities in the L1157 and L1448 molecular outflows. The inferred temperature, opacity, and SiO column density in the emission region differ significantly from those estimated by means of LVG ``slab'' models. The fractional abundance of SiO is deduced and found to be in the range 4 $\times $ $10^{-8} \la n({\rm SiO})/n_{\rm H} \la 3$ $\times $ 10-7. Observed line profiles are wider than predicted and imply multiple, unresolved shock regions within the beam.

Key words: astrochemistry - atomic processes - magnetohydrodynamics (MHD) - molecular processes - radiative transfer - shock waves

   
1 Introduction

Unlike CO, which is observed extensively in the interstellar medium of our own and other galaxies, its homologue SiO is observed principally in outflows associated with regions of star formation. This striking difference in behaviour is a consequence of the lower elemental abundance and the more complete depletion of silicon from the gas phase during grain formation. Both carbon and silicon form refractory materials - graphite and silicates - of which the cores of interstellar grains are believed to be composed; but the much higher elemental abundance of carbon, $n_{\rm C}/n_{\rm H} = 3.55$ $\times $ 10-4, compared with silicon, $n_{\rm Si}/n_{\rm H} = 3.37$ $\times $ 10-5 (Anders & Grevesse 1989), leads to some of the carbon remaining in the gas phase.

The SiO molecule was first detected in the Galactic centre (Sgr B2) by Wilson et al. (1971) and subsequently in Ori A by Dickinson (1972). More recent observations of SiO in jets (see, for example, Bachiller et al. 1991; Martín-Pintado et al. 1992; Codella et al. 1999; Nisini et al. 2007) imply that, in these objects, at least some of the silicon has been restored to the gas phase; this can be achieved through sputtering of the grain material, probably in shock waves, which are features of the outflows. It has been known since the study by Draine et al. (1983) that grain ice-mantles can be eroded in C-type shock waves, owing to impact of neutral particles on charged grains at the ion-neutral drift speed, which is the speed of ambipolar diffusion. More recent work (Flower & Pineau des Forêts 1995; Flower et al. 1996; Field et al. 1997; May et al. 2000) has shown that this process might result also in the partial erosion of the refractory grain cores. The simulations undertaken by May et al. of the sputtering of various silicates (forsterite, fayalite and olivine) by neutral atoms showed that C-shock speeds in excess of 30 km s-1 are necessary to erode a significant fraction (more than a few per cent) of these materials. On the other hand, it has since been recognized (Le Bourlot et al. 2002; Ciolek et al. 2004) that the speeds at which C-type shock waves can propagate are limited, both by the collisional dissociation of molecular hydrogen, which leads to a thermal runaway and the formation of a J-type shock wave, and by the ion magnetosonic speed, whose value is constrained by the inertia of the charged grains. Consequently, the efficacy of the erosion of silicates in C-type shock waves is restricted by the maximum C-shock speed, which depends on the preshock density and transverse magnetic field strength, and on the fraction of charged grains (Flower & Pineau des Forêts 2003).

In a previous study of SiO production in the interstellar medium, Schilke et al. (1997, henceforth Sch97) considered in detail the erosion of silicon from grain cores and from their mantles by C-type shock waves and the resulting SiO emission spectrum; this study remains the only one of its kind that has been published to date. In the intervening decade, there has been progress in our understanding of the sputtering process (May et al. 2000), the gas-phase chemistry of silicon (Le Picard et al. 2001), and the grain dynamics (Flower & Pineau des Forêts 2003), and so it seems timely to revisit this topic. Recent observations of SiO in jets (Nisini et al. 2007) extend to higher rotational levels than previously and provide a further motivation for such an study. Accordingly, we have reconsidered the chemistry of silicon in steady-state, planar MHD shock waves, with a view to providing a grid of models which may serve in the analysis of current and future observations of rotational transitions of SiO in outflows. These models yield additional results relating to rovibrational transitions of H2 and rotational lines of CO, as well as forbidden lines of atoms and atomic ions, including [Fe II] (cf. Giannini et al. 2006).

Grain-grain collisions and the sputtering of silicon in oblique C-type shock waves, in which the preshock magnetic field direction is inclined to the plane of the shock front, have been considered by Caselli et al. (1997). However, such simulations have yet to incorporate the gas-phase chemistry and the radiative cooling processes, in order to enable quantitative analyses of the spectral line observations to be made.

   
2 The model

We summarize below those developments, in both the MHD code and the data used and produced by the code, which are relevant to the modelling of the intensities and the profiles of the rotational lines of SiO. We take as our baseline the study by Sch97. The reader is referred to the more recent papers cited below for details of the modifications.

2.1 The dynamics of charged grains

The main revision of the code itself concerns the treatment of the charge and of the dynamical effects of the grains. The charge distribution of both the grains and the ``very small grains'' (VSG), represented by polycyclic aromatic hydrocarbons in our model, are calculated, assuming that collisions with gas-phase particles (electrons, ions and neutrals) determine the charge; see Flower & Pineau des Forêts (2003). As was mentioned in the Introduction, allowance for the mass density of the (mainly negatively) charged grains can reduce significantly the magnetosonic speed in the charged fluid

\begin{displaymath}c_{\rm m}^2 = \frac {5 k_{\rm B}(T_+ + T_-)} {3(\mu _+ + \mu _-)} + \frac {B^2}{4\pi (\rho _+ + \rho _-)} \end{displaymath}

where $T_{\pm }$ are the temperatures, $\mu _{\pm }$ are the mean masses and $\rho _{\pm }$ are the mass densities of the positively and negatively charged fluids; B is the magnetic field strength in the preshock gas. The magnetosonic speed is the maximum speed at which a C-type shock can propagate in the medium. Furthermore, the momentum transfer between the charged grains and the neutral fluid affects the ion-neutral drift speed and has consequences for the degree of sputtering of the grains within a C-type shock wave.

2.2 Radiative cooling by H2

The thermal balance of the medium, particularly the radiative cooling by H2, is treated more exactly in the current model. Rovibrational excitation of H2, principally by H, H2 itself, and He, followed by radiative decay, is the principal mechanism of cooling the shock-heated gas. Following Le Bourlot et al. (2002), the populations of the rovibrational levels of H2 are calculated in parallel with the hydrodynamical and chemical rate equations, yielding the most accurate determination of the rate of cooling by H2 that is achievable within the framework of a time-independent (steady-state) model of the shock structure. The rate coefficients for the collisional excitation by H of rovibrational transitions of H2 are from the recent work of Wrathmall et al. (2007). Collisional dissociation and ionization of H2, as well as ionization of H, are taken into account. Collisional dissociation of H2 is a particularly important process, as the removal of H2 can lead to a thermal runaway. The associated increase in the kinetic temperature, $T_{\rm n}$, of the neutral fluid, and hence of the adiabatic sound speed

\begin{displaymath}c_{\rm s}^2 = \frac {5 k_{\rm B}T_{\rm n}}{3\mu _{\rm n}} \end{displaymath}

where $\mu _{\rm n}$ is the mean mass of the neutral fluid, can give rise to a sonic point in the flow and hence to a J-type discontinuity. This phenomenon imposes an additional constraint on the maximum speed of a C-type shock wave in a molecular medium.

In the compressed gas of the postshock region where most of SiO forms, cooling by 12CO, 13CO, and H2O starts to dominate that by H2. In order to predict accurately the emitted SiO emission spectrum, an LVG treatment of the cooling by these species has been introduced, following the procedures of Neufeld & Kaufman (1993). A ``thermal gradient'' $c_{\rm s}/z^{\prime }$, where $z^{\prime } = z + 1.0$ $\times $ 1013 cm and z is the independent integration variable, is added quadratically to the macroscopic velocity gradient, in order to simulate photon escape through thermal line broadening.

2.3 The sputtering of grains

The sputtering probabilities computed by May et al. (2000) for olivine (MgFeSiO4) are used to determine the rate of erosion of Si from (charged) silicate grains by neutral particles. We include also the sputtering of carbonaceous (amorphous carbon) grains, using the sputtering yields of Field et al. (1997). Impacts (at the ion-neutral drift speed) of abundant heavy neutral species are the most effective in eroding the grain cores. For collisions with CO, for example, we adopted the sputtering probabilities calculated for impacts of Si atoms, which have the same mass as CO and hence the same impact energy. As May et al. (2000) have shown, similar yields of Si are obtained from the three types of silicate: fayalite, Fe2SiO4; forsterite, Mg2SiO4; and olivine, MgFeSiO4.

The grain mantles are eroded first, through impacts with the most abundant neutral species, H, H2 and He, at ion-neutral drift speeds which are below the threshold for erosion of the cores (Draine et al. 1983; Flower & Pineau des Forêts 1994). The initial composition of the gas is calculated in chemical equilibrium, whilst that of the grain mantles, and the elemental depletion into the grain cores, is based on observations (cf. Flower & Pineau des Forêts 2003, Table 1). We incorporated a representative polycyclic aromatic hydrocarbon (PAH), with a fractional abundance $n_{\rm PAH}/n_{\rm H} = 10^{-6}$, an upper limit which is consistent with estimates of the fraction of elemental carbon likely to be present in the form of very small grains (Li & Draine 2001; Weingartner & Draine 2001). The state of charge of this species was computed following Flower & Pineau des Forêts (2003), who showed that increasing the fractional abundance of the PAH results in higher values of the magnetosonic speed in the preshock gas, thereby enabling C-type shock waves to propagate at higher speeds. In their turn, the higher speed shocks erode the silicate grains more effectively. However, the adopted PAH abundance has little effect on the internal structure of the shock wave, as the grains become negatively charged in the shock wave and dominate grain-neutral momentum transfer.

The total number density of the grains was computed

The thickness of the grain mantles was determined from their molecular composition (Flower & Pineau des Forêts 2003, Table 2), assuming a mean number of 5 $\times $ 104 molecular binding sites per layer of the mantle and a thickness of 2 $\times $ 10-4 $\mu $m for each layer, independent of the size of the grain core; there are no Si-containing species in the mantles. This procedure yields an initial mantle thickness of 0.015 $\mu $m on the grain cores, whose mean radius is $a_{\rm g} = 0.02$ $\mu $m. However, the erosion of the grain mantles occurs sufficiently rapidly, as the ion and neutral flow velocities begin to decouple in the shock wave, that the existence of thick ice-mantles on the grains in the preshock gas has little effect on the shock dynamics (see Fig. 6 of Guillet et al. 2007).

   
2.4 Chemistry

The chemical reaction network comprises over 900 reactions connecting the abundances of more than 100 species, in both the gas and the solid phases. The complete list of species and reactions is available from http://massey.dur.ac.uk/drf/outflows_test/species_chemistry_shock/. In the context of the present study, we emphasize the gas-phase chemistry of silicon and the formation of SiO, in particular. The total rate of cosmic ray ionization of hydrogen, $\zeta$, was taken to be 5 $\times $ 10-17 s-1.

Two reactions are important for the formation of SiO from Si, which is released into the gas phase by the sputtering of the grains, namely

 \begin{displaymath}% \begin{array}{ccc} \rm {Si} + \rm {O}_{2} & \longrightarrow & \rm {SiO} + \rm {O}\\ \end{array}\end{displaymath} (1)


 \begin{displaymath}% \begin{array}{ccc} \rm {Si} + \rm {OH} & \longrightarrow & \rm {SiO} + \rm {H} \\ \end{array}\end{displaymath} (2)

for which the following rate coefficients (cm3 s-1)

 \begin{displaymath}% \begin{array}{ccc} k_1 & = & 2.7\times 10^{-10} \exp({-111/T}) \\ \end{array}\end{displaymath} (3)


 \begin{displaymath}% \begin{array}{ccc} k_2 & = & 1.0\times 10^{-10} \exp({-111/T}) \\ \end{array}\end{displaymath} (4)

were adopted by Sch97 from the compilation of Langer & Glassgold (1990). The exponential factor in Eq. (3) derives from the argument (Graff 1989) that the reactions proceed only with the fine-structure states of Si (3p2 3PJ) with J > 0, of which J = 1, which lies 111 K above the J = 0 ground state, is the more significantly populated at low temperatures. More recently, the rate coefficient for reaction (1) has been measured at low temperatures ( $15 \le T \le 300$ K) by Le Picard et al. (2001) and found to be given by

 \begin{displaymath}% k_1 = 1.72\times 10^{-10} (T/300)^{-0.53} \exp(-17/T). \end{displaymath} (5)

We have adopted (5) and the same expression for k2. Evidently, the differences between the present and previous values of these rate coefficients are most significant for temperatures $T \la 100$ K, i.e. in the cooling flow of the shocked gas.

The abundance of SiO is limited by its conversion to SiO2 in the reaction with OH

 \begin{displaymath}% \begin{array}{ccc} \rm {SiO} + \rm {OH} & \longrightarrow & \rm {SiO_2} + \rm {H} \\ \end{array}\end{displaymath} (6)

whose rate coefficient remains subject to considerable uncertainty. We adopt the same expression as Sch97, viz.

 \begin{displaymath}% k_6 = 1.0\times 10^{-11} (T/300)^{-0.7} \end{displaymath} (7)

in units of cm3 s-1. However, we note that Zachariah & Tsang (1995) calculated a barrier of 433 K to reaction (6), and a rate coefficient

\begin{eqnarray*}k_6 = 2.5\times 10^{-12} (T/300)^{0.78}\exp(-613/T); \end{eqnarray*}


see the discussion of Le Picard et al. (2001). At T = 300 K, the latter rate coefficient is 30 times smaller than the former. In the ambient (preshock) and the postshock gas, where $T \approx 10$ K, the existence of an activation energy of several hundred kelvin would prevent the oxidation of SiO in reaction (6) from occurring. The rate coefficient for reaction (6) in the UMIST data base (Le Teuff et al. 2000) is

\begin{eqnarray*}k_6 = 2.0\times 10^{-12} \end{eqnarray*}


in cm3 s-1, which is 50 times smaller than (7) at T = 10 K. Fortunately, the conversion of SiO into SiO2 occurs in a region which is too cold and optically thick to contribute to the SiO line intensities, and so the uncertainty in the rate coefficient for reaction (6) is not significant in the present context.

Re-adsorption of molecules on to grains (``freeze-out'') in the postshock gas is included, as in Sch97. The effects of freeze-out on the abundance of SiO and its spectrum are considered below.

   
2.5 Line transfer

The physical and chemical profiles which derive from the shock model summarized above are used in a large velocity gradient (LVG) calculation of the intensities of the emission lines of SiO and of CO. Our implementation of this technique is described in Appendix A, where we note that an expression for the escape probability which differs from that of Sch97 has been adopted.

   
3 Results

   
3.1 Comparison with the calculations of Schilke et al. (1997)

First, we compare the computed shock structure with that of Sch97, for a reference C-type shock model, in which the preshock density $n_{\rm H} = 10^{5}$ cm-3 and the magnetic field strength B = 200 $\mu $G, and the shock velocity $v_{\rm s} = 30$ km s-1. The most striking difference between the current and the previous model is that the width of the shock wave decreases by a factor of approximately 4, to 5 $\times $ 1015 cm, from 2 $\times $ 1016 cm in the study of Sch97; see Fig. 1a. This difference is attributable to the more accurate treatment of the coupling between the neutral fluid and the charged grains in the current model and is an indication of the significance of the inertia of the (negatively) charged grains in dark clouds, in which the degree of ionization is low. With the narrower shock wave is associated a higher maximum temperature of the neutral fluid, as there is less time for the initial energy flux, $\rho_{\rm n}v_{\rm s}^3/2$, associated with the bulk flow, to be converted into internal energy of the H2 molecules or to be radiated away. Thus, $T_{\rm n} \approx 4000$ K here, compared with $T_{\rm n} \approx 2000$ K in the model of Sch97.


  \begin{figure} \par\includegraphics[width=6.5cm,clip]{Fig.1.eps}\end{figure} Figure 1: a) Temperature of the neutral fluid and velocity profiles of the neutral and charged fluids, predicted by the present C-type shock model. The shock parameters are $n_{\rm H} = 10^{5}$ cm-3 and B = 200 $\mu $G in the preshock gas, and $v_{\rm s} = 30$ km s-1, $\zeta = 5$ $\times $ 10-17 s-1. The fractional gas-phase abundances of selected Si-and O-bearing species are plotted in panels  b) and  c); cf. Sch97, Fig. 2. The broken curves in panel  b) are the corresponding results obtained when re-adsorption on to grains is neglected. The discontinuous curves in panel  c) show the effects of assuming the initial abundance of O2 ice to be negligible, i.e. the second of the two scenarios described in Sect. 3.1.
Open with DEXTER

There are chemical differences between the models also, which are related in part to the changes in the shock structure; these differences may be summarized as follows.

(i)
The fraction of the Si in the grain cores which is released into the gas phase by sputtering is approximately ten times smaller in the current model than in the model of Sch97. This reduction is attributable partly to the sputtering yields, which have higher thresholds and are smaller for olivine (MgFeSiO4) than for the amorphous silica (SiO2) considered by Sch97; but the main reason for the decrease is the enhanced coupling between the neutral fluid and the charged grains, which reduces the shock width and hence the time available to erode the grains. On the other hand, the magnitude of the ion-neutral drift speed is similar in both calculations. As a consequence of the reduction in the shock width, the integrated SiO line intensities predicted by the current model are smaller, in general, than calculated by Sch97; see Sect. 4.

(ii)
The displacement of the maximum fractional abundance of SiO, which forms in the gas-phase reactions (1) and (2), from that of Si, which is eroded from the grains, is a more significant fraction of the shock width in the current model; cf. Fig. 1b. The initial fractional abundance of O2 in the preshock medium is lower here, by a factor of approximately 10, than in the model of Sch97, delaying the initial formation of SiO. The O2 is assumed to be predominantly in the form of ice, which is sputtered rapidly from the grains in the early stages of development of the shock wave, as may be seen from the two orders of magnitude increase in the fractional abundance of gas-phase O2, apparent in Fig. 1c. The fractional abundances of O2 and OH decrease subsequently, at high kinetic temperatures, owing to their dissociation by H in the chemical reactions O2(H, O)OH and OH(H, O)H2. The former reaction, which is endothermic by over 8000 K, proves to be less effective in destroying O2 over the smaller width of the current shock model (see Fig. 1c) than was the case in the calculations of Sch97. On the other hand, the lower energy threshold of 17 K in reaction (5) allows oxidation reactions to proceed further into the postshock region, compared with Sch97, whose adopted threshold was 111 K. As a consequence, conversion of Si into SiO is slower initially but more complete eventually than predicted by Sch97.

(iii)
SiO2 is removed more rapidly from the gas phase in the current model. The compression is more rapid than in the model of Sch97, and so the rate of adsorption of molecules to grains (``freeze-out'') is higher. If the oxidation of SiO in the reaction SiO(OH, H)SiO2 has an activation energy of several hundred kelvin (Zachariah & Tsang 1995; see Sect. 2.4), the maximum fractional abundance of SiO2 would be reduced still further.
Also shown in Fig. 1b are the fractional abundances which are obtained neglecting re-adsorption on to the grains. The timescale for freeze-out is sufficiently large that the chemical profiles are modified only in the cold postshock gas, where the flow speed is practically constant and the optical depths in the SiO lines are large. Consequently, whilst the effects re-adsorption on the composition of the postshock gas are dramatic, the freeze-out of SiO has no effect on the predicted line intensities.

In chemical equilibrium, the initial fractional abundance of O2 is approximately 10-5, whereas observations with the Odin satellite (Pagani et al. 2003; Larsson et al. 2007) have placed upper limits of $n({\rm O}_2)/n({\rm H}_2) \la 10^{-7}$. In view of these measurements, we have considered two scenarios, both with an initial gas-phase fractional abundance $n({\rm O}_2)/n_{\rm H} = 1.0$ $\times $ 10-7, as a consequence of the freeze-out of oxygen on to grains, but with differing assumptions regarding its chemical form in the grain mantles.

It turns out that the first scenario is practically equivalent to assuming that the O2 is initially in the gas-phase (see Fig. 1c), as its release from the grain mantles occurs early and rapidly (on a timescale of a few years for the model in Fig. 1) in the shock wave. On the other hand, the second scenario can result in reduced levels of oxidation of Si to SiO in the gas-phase (cf. Fig. 1), depending on the relative importance of reactions 1 and 2 in the oxidation process. In what follows, we present results corresponding principally to the first scenario, with the second being considered mainly in Appendix C.

   
3.2 A grid of models

We have computed a grid of shock models with the following parameters:

In all of these models, we characterized the preshock magnetic field strength by

\begin{eqnarray*}B = bn_{\rm H}^{0.5}, \end{eqnarray*}


where $n_{\rm H}$ is in cm-3 and B is in $\mu $G; the scaling parameter b was taken equal to 1. (The effect on Si sputtering of varying b is discussed in Sect. 3.3.) The maximum shock speed for $n_{\rm H} \ga 10^4$ cm-3 is determined by the collisional dissociation of H2, the main coolant, which leads to a thermal runaway and a J-discontinuity (Le Bourlot et al. 2002; Flower & Pineau des Forêts 2003).

In fact, we computed two grids, one for each of the scenarios concerning the initial distribution of oxygen between O2 and H2O ices, as specified towards the end of the previous Sect. 3.1. We concentrate on the first of these two scenarios, but some additional figures for the second scenario are given in Appendix COnline Material. (Our results are available in digital tabular format on request to the authors.)

Because of the sharply defined sputtering threshold energy of approximately 50 eV, there is negligible sputtering of Si from the olivine (MgFeSiO4) for shock speeds of 20 km s-1 or less. The fractions of the Mg, Si and Fe which are released from the olivine into the gas phase are shown in Fig. 2. Comparing Fig. 2 with the corresponding Fig. 4 of May et al. (2000), whose sputtering yields are used in the present calculations, shows that the fractions of Mg, Si and Fe which are sputtered from olivine have decreased by an order of magnitude. As the same sputtering yields have been used in both studies, this change is attributable to the reduction in the shock width, resulting from the improved treatment of grain-neutral coupling. We note that CO is the principal eroding partner (cf. May et al. 2000).


  \begin{figure} \par\includegraphics[width=7cm,clip]{Fig.2.eps}\end{figure} Figure 2: The fractions of Mg, Si and Fe, initially in the form of olivine (MgFeSiO4), which are released into the gas phase by sputtering within a steady-state C-type shock wave.
Open with DEXTER

Figure 2 shows that the degree of sputtering is, in fact, insensitive to the preshock gas density (cf. Caselli et al. 1997); it depends essentially on the shock speed. Polynomial fits of the sputtered fractions of Fe, Si and Mg, as functions of the shock speed, are given in Appendix B.

In Fig. 3, we display the fractional gas-phase abundances of Si and SiO, as functions of the relevant flow time. Silicon is produced by erosion of the charged grains by collisions, principally with molecules, at the ion-neutral drift speed. Once the drift speed exceeds the sputtering threshold velocity, the erosion of Si occurs rapidly, as Fig. 3 shows. Thus, the flow time which is directly relevant to the release of Si into the gas phase is that of the charged fluid, rather than that of the neutrals, which is the appropriate measure of the total time for formation of SiO. As noted in item (ii) of Sect. 3.1, there is an additional, chemical delay to the conversion of Si into SiO, in reactions (1) and (2), which is apparent in our Fig. 3, owing to the low abundance and partial destruction of O2. The magnitude of this delay depends on the parameters of the model, notably the shock speed, $v_{\rm s}$, and the preshock gas density, $n_{\rm H}$. Conversion is almost instantaneous for $v_{\rm s} \ge 30$ km s-1, $n_{\rm H} = 10^{6}$ cm-3, when OH is formed abundantly at the start of the shock and reaction (2) dominates the oxidation process.


  \begin{figure} \par\includegraphics[width=12cm,clip]{Fig.3.eps}\end{figure} Figure 3: The fractional abundances of Si, released into the gas phase by the sputtering of olivine (MgFeSiO4), and of SiO, which subsequently forms in reactions (1) and (2). In the left-hand panels, the independent variable (abscissa) is the flow time of the charged fluid: see text, Sect. 3.2.
Open with DEXTER


  \begin{figure} \par\includegraphics[width=12cm,clip]{Fig.4.eps}\end{figure} Figure 4: The fractional abundance of SiO, n(SiO)/$n_{\rm H}$, computed for the grid of shock models and plotted as a function of distance ( left) and neutral flow time ( right); $n({\rm SiO})/n_{\rm H}$ is negligible for $v_{\rm s} \protect\la 20$ km s-1. In addition, the temperature of the neutral fluid is plotted. (See also Fig. C.1 of Appendix C.)
Open with DEXTER

Figure 4 shows the variation with shock speed and preshock gas density of the fractional abundance of SiO, computed through the entire shock wave. It is evident from Fig. 4 that the duration of the C-type shock wave, as measured by the temperature profile, is of the order of 104, 103 and 102 yr for preshock gas densities $n_{\rm H} = 10^{4}$, 105 and 106 cm-3, respectively. The peak SiO abundance is reached over similar timescales. It may be seen that the highest fractional abundances of SiO are attained for the lowest preshock density, $n_{\rm H} = 10^{4}$ cm-3. At higher densities, both O2 and OH, which are reactants in (1) and (2), are destroyed by the atomic hydrogen which is produced in the shock wave. Thus, the conversion of Si into SiO becomes incomplete at high density, and the gas-phase SiO abundance depends on $n_{\rm H}$ even though the Si sputtered fraction does not (cf. Fig. 2).

   
3.3 Dependence on the transverse magnetic field strength

The existence of a magnetic field transverse to the direction of propagation is a necessary condition for a C-type shock wave to form, and it is instructive to consider the variation of the structure of the shock wave with the strength of the magnetic field. Energy equipartition arguments, applied to the magnetic and thermal energy densities in the preshock molecular gas, of particle density $n({\rm H}_2) + n({\rm He}) = 0.6n_{\rm H}$ and kinetic temperature T, suggest that $B^2/(8\pi ) \approx n_{\rm H}k_{\rm B}T$ and hence that $B = bn_{\rm H}^{0.5}$, where b is a scaling parameter (cf. Sect. 3.2) such that B is in $\mu $G when $n_{\rm H}$ is in cm-3; this is the proportionality adopted in the grid of models presented in Sect. 3.2. In gas of T = 10 K, equipartition with the thermal energy implies b = 0.18. However, we note that such a low value of b is inconsistent with the existence of a steady-state C-type shock wave when $n_{\rm H} = 10^{5}$ cm-3 and $v_{\rm s} \ge 10$ km s-1, as the corresponding ion magnetosonic speed in the preshock gas (9.7 km s-1) is lower than the shock speed.


  \begin{figure} \par\includegraphics[width=6.5cm,clip]{Fig.5.eps}\end{figure} Figure 5: a) The ion-neutral velocity difference, $\Delta v = \vert v_{\rm i} - v_{\rm n}\vert$, and b) the fractions of Mg, Si and Fe eroded from olivine (MgFeSiO4) grains, computed as functions of the transverse magnetic field strength, B, in the preshock gas; $n_{\rm H} = 10^{5}$ cm-3 and $v_{\rm s} = 30$ km s-1. Note that $B = bn_{\rm H}^{0.5}$, where b is a scaling parameter (cf. Sect. 3.2) such that B is in $\mu $G when $n_{\rm H}$ is in cm-3.
Open with DEXTER

In Fig. 5, we present results as a function of the transverse magnetic field strength, for the parameters of the model of Sch97: $n_{\rm H} = 10^{5}$ cm-3, $v_{\rm s} = 30$ km s-1, and a magnetic field scaling parameter b in the range $0.5 \le b \le 5$. We recall that Sch97 adopted B = 200 $\mu $G, corresponding to b = 0.63. This value of b is consistent with the analysis of Zeeman measurements by Crutcher (1999), who concluded that there was an approximate equipartition of the magnetic and kinetic energy densities in the molecular clouds that he had observed. It may be seen from Fig. 5 that increasing the magnetic field inhibits the release of Si from refractory grain cores in the shock wave, owing to the reduction in the maximum ion-neutral velocity difference, $\Delta v$. In Sect. 4.5, we consider how the magnetic field strength affects the relative intensities of the rotational transitions of SiO.

   
4 SiO rotational emission lines

The observable quantities are the intensities of the rotational transitions of SiO and the velocity-profiles of these emission lines. Having computed the shock structure, we evaluate the line intensities and profiles as described in Appendix A, assuming that the shock is viewed face-on.

   
4.1 Physical conditions in the SiO emission region

Figure 6 illustrates the variation of physical conditions throughout the formation region of the SiO 5-4 rotational line for our reference model with $n_{\rm H} = 10^{5}$ cm-3, $v_{\rm s} = 30$ km s-1, and b = 0.63. It may be seen that the line is optically thin through most of the hot precursor (where the flow speeds of the charged and neutral fluids differ), due to both the low SiO abundance and the large velocity gradient there. Therefore the line intensity is low despite a high kinetic temperature. At the rear of the shock wave, approaching maximum compression, the synthesis of SiO and the steady decrease in velocity gradient eventually raise the optical depth in the line, and the 5-4 intensity peaks, with the line temperature attaining values close to the local kinetic temperature of the neutral fluid, $T_{\rm n}$; that is, the line approaches LTE. The intensity then declines rapidly as the gas cools; the decline occurs in 500 yr for the model shown here. This behaviour is insensitive to the rate of re-adsorption of SiO on to the grains, which occurs over much longer timescales.


  \begin{figure} \par\includegraphics[width=7cm,clip]{Fig.6.eps}\end{figure} Figure 6: a) The temperature of the neutral fluid, $T_{\rm n}$, the brightness temperature, T(5-4), in the j = 5-4 line, and the compression factor, $n_{\rm H}/n_{\rm H}({\rm initial})$; b) the optical depth, $\tau _{5-4}$ in the 5-4 transition and the fractional abundance of SiO, $n({\rm SiO})/n_{\rm H}$, as functions of the flow time of the neutral fluid, $t_{\rm n}$. The model parameters are $n_{\rm H} = 10^{5}$ cm-3, $v_{\rm s} = 30$ km s-1, and b = 0.63.
Open with DEXTER


  \begin{figure} \par\includegraphics[width=9cm,clip]{Fig.7.eps}\end{figure} Figure 7: Physical conditions at the position of the peak in the SiO 5-4 line intensity [ $T_{\rm peak}(5{-}4)$] as functions of the shock speed, $v_{\rm s}$, for all models of the grid: a) neutral temperature, $T_{\rm n}$; b) preshock density, $n_{\rm H}$; c) the LVG parameter, $n({\rm SiO})/({\rm d}v_z/{\rm d}z)$, and the fractional abundance of SiO, $x({\rm SiO}) \equiv n({\rm SiO})/n_{\rm H}$. (See also Fig. C.2 of Appendix C.)
Open with DEXTER


  \begin{figure} \par\includegraphics[width=7cm,clip]{Fig.8.eps}\end{figure} Figure 8: The velocity profiles of transitions from rotational levels $j \rightarrow $ j-1 of SiO, computed for our reference model, viewed face-on. Only those lines detectable from the ground are shown. The model parameters are $n_{\rm H} = 10^{5}$ cm-3, B = 200 $\mu $G, and $v_{\rm s} = 30$ km s-1. The flow speed is in the reference frame of the preshock gas. The neutral temperature profile, $T_{\rm n}$, is shown also, as are indicative values of the flow time of the neutral fluid. (See also Fig. C.3 of Appendix C.)
Open with DEXTER

In order to illustrate the dependence of conditions in the SiO emission region on the shock parameters, we present, in Fig. 7, the important physical quantities, evaluated at the peak of the SiO 5-4 line, for all models in our grid. We have verified that this is equivalent to computing intensity-weighted geometric means of the same quantities over the region where the 5-4 line intensity is more than 50% of its maximum value[*]. Hence, it provides a good indication of the mean characteristics of the region producing the peak of the emission.

Figure 7 shows that the SiO emission peak always occurs in the cool and dense postshock region: $T_{\rm n} \approx 50$ K (Fig. 7a) and a density of 10-40 times the preshock value, close to maximum compression (Fig. 7b). The SiO fractional abundance, $n({\rm SiO})/n_{\rm H}$, is also approximately equal to its maximum value in the shock wave (compare Figs. 7c with 4). The trend to lower SiO abundance at higher $n_{\rm H}$, noted in Sect. 3.2, is clearly visible here. Finally, the ``LVG parameter'', $n({\rm SiO})/({\rm d}v_z/{\rm d}z)$, lies typically in the range 1014-1016 cm-2 km-1 s, implying that the 5-4 line is optically thick at its peak for most models of our grid. In Sect. 4.5, we compare these physical parameters to values inferred previously from LVG analyses of observations, assuming a uniform slab which fills the beam.

   
4.2 Line profiles and peak line temperatures

In Fig. 8, we compare the intensity profiles of various rotational lines, as functions of the flow speed of the neutral fluid, expressed in the frame of the preshock gas, for our reference model: $n_{\rm H} = 10^{5}$ cm-3, $v_{\rm s} = 30$ km s-1, and b = 0.63. The profiles are seen to be narrow (widths of 1-2 km s-1), with similar shapes and peaking within 2 km s-1 of $v_{\rm s}$, as expected for compressed material at the rear of the shock wave[*]. Note that the transition 2-1 peaks further into the cooling flow, because the lower j-levels are repopulated from the higher levels as the temperature falls.

The general shape and centroid velocities are globally similar to those found by Sch97 (their Fig. 3b where the profiles were plotted in the shock frame), although the differences between the various lines are less significant in the present calculations. Also, the emission wing from the fast precursor, at the start of the shock wave, is weaker in the current models, owing to the delay in SiO formation (see Sect. 3.1), making our line profiles narrower than in Sch97. Note that including local thermal broadening in our profile calculations would not significantly change our predicted SiO line width of 1-2 km s-1, because the SiO emission peaks at low temperatures, $T_n \la 100$ K, where the Doppler width is $\sqrt{kT/44 m_{\rm H}} \le 0.1$ km s-1.

In the left column of Fig. 9, we show the predicted peak temperature of the SiO 5-4 line for the grid of models considered in Sect. 3.2, as well as the variation with $j_{\rm up}$ of the peak brightness temperatures of various lines, relative to that of 5-4. The relative intensities have the advantage of being independent of the beam filling factor, and thus they are comparable directly to observations, without prior knowledge of the source size.

The relative peak intensities are within 20% of unity for $j_{\rm up} \le 7$ over a broad range of model parameters ( $v_{\rm s} \ge 30$ km s-1 and $n_{\rm H} < 10^6$ cm-3), owing to the large opacity and near-LTE excitation conditions. For larger values of $j_{\rm up}$, the relative intensities are more dependent on the shock speed and $\approx$1 only when the limiting speed is approached. The absolute peak brightness temperature in the 5-4 line is typically 10-50 K for $v_{\rm s} \ge 30$ km s-1, similar to the kinetic temperature in the emission region, but drops sharply at lower shock speeds, for which the SiO abundance (and opacity) is small. The broken curves in Figs. 9d and h are the results obtained assuming that the initial abundance of O2 ice is negligible, i.e. the second of the two scenarios described in Sect. 3.1.

   
4.3 Integrated line intensities

In the right column of Fig. 9 are presented the integrated intensities (denoted $T{\rm d}V$) of the rotational emission lines of SiO, relative to the 5-4 line, computed for the grid of models considered in Sect. 3.2. There are significant differences between the relative integrated and peak ( $T_{\rm peak}$) line temperatures (right and left columns, respectively, of Fig. 9), owing to systematic variations in linewidth with $j_{\rm up}$, i.e. in the extent of the region where the line is significantly excited. The relative integrated intensities of lines with $j_{\rm up} \ge 7$ remain the most sensitive to the shock speed. As may be seen in Sch97, the maximum value of $T{\rm d}V$ occurs at higher values of the rotational quantum number, $j_{\rm up}$, for higher $n_{\rm H}$. Although the shock temperature varies only weakly with $n_{\rm H}$ (see Fig. 4), a higher density enhances the rates of collisional excitation of the high-j levels, for any given value of the shock speed, $v_{\rm s}$. In Sect. 5, we explore the usefulness of this effect for constraining the preshock density, based on a comparison with actual observations.

The variation of the absolute integrated intensity of the 5-4 rotational emission line with the shock parameters is shown also in Fig. 9h. The bump in $T{\rm d}V$ of the 5-4 transition, for $30 \leq v_{\rm s} \leq 40$ km s-1 and $10^4 \leq n_{\rm H} \leq 10^5$ cm-3 is caused by incomplete O2 destruction in the shock wave, resulting in more rapid SiO formation and warmer emission zones (cf. Fig. 7). As the broken curves in Figs. 9d and h show, this ``bump'' is absent in our second scenario, where O2 is never abundant in the gas phase. At higher shock speeds, the results from the two scenarios become identical, as OH dominates the oxidation of Si in both cases.


  \begin{figure} \par\includegraphics[width=12cm,clip]{Fig.9.eps}\end{figure} Figure 9: a)- c) The peak line temperatures, $T_{\rm peak}$, of the rotational emission lines of SiO, relative to the 5-4 line, as functions of the rotational quantum number of the upper level of the transition, $j_{\rm up}$, for the grid of models in Sect. 3.2. The value of density, $n_{\rm H}$, of the preshock gas is indicated in each panel. d) The absolute peak brightness temperature of the 5-4 line, $T_{\rm peak}(5{-}4)$, as a function of shock speed, $v_{\rm s}$, for all three values of the preshock gas density, $n_{\rm H}$. Shock speeds in excess of 34 km s-1 are absent when $n_{\rm H} = 10^{6}$ cm-3, as they give rise to a J-type discontinuity, and the shock wave is no longer C-type; see Sect. 3.2. The right-hand panels e)-h) show the corresponding values of the integrated line intensities, $T{\rm d}V$. The values of $T{\rm d}V(5{-}4)$ observed in L1157 and L1448 are indicated. In panels  d) and  h), the broken curves show the results obtained assuming that the initial abundance of O2 ice is negligible, i.e. the second of the two scenarios described in Sect. 3.1.
Open with DEXTER

   
4.4 Influence of viewing angle


  \begin{figure} \par\includegraphics[width=7.5cm,clip]{Fig.10.eps}\end{figure} Figure 10: Effect of the inclination angle, $\theta $, on a) the peak, and b) the integrated intensities of the rotational emission lines of SiO, relative to the 5-4 line; c) the integrated and peak intensities of the 5-4 line.
Open with DEXTER

Statistically, there is a low probability that a planar shock should happen to be viewed face-on. Accordingly, we have explored the effects on the SiO rotational line intensities of varying the viewing angle, for the case of our reference model, using the formula (A.16), derived in Appendix A.

The variations of $T_{\rm peak}(5{-}4)$ and $T{\rm d}V(5{-}4)$ with viewing angle, $\theta $, are plotted in the bottom panel of Fig. 10. The peak intensity is almost unaffected by the inclination, as the line is already optically thick for a face-on view (see Eq. (A.16)). On the other hand, the velocity projection reduces the line width, and the integrated intensity, $T{\rm d}V$, declines steadily with increasing viewing angle - by up to a factor of 2 at 75$^\circ$. However, such a variation would be difficult to deduce from observations, given the typical uncertainties in beam filling factors.

Panels (a) and (b) of Fig. 10 illustrate the changes in the (relative to 5-4) peak and integrated SiO line temperatures. Significant changes, compared with viewing face-on, are seen only for inclinations greater than 60$^\circ$ from the normal, and they affect only the optically thin lines $j_{\rm up} \le 3$ and $j_{\rm up} \ge 12$. The main change in the curves for the relative integrated line temperatures (panel (b) of Fig. 10) is that the maximum occurs at lower  $j_{\rm up}$ as $\theta $ increases; this effect could be easily confused with a face-on shock of slightly lower preshock density (cf. Fig. 9). The relative peak temperature (panel (a)) is even more strongly modified, with an upward turn of the curve at $j_{\rm up} \le 4$. The latter characteristic appears to be the only unambiguous signature of a viewing angle >$60^\circ$ in the context of our one-dimensional models.

   
4.5 Influence of the transverse magnetic field strength

In Sect. 3.3, we have seen that the efficiency of sputtering Si from grains decreases monotonically with increasing magnetic field strength. The effect of varying the scaling parameter, b, on the emergent SiO line intensities is shown in Fig. 11, for our reference model.

In panels a and b of Fig. 11 are plotted the predicted line profiles, peak temperatures, and integrated intensities of the SiO 5-4 line, for various values of b. It may be seen that the maximum intensity is reached for intermediate values of $0.63 \la b \la 1$. At smaller b, SiO is less abundant: the shock wave is narrower and hotter, and so O2 is more readily destroyed by H, resulting in incomplete oxidation of Si into SiO. At larger b, the SiO emission decreases owing to less efficient sputtering of Si (see Fig. 5b) at the lower ion-neutral drift speeds. In particular, the predicted intensity drops from $T_{\rm peak} = 10$ K at b = 2 to practically zero at b = 3.

Panels c and d of Fig. 11 show the relative peak and integrated line temperatures, as functions of $j_{\rm up}$, for various values of $b \le 2$(curves for b > 2 are not shown as they lead to negligible SiO emission). It may be seen that the values $0.63 \la b \la 1$, which give rise to the strongest 5-4 emission, also yield the highest relative intensities of lines from $j_{\rm up} \ge 7$. For example, the intensity of the 11-10 line, relative to 5-4, is 3 to 4 times larger than when b = 0.5 or b = 1.5. Comparison with Fig. 9 suggests that this dependence on b may be difficult to distinguish observationally from variations in shock speed. A less ambiguous indication of the value of b might be obtained from the width of the cooling zone, which increases as b2, from 1015 cm for b = 0.4 to 2 $\times $ 1016 cm for b = 2 and for the parameters of our reference model (see Fig. 5a).


  \begin{figure} \par\includegraphics[width=12cm,clip]{Fig.11.eps}\end{figure} Figure 11: a) The SiO 5-4 line temperature, as a function of the flow speed of the neutral fluid, in the reference frame of the preshock gas, for the specified values of the magnetic field parameter, b; b) the peak and integrated intensities of the 5-4 line, as functions of b; c) the integrated and d) the peak intensities of rotational emission lines $j_{\rm up} \rightarrow j_{\rm up}-1$ of SiO, relative to the 5-4 transition, for the specified values of b. All calculations for $v_{\rm s} = 30$ km s-1 and $n_{\rm H} = 10^{5}$ cm-3.
Open with DEXTER


  \begin{figure} \par\includegraphics[width=12cm,clip]{Fig.12.eps}\end{figure} Figure 12: The relative intensities of the rotational emission lines of SiO observed in the outflow sources L1157 and L1448 (Nisini et al. 2007: points with error bars) and predicted by the C-type shock models (curves) with the parameters ($n_{\rm H}$, $v_{\rm s}$) indicated; see text, Sect. 5. The data points which are plotted include a correction for differing beam sizes, which is significant for low-j lines. Full curves denote the best-fitting models of our grid. Broken curves show ``near-miss'' models with a different value of $n_{\rm H}$, which either yield a worse fit to the data points or do not reproduce the absolute $T{\rm d}V$ of the 5-4 line. (See also Fig C.4 of Appendix C.)
Open with DEXTER

   
5 Comparisons with observations

5.1 Rotational line profiles

As Sch97 first noted, the generic SiO line profile predicted by steady planar C-type shock waves, with a peak at high velocity (in the postshock gas) and a tail at lower velocity (in the accelerating precursor), is reminiscent of the SiO line profiles in the L1448 molecular jet (Bachiller et al. 1991). Similarly, we note that the reversed shape of SiO profiles in the L1157 bowshocks, with a peak at low velocity and a high velocity tail (Zhang et al. 1995), could arise if the postshock gas is stationary in the cloud frame, i.e. if one observes the reverse shock, in which the jet is being decelerated. However, in either case, the SiO profiles predicted by our models remain narrower than those observed, with widths of 0.5-2 km s-1, as compared to the observed widths of 5-20 km s-1 in 3 $^{\prime\prime}$-10 $^{\prime\prime}$ beams.

Broader line profiles from steady C-type shocks could arise if Si was sputtered not only from grain cores, as assumed here, but also from SiO-containing grain mantles, with lower binding energy. Then, the SiO abundance would be much enhanced at intermediate velocities, in the precursor (see Sch97). However, owing to the steep temperature decline across the shock wave, this situation results in large variations of the line widths and velocity centroids with the emitting rotational level, $j_{\rm up}$ (cf. Fig. 5 of Sch97). In fact, the observed profiles are very similar from line to line, with the (8-7)/(2-1) intensity ratio showing only modest variations with velocity (see, for example, Fig. 9 of Nisini et al. 2007). These observations suggest that the broad SiO lines are not attributable entirely to intrinsic velocity gradients through a single, planar C-type shock wave. There may be several shock-cooling zones inside the beam, each with a narrow intrinsic profile, which appear spread out in radial velocity owing to a range of inclination angles or propagation speeds, in the observer's frame. This conclusion is supported by interferometric observations of L1157 and L1448 (Guilloteau et al. 1992; Gueth et al. 1998; Benedettini et al. 2007), which reveal systematic velocity gradients across the SiO emitting knots (reminiscent, in some cases, of a bowshock geometry; Dutrey et al. 1997) down to 2 $^{\prime\prime}$-3 $^{\prime\prime}$ resolution; at the distances of L1157 and L1448, the angular dimensions of the SiO emitting regions in Fig. 8, for example, are a few tenths of an arcsec. Such complex two-dimensional modelling lies outside the scope of the present paper. However, we argue in Sect. 5.3 that we may still perform a meaningful comparison of our predicted SiO line intensities with observations of knots in outflows, without reproducing in detail the line profiles, provided that the shock conditions do not vary too much across the beam.

5.2 Narrow SiO lines near ambient velocity

In addition to the typically broad SiO line profiles mentioned above, some outflow regions such as NGC 1333 and L1448 exhibit extremely narrow SiO emission lines, with $\Delta v \approx 0.5$ km s-1, near rest velocity (Lefloch et al. 1998; Codella et al. 1999; Jiménez-Serra et al. 2004, 2005). The corresponding SiO abundance of 10-11-10-10 is two to three orders of magnitude smaller than in the broad SiO components (Codella et al. 1999; Jiménez-Serra et al. 2005). Jiménez-Serra et al. proposed that this feature in L1448 traces a magnetic precursor, where neutral gas is just beginning to accelerate and grain species are starting to be released into the gas phase. However, as we now explain, detailed multifluid shock models do not support this interpretation.

Our SiO line profile calculations show that emission from a magnetic precursor does not give rise to a narrow feature near the speed of the preshock gas. The line intensity increases as the neutral fluid is accelerated, heated, and enriched in SiO - by orders of magnitude by the time that grain-sputtering is complete. Indeed, this deduction could have been made already, on the basis of Figs. 3 and 5 of Sch97, which cover the entire velocity range relevant to predicting the SiO line profiles. Truncation of the precursor when the neutral fluid has been accelerated to only $v_{\rm n} = 0.5$ km s-1 would imply a very finely-tuned shock age, a circumstance which appears to us to be improbable. Furthermore, an ion-neutral drift speed of at least 5 km s-1 is needed to start releasing species from grain mantles, where binding energies are a few tenths of an eV (Flower & Pineau des Forêts 1994), and of at least 20 km s-1 to start sputtering grain cores (May et al. 2000). However, the H13CO+ line does not show evidence of this predicted acceleration: its emission peak is shifted by only +0.5 km s-1from the velocity of the ambient gas, like the narrow SiO feature (Jiménez-Serra et al. 2004).

We believe that a more likely explanation of the narrow SiO feature in L1448 is that it traces Si-enriched postshock material that has been decelerated by and mixed with the ambient gas, as proposed originally by Lefloch et al. (1998) and Codella et al. (1999) in connection with other regions. Given a shock speed $v_{\rm s} \leq 30$ km s-1 and the high ambient density characteristic of Class 0 protostellar envelopes, deceleration could be achieved readily within the L1448 flow age of approximately 3500 yr. The low SiO fractional abundance would then be a consequence of mixing with SiO-poor ambient gas. Alternatively, the narrow feature might arise in a reverse C-type shock, where outflow material at v < 20 km s-1is brought almost to rest by the much denser ambient medium, and the shock speed is too low to produce abundant SiO. Both interpretations are consistent with NH3 observations of dense gas in the envelope of the L1448 protostar, with radial velocity and spatial extent similar to that of the narrow SiO feature and signs of heating near the path of the fast L1448 jet (Curiel et al. 1999).

   
5.3 SiO line intensities

If our explanation of SiO profile broadening is correct, one could in principle recover the parameters of each individual emission zone in the beam by analysing the relative intensity ratios as functions of velocity. Unfortunately, such data are currently quite noisy and not yet available for a wide range of values of $j_{\rm up}$. Furthermore, knowledge of the beam-filling factor as a function of velocity would be necessary to obtain absolute intensities and remove ambiguity in the shock parameters; but this would require sub-arcsecond angular resolution, which is not yet available. Nevertheless, one may still derive some approximate beam-averaged shock properties, if all of the shock components have similar excitation conditions, as is suggested by single-dish data, which show the line ratios to be insensitive to velocity, v. In this case, the observed profile will be simply a convolution of the individual, narrow shock profiles with the (unknown) filling-factor, $\phi(v)$. The observed absolute $T{\rm d}V$ is simply that for a single shock, multiplied by the total beam filling factor of the SiO-emitting region in the beam, $f = \int{\phi(v) {\rm d}v}$, as inferred from its overall size in single-dish maps. The values of $T{\rm d}V$ for different $j_{\rm up}$, relative to the 5-4 transition, remain unchanged compared to a single shock, because f cancels out in the ratios, thereby enabling direct comparison with our models. In the following, we assume that this situation prevails.

By way of illustration of the applicability of the shock models, we show in Fig. 12 the relative integrated intensities of the rotational transitions of SiO observed in the outflow sources L1157 and L1448 (Nisini et al. 2007) and predicted by the grid of models, whose parameters are specified; in all cases, the magnetic field scaling parameter b = 1. As noted in Sect. 4.5, this value of b yields the largest relative intensities of the high-j lines and hence will yield a lower limit to the shock speed required to reproduce the observations (except at $n_{\rm H} = 10^{5}$ cm-3, where high-j excitation is a non-monotonic function of $v_{\rm s}$; see Fig. 9f). We assume also that the shock wave is viewed face-on; if the true inclination exceeds $60^\circ$, this assumption results in the preshock density being slightly underestimated (see Sect. 4.4). The models shown as the full curves are those which provide the best fits to the observations. In order to illustrate how well the shock parameters are constrained, we plot as dashed curves ``near-miss'' models that fit most of the data points, or fit all points but do not reproduce the absolute intensity of the 5-4 line (see below).

Figure 12 demonstrates that steady-state C-type shocks with Si-sputtering from grain cores can reproduce successfully the relative integrated intensities of SiO lines in these molecular outflows. Furthermore, Fig. 9h shows that the models can reproduce also the absolute integrated intensity of SiO 5-4 with the estimated beam filling factors $f \approx 1$ in L1157 and L1448-R4 and $f \simeq 1/4$ in L1448 R1 and B1 (cf. Nisini et al. 2007). Alternative fits with higher density and lower shock speeds ( $n_{\rm H} = 10^{6}$ cm-3 and $27 \la v_{\rm s} \la 30$ km s-1; $n_{\rm H} = 10^{5}$ cm-3 and $v_{\rm s} = 25$ km s-1) underestimate $T{\rm d}V(5{-}4)$ and are thus ruled out. The shock speed, $v_{\rm s}$, is constrained to within 15 km s-1 and the preshock density to within a factor of 10, with one of our grid values of $n_{\rm H}$ yielding a clear best fit in all cases.

In Appendix C, we present, in Fig. C.4 results equivalent to those in Fig. 12, but for the secondary grid of models, in which oxygen is initially in the form of H2O ice rather than O2 ice. Again, the relative intensities can be well reproduced by steady-state C-type shocks. The absolute SiO  $T{\rm d}V(5{-}4)$ favour the low $n_{\rm H}$, high $v_{\rm s}$ cases (cf. the dashed curved in Fig. 9h). The best-fit shock parameters and the inferred physical conditions at the (5-4) line peak are almost unchanged, as the range of shock speeds is such that oxidation of Si by OH is dominant.

Table 1: Properties of SiO-emission regions deduced from face-on C-type shock models or homogeneous-slab LVG models.

It is instructive to compare the physical parameters at the SiO peak of our best-fit, steady-state C-type shock models to those previously inferred from an LVG analysis, assuming a slab of constant density, temperature, and velocity gradient along the line of sight (Nisini et al. 2007). From Table 1, it may be seen that the ``slab LVG'' approach yields similar values of the density to our shock models but overestimates by a factor 5-10 the kinetic temperature and underestimates by several orders of magnitude the Sobolev LVG opacity parameter. The cool postshock layer emits over a narrow velocity range and needs to be more optically thick in SiO to produce the same $T{\rm d}V$ as a hot slab with a large velocity gradient along the line of sight. On the other hand, similar SiO abundances are deduced using both approaches, to within typically a factor of 3.

We note that the models should be able to simulate also the other spectral observations of the sources, notably the H2 line intensities. In a forthcoming publication, we shall consider in detail the outflow L1157 and make a more comprehensive comparison of its observed spectrum with the predictions of shock models.

   
6 Concluding remarks

We have considered the structure of C-type shock waves propagating into molecular gas containing amorphous carbon and silicate grains; olivine (MgFeSiO4) was chosen as the representative silicate-grain material. We find that the degree of sputtering of silicon from the grains is smaller, by about an order of magnitude, than was predicted by the calculations of Sch97, owing partly to higher sputtering thresholds and lower sputtering yields but principally to the reduced width of the shock wave in the present calculations. The reduction in the shock width is a consequence of a more accurate treatment of the coupling between the neutral fluid and the charged grains in the current model.

A grid of C-type shock models has been computed, for values of the preshock gas density and the shock speed which are believed to span the ranges of these parameters in molecular outflows; two scenarios were considered regarding the initial distribution of oxygen in the gas and solid phases. The maximum speed of a C-type shock is limited by the inertia of the (charged) grains in the preshock gas and by the collisional dissociation of molecular hydrogen within the shock wave, which can lead to a J-type discontinuity (Le Bourlot et al. 2002; Flower & Pineau des Forêts 2003). We find that significant sputtering of the grain cores occurs only for shock speeds $v_{\rm s} \ge 25$ km s-1 and moderate magnetic field strengths close to equipartition with the cloud kinetic energy (magnetic field parameter $0.5 \le b \le 2$). However, we note that the sputtering threshold energy is determined principally by the so-called ``displacement energy'', $E_{\rm D}$, of the material, whose value for the silicates of relevance here remains uncertain (May et al. 2000). The sputtering yields increase rapidly from threshold, and a reduction in the threshold energy would enable significant sputtering to occur at lower shock speeds or higher magnetic field strengths.

We find that, in the absence of silicon-containing grain mantles, SiO line emission in steady-state planar C-type shock waves arises predominantly from cool postshock gas, close to maximum compression, with negligible emission from the precursor. Except at the lowest shock speeds, the SiO emission is optically thick and close to LTE for $4 \le j_{\rm up} \le 7$. The relative line intensities, as functions of $j_{\rm up}$, together with the absolute 5-4 line intensity, provide good diagnostics of the shock parameters, $n_{\rm H}$ and $v_{\rm s}$. The influence of the viewing angle and transverse magnetic field strength is found to be relatively minor over the typical ranges of their values.

Our shock models provide good fits to both the relative and absolute SiO intensities in the molecular outflows L1157 and L1448; the SiO fractional abundance is deduced to be in the range 4 $\times $ $10^{-8} \la n({\rm SiO})/n_{\rm H} \la 3$ $\times $ 10-7. The emission regions of the shock wave are much colder, more optically thick, and have 10 to 100 times greater SiO column density than estimated previously from optically thin LVG slab models (Nisini et al. 2007). Our results are in line with a recent analysis of interferometric maps of the HH212 jet (Cabrit et al. 2007), demonstrating that the SiO emission is optically thick and close to LTE, with an intrinsic peak brightness temperature of approximately 50 K. On the other hand, the line profiles predicted by our planar C-type shocks are typically 10 times narrower than observed in L1157 and L1448, suggesting that the single-dish beam includes shocks with various inclinations and speeds, and/or mixing layers. Detailed modelling of the line ratios as functions of velocity, with a careful correction for differing beam widths, would be needed to clarify this issue.

Whilst the grid of models presented here is intended to provide a guide to interpreting observations of outflow sources, it should be recalled that the dynamical timescales which characterize these regions are often too short to enable C-type shock waves to attain their steady-state structure (Chièze et al. 1998; Lesaffre et al. 2004). The presence of an embedded J-type discontinuity has then to be considered when modelling specific sources. Studies of the outflow source in Orion (Le Bourlot et al. 2002), of jets associated with low-mass star formation (Giannini et al. 2004, 2006; McCoey et al. 2004), and of the supernova remnant IC 443 (Cesarsky et al. 1999) have shown that the rovibrational line spectrum of H2 can be reproduced successfully by hybrid shock waves, but not by pure C- or J-type shocks. The influence of the discontinuity on the C-component (magnetic precursor), and hence on the formation of SiO and its emission line spectrum, becomes important for ages smaller than the flow time to maximum compression.

The existence of Si-containing grain mantles is another circumstance that would significantly modify the intensities of SiO lines and their profiles. The fractional abundance of SiO in the warm shock precursor would be enhanced, compared with the models considered here, in which Si is sputtered exclusively from olivine grain cores. The effects of non-steady C-type shocks and Si-containing mantles will be considered in a forthcoming paper.

Acknowledgements
Antoine Gusdorf and the University of Durham acknowledge the support of the European Commission under the Marie Curie Research Training Network ``The Molecular Universe'' MRTN-CT-2004-512302. We thank Brunella Nisini for helpful correspondence relating to the SiO observations of Nisini et al. (2007). We thank also Paul Goldsmith and Laurent Pagani for information regarding SWAS and Odin observations of O2.

References

 

   
Appendix A: SiO radiative transfer

A.1 Photon escape probabilities

The SiO rotational level populations and excitation temperatures in our planar shock models are calculated by means of a large velocity gradient (LVG) method. This approximate treatment assumes that, owing to the macroscopic velocity field and resulting Doppler shifts, emitted line photons are either re-absorbed locally or escape to infinity. The escape probability of a line photon in a given direction  $\vec{\hat{s}}$ is then (see Surdej 1977, for a pedagogical derivation):

\begin{displaymath}% \beta_s = \frac{1 - {\rm e}^{-\tau_s}}{\tau_s} \end{displaymath} (A.1)

where the ``LVG optical depth'' $\tau_s$ is defined as

 \begin{displaymath}% \tau_s = \frac{hc}{4 \pi} \frac{n_{\rm l}}{\partial (\vec{v... ...eft( 1 - \frac{g_{\rm l}n_{\rm u}}{g_{\rm u}n_{\rm l}}\right), \end{displaymath} (A.2)

with ${\partial (\vec{v} \cdot \vec{\hat{s}}) /\partial{s}}$ the radial velocity gradient along direction  $\vec{\hat{s}}$, nl and nu the number density of molecules in the lower and upper levels of the transition, respectively, $g_{\rm l}$ and $g_{\rm u}$ the corresponding statistical weights, and $B_{\rm lu}$ the Einstein coefficient for stimulated absorption. The mean intensity of the radiation field at the local frequency $\nu$ of the transition, averaged over all angles, may then be expressed in terms of local quantities only as

 \begin{displaymath}% \bar{J_\nu} = S_\nu(1-\bar\beta)+I_{\rm c}\bar\beta, \end{displaymath} (A.3)

where $\bar\beta = \int{\beta_s {\rm d}\Omega}/4\pi$ is the escape probability averaged over all solid angles; $I_{\rm c}$ is the mean intensity of the continuum radiation field, taken to be a blackbody (Planck) function $B_\nu(T)$ at the temperature of the cosmic background, $T_{\rm bg} = 2.73$ K; and $S_\nu = B_\nu(T_{\rm ex})$ is the local source function, where the excitation temperature  $T_{\rm ex}$ of the transition is defined through $n_{\rm u}/n_{\rm l} \equiv g_{\rm u}/g_{\rm l} \exp (-h\nu/k_{\rm B} T_{\rm ex})$. Two expressions have been used for the average escape probability $\bar\beta$: The first expression was ultimately adopted in the present study for consistency with our one-dimensional shock geometry, whereas the second was used by Sch97. Since $\bar\beta$ is always smaller in the plane-parallel case (owing to the reduced velocity gradients at small $\mu $), photon trapping is more efficient and the excitation temperatures are increased compared to the isotropic approximation. Thus, adopting the expression of Neufeld & Kaufman (1993) leads to larger integrated SiO line intensities; this effect is illustrated in Fig. A.1 for our reference shock model. It is marginally significant for small and large values of the rotational quantum number but reaches a factor 3 for $j_{\rm up} \approx 12$.

  \begin{figure} \par\includegraphics[width=9cm,clip]{Fig.A1.eps}\end{figure} Figure A.1: The integrated intensities of the rotational transitions of SiO. Full curves correspond to the Neufeld & Kaufman (1993) expression (A.4) for the escape probability in a plane-parallel flow, dashed curves to the isotropic case, Eq. (A.5). The model parameters are $n_{\rm H} = 10^{5}$ cm-3, B = 200 $\mu $G, and $v_{\rm s} = 30$ km s-1.
Open with DEXTER

  \begin{figure} \par\includegraphics[width=9cm,clip]{Fig.A2.eps}\end{figure} Figure A.2: Consistency of the LVG method: velocity gradient criterion (A.6). The model parameters are $n_{\rm H} = 10^{5}$ cm-3, B = 200 $\mu $G, and $v_{\rm s} = 30$ km s-1.
Open with DEXTER

A necessary condition for the local LVG approximation to be valid is that, over the characteristic distance L where physical conditions vary, the velocity shift $\Delta v$ arising from the velocity gradient should be larger than the local line thermal width  $v_{\rm th}$. Then, line photons are re-absorbed only within a region of size <L, where physical conditions and SiO excitation are uniform. This criterion can be rewritten as

 \begin{displaymath}% \vert{\rm d}v_z/{\rm d}z\vert > v_{\rm th}/L, \end{displaymath} (A.6)

where

\begin{eqnarray*}v_{\rm th} = \sqrt{8k_{\rm B} T_{\rm n} \over \pi m}, \end{eqnarray*}


$T_{\rm n}$ is the temperature of the neutral fluid, and m is the mass of the molecule.

In Fig. A.2 we compare the left and right hand sides of (A.6) for our reference shock model, using $L \approx T_{\rm n}/\vert{\rm d}T_{\rm n}/{\rm d}z\vert$ as a characteristic distance. The LVG criterion may be seen to be verified throughout the cooling flow of the shock wave, where the bulk of SiO emission arises. It is not verified in the far postshock region, where the computed velocity gradient tends to zero; but this region makes a negligible contribution to the line flux, owing to the low escape probabilities.

   
A.2 Numerical implementation

In the diatomic SiO molecule, electric dipole transitions take place only between adjacent rotational levels of the ground vibrational state ( $\Delta j = \pm 1$ with j the rotational quantum number), while collisionally induced transitions can, in principle, connect any pair of levels (in practice they become less probable with increasing level separation). The evolution of the population density ni of level i may thus be written in the following matrix form:

 
                      $\displaystyle % \frac{{\rm d}n_{i}}{{\rm d}t}$ = $\displaystyle n_{i+1} (A_{i+1,i} + B_{i+1,i}{\bar{J}_{i,i+1}})$  
    $\displaystyle - n_{i} (A_{i,i-1} + B_{i,i+1}\bar{J}_{i,i+1} + B_{i,i-1}\bar{J}_{i-1,i})$  
    $\displaystyle + n_{i-1}B_{i-1,i}\bar{J}_{i-1,i}$  
    $\displaystyle + \sum_{\rm coll} \sum_{j \neq i} n_{\rm coll} (n_{j}C_{ji} - n_{i}C_{ij}),$ (A.7)

where $n_{\rm coll}$ is the number density of each collisional partner (here H2 and He), Cij is the collisional rate coefficient from level i to level j, and A and B denote the Einstein coefficients for spontaneous and induced radiative transitions. Two approaches to solving for the level populations may be taken: The former approach is preferable, particularly when the molecule under consideration is an important shock coolant (e.g. H2) - which is not the case of SiO. Furthermore, the flow time in the SiO emission zone is sufficiently long that the assumption of local statistical equilibrium is justified. Accordingly, we adopted the latter approach in the present work. It is common to most applications of the LVG method to astrophysical problems (e.g. Sch97, Neufeld & Kaufman 1993) and has the advantage of being less demanding in CPU time. For more significant coolants, such as CO, 13CO, and H2O, the rate of cooling was computed in parallel with the shock dynamics, using the cooling functions calculated by Neufeld & Kaufman (1993), by means of the LVG method.

Note that Eqs. (A.7) may be put in a different (matrix) form by expressing explicitly ${\bar J}$ as a function of $B_\nu(T_{\rm ex})$ and $I_{\rm c}$. Using the standard relationships between Einstein coefficients, and the definition of $T_{\rm ex}$ in terms of nj and nj+1, all terms involving $T_{\rm ex}$ cancel from the equations, leaving radiative terms which are proportional to the $\bar\beta$'s (see Goldreich & Kwan 1974):

 
                             $\displaystyle % \frac{{\rm d}n_{i}}{{\rm d}t} =$ + $\displaystyle n_{i+1} \bar\beta_{i,i+1} \left (A_{i+1,i} + B_{i+1,i}I_{\rm c} \right )$  
  - $\displaystyle n_{i}\left(\bar\beta_{i-1,i}A_{i,i-1}+\bar\beta_{i,i+1}B_{i,i+1}I_{\rm c}+\bar\beta_{i-1,i}B_{i,i-1}I_{\rm c} \right )$  
  + $\displaystyle n_{i-1} \bar\beta_{i-1,i} B_{i-1,i}I_{\rm c}$  
  + $\displaystyle \sum_{\rm coll} \sum_{j \neq i} n_{\rm coll} (n_{j}C_{ji} - n_{i}C_{ij}).$ (A.8)

Although the two formulations are equivalent and converge to the same solution, we have found that inversion of the latter matrix encounters numerical instabilities and convergence problems at high optical depths, possibly due to round-off errors in the (vanishingly small) $\bar\beta$ terms. In the first formulation, the radiative elements of the matrix are never zero, even at high opacity, and we obtain much better convergence and accuracy. Hence we have adopted (A.7) in the present calculations. The ``lambda-iteration'' is terminated when we reach convergence of the level populations to 1 part in 104.

Our code was tested thoroughly against the routine used by Sch97, in the case $\bar\beta$ = $\beta_\perp$, and against the web-based online version of RADEX (www.sron.rug.nl/vdtak/radex/radex.php), which uses yet another expression for the escape probability, valid for a turbulent homogeneous sphere (Van der Tak et al. 2007). Discrepancies never exceeded a few percent.

The MHD shock code provides the physical and chemical profiles (of the temperatures, densities, velocities, and abundances) which are required in order to apply the LVG technique. For SiO, we used the rate coefficients for collisional de-excitation published by Turner et al. (1992), for which the collision partner is ground state para-H2. These data are available for rotational quantum numbers $0 \le J \le 20$ and for kinetic temperatures T = 20, 40, 70, 100, 150, 200, 250, 300 K and are interpolated to intermediate values of T. We made use also of the extrapolated data of Schöier et al. (2005), which extend to J = 40 and T = 2000 K. Subsequent calculations, using the rate coefficients of Dayou & Balança (2006) for collisions of SiO with para-H2, have shown that the rotational line intensities are insensitive to these collision rates because the lines are formed under conditions which approach LTE. For collisions of SiO with He, we used the rate coefficients of Dayou & Balança (2006), for $J \le 26$ and kinetic temperatures in the range $10 \le T \le 300$ K. At temperatures higher than the maximum for which the rate coefficients were calculated, we assumed that they remain constant. Upwards (excitation) rate coefficients were obtained from detailed balance.

The Einstein A-values, and the rotational constant B0 = 21711.967 MHz (corresponding to $hB_0/k_{\rm B} = 1.042$ K) for the ground vibrational state of 28Si16O were taken from the NIST database (www.nist.gov/data/).

   
A.3 Emergent line intensities and radiation temperatures

Each layer of the shock wave, of thickness $\vert\Delta z\vert$, elementary surface $\Delta S$, and velocity vz, emits the following luminosity (in erg s-1) over all directions in a given transition $j+1 \rightarrow j$ of frequency $\nu = 2B_0(j+1)$:

 \begin{displaymath}% F_\nu(v_z) = \bar{\beta} ~ {h\nu}~ n_{j+1}A_{j+1,j} ~ \mid \Delta z \mid ~ \Delta S. \end{displaymath} (A.9)

In order to compute line profiles, however, we need to evaluate the specific intensity per unit solid angle, frequency interval, and projected emitting area (in erg cm-2 s-1 Hz-1 sr-1). Following Sch97, we assume that the shock is viewed face-on, i.e. in the z-direction, normal to the shock front. The intensity is

 \begin{displaymath}% I_\perp(v_z) = \frac{\beta_\perp }{4\pi} ~ h\nu ~ n_{j+1}A_{j+1,j} \frac{\mid \Delta z\mid}{\vert\Delta \nu_z\vert}, \end{displaymath} (A.10)

where ${\vert\Delta \nu_z\vert}$ is the Doppler width of the layer, viewed along the z-axis:

 \begin{displaymath}% {\vert\Delta \nu}_z\vert = \frac{\nu}{c} \vert \Delta{v_z}\... ...u}{c} \vert\Delta z\vert \times \vert{\rm d}v_z/{\rm d}z\vert. \end{displaymath} (A.11)

Consequently the intensity becomes

 \begin{displaymath}% I_\perp(v_z) = \frac{hc}{4\pi} \frac{\beta_\perp}{\mid {\rm d}v_z/{\rm d}z \mid} n_{j+1}A_{j+1,j}. \end{displaymath} (A.12)

Noting that $\beta_\perp$= $(1-{\rm e}^{-\tau_\perp})$/ $\tau_\perp$ and using the definition of $\tau$ in Eq. (A.2) as well as the relationships between A and B-Einstein coefficients, one obtains the alternative, simpler expression

 \begin{displaymath}% I_\perp(v_z) = B_\nu(T_{\rm ex}) \left (1 - {\rm e}^{-\tau_\perp} \right ) \end{displaymath} (A.13)

which corresponds to the well-known radiative transfer result for a uniform slab of excitation temperature  $T_{\rm ex}$ and opacity  $\tau_\perp$.

In order to compare with actual observed SiO spectra, one needs to add the cosmic background, and take into account the ON-OFF subtraction applied to radio spectra (to remove atmospheric noise). The cosmic background at the ON position and velocity vz is $B_\nu(T_{\rm bg})$ exp(- $\tau_\perp$), due to attenuation by the layer, while at the OFF position it is simply $B_\nu(T_{\rm bg})$. The observed intensity is thus

 
$\displaystyle % I_\perp^{\rm obs}(v_z) = {\tt ON - OFF}$ = $\displaystyle \left[ B_\nu(T_{\rm ex}) - B_\nu(T_{\rm bg})\right ] \left (1 - {\rm e}^{-\tau_\perp} \right )$  
  = $\displaystyle I_\perp(v_z) \times \left[ 1 - B_\nu(T_{\rm bg})/B_\nu(T_{\rm ex})\right].$ (A.14)

Finally, we convert the observed specific intensity of the layer to a line radiation temperature $T_{\rm R}$ (in kelvin) using the definition

 \begin{displaymath}% T_{\rm R} \equiv \frac{I^{\rm obs} c^2}{2k_B \nu^2} \end{displaymath} (A.15)

which is standard in radioastronomy. Once radiation temperatures are calculated for each shock layer, and thus each vz, the line profile is integrated over vz to yield integrated line intensities, $T{\rm d}V$, in K km s-1. Note that the specific intensities, radiation temperatures and $T{\rm d}V$ are valid for a shock that entirely fills the observing beam. Otherwise, the observed brightness temperature has to be reduced by an appropriate ``surface filling factor'', $f \approx \Delta S$/(beam area).

A.4 Line profile at arbitrary inclinations

Table A.1: The numerical values of the coefficients, ai(X), of the polynomial fits to the fractions of X $\equiv $ Fe, Si and Mg sputtered from olivine grains. Numbers in parentheses are powers of 10.


  \begin{figure} \par\includegraphics[width=9.5cm,clip]{Fig.B1.eps} % \end{figure} Figure A.3: Fits of the fractions of Fe, Si and Mg sputtered from olivine grains, as functions of the shock speed. The points on the curves indicate the limiting speeds of steady-state C-type shock waves, for preshock gas densities $n_{\rm H} = 10^{5}$ and 106 cm-3, with the lower speed corresponding to the higher density. The parameters of the fits to the curves are given in Table B.1.
Open with DEXTER

In the general case of an arbitrary viewing angle, $\mu = \cos\theta$, the term  $\tau_\perp$ in Eq. (A.13) is replaced by the LVG opacity in the chosen direction, $\tau(\mu) = \tau_\perp/\mu^2$. Therefore

 \begin{displaymath}% T_{\rm R}(\mu) = {T_{\rm R}}_\perp \frac{\left (1 - {\rm e}... ...^2} \right )} {\left (1 - {\rm e}^{-\tau_\perp} \right )}\cdot \end{displaymath} (A.16)

Thus, $T_{\rm R}$ will be multiplied by $1/\mu^2$ in the optically thin regime and will be unchanged in the optically thick regime. At the same time, the velocity of the layer, vz, is replaced by its projection along the photon path, $\mu v_z$, so the line profile becomes narrower. A view which is not face-on will leads to a larger $T{\rm d}V$ for optically thin lines, and to a smaller $T{\rm d}V$ for optically thick lines (assuming that the source still fills the beam entirely). A quantitative evaluation of this effect, for our reference shock model, is presented in Sect. 4.4.

Appendix B: Fractions of Fe, Si and Mg sputtered from olivine

 

In Sect. 3.2, we presented and discussed the sputtering of Fe, Si and Mg from grains composed of olivine (MgFeSiO4). In Fig. B.1, we show the numerical fits to these results, as functions of the shock speed; they are practically independent of the preshock gas density. Also shown in this figure are the limiting speeds of steady-state C-type shock waves for preshock densities $n_{\rm H} = 10^{5}$ and 106 cm-3; the limiting speed is lower for the higher density. The limit is associated with the thermal runaway which occurs, owing to the collisional dissociation of H2 within the shock wave. The numerical values of the coefficients of the polynomial fits, of the form

\begin{eqnarray*}y({\rm X}) = \sum _{i=0}^5 a_i({\rm X})v_{\rm s}^i, \end{eqnarray*}


are given in Table B.1; the shock speed, $v_{\rm s}$, is in km s-1. We have assumed that the magnetic field parameter b = 1.

   
Appendix C: Initial gas-phase abundance of O2


  \begin{figure} \par\includegraphics[width=12cm,clip]{Fig.C1.eps}\end{figure} Figure C.1: As Fig. 4, but assuming that the initial abundance of O2 ice is negligible (the second of the two scenarios described in Sect. 3.1).
Open with DEXTER

We present, in Figs. C.1-C.4, the results corresponding to those in Figs. 4, 7, 8 and 12, respectively, of the main text, but derived from our secondary grid of models, in which $n({\rm O}_2)/n_{\rm H} = 1.0$ $\times $ 10-7 initially in the gas phase but the excess oxygen is in form of H2O ice, rather than O2 ice as in our primary grid. We recall that the initial fractional abundance of O2 in chemical equilibrium is $n({\rm O}_2)/n_{\rm H} \approx 10^{-5}$.

In these models, O2 is never abundant in the gas phase, and Si oxidation occurs only in reaction 2, with OH. Thus, SiO formation is less efficient at intermediate shock speeds and not significant at $v_{\rm s} = 25$ km s-1.


  \begin{figure} \par\includegraphics[width=8.5cm,clip]{Fig.C2.eps}\end{figure} Figure C.2: As Fig. 7, but assuming that the initial abundance of O2 ice is negligible (the second of the two scenarios described in Sect. 3.1).
Open with DEXTER


  \begin{figure} \par\includegraphics[width=9cm,clip]{Fig.C4.eps}\end{figure} Figure C.3: As Fig. 8, but assuming that the initial abundance of O2 ice is negligible (the second of the two scenarios described in Sect. 3.1).
Open with DEXTER


  \begin{figure} \par\includegraphics[width=12cm,clip]{Fig.C3.eps}\end{figure} Figure C.4: As Fig. 12, but assuming that the initial abundance of O2 ice is negligible (the second of the two scenarios described in Sect. 3.1).
Open with DEXTER


Copyright ESO 2008