F. F. S. van der Tak ¹ - A. M. S. Boonman ² - R. Braakman ² - E. F. van Dishoeck²

1 - Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany
2 - Sterrewacht, Postbus 9513, 2300 RA Leiden, The Netherlands

Abstract
The sulphur chemistry in nine regions in the earliest stages of high-mass star formation is studied through single-dish submillimeter spectroscopy. The line profiles indicate that 10-50% of the SO and SO₂ emission arises in high-velocity gas, either infalling or outflowing. For the low-velocity gas, excitation temperatures are 25 K for H₂S, 50 K for SO, H₂CS, NS and HCS⁺, and 100 K for OCS and SO₂, indicating that most observed emission traces the outer parts (T<100 K) of the molecular envelopes, except high-excitation OCS and SO₂ lines. Abundances in the outer envelopes, calculated with a Monte Carlo program, using the physical structures of the sources derived from previous submillimeter continuum and CS line data, are $\sim$ 10^-8 for OCS, $\sim$ 10^-9 for H₂S, H₂CS, SO and SO₂, and $\sim$ 10^-10 for HCS⁺ and NS. In the inner envelopes (T>100 K) of six sources, the SO₂ abundance is enhanced by a factor of $\sim$ 100-1000. This region of hot, abundant SO₂ has been seen before in infrared absorption, and must be small, $\la$ 0 $.\!\!^{\prime\prime}$ 2 (180 AU radius). The derived abundance profiles are consistent with models of envelope chemistry which invoke ice evaporation at $T\sim 100$ K. Shock chemistry is unlikely to contribute. A major sulphur carrier in the ices is probably OCS, not H₂S as most models assume. The source-to-source abundance variations of most molecules by factors of $\sim$ 10 do not correlate with previous systematic tracers of envelope heating. Without observations of H₂S and SO lines probing warm ( $\ga$ 100 K) gas, sulphur-bearing molecules cannot be used as evolutionary tracers during star formation.

Key words: ISM: molecules - molecular processes - stars: circumstellar matter - stars: formation

1 Introduction

Spectral line surveys at submillimeter wavelengths have revealed considerable chemical differences between star-forming regions (e.g., Blake et al. 1987; Nummelin et al. 2000; Helmich & van Dishoeck 1997; Schilke et al. 1997; Hatchell et al. 1998a; Sutton et al. 1995; see van Dishoeck 2001 for a complete list). While these differences indicate activity, the dependence of molecular abundances on evolutionary state and physical parameters is poorly understood. Better insight into this relation would be valuable for probing the earliest, deeply embedded phases of star formation where diagnostics at optical and near-infrared wavelengths are unavailable. This is especially true for the formation of high-mass stars, for which the order in which phenomena occur is much less well understood than in the low-mass case, and for which the embedded phase is a significant fraction ( $\approx$ 10%) of the total lifetime of $\sim$ 10⁶ yr.

Sulphur-bearing molecules are attractive as candidate tracers of early protostellar evolution (Buckle & Fuller 2003; Hatchell et al. 1998b). In models of "hot cores'' (Charnley 1997), where the chemistry is driven by the evaporation of icy grain mantles due to protostellar heating, the abundances of sulphur-bearing molecules exhibit a characteristic variation with time. After its release from the grain mantles, H₂S is converted into SO and SO₂ in $\sim$ 10³ yr, and further into CS, H₂CS and OCS after $\sim$ 10⁵ yr. This type of model has had some success in reproducing the chemistry of star-forming regions (Langer et al. 2000).

Besides thermal heating, shocks are potentially relevant to the chemistry in star-forming regions. Young stars drive outflows which expose their surroundings to the effects of shocks. The submillimeter lines of SO₂ in Orion are very broad, and likely arise in the outflow (Schilke et al. 1997). Models of chemistry behind interstellar shocks (Pineau des Forêts et al. 1993; Mitchell 1984; Leen & Graff 1988) predict enhancements of H₂S, SO and SO₂ on time scales of $\sim$ 10⁴ yr. Thus, both by ice evaporation and by shock interaction, sulphur could act as a clock on the time scale of 10⁴ yr relevant for the embedded phase of star formation.

$\begin{figure} \par\resizebox{16.2cm}{!}{\includegraphics[angle=-90,clip]{ms3386.f1}} \end{figure}$

Figure 1: Two examples of spectral settings toward eight of our sources, showing differences in the relative and in the absolute strengths of lines of sulphur-bearing molecules. The top and bottom frequency scales are the two receiver sidebands. The data for NGC 6334 IRS1 have been divided by 10 and 3 respectively, for clarity. The non-labeled lines are due to non-sulphur-bearing molecules, of which the strongest are labeled in the AFGL 2591 panels.

