© ESO, 2014

1. Introduction

The formation mechanism of high-mass stars (M> 8 M_⊙) has been an open question for the past few decades, the main reason being that the strong radiation pressure exerted by the young massive star overcomes the star’s gravitational attraction (Kahn 1974). Different mechanisms have been proposed to lead to the formation of massive stars, some based on disc-mediated accretion (e.g., McKee & Tan 2003; Krumholz et al. 2005), others on alternative scenarios where massive cores accrete mass from their surroundings while collapsing (e.g., Bonnell & Bate 2006).

Molecular outflows are predicted in both scenarios and can help distinguish between models of massive star-formation. In particular, a study of their properties and their evolution with time would allow us to understand the different phases involved in the formation of massive stars. A full understanding of the properties of outflows as function of mass of the powering source would also help to evaluate whether there are similarities with low-mass star formation because the relations derived for low-mass young stellar objects (YSOs; e.g., Bontemps et al. 1996) should hold over the whole range of masses if massive stars form in the same way as low-mass stars. Given their relatively large scales (~1 pc), bipolar molecular outflows often represent the easiest, if not the only, means to observationally investigate high-mass star formation in a statistically significant number of objects.

SiO thermal emission is the best tool to investigate molecular outflows. Its appearance in the gas phase is attributed to sputtering of silicon-bearing material from grains (Gusdorf et al. 2008a,b), supplemented by grain-grain interactions (e.g. vaporisation and shattering, Guillet et al. 2007, 2009, 2011; Anderl et al. 2013) in higher-density regions. As a result, SiO suffers minimal contamination from the more quiescent gas. Observations towards low-mass protostars point to an evolutionary picture in which bipolar jets and outflows are always associated with young Class 0 objects (e.g., Gueth & Guilloteau 1999; Codella et al. 2007), while Class I/II objects show less evidence of powerful jets. Finding an analogous trend in high-mass YSOs would be strongly suggestive of a similar formation scenario for low- and high-mass stars.

In an attempt to investigate how the properties of molecular outflows change with the evolution of the source in massive YSOs, López-Sepulcre et al. (2011) and Sánchez-Monge et al. (2013) investigated massive parsec-scale molecular clumps in SiO transitions with the IRAM 30 m telescope. They concluded that the SiO abundance decreases with time, favouring a scenario in which the amount of SiO is mainly enhanced in the first evolutionary stages, probably by strong shocks produced by the protostellar jet, and decreases as the object evolves, most likely because of a decaying jet activity.

A limitation of these studies is that they used low -excitation SiO lines with upper level energies ranging from 6 K (J = 2−1)¹ to 30 K (J = 5−4). Models of SiO emission show that a more exhaustive study of the excitation of SiO needs higher J transitions: Nisini et al. (2007) demonstrated that the SiO(8−7) and SiO(3–2) lines can be used to distinguish between different excitations in shocks as they cover a wider energy range (E_{lowSiO(3−2)} ~ 12 K, E_{lowSiO(8−7)} ~ 73 K). Higher J SiO lines are also needed to distinguish between different shock conditions (Gusdorf et al. 2008a). In the current study, we present observations of a sub-sample of sources from López-Sepulcre et al. (2011) in the CO(4−3) and SiO(8−7) transitions with the APEX telescope with a resolution that well matches that of the SiO(3−2) data presented by López-Sepulcre et al. (2011). We aim to study the SiO(3−2) and SiO(8−7) data (1) to verify whether the decrease in SiO luminosity with evolution seen in previous studies is confirmed when excitation effects are taken into account; (2) to derive a direct measurement of the SiO abundance using CO(4−3); (3) to determine whether the SiO abundance changes with evolution; (4) to discuss the properties of SiO as function of evolution; and (5) to compare for the first time SiO observations with shock models that include grain-grain processing at the high densities typical of high-mass star-forming clumps. In Sect. 2 we discuss our selection criteria and present the properties of the observed sources. In Sect. 3 we present the observations, while in Sect. 4 we discuss the observational results. Finally in Sect. 5 we compute the SiO line luminosities and abundances for our sources, study their behaviour as function of evolution, and interprete our observations in the framework of shock models and a non-local thermodynamic equilibrium (non-LTE) analysis.

2. Selection of the sample

From the sample of López-Sepulcre et al. (2011), we selected fifteen sources with a line-width in SiO(3−2) broader than 40 km s^-1 and a peak intensity higher than 0.2 K. The third criterion was that the sources were observable with the APEX telescope². The list of observed sources and their properties is given in Table 1. Luminosities and masses come from López-Sepulcre et al. (2011) and, where available, from Sánchez-Monge et al. (2013), who improved the previous spectral energy distribution fits with data from the Herschel infrared Galactic Plane Survey (Molinari et al. 2010). For IRAS 18151-1208_2, there were not enough continuum data available to model the spectral energy distribution, which resulted in a lack of mass and luminosity estimates for this source. The source G19.61–0.24A was assumed to be at the near distance of 3.8 kpc by López-Sepulcre et al. (2011). However, Kolpak et al. (2003) located it at the far distance of 12.6 kpc thanks to HI interferometric observations, which show absorption features at velocities higher than the systemic velocity. Therefore, we here adopt the far distance for the source, and correct the mass and luminosity from López-Sepulcre et al. (2011) for the new value. We note, however, that the mass and luminosity of G19.61−0.24A have large uncertanties because of the complexity of the source and the coarse spatial resolution of the data used to fit its spectral energy distribution. For example, Furuya et al. (2011) estimated a luminosity of ~10⁵ L_⊙ for the powering source of the main outflow in the G19.61−0.24A cluster, while the total mass of the clumps in the region is 1580 M_⊙ from interferometric measurements.

Fig. 1 SiO(2−1) luminosity against Lbol/M for the full sample of sources studied by López-Sepulcre et al. (2011). Full circles depict the sources observed in SiO(8−7) that we studied here.

The variation of the SiO(2−1) luminosity as function of the bolometric luminosity-to-mass ratio, L_bol/M (a rough estimator of evolution in the star formation process, e.g., Molinari et al. 2008; Ma et al. 2013), for the selected sources is shown in Fig. 1. Although our selection criteria bias us towards the more luminous sources in SiO(2−1) of the sample of López-Sepulcre et al. (2011), we cover a wide range of L_bol/M values up to ~350 L_⊙/M_⊙.

To probe the evolutionary phase of the clumps, we investigated the emission at 22 and 8 μm in the WISE and Midcourse Space eXperiment (MSX) missions. We expect the youngest sources in the sample to be still dark at 22 μm and the most evolved to be detected at 22 and 8 μm. The association of the sources with mid-IR emission at 21−22 μm is described by Csengeri et al. (2014). Sources are classified as 22 μm-quiet if there is no association with a WISE source and as 22 μm-loud in case of association, and correspondingly infrared loud at 8 μm (8 μm-loud) and infrared quiet at 8 μm (8 μm-quiet) if associated with 8 μm emission or not. Table 1 summarises the IR properties of the sources.

3. Observations

All observations were performed with the APEX telescope. The majority of the data was observed under Max Planck Society (MPS) time in April and November 2011 in wobbler-switching mode, with a wobbler throw of 120′′. Additional data were taken under ESO programme 089.C-0203 in September 2012. For the observations performed in April and November 2011, the FLASH⁺ receiver was used (Klein et al. 2014). This is a dual-frequency MPIfR principal investigator receiver that operates simultaneously in the 345 GHz and the 460 GHz atmospheric windows. The 345 GHz channel is a sideband-separating detector with 4 GHz bandwidth per sideband. The 345 GHz channel was tuned to 347.33 GHz in upper sideband to observe SiO(8−7) and CO(3−2), and C¹⁷O(3−2) in lower side band. The 460 GHz channel was only available for the observing run in October 2011 and it was tuned to 461.04 GHz in lower sideband to observe CO(4−3). G19.27+0.1M1, IRAS 18316−0602, and G34.43+0.2M3 were not observed in CO(4−3).The ESO observations in SiO(8−7) and CO(3−2) were performed with the APEX-2 receiver (Vassilev et al. 2008).

After comparison of the ESO and MPS data at 345 GHz source by source for the objects observed in both programmes, observations were averaged to increase the signal-to-noise ratio. Data were converted into T_MB units assuming a forward efficiency of 0.95 for both receivers, and a beam efficiency of 0.75 for SiO(8−7), CO(3−2) and C¹⁷O(3−2), and of 0.60 for CO(4−3). The original velocity resolution of the data was 0.07 km s^-1 for SiO(8−7), CO(3−2), and C¹⁷O(3−2), and 0.1 km s^-1 for CO(4−3). The spectra were smoothed to a final velocity resolution of 2 km s^-1 in the low-frequency data, and of 1 km s^-1 for CO(4−3). The typical r.m.s. of the data at these spectral resolutions is 0.01 K at 335 and 345 GHz, and 0.1 K at 460 GHz. The observations are summarised in Table 2. Data reduction and analysis were performed with the GILDAS software³.

4. Observational results

4.1. SiO(8–7)