Open with DEXTER

This paper uses submillimeter observations of sulphur-bearing molecules to test models of chemistry during star formation, and attempts to order these regions chronologically. The sources are nine regions of high-mass star formation with luminosities $1\times 10^4 {-} 2\times 10^5$ $L_{\odot}$ at distances of 1-4 kpc (Table 1). Originally selected for their mid-infrared brightness, the sources were found to have large envelope masses (40-1100 $M_{\odot}$ within r=0.15-0.36 pc) based on submillimeter data (van der Tak et al. 2000b). Together with their weak radio emission and their strong molecular outflows, these masses indicate very early evolutionary stages, probably preceding the formation of hot cores. Two of our sources may be more evolved than the others: NGC 6334 IRS1, which has a very rich, "hot core''-type submillimeter spectrum (McCutcheon et al. 2000; Thorwirth et al., in prep.), and NGC 7538 IRS1, which has strong free-free emission.

Observations of CH₃OH lines towards some of these sources (van der Tak et al. 2000a) indicate an increase in the CH₃OH abundance by a factor of 100 in the warm inner envelope, which is likely due to evaporation of CH₃OH-rich ices. The excitation of CH₃OH, as well as that of CO and C₂H₂, correlates well with other temperature tracers, such as the gas/solid ratios of CO₂ and H₂O, the 45/100 $\mu$ m colour and the fraction of heated solid ¹³CO₂ (Boogert et al. 2000; Boonman et al. 2003; van der Tak et al. 2000b). Since these quantities also correlate with the ratio of envelope mass to stellar mass, their variation likely reflects evolution through the dispersal of the protostellar envelopes. These findings make the sample a good testbed for theories of chemical evolution during the earliest deeply embedded phases of high-mass star formation.

2 Observations

The data presented in this paper are part of a targeted spectral line survey of embedded massive stars, aimed to track chemical evolution during the earliest stages of high-mass star formation. Fifteen frequency settings have been observed, which also contain many lines of sulphurless molecules. Some of these data have been published by van der Tak et al. (1999,2000b,a). These observations were performed with the 15-m James Clerk Maxwell Telescope (JCMT) on Mauna Kea, Hawaii between 1995 and 1998. The beam size (FWHM) and main beam efficiency of the antenna were 18'' and 64-69% at 230 GHz and 14'' and 58-63% at 345 GHz. The frontends were the receivers A2, B3 and B3i; the backend was the Digital Autocorrelation Spectrometer, covering 500 MHz instantaneous bandwidth. Pointing was checked every 2 hours during the observing and was always found to be within 5''. To subtract the atmospheric and instrumental background, the chopping secondary was used with 180''offsets. Total integration times are 30-40 min for each frequency setting. Figure 1 shows examples of the data towards all but one of our sources; the data on W3 IRS5 were presented by Helmich & van Dishoeck (1997).

Additional observations of H₂S and OCS in the 130-180 GHz window were carried out with the 30-m telescope of the Institut de Radio Astronomie Millimétrique (IRAM) on Pico Veleta, Spain, in August 2002. The frontend was the facility receiver C150 and the backend the Versatile Spectral Assembly (VESPA) autocorrelator. Integration times were 10-20 min using frequency switching with a throw of 7.3 MHz. Data were calibrated onto $T_{\rm MB}$ scale by multiplying by 1.43, the ratio of forward and main beam efficiencies. The beam size is 15'' at these frequencies.

The H₂S 393 GHz line was observed in May 1995 at the 10.4-m Caltech Submillimeter Observatory (CSO), with the facility receiver as frontend and the 50 MHz AOS as backend. The telescope has a beam size of 25'' and a main beam efficiency of 60% at this frequency.

Reduction was carried out using the CLASS package developed at IRAM. Linear baselines were subtracted and the spectra were smoothed once and calibrated onto $T_{\rm MB}$ scale. Only the frequency switched 30 m data required higher order polynomials to fit the baseline. The final spectra have a resolution of 0.3-1.5 km s^-1 and rms noise levels (in $T_{\rm MB}$ ) of 20-30 mK for the 150 and 230 GHz bands, and 30-50 mK for the 345 GHz band. Multiplying these values by the line width (Table 4) yields the uncertainties of the line fluxes in Table 3. Although the absolute calibration is only correct to $\approx$ 30%, the relative strength of lines within one frequency setting is much more accurate.

Line identification was performed using the JPL catalog (Pickett et al. 1998). A matching frequency is the main criterion for assignment, but the upper energy level of the transition and the complexity of the molecule are also considered. Table 2 has details about the observed transitions.

For each detected feature, a satisfying match could be found in the catalog within 0.5 MHz, except for the H₂CS 10₁₉-9₁₈ line at 348.5 GHz, whose residual is 2 MHz. Laboratory spectroscopy of H₂CS only exists up to 250 GHz (Beers et al. 1972), and the analysis in the JPL catalog only includes pure rotation and centrifugal distortion terms. At higher frequencies and for higher J-values, higher order terms become important, which may be why the catalog prediction is inaccurate for this line. Our data indicate a frequency of 348 534.2 MHz.

3 Results

3.1 Line profiles

We have fitted the detected lines with Gaussian profiles. Table 3 reports the observed fluxes and Table 4 the widths of the lines. While most lines could be satisfactorily fitted with a single Gaussian, the lines of SO, which are the strongest lines that we have observed, required two Gaussians for a good fit to their profiles (Fig. 2; Table 4). One component agrees in position and width with the values measured for C¹⁷O and C³⁴S by van der Tak et al. (2000b) and is attributed to a static envelope. Table 5 reports similar measurements for the source MonR2 IRS3 which van der Tak et al. did not study. The other component has 2.8 times as large a width on average. This high-velocity emission may trace outflowing or infalling motions. The contributions from high- and low-velocity gas to the SO line fluxes are about equal, unlike CS, where the envelope carries 90% of the fluxes (van der Tak et al. 2000b; Table 5). The infrared absorption and submillimeter emission lines of CO towards these sources also show multiple velocity components (Mitchell et al. 1992; Mitchell & Hasegawa 1991; Mitchell et al. 1991). However, the velocities do not match in detail, as expected because of differences in spectral and angular resolution between infrared and submillimeter data on one hand, and excitation differences between CO and SO on the other.

The fit of two Gaussians to the SO lines is excellent, which is not always the case for CO. An alternative method is to use area integrals, which do not attribute any low-velocity gas to the high-velocity component. Our approach may overestimate the high-velocity contribution to the line flux by factors up to 2.

The ratios of the high- to low-velocity contributions to the SO lines are $\la$ 1 for the lines in the 230 GHz band, and $\ga$ 1 for the lines in the 345 GHz band. This difference could be due to the higher temperatures and densities required for excitation of the 345 GHz lines, or to the smaller beams in which they are observed. Since the high-velocity wings on the CS lines in these sources also become more pronounced in smaller beams, independent of transition (van der Tak et al. 2000b), the broad component probably is more compact than the envelope, but not necessarily warmer or denser. The high-velocity gas in these sources thus seems to be confined to small radii. This situation differs from that in molecular clouds, where line width increases with size, and from that in many molecular outflows, where velocities usually increase with radius (e.g., Lee et al. 2002). However, while CO data show arcminute-sized outflows in our sources, higher velocities naturally occur at smaller radii in the case of infall. For example, gas at 1000 AU (1'' at 1 kpc) in free fall onto a 10 $M_{\odot}$ object would have $V_{\rm inf}=4$ km s^-1, consistent with the observed width of the broad component (Table 4).

The broad component peaks at slightly more blueshifted velocities (by 0.17 km s^-1 on average) than the envelope, a situation found before for CS and H₂CO (van der Tak et al. 2000b). In the case of outflow, obscuration by A_V=10⁴ mag of dust could explain the absence of redshifted wings. Since some redshifted SO is seen (Fig. 2), the requirement relaxes to several 1000 mag. Of the sources discussed here, only NGC 6334 IRS1 may have such a large amount of dust in its envelope (Table 1); the other sources may need circumstellar disks seen face-on. In the case of infall, a mixture of absorption and emission is expected at redshifted velocities (Choi 2002), which may partly cancel each other if unresolved, to give the impression of predominantly blueshifted emission. We conclude that both infall and outflow are plausible explanations for the high-velocity gas, and recommend interferometric observations of SO to decide between these options.

$\begin{figure} \par\resizebox{8.8cm}{!}{\includegraphics[angle=-90,clip]{ms3386.f2}} \end{figure}$	Figure 2: Observed velocity profiles of the SO 6₆-5₅ and 8₇-7₆ lines in the sources AFGL 2591 and W3 IRS5, showing the low- and high-velocity components. The dashed lines are our fits to the individual components and the dotted line is the sum of the two components.
Open with DEXTER

For most sources in our sample, the width of the SO₂ lines is between the values for the low- and high-velocity components. The only source where the SO₂ lines are strong enough to fit two components is W3 IRS5, where the line profiles do not show wings. However, most of these data, taken from Helmich & van Dishoeck (1997), were observed with an AOS as backend, and have lower spectral resolution than our data. It seems therefore plausible that, like for SO, about half of the SO₂ emission arises in high-velocity gas.