Figures A.1 and A.2 show all SiO(8−7) detected spectra of our sample together with the SiO(3−2) data of López-Sepulcre et al. (2011, beam⁴=19″ and with the SiO(5−4) line from Sánchez-Monge et al. (2013), when available, smoothed to the SiO(8−7) resolution. The SiO(8−7) line is detected above 4σ towards all sources except three (G25.04−0.2M1, G23.60+0.0M2, and G19.27+0.1M2, Fig. A.3). These results give us a detection rate of ~80%, lower than the 88% rate found in López-Sepulcre et al. (2011) using lower-J SiO emission, but higher than the 60% detection rate found by Klaassen et al. (2012) in a sample of IR-loud high-mass YSOs in SiO(8−7) with a lower sensitivity limit than ours. We note, however, that the detection rate for our sample would decrease to 66% with the sensitivity of their observations.

The three undetected sources are classified as dark at 8 μm and 22 μm, which suggests that they are in a very early evolutionary phase. The non-detection of SiO(8−7) in these very early evolutionary phases seems to contradict the findings of López-Sepulcre et al. (2011), who reported that the SiO luminosity decreases with the evolution of the source and is highest in the earliest stages of star formation. In Sect. 5 we investigate how the luminosity of the SiO(8−7) transition and the SiO abundance vary with time to determine the effects of excitation on the results of López-Sepulcre et al. (2011) and understand our low detection rate in IR-quiet sources.

SiO(8−7) emission is typically very broad, with an average full width at zero power (FWZP) of 45 km s^-1. Table 3 lists the measured SiO(8−7) FWZPs (and for comparison the SiO(3−2) parameters as well) and the velocity-integrated intensities over the FWZP velocity range for all detected sources. We note that the blue-shifted profile of SiO(8−7) in G24.33+0.1M1 shows a feature at 347.35 GHz (full width at half maximum of 10 km s^-1) that is neither detected in SiO(3−2), nor in the (5−4) line (Fig. A.1). According to the JPL catalogue, this feature might be acetaldehyde (CH₃CHO) (18_5,13 − 17_5,12) and (18_5,14 − 17_5,13) at 347.346 GHz and 347.349 GHz, respectively, with an upper level energy of about 214 K. Since there is no other lower energy CH₃CHO line in our spectra we cannot unambiguously identify this feature with CH₃CHO. To determine the level of possible contamination from this feature for the whole sample, we compared the integrated intensity over a velocity range of 10 km s^-1 centred at 347.35 GHz with the integrated intensity over the whole FWZP. Contamination is of the order of 7−8% of the total emission for all sources, except for G24.33+0.1M1 (33%), IRAS 18316−0602 (30%), G34.43+0.2M1 (21%), and G19.61−0.24A (20%).

4.2. CO and C¹⁷O lines

Figures A.1 and A.2 show the spectra of the CO isotopologues for the sources where SiO(8−7) was detected. Table 4 lists line-widths of the CO(4−3) and C¹⁷O(3−2) lines. CO(3−2) profiles are not analysed here because they are often severely affected by absorptions caused by the off position used for the observations. As expected, the CO profiles show very prominent high-velocity wings in all cases. Extreme sources are IRAS 18264−1152 and IRAS 18151−1208_2 with FWZP line-widths in excess of 100 km s^-1. On the other hand, the C¹⁷O(3−2) line usually traces gas associated with the dense envelope. Exceptions are IRAS 18264−1152, G20.08−0.14, G19.61−0.24A, and IRAS 19095+0930, where the C¹⁷O(3−2) transition also shows non-Gaussian wings.

The line profiles of SiO(8−7) and CO(4−3) are very similar in particular for IRAS 18264-1152, with emission at high velocities up to | ν− ν_LSR | = 63 km s^-1 in the red-shifted wing. For IRAS 18151−1208_2, the SiO(8−7) and CO(4−3) transitions have the same behaviour at high velocities, with emission up to | ν− ν_LSR | = 40−45 km s^-1. Emission from CO and SiO often arises from different components in outflows: the former traces swept-up material, the latter traces molecular gas close to the primary jet (Palau et al. 2006). However, the difference between the two species strongly depends on the considered CO transition and velocity range. Recent observations show that high-velocity CO emission may coincide with SiO even in the J = 2−1 line (Fig. 1 of Cabrit et al. 2012). In addition, the probability that CO emission traces the same component as SiO increases with J (e.g., Lee et al. 2006). Therefore, in Sect. 5 we assume that the SiO(8−7) and CO(4−3) emissions at high velocity originate from the same gas and derive the abundance of SiO by means of the SiO(8−7) to CO(4−3) line ratio at high velocities.

Fig. 2 Variation of the SiO(8−7) to SiO(3−2) luminosity ratio (bottom left), of the SiO(8−7) luminosity (bottom right), of the SiO(3−2) to SiO(2−1) luminosity ratio (top left), and of the SiO(3−2) luminosity (top right) as a function of Lbol/M. 22 μ-quiet and 22 μ-loud sources are marked as black and red filled circles. Red stars are 8 μ-loud sources.

5. Analysis

5.1. SiO luminosity

Following the same approach as López-Sepulcre et al. (2011), we investigated the dependence of the SiO(8−7) luminosity (L_SiO87 = ∫SiO(8 − 7)dv × 4πd²) on the evolutionary phase of the source described by L_bol/M. Although Faúndez et al. (2004) suggested that this ratio is mainly determined by the most luminous member of the cluster in high-mass star-forming regions, several authors proposed L_bol/M to be a measure of the evolutionary phase in the star formation process (e.g., Molinari et al. 2008; Ma et al. 2013). Since the APEX telescope has a beam size of 18′′at the frequency of SiO(8−7) and the beam size of the IRAM 30 m telescope is 19′′ at the frequency of SiO(3−2), the two datasets are similarly affected by beam dilution if the two lines are emitted by the same gas. Therefore, we can investigate how excitation evolves with time through the SiO(8−7) to SiO(3−2) integrated intensity ratio without any assumption on the size of the emitting gas. In the following, we use the values reported in Table 3, which are integrated over the whole line profile, including the blue-shifted emission. This is because possible contamination is only up to ~20−30% of the total integrated intensity for four sources, while for the others is negligible. The following results are also valid when only the red-shifted wing is considered, however.

Figure 2 shows L_SiO87 and the SiO(8−7) to SiO(3−2) integrated intensity ratio as function of L_bol/M. For comparison, we also show the SiO(3−2) luminosity and the SiO(3−2) to SiO(2−1) ratio taken from López-Sepulcre et al. (2011) and without any correction for the different beam sizes. Given the large uncertanties associated with the mass and the luminosity estimates of G19.61−0.24A (Sect. 2), we conservatively decided to exclude the source from the analysis.

In all cases, the SiO luminosities and their ratios seem to weakly increase with L_bol/M. The SiO(8−7) luminosity, for example, changes by more than a factor 10 in the sample with an average value of ~300 K km s^-1 kpc². On the other hand, the SiO(3−2) luminosity increases by only a factor 6 (with an average value of 1700 K km s^-1 kpc²). A similar comparison with the SiO line ratios is complicated by the different dilution factors that affect the SiO(3−2) to SiO(2−1) ratio. Without correcting for the different beams, the SiO(3−2) to SiO(2−1) ratio varies within a factor of less than three in our sample of sources. However, we stress that not taking into account the different beam sizes may strongly affect these findings. On the other hand, the corresponding SiO(8−7) to SiO(3−2) ratio is not affected by this problem as long as the two lines are emitted by the same gas; this ratio changes by a factor 5 in the sample.

These trends indicate an increase of excitation conditions with L_bol/M. This interpretation is strengthened by the fact that the only sources not detected in SiO(8−7) are 22 μm-dark and therefore probably are in a very early evolutionary phase (see Sect. 3). The increase of SiO(8−7) luminosity with evolution agrees with the result of Klaassen et al. (2012), who found that the integrated intensity of the SiO(8−7) line increases with time for a sample of IR-loud high-mass star-forming regions (but three sources, IRAS 18264−1152, G19.61−0.24A, and G20.08−0.14 are common to both papers). Unfortunately, our statistics is too small to conclude about this.

5.2. SiO abundance

Assuming that the SiO(8−7), and CO(4−3) lines are emitted by the same gas at high-velocities (see Sect. 4.2), we computed the SiO abundance relative to CO in the wings of the SiO(8−7) and CO(4−3) lines and then converted it into abundance relative to H₂, X_SiO = N_SiO/N_H2, assuming and abundance of CO to H₂ of 10^-4 (Lacy et al. 1994). To estimate the CO and SiO column densities, we used their integrated intensities over the same velocity range (Table 5). Only the red-shifted emission is considered. We excluded G23.60+0.0M1 from the analysis because the SiO(8−7) transition is too narrow compared with CO(4−3). We calculated the rotation temperature of SiO for each source from the intensity ratio of the SiO (8−7) and (3−2) lines in the velocity range defined on the CO(4−3) line. We checked results in the optically thin and thick limit and found that apart for two sources (G19.61−0.24A and IRAS 18264−1152) the estimated temperature does not significantly depend on the optical depth. Typical rotation temperatures, T_rot, are approximately 5 K, and 11 K for G19.61−0.24A and IRAS 18264−1152 under the assumption of optically thin emission. For the same sources, Sánchez-Monge et al. (2013) found an average temperature of 10 K based on lower J lines.

For the estimate of X_SiO, we assumed a rotation temperature of 10 K for SiO and CO. Finally, since the CO(4−3) transition observed with APEX has a slightly smaller beam than SiO(8−7), we also assumed that all the SiO(8−7) emission comes from a region as large as the APEX beam at 461 GHz and corrected the SiO(8−7) emission for dilution. Under these assumptions, the abundance of SiO varies between 1.2 × 10^-8 and 6.8 × 10^-8 with an average value of 2.6 × 10^-8. The average value of X_SiO in our sample increased to 1.8 × 10^-7 when we assume T_rot = 5 K. The derived abundances are plotted against L_bol/M in Fig. 3, where G19.61−0.24A is not considered because of the uncertainties on its mass and luminosity. Although this analysis can be performed for only a few sources, they cover the full range of L_bol/M values of our whole sample. We found no decrease of SiO abundance with L_bol/M, as reported by Sánchez-Monge et al. (2013) on a similar sample of sources and based on the SiO(2−1) and (5−4) transitions. On the other hand, Fig. 3 suggests a rather constant SiO abundance with time given the large uncertainties in the determination of X_SiO (which arise because of the LTE assumption and the selected velocity range, for example). We note that the trend seen by Sánchez-Monge et al. (2013) is mostly driven by two sources of their sample (IRAS 18151−1208_1 and IRAS 19095+0930), while the other sources show a small variation in X_SiO compatible with our results. IRAS 19095+0930 is also part of our sample and has the second smallest abundance of SiO reported in this study.

Fig. 3 Distribution of SiO abundance derived from the ratio of SiO(8−7) to CO(4−3) as function of Lbol/M.

To further test our results, we also investigated the L_SiO87/L_bol ratio as function of L_bol/M, because López-Sepulcre et al. (2011) suggested that the variation of SiO luminosity, measured through the (2−1) line, with L_bol/M might reflect a change of SiO abundance. Contrary to the result of López-Sepulcre et al. (2011), we found no correlation between L_SiO87/L_bol and L_bol/M. This is probably due to our selection criterion of broad profiles, however, which excluded weak SiO(2−1) sources; a similar behaviour was also found in the L_SiO21/L_bol for the sources we analysed (Fig. 1).

5.3. Large velocity gradient calculations

We analysed the observed SiO(8−7)/SiO(3−2) (and SiO(8−7)/SiO(5−4) when available) line ratios using the RADEX⁵ non-LTE code (van der Tak et al. 2007) with the rate coefficients for collisions with H₂ (Dayou & Balança 2006) using a plane parallel geometry. The SiO(5−4) was taken from Sánchez-Monge et al. (2013) and was smoothed over the SiO(8−7) beam to avoid different dilution effects. We explored H₂ densities from 10² to 10⁸ cm^-3, kinetic temperatures, T_kin, from 50 to 500 K, and a specific column density N(SiO)/(dV/dz) ranging from 10¹⁰ to 10¹⁶ cm^-2 (km s^-1)^-1, that is, from the fully optically thin to the optically thick regime. No correction for beam dilution was applied, which means that we implicitly assumed the same emitting size for the lines. The solutions therefore are representative of the beam-averaged emission (~18″ for both lines).

Although three transitions are not enough to properly analyse the excitation conditions, useful constraints can be still placed. For the sources that were also observed in SiO(5−4), the SiO(8−7)/SiO(3−2) and SiO(8−7)/SiO(5−4) line ratios behave similarly. As an example, Fig. 4 shows the solutions for the observed SiO(8−7)/SiO(3−2) intensity ratio towards two sources of the sample, IRAS 18264−1152 and IRAS 18151−1208_2, for typical red-, blue-shifted, and systemic emission in the n_H2–T_kin plane and for both the optically thin (N(SiO) = 10¹¹ cm^-2) and optically thick (N(SiO) = 10¹⁶ cm^-2) cases. As already learnt from previous studies of low-mass protostellar systems (e.g., Cabrit et al. 2007, and references therein), the SiO line ratios are not very sensitive to kinetic temperature. On the other hand, densities always higher than 10³ cm^-3 can be inferred. The average SiO(8−7)/SiO(3−2) line ratios measured at typical red- (0.48±0.34) and blue-shifted (0.66±0.36) outflow velocities (|V − V_sys|~ 10–15 km s^-1) are higher by a factor ~2–3 than those measured at the systemic velocity (0.15±0.10). This trend is better visible in Fig. 5, where the SiO(8−7)/SiO(3−2) line ratio is shown as function of velocity for IRAS 18151−1208_2. These results show that the excitation conditions of SiO increase with the velocity of the emitting gas, and that emission at high velocities traces a gas closer to the primary jet, as found in low-mass protostars by Nisini et al. (2007). In the red- and blue-shifted outflow velocity ranges we find a lower limit to the density of 10⁴ cm^-3, while at ambient velocities the SiO(8−7)/SiO(3−2) line ratio constrains the density to values higher than 10³ cm^-3. Interestingly, this is consistent with the findings of Nisini et al. (2007) of volume densities between 10⁵ and 10⁶ cm^-3 towards the nearby prototypical Class 0 low-mass objects L1148 and L1157 on angular scales of 10″–30″. Finally, we can obtain lower limits on column densities by using the SiO brightness temperatures and assuming the emitting size >3″, consistently with the SiO(5−4) maps by Sánchez-Monge et al. (2013), which shows outflow lobes barely resolved with a beam of 11″. In this way, we obtain N(SiO) ≥ 10¹³–10¹⁴ cm^-2, in agreement with the values derived in Sect. 5.2.

Fig. 4 Analysis of the SiO line emission in the IRAS 18264−1152 a) and IRAS 18151−1208_2 b) outflows observed with the APEX (J = 8–7, present work) and IRAM 30 m (J = 3–2, from López-Sepulcre et al. 2011) antennas. Both lines have been observed with a well-matching beam (~18″): the solutions for the observed SiO(8−7)/SiO(3−2) intensity ratio are shown in the nH2–Tkin-plane for non-LTE (RADEX) plane-parallel models for both the optically thin (N(SiO) = 1011 cm-2) and optically thick (N(SiO) = 1016 cm-2) cases. Contours are for the measured ratios at the typical red- (+54 km s-1, for IRAS 18264−1152; +50 km s-1, for IRAS 18151−1208_2) and blue-shifted outflow velocities (+15 km s-1 for IRAS 18151−1208_2), as well as at the systemic velocity (IRAS 18264−1152: +43.9 km s-1, IRAS 18151−1208_2: +30 km s-1). Dashed contours are for uncertainties.