Table 3: Observed line fluxes $\int T_{\rm{MB}} {\rm d}V$ (K km s^-1) from Gaussian fits, and 1 $\sigma$ upper limits^a.

Table 4: Widths $\Delta V$ (km s^-1) of the emission line components, averaged over transition. Numbers in brackets denote uncertainties in units of the last decimal, except when only one transition was observed. Dots denote missing data. For SO, values for the narrow and the broad components are reported, whose positions are also given.

3.2 Excitation temperatures

For SO and SO₂, the number of detected lines is large enough to construct rotation diagrams. This method, described in detail by Blake et al. (1987) and Helmich et al. (1994), assumes that the lines are optically thin and that the molecular excitation can be described by a single temperature, the rotation temperature. If radiative decay competes with collisional excitation, the rotation temperature lies below the kinetic temperature. If the lines are optically thick, the rotation diagram underestimates the column density and overestimates the excitation temperature. Another complication is that the lines are measured in beams of unequal size, so that somewhat different volumes of gas are probed. Despite these caveats, the rotation diagram provides a useful first estimate of excitation conditions and molecular column densities, and a stepping stone towards a more sophisticated analysis (Sect. 4).

Table 6: Rotation temperatures (K). Numbers in brackets denote uncertainties.

Table 6 presents the results. In the case of SO, separate fits were made to the low- and high-velocity components. For NGC 7538 IRS9, the detected SO₂ lines do not cover a large enough range of energy levels to constrain $T_{\rm rot}$ . The derived temperatures are 30...70 K for SO and 50...190 K for SO₂, reflecting the available range of energy levels. Lower limits to $T_{\rm rot}$ are poor fits to the data, probably caused by optically thick lines. Indeed, the SO/³⁴SO and SO₂/³⁴SO₂ line ratios of 1.3... $~\ga$ 10 are below the isotopic abundance ratio of ³²S/ $^{34}\rm S=22$ , which indicates optical depths of a few in the SO and SO₂ lines.

The last column of Table 6 lists rotation temperatures of CH₃OH, measured previously. Strictly, these numbers are only lower limits to the kinetic temperatures in a 15'' beam, but they show good correlation with $T_{\rm ex}$ (C₂H₂) which does directly trace the kinetic temperature (van der Tak et al. 2000a). The rotation temperatures of SO and SO₂ are not clearly correlated with $T_{\rm rot}$ (CH₃OH) within the uncertainties (Fig. 3). Therefore we cannot conclude at this point whether these molecules trace the same gas within the envelopes. In any case, differences in beam sizes and optical depth effects in the lines make it hard to draw firm conclusions before a more detailed radiative transfer analysis is presented in Sect. 4.

For the other molecules, for which only few lines are available, the line ratios have been used to constrain $T_{\rm ex}$ . Source-averaged line ratios indicate $T_{\rm ex}$ $~~\approx~~$ 50 K for H₂CS, NS and HCS⁺, $T_{\rm ex}$ $~~\approx~~$ 25 K for H₂S and $T_{\rm ex}$ $~~\approx~~$ 100 K for OCS. These values are only rough estimates, especially if based on few sources and/or lines.

3.3 Column densities

Table 7 presents beam-averaged molecular column densities, estimated from the observed line fluxes using the excitation temperatures found above. The column densities for H₂CS and H₂S are the sums of the ortho and para species, assuming an ortho to para ratio of 3. After calibration, the adopted $T_{\rm ex}$ is the main source of uncertainty in the column densities, especially if only one or two lines have been observed. For example, decreasing $T_{\rm ex}$ from 50 to 25 K increases column density estimates by factors of $\approx$ 2, while increasing $T_{\rm ex}$ from 50 to 75 K decreases column density estimates by $\la$ 20%. In addition, column densities may have been underestimated by factors of a few due to nonnegligible optical depth.

To search for trends of column densities with temperature, Table 7 lists the sources in order of increasing T(CH₃OH). The only trend seems to be that the warmest source has the largest column density of all molecules, except SO and SO₂. In order to investigate whether these trends in the column density reflect chemical differences between sources or are due to optical depth or beam size effects, a more detailed radiative transfer analysis is presented in Sect. 4.

3.4 Comparison with infrared data

Keane et al. (2001) observed the $\nu_3$ band of SO₂ around 7.3 $\mu$ m in absorption towards five of the sources studied here, and found column densities of SO₂ of a few 10¹⁶ cm^-2, a factor of $\sim$ 100 higher than the submillimeter data indicate. In contrast, infrared observations of CO and dust yield 3-5 times lower column densities than submillimeter data, as a result of non-spherical geometry on scales ${\la} 1''$ (van der Tak et al. 2000b). Correcting for optical depth in the SO₂ submillimeter lines increases the column density estimate, but only by factors of a few (Sect. 3.2), which does not explain the discrepancy with the infrared values. More likely to be important are chemical effects on small scales. The sources are centrally concentrated, n(H₂) $\propto r^{-\alpha}$ with $\alpha=1.0{-}1.5$ , so that absorption data are more sensitive to warm, dense gas at small radii, while emission data probe more extended, cooler and less dense gas. The infrared excitation temperatures of 225-750 K also suggest that the SO₂ absorption arises in warm gas.

The infrared estimates of N(SO₂) imply optically thick submillimeter lines of SO₂, but their brightness measured with the JCMT is much lower than the infrared $T_{\rm ex}$ . The implied limit on the source size is $\theta_{\rm S} = 15'' (T_{\rm ex} / T_0)^{-0.5}$ , with T₀the typical $T_{\rm MB}$ of 0.1 K. Inferred values of $\theta_{\rm S}$ range from 0 $.\!\!^{\prime\prime}$ 17 for AFGL 2591 to 0 $.\!\!^{\prime\prime}$ 32 for MonR2 IRS3, corresponding to linear radii of 100-250 AU. The limited spectral resolution of the infrared data prohibits assignment of the absorbers to either velocity component in emission, which however does not affect our conclusion.

4 Radiative transfer analysis

4.1 Model description

To estimate molecular abundances we have used the Monte Carlo radiative transfer program by Hogerheijde & van der Tak (2000). Starting from radial profiles of the density and kinetic temperature, this program solves for the molecular excitation as a function of radius. Besides collisional excitation, radiation from the cosmic microwave background and thermal radiation from local dust are taken into account. The result is integrated over the line of sight and convolved with the appropriate telescope beam. Observed and synthetic line fluxes are compared with a $\chi^2$ statistic to find the best-fit abundance. Initially, abundances were assumed to be constant with radius.

Molecular data were mostly taken from the database by Schöier et al. (in prep.), which paper also describes our extension of the SO₂ rate coefficients by Green (1995) to higher energy levels. No rate coefficients are available for NS, so its excitation was assumed to be thermalized at the temperature of each grid point of the model. The population of the $^2\Pi_{3/2}$ state of NS, which lies 225 K above the $^2\Pi_{1/2}$ state, was assumed to be negligible. In the case of H₂S, Ball et al. (1999) measured the 1₁₀-1₀₁ rate coefficient in collisions with He at T=1.36-35.3 K. Our model uses the rate coefficients by Green et al. (1993) for H₂O-He, multiplied by 2.6 based on comparison with the Ball et al. data, and scaled for the different reduced mass of the H₂S-H₂ system. Turner (1996) lists reasons why H₂S rate coefficients are likely to be larger than H₂O values.

4.2 Physical structure

The temperature and density structures of all our sources but MonR2 IRS3 were derived by van der Tak et al. (2000b). We have followed the same procedure to infer the structure of MonR2 IRS3, assuming a power-law density structure $n=n_0(r/r_0)^{-\alpha }$ . The temperature profile is derived from the observed luminosity of $1.3\times 10^4$ $L_{\odot}$ (Henning et al. 1992), while n₀ is constrained by the submillimeter photometry of Giannakopoulou et al. (1997). These dust models solve the grain heating and cooling self-consistenly, using opacity model 5 from Ossenkopf & Henning (1994). These "OH5'' opacities are the only ones that yield dust masses consistent with C¹⁷O measurements of AFGL 2591, where CO depletion is known to be negligible (van der Tak et al. 1999). To constrain $\alpha$ , the CS and C³⁴S line spectrum of MonR2 IRS3 was modeled with the Monte Carlo program. Besides our own JCMT measurements (Table 5), data from Tafalla et al. (1997) and Choi et al. (2000) were used. The best fit was obtained for $\alpha=1.25$ (Table 8), with a reduced $\chi^2=2.33$ over 11 degrees of freedom, although $\alpha=1.0$ ( $\chi^2=3.1$ ) and $\alpha=1.5$ ( $\chi^2=3.8$ ) also give acceptable fits. Radial profiles of 350 $\mu$ m dust emission for our sources indicate larger values of $\alpha$ , by up to 0.5, than CS line emission (Mueller et al. 2002). If this difference holds more generally, it may reflect a reduced CS abundance at small radii, although the models of Sect. 5.2 do not predict that.

Table 8: Power law model $n=n_0(r/r_0)^{-\alpha }$ for the density structure of MonR2 IRS3.

4.3 Molecular abundances