Fig. 5 Upper panel: SiO(3−2) (red) and (8−7) (black) spectra of IRAS 18151-1208_2. Lower panel: SiO(8−7)/SiO(3−2) line ratio as function of velocity (red line). The black line is the SiO(8−7) spectrum. The dashed black line marks the zero TMB level of the SiO(8−7) spectra. The dashed red lines label the velocity range used in Sect. 5.4 for the shock analysis and reported in Table 6.

5.4. Shock analysis

We compared the SiO data with state-of-the-art one-dimensional models of shock propagation in the interstellar medium. We extracted the observed integrated intensity of the shocked gas for both SiO (3−2) and (8−7) lines, and we computed and used the (8−7)/(3−2) integrated intensity ratio. Owing to contamination (Sect. 4.1) of the blue lobe in the SiO (8−7) line for several objects in our sample, we chose to extract these three quantities only from the red-shifted shocked component. Table 6 provides the details relevant to this operation for each object, including the velocity range on which the ratio was extracted. The velocity ranges used here differ from those discussed in Sect. 5.2 because CO(4−3) is more contaminated by ambient gas emission than SiO at low-velocities.

Figure 6 compares the observed ratios with the values for three classes of 1D stationary shock models: the observations are displayed in colours, while the model points are indicated in shades of grey. We adopted two different x-axes: for the models, the x-axis is the shock velocity, for the observations it is the width of the velocity range used for the integration. When trying to fit observations with 1D models, these two quantities should be of the same order. We subsequently chose to restrict the analysis to observations for which the velocity range used for the observations could be compared with the shock velocity of our models, namely between 20 and 50 km s^-1. The three model classes shown in Fig. 6 are:

• class 1 models from Gusdorf et al. (2008a): pre-shock densities of 10⁴, 10⁵, 10⁶ cm^-3, magnetic field parameter⁶b = 1, and Si-bearing material only distributed in the grain cores in the pre-shock phase;

• class 2 models from Anderl et al. (2013): they include grain-grain interactions. Three different shock velocities (20, 30 and 40 km s^-1) are covered, with three different values of the magnetic field for each velocity (b = 1.0,1.5,2.0 for 20 km s${-1}^{}$, b = 1.5,2.0,2.5 for 30 km s${-1}^{}$, and b = 2.0,2.5,3.0 for 40 km s${-1}^{}$) at a pre-shock density of 10⁵cm^-3. Compared with models that do not incorporate grain-grain processing, higher values of the magnetic field parameter b are required for continuous magneto-hydrodynamic shocks to propagate at higher velocities because of the different coupling of small charged grains to the magnetic field;

• class 3 models are of unprecedented kind: they are similar to class 2 models, but additionally, 10% of the elemental Si is assumed to be in the form of SiO in the grain mantles in the pre-shock phase. One of the important conclusions in Anderl et al. (2013) is that these models might be the only way to fit SiO observations of outflows from low-mass YSOs. For this class, the range of covered input parameters is the same as for class 2.

In a first step, we compared the observations with the class 1 models, in which SiO cannot be efficiently produced with shock velocities lower than 25 km s^-1 because this is the threshold for core erosion to efficiently release Si in the gas phase. These models unambiguously indicate that only a pre-shock density greater than 10⁴ cm^-3 could yield solutions compatible with the observed ratios. At such high densities, however, the effects of grain-grain interactions (mostly shattering and vaporisation, Guillet et al. 2011) must be taken into account (Anderl et al. 2013). These effects are not considered in the class 1 models.

In class 2 and 3 models, SiO can be released into the gas phase by gas-grain interactions, leading to the sputtering of the grain cores (as in Gusdorf et al. 2008a), but also through grain-grain collisions leading to the vaporisation and shattering of the grains. The latter effect results in a production of small grain fragments in large numbers, which increases the total dust grain surface area and thereby changes the coupling between the neutral and the charged fluids within the shocked layer. The corresponding shocks become much hotter and subsequently thinner, because the kinetic energy of the shock is transformed into heat within a shorter time. In addition to this release of SiO by vaporisation, this strong effect of shattering on the structure of magneto-hydrodynamic shocks is the second reason why it is necessary to include grain-grain processing in molecular shocks at densities higher than ~10⁵cm^-3. In class 3 models, part of the SiO can additionally originate from the sputtering of the grain mantles. Two behaviours can be identified:

• all class 2 models can provide an overall decent fit to most observations in terms of integrated intensity ratios (except for G19.61−0.24A, which might be associated to a higher pre-shock density). However, for a shock velocity equal to 20 km s^-1, or for the shock at 30 km s^-1 and b = 2.5, the SiO absolute integrated intensities are simply too low to match the observed values assuming an emitting size >3′′(see the end of Sect. 5.3). This is because at low shock velocities, the grain-grain interactions are not efficient, and because the combination of a moderate (30 km s^-1) shock velocity and high magnetic field strength results in a critically narrow SiO emission layer;

• the class 3 models do provide convincing fits in terms of absolute integrated intensities, but they seem to generally overestimate the integrated intensity ratios as displayed in Fig. 6 (except for G19.61-0.24A and G23.60+0.0M1). This is because we chose to put 10% of SiO in the grain mantles, similar to Gusdorf et al. (2008b). Figure 6 and the absolute integrated intensity criterion suggest that models with fewer SiO in the grain mantles would provide the best fit.

Fig. 6 Comparisons of observations and models of SiO (8−7) to (3−2) integrated intensity ratios, extracted over the red-shifted wing of the observed line profiles. For observations, the ratios are displayed against the width of the velocity range over which they are extracted. For models, the ratios are plotted as a function of the shock velocity of the model. In the context of 1D shock modelling, the width of the velocity range used for observations should be directly comparable with the shock velocity of the models for each observed outflow. The observations are shown in coloured squares for various objects in our sample. The empty circles are the models of Gusdorf et al. (2008a) for pre-shock densities of 104, 105, 106 cm-3, the filled diamonds are the models (nH = 105 cm-3) of Anderl et al. (2013) without SiO initially in the grain mantles, and the filled flat hexagons are the same, including an initial amount of SiO in the grain mantles.

6. Discussion

Under the assumption that the SiO emission arises from gas closely associated with the primary jet, as found in molecular jets from low-mass YSOs, our observations confirm the presence of jets up to L_bol/M values at least of 50 L_⊙/M_⊙. The results presented in Sect. 5.1 suggest that excitation conditions of the gas responsible for the SiO emission increase with time. The trend of SiO luminosity and luminosity ratios as a function of L_bol/M are weak and similar observations on a larger sample of sources need to be performed to confirm our findings. However, these trends might explain the non-detection of SiO(8−7) in the youngest sources of our sample. Large velocity gradient calculations showed that the SiO(8−7) to SiO(3−2) line ratio depends more strongly on density than on temperature. This suggests that the increase of SiO(8−7) to SiO(3−2) seen in our sources reflects an increase of density of the gas in the jet. We speculate that this increase in jet density might be linked to an increased mass-accretion rate at very early stages of star formation (as noted for HH212 by Cabrit et al. 2007). In this case, the increase of density seen in our sample of sources might be interpreted as an increase of the accretion rate during the evolution of a massive YSO or, alternatively, as an increase of the accretion rate with the mass of the central object. Indeed, an increase of L_bol/M is always linked to an increase of mass of the central object when accretion is still significant, when L_bol/M is either interpreted as a time estimator, or when it mostly reflects the mass of the most massive object in the cluster. Recently, Kuiper & Yorke (2013) modelled the simultaneous evolution of a massive protostar and its host molecular core. Their calculations showed that the accretion rate increases with time and with the mass of the protostellar mass in the very early stages of evolution until the beginning of the Kelvin-Helmholtz contraction. A note of caution here is that effects of radiative pumping processes are not considered in our analysis, while Godard & Cernicharo (2013) showed that SiO J ≥ 3 lines are sensitive to optical/UV radiation fields. This effect becomes more relevant with time, especially for high J transitions, as a consequence of the stronger radiation field from the central object. However, this mechanism probably mostly affects gas close to the YSO due to extinction, unless the SiO emission comes from the walls of the outflow cavities that allow the radiation to reach longer distances from the protostar.