Table 9 lists the results of the radiative transfer analysis, with the same source ordering as in Table 7. Abundances, if assumed constant, are $\sim$ 10^-8 for OCS, $\sim$ 10^-9 for H₂S, H₂CS, SO and SO₂, and $\sim$ 10^-10 for HCS⁺ and NS. These trends are the same as found above for beam-averaged column densities (Sect. 3.3); the Monte Carlo analysis does not change results by factors of more than a few. Except for SO, SO₂ and OCS, the source-to-source spread in the abundances is less than a factor of 10, and not related to temperature. These molecules (CS, H₂S, H₂CS, NS and HCS⁺) appear to trace the chemically inactive outer envelope (T<100 K), at least in the observed transitions. However, the high OCS abundances occur all in warm sources.

$\begin{figure} \par\resizebox{13.1cm}{!}{\includegraphics[clip]{ms3386.f4}} \end{figure}$	Figure 4: Observations of SO (open circles; envelope component only), and predictions from Monte Carlo models (filled dots). Line fluxes have been converted to upper state column densities following Helmich et al. (1994). The estimated error on the line fluxes is a factor of two, due to the uncertain contribution from high-velocity gas.
Open with DEXTER

For SO and SO₂, the number of detected lines is large enough to investigate their abundances as a function of radius. Figure 4 shows that the Monte Carlo models fit the SO lines uniformly well over the observed range of energy levels for all sources except W 33 A: there is no indication for changes in the SO abundance between T=30 and 90 K. In the case of W 33 A, some of the observed values appear to be higher than the models, which may be due to an uncertain high-velocity contribution for these lines. Observations of higher-excitation SO lines ( $N=10\to9$ near 430 GHz or $N=11\to10$ near 473 GHz) would be valuable to investigate the effect of evaporating ice mantles on the abundance of SO.

$\begin{figure} \par\resizebox{13.16cm}{!}{\includegraphics[clip]{ms3386.f5}} \end{figure}$	Figure 5: As previous figure, for SO₂. Triangles indicate models with a temperature-dependent SO₂ abundance.
Open with DEXTER