Alternatively, when no pre-shock SiO is placed in the grain mantles, and all other shock model parameters are equal, Gusdorf et al. (2008b) predicted that the integrated intensity of SiO transition tends to increase with the age of the shock in outflows from low-mass Class 0 objects, owing to the greater velocity extent of the emitting region. To verify this, we looked for correlations between the luminosity of the SiO(8−7), (3−2) and (2−1) transitions and their line width (Fig. A.4). The luminosity of the SiO(8−7) increases for line widths between 30 and 60 km s^-1, although two objects (IRAS 1824−152 and G23.60+0.0M1 ) do not follow the general trend. The SiO(3−2) luminosity (and also that of the SiO(2−1) line) neither shows any strong correlation with its corresponding line width, nor is there an obvious correlation between the line width of the SiO(3−2) and (8−7) transitions and L_bol/M because the line widths are relatively constant in our sample (see Table 3), which probably is a consequence of our selection criteria. On the other hand, the SiO(2−1) line width seems to decrease with increasing values of L_bol/M. We note, however, that the observed line width strongly depends on the inclination angle of the outflow with respect to the sky, which is unknown for the sources of our sample.

The second result of our study is that we do not find any significant variation of SiO abundance in our sources. This is estimated through comparisons with CO column densities obtained by mean of the CO(4−3) integrated intensity when available, but also through L_SiO87/L_bol for all the sources in the sample.This result contradicts the findings of other authors (Miettinen et al. 2006; López-Sepulcre et al. 2011; Sánchez-Monge et al. 2013), who reported a decrease of SiO abundance with time. However, there are two differences between our and the previous studies. First, previous results from the literature were obtained through observations of the SiO(2−1) and (3−2) lines, which are not sensitive to changes of excitations. Second, our study does not sample the most evolved sources where the SiO abundance is expected to drop more dramatically because of a decay of jet activity. To investigate the variation of SiO abundance as function of time and of excitation in more detail, we also need observations of high J SiO transitions towards more evolved objects than those presented here. Moreover, a more accurate analysis of the SiO emission requires high angular resolution observations to resolve the emitting gas and verify that the SiO emission is dominated by one object and does not arise from a population of lower mass YSOs.

Finally, our analysis (Fig. 6) shows that shocks with a rather high pre-shock density (n_H = 10⁵ cm^-3) provide the best fit to the observations. The analysis shows that grain-grain interactions must be taken into account, and that part of the SiO emission arises from the sputtering of a fraction of silicon-bearing material (corresponding to less than 10% of the total silicon abundance) from the grain mantles. In this sense, we corroborate the findings of Anderl et al. (2013). Silicon-bearing material in the mantles of grains around massive YSOs might be due to previous shock processing of the regions, either owing to previous episodes of ejection in the targeted regions (as for the B1 knot of the L1157 outflow, invoked by Gusdorf et al. 2008b), or to converging flows that might lead to the formation of massive stars in filamentary structures (‘ridges’, like in the W43 region, e.g., Nguyen-Luong et al. 2013). Unfortunately, this shock modelling cannot be used for evolutionary trend studies because the shocks are stationary. However, dense shocks can account for the observations of SiO in massive star-forming regions, as was also reported by Leurini et al. (2013).

7. Conclusions

We presented SiO(8−7) and CO(4−3) APEX observations of a sample of massive clumps in different evolutionary phases of star formation. The data were complemented with SiO(3−2) (and SiO(5−4) when available) observations from the literature with a similar angular resolution. The main results can be summarised as follows:

• 1.

  We detected SiO(8−7) emission in all sources of our sample except in three objects that are probably in very early phases of star formation. The line-profiles are very broad, with an average FWZP of 45 km s^-1.

• 2.

  The SiO(8−7) luminosity and the SiO(8−7) to SiO(3−2) line ratio increase weakly with the evolutionary phase of the source described by L_bol/M. This might be explained with an increase of density in the gas traced by SiO.

• 3.

  We analysed the SiO(8−7) to SiO(3−2) line ratio with a non-LTE radiative transfer code. At high outflow velocities we found a lower limit to the density of 10⁴ cm^-3, in agreement with measurements towards nearby prototypical Class 0 low-mass objects. We also found that the excitation conditions of SiO increase with the velocity of the emitting gas. This suggests that emission at high velocities traces a gas close to the primary jet, as found in low-mass protostars.

• 4.

  We estimated the abundance of SiO at high velocities through comparison with CO(4−3) emission in the same velocity range as SiO. Although the estimates are affected by large uncertainties, the SiO abundance is realtively constant, with an average value of 2.6 × 10^-8 and no clear change of X_SiO versus L_bol/M.

• 5.

  Shock modelling reveals high pre-shock densities (n ≥ 10⁵cm^-3). At these densities the effects of grain-grain collisions cannot be ignored anymore. For the first time, the models of Anderl et al. (2013) of magnetohydrodynamical shocks including the effects of shattering and vaporisation were used to interpret observations. The observations are globally compatible with these models at a pre-shock density of 10⁵cm^-3, in which sputtering of silicon-bearing material (corresponding to less than 10% of the total silicon abundance) from the grain mantles occurs.

Online material

Appendix A: Figures

Fig. A.1 Spectra of CO(4−3) and (3−2), C17O(3−2), SiO(8−7) (multiplied by a factor of five for clarity), SiO(5−4) (blue, where available) and SiO(3−2) (red). The CH3CHO features at 347.35 GHz (Sect. 3) are labelled.

Fig. A.2 Spectra of CO(3−2), C17O(3−2), SiO(8−7) (multiplied by a factor of five for clarity), SiO(5−4) (blue, where available) and SiO(3−2) (red). The CH3CHO features at 347.35 GHz (Sect. 3) are labelled.

Fig. A.3 Spectra of CO(4−3) and (3−2), C17O(3−2), SiO(8−7) (black), SiO(5−4) (blue, where available) and SiO(3−2) (red) for the three sources not detected in SiO(8−7) (see Table 3).

Fig. A.4 Luminosities of SiO(8−7) a), SiO(3−2) b) and SiO(2−1) c) as a function of their corresponding line width. 22 μ-quiet and 22 μ-loud sources are marked as black and red filled circles. Red stars are 8 μ-loud sources.

Acknowledgments

The authors thank the anonymous referee for helpful comments that improved the clarity and content of the paper.

References

All Tables

All Figures

	Fig. 1 SiO(2−1) luminosity against Lbol/M for the full sample of sources studied by López-Sepulcre et al. (2011). Full circles depict the sources observed in SiO(8−7) that we studied here.
In the text

	Fig. 2 Variation of the SiO(8−7) to SiO(3−2) luminosity ratio (bottom left), of the SiO(8−7) luminosity (bottom right), of the SiO(3−2) to SiO(2−1) luminosity ratio (top left), and of the SiO(3−2) luminosity (top right) as a function of Lbol/M. 22 μ-quiet and 22 μ-loud sources are marked as black and red filled circles. Red stars are 8 μ-loud sources.
In the text

	Fig. 3 Distribution of SiO abundance derived from the ratio of SiO(8−7) to CO(4−3) as function of Lbol/M.
In the text

Fig. 4 Analysis of the SiO line emission in the IRAS 18264−1152 a) and IRAS 18151−1208_2 b) outflows observed with the APEX (J = 8–7, present work) and IRAM 30 m (J = 3–2, from López-Sepulcre et al. 2011) antennas. Both lines have been observed with a well-matching beam (~18″): the solutions for the observed SiO(8−7)/SiO(3−2) intensity ratio are shown in the nH2–Tkin-plane for non-LTE (RADEX) plane-parallel models for both the optically thin (N(SiO) = 1011 cm-2) and optically thick (N(SiO) = 1016 cm-2) cases. Contours are for the measured ratios at the typical red- (+54 km s-1, for IRAS 18264−1152; +50 km s-1, for IRAS 18151−1208_2) and blue-shifted outflow velocities (+15 km s-1 for IRAS 18151−1208_2), as well as at the systemic velocity (IRAS 18264−1152: +43.9 km s-1, IRAS 18151−1208_2: +30 km s-1). Dashed contours are for uncertainties.

In the text

	Fig. 5 Upper panel: SiO(3−2) (red) and (8−7) (black) spectra of IRAS 18151-1208_2. Lower panel: SiO(8−7)/SiO(3−2) line ratio as function of velocity (red line). The black line is the SiO(8−7) spectrum. The dashed black line marks the zero TMB level of the SiO(8−7) spectra. The dashed red lines label the velocity range used in Sect. 5.4 for the shock analysis and reported in Table 6.
In the text

Fig. 6 Comparisons of observations and models of SiO (8−7) to (3−2) integrated intensity ratios, extracted over the red-shifted wing of the observed line profiles. For observations, the ratios are displayed against the width of the velocity range over which they are extracted. For models, the ratios are plotted as a function of the shock velocity of the model. In the context of 1D shock modelling, the width of the velocity range used for observations should be directly comparable with the shock velocity of the models for each observed outflow. The observations are shown in coloured squares for various objects in our sample. The empty circles are the models of Gusdorf et al. (2008a) for pre-shock densities of 104, 105, 106 cm-3, the filled diamonds are the models (nH = 105 cm-3) of Anderl et al. (2013) without SiO initially in the grain mantles, and the filled flat hexagons are the same, including an initial amount of SiO in the grain mantles.

In the text

	Fig. A.1 Spectra of CO(4−3) and (3−2), C17O(3−2), SiO(8−7) (multiplied by a factor of five for clarity), SiO(5−4) (blue, where available) and SiO(3−2) (red). The CH3CHO features at 347.35 GHz (Sect. 3) are labelled.
In the text

	Fig. A.2 Spectra of CO(3−2), C17O(3−2), SiO(8−7) (multiplied by a factor of five for clarity), SiO(5−4) (blue, where available) and SiO(3−2) (red). The CH3CHO features at 347.35 GHz (Sect. 3) are labelled.
In the text

	Fig. A.3 Spectra of CO(4−3) and (3−2), C17O(3−2), SiO(8−7) (black), SiO(5−4) (blue, where available) and SiO(3−2) (red) for the three sources not detected in SiO(8−7) (see Table 3).
In the text

	Fig. A.4 Luminosities of SiO(8−7) a), SiO(3−2) b) and SiO(2−1) c) as a function of their corresponding line width. 22 μ-quiet and 22 μ-loud sources are marked as black and red filled circles. Red stars are 8 μ-loud sources.
In the text