Figure 5 shows the results of the models for SO₂. Models with SO₂/H₂ $\sim 10^{-9}$ fit the lines with $E_u \la 100$ K within factors of a few. However, these models underproduce the lines with $E_u \ga100$ K by factors of $\ga$ 100. This discrepancy is best seen for W3 IRS5, W33A and AFGL 2591, where the most SO₂ lines are detected. For these sources, we have run models where the SO₂ abundance increases sharply ("jumps'') when T>100 K. These models are shown as triangles in Fig. 5 and quantified in Table 10. The SO₂ abundance increases by a factor of 100 in the case of W3 IRS5 and W33A, and by a factor of 1000 for AFGL 2591. These models reproduce all observed SO₂ lines to factors $\sim$ 10. The residuals are considerably bigger than the expected error bars of factors of 2...3, mainly due to the uncertain high-velocity contribution, which suggests that more detailed modeling is needed. However, these "jump'' models fit the data much better than models where the SO₂ abundance is constant. The SO₂ abundances derived from the high-excitation lines are similar to those based on the infrared data, indicating that the same gas is probed. The data for the other sources are consistent with similar abundance jumps, but not enough high-excitation SO₂ lines have been detected to constrain such models. However, the data of AFGL 2136, NGC 7538 IRS1 and MonR2 IRS3 allow estimates of SO₂ abundances in their inner envelopes, which Table 10 also lists. Except NGC 6334 IRS1, the jumps in the SO₂ abundance occur all in sources with high T(CH₃OH), confirming that the abundance enhancement is related to high temperatures, probably through ice evaporation.

The models discussed so far assumed a static envelope and aimed to describe the low-velocity gas. To investigate whether the high-velocity gas could be due to infalling motions, we have modeled the SO emission from W3 IRS5 under the assumption of gravitational infall. In this model, the infall velocity $V_{\rm inf}$ varies with radius R as $\sqrt{GM/R}$ , where G is the gravitational constant. The luminosity of W3 IRS5, if due to a single star, suggests a stellar mass of $\approx$ 30 $M_{\odot}$ (Vacca et al. 1996), so that $V_{\rm inf}$ ranges from 0.45 to 10.66 km s^-1 in the model. However, the synthetic SO line profiles in 14-18'' beams have widths of 5-6 km s^-1 FWHM, much less than the observed value of 13.5 km s^-1. Since the other sources have lower luminosities (hence infall speeds), it seems unlikely that the high-velocity SO and SO₂ emission is solely due to infall.

5 Discussion

5.1 Comparison to other objects

Tables 9 and 10 also list molecular abundances in other regions. Only the case of the low-mass protostar IRAS 16293 has been analyzed to the same level of detail as done here. To reduce uncertainties introduced by methodological differences, the table contains abundances scaled to the value of the commonly observed CS molecule. Using ratios also facilitates comparison with chemical models, because these are independent of the total amount of sulphur in the gas phase, which is not well-understood. Ultraviolet observations of diffuse clouds toward $\zeta$ Oph indicate that sulphur is undepleted from its solar abundance of ${\sim}2\times 10^{-5}$ (Savage & Sembach 1996). However, for dense clouds, it is unclear what fraction of sulphur is locked up in dust grains.

The ratios observed in our sources differ considerably from those in PDRs, dark clouds and shocks, indicating that in those regions, other chemical processes dominate. The ratios in IRAS 16293 are at the low end of the values found here, suggesting that in low-mass objects, depletion of molecules onto grains is more important.

In the case of hot cores, both abundance ratios and absolute abundances are similar to those observed here. The chemical similarity of our sources with hot cores prompts us to compare our results with chemical models invoking ice evaporation. Our assumption is that source-to-source differences in abundances are due to age differences, rather than to differences in initial physical or chemical conditions.

5.2 Comparison with chemical models

The sulphur chemistry of hot cores, where ice mantles are evaporating off dust grains and start a high-temperature chemistry, was studied by Charnley (1997). The model assumes that the dominant form of sulphur in the ices is H₂S, which after evaporation is destroyed by H and H₃O⁺ to form S and H₃S⁺. Reactions with OH and O₂ make first SO and then SO₂. Some H₂S reacts with C⁺ into CS, which reacts with OH to make OCS. The sulphur chemistry strongly depends on temperature through its connection with the oxygen chemistry. At $T\ga300$ K, the SO₂ abundance remains factors of 10-100 below that of SO because O and OH are driven into H₂O. Results also depend on the presence of O₂ in the ice mantles. Observational limits on solid and gaseous O₂ (Goldsmith et al. 2000; Pagani et al. 2003; Vandenbussche et al. 1999) do not firmly rule out either assumption, but do make the model without O₂ ice more plausible. However, this model cannot directly be compared with our observations, which probe a mixture of warm and cold gas.

The chemistry of protostellar envelopes was modeled by Doty et al. (2002). In the inner warm region, a hot core chemistry similar to that of Charnley (1997) is used, whereas in the cold outer parts, a low-temperature chemistry including freeze-out is adopted. The initial conditions of the model depend on temperature to mimic the effect of ice evaporation. Although Doty et al. adopt the specific temperature and density structure of AFGL 2591 as modeled by van der Tak et al. (2000b), their results do not change by factors of more than a few for our other sources, which have temperature and density distributions similar to those of AFGL 2591.

The Doty et al. (2002) model has most gas-phase sulphur initially in S at T<100 K (S/H₂ $=6\times10^{-9}$ ) and in H₂S at T>100 K (H₂S/H₂ $=1.6\times10^{-6}$ ) to mimic the effect of ice evaporation. This choice of initial conditions is consistent with the nondetection of the [S I] line at 25.249 $\mu$ m toward our sources with ISO. The upper limit of $\tau<0.01$ in absorption implies N(S) $<5\times 10^{15}$ cm^-2, so that for N(H₂) $\sim 10^{23}$ cm^-2, S/H₂ $~~<5\times 10^{-8}$ , indicating that not all gas-phase sulphur in the inner envelope can be in atomic form. The inferred abundances imply that SO₂ is one of the dominant sulphur-bearing gas-phase molecules in the warm inner envelope. However, the sum of all sulphur-bearing molecules in the gas phase is still much lower than the value derived for diffuse clouds (Savage & Sembach 1996). Unless a major sulphur-bearing gas-phase species has been missed in this survey, $\sim$ 90% of the sulphur must be in a solid form which has not yet been detected.

Infrared observations of our sources with ISO-SWS do not support the assumption that H₂S is the main sulphur reservoir of the grain mantles. Toward W33A, a source with large ice abundances, the 3.98 $\mu$ m band of solid H₂S is not detected to $\tau < 0.02$ , which using a band strength of $5\times 10^{-18}$ cm mol^-1 gives $N_{\rm s}$ (H₂S) $<3.6\times 10^{17}$ cm^-2, or $N_{\rm s}$ (H₂S)/ $N_{\rm s}$ (H₂O) <0.03 (W. Schutte, priv. comm.). However, solid OCS was detected with an abundance of $N_{\rm s}$ (OCS)/ $N_{\rm s}$ (H₂O) $=(0.4-1)\times 10^{-3}$ or OCS/H₂ $=(4-14)\times 10^{-8}$ (Palumbo et al. 1997). Evaporation of solid OCS may explain our measured gas-phase OCS abundances which are otherwise unaccounted for. W 33 A is the only source in our sample for which both gas-phase and solid OCS have been detected. The results imply a gas/solid ratio of $\sim$ 0.5, much higher than that found for e.g. H₂O and CO₂ in this source (Boonman & van Dishoeck 2003; Boonman et al. 2003). Observations of both gas-phase and solid-state OCS toward other sources are needed to determine the role of grain-mantle evaporation. Observations of OCS lines from a wide range of energy levels may reveal an abundance increase at $T\ga 100$ K, as indeed found for IRAS 16293 (Schöier et al. 2002). The low $T_{\rm ex}$ of H₂S and the high value of OCS (Sect. 3.2) support these conclusions.

The effect of evaporating OCS-rich ice mantles on hot core chemistry is explored by Doty et al. (2003). Briefly, the gas-phase abundances of CS, HCS⁺ and OCS increase about linearly with the fraction of solid OCS, while those of H₂S, SO and SO₂ decrease about linearly. The H₂CS molecule is not affected much. Since the gas-phase chemistry of CS, HCS⁺ and OCS depends only weakly on time, the use of sulphur molecules as chemical clock diminishes when H₂S is only a minor sulphur carrier in the ice mantles. However, the "jump'' models for SO₂ (Table 10) indicate that either H₂S or SO₂ is present in the grain mantles. The limits on solid SO₂ from infrared observations are not very stringent (Boogert et al. 1998).

We have compared our observations to the models of Doty et al. (2002), modified to have OCS in the evaporating ice mantles (S. Doty, priv. comm.). This model reproduces our observed abundances of CS, SO and H₂CS in the outer envelopes (Table 9) to within factors of a few for chemical ages of $\sim$ 10⁵ yr. The abundances of SO₂ in the inner envelopes (Table 10) are also reproduced for the same chemical age. The model abundances of HCS⁺ and SO₂ in the outer envelopes are factors of $\sim$ 10 below the observed values, which may be due to the adopted initial conditions, or the chemical network, which is based on the UMIST reaction rate database (Le Teuff et al. 2000). The outer envelope abundances of H₂S and OCS are underproduced even more (factors of $\sim$ 100). However, the small number of observed lines forced us to assume constant abundances for these species, while they may in fact be confined to the inner envelopes.

5.3 The role of shock chemistry

As an alternative to ice evaporation models, we consider the effect of non-dissociative shocks on the molecular envelopes surrounding the protostars. Such shocks occur by the interaction of the molecular outflows of protostars with their envelopes. Shocks fast enough to dissociate H₂ ( $v\ga40$ km s^-1) would have a major impact on their molecular composition which our data do not indicate. Although our sources show outflows faster than that, both in infrared absorption (Mitchell et al. 1991) and in submillimeter emission (Mitchell et al. 1992; Mitchell & Hasegawa 1991) of CO, their filling factor must be small.

Our observed SO, SO₂ and H₂S abundances in the outer envelopes are consistent with those in postshock gas if sulphur was initially in molecular form (Leen & Graff 1988). The observed similarity of these abundances is not reproduced by the Pineau des Forêts et al. (1993) model which has most sulphur initially in atomic form. However, the two models are hard to compare because one is a two-fluid treatment, while the other treats neutral and charged particles as one fluid. Before concluding that postshock gas reproduces our data, we recommend that a two-fluid calculation is performed with most sulphur initially in molecular form.

Two other molecules suggest that the envelopes of young high-mass stars may have been processed by shocks, although neither one conclusively. First, ISO-SWS data show low ratios of gas-phase to solid-state CO₂ ( $\la$ 0.1; van Dishoeck 1998) and the CO₂ gas phase abundance remains low through the 100-300 K temperature regime. After evaporating off grains, CO₂ must be promptly destroyed, which shocks can do in reactions with H (Charnley & Kaufman 2000) or perhaps H₂ (Doty et al. 2002; Talbi & Herbst 2002). Second, Hatchell & Viti (2002) measured NS/CS ratios of 0.02-0.05 in a sample of six hot cores and interpreted these as evidence for shocks. The main reactions to form NS require SH and NH which are produced from OH. Shocks use OH to form H₂O and suppress the production of NS. The values of NS/CS = 0.001-0.01 measured here may indicate that shocks play a role too. However, for $t=3\times 10^4$ yr, the Doty et al. (2002) model also predicts NS/CS $\sim 10^{-3}$ , so this ratio cannot be used to demonstrate the influence of shocks.

6 Conclusions

Submillimeter observations of SO, SO₂, H₂S, H₂CS, OCS, NS and HCS⁺ toward nine embedded massive stars show that:

Sulphur chemistry in the envelopes of massive young stars