Issue 
A&A
Volume 588, April 2016



Article Number  A119  
Number of page(s)  24  
Section  Stellar structure and evolution  
DOI  https://doi.org/10.1051/00046361/201527943  
Published online  30 March 2016 
Sulphur molecules in the circumstellar envelopes of Mtype AGB stars^{⋆}
Onsala Space Observatory, Department of Earth and Space Sciences, Chalmers University of Technology, 439 92 Onsala, Sweden
email: taissa@chalmers.se
Received: 10 December 2015
Accepted: 25 January 2016
Aims. The sulphur compounds SO and SO_{2} have not been widely studied in the circumstellar envelopes of asymptotic giant branch (AGB) stars. By presenting and modelling a large number of SO and SO_{2} lines in the low massloss rate Mtype AGB star R Dor, and modelling the available lines of those molecules in a further four Mtype AGB stars, we aim to determine their circumstellar abundances and distributions.
Methods. We use a detailed radiative transfer analysis based on the accelerated lambda iteration method to model circumstellar SO and SO_{2} line emission. We use molecular data files for both SO and SO_{2} that are more extensive than those previously available.
Results. Using 17 SO lines and 98 SO_{2} lines to constrain our models for R Dor, we find an SO abundance of (6.7 ± 0.9) × 10^{6} and an SO_{2} abundance of 5 × 10^{6} with both species having high abundances close to the star. We also modelled ^{34}SO and found an abundance of (3.1 ± 0.8) × 10^{7}, giving an ^{32}SO/^{34}SO ratio of 21.6 ± 8.5. We derive similar results for the circumstellar SO and SO_{2} abundances and their distributions for the low massloss rate object W Hya. For the higher massloss rate stars, we find shelllike SO distributions with peak abundances that decrease and peak abundance radii that increase with increasing massloss rate. The positions of the peak SO abundance agree very well with the photodissociation radii of H_{2}O. We also modelled SO_{2} in two higher massloss rate stars but our models for these were less conclusive.
Conclusions. We conclude that for the low massloss rate stars, the circumstellar SO and SO_{2} abundances are much higher than predicted by chemical models of the extended stellar atmosphere. These two species may also account for all the available sulphur. For the higher massloss rate stars we find evidence that SO is most efficiently formed in the circumstellar envelope, most likely through the photodissociation of H_{2}O and the subsequent reaction between S and OH. The Sbearing parent molecule does not appear to be H_{2}S. The SO_{2} models for the higher massloss rate stars are less conclusive, but suggest an origin close to the star for this species. This is not consistent with current chemical models. The combined circumstellar SO and SO_{2} abundances are significantly lower than that of sulphur for these higher massloss rate objects.
Key words: stars: massloss / stars: AGB and postAGB / circumstellar matter / stars: evolution
© ESO, 2016
1. Introduction
Low to intermediatemass stars eventually evolve from the main sequence to the asymptotic giant branch (AGB). AGB stars lose mass rapidly, producing a circumstellar envelope (CSE) of atomic and molecular matter and dust, rich in chemical diversity. The variety and abundances of molecules that can be found in the CSEs of AGB stars depend on the chemistry of the individual star. For example, carbon stars, which have carbontooxygen ratio C/O > 1, are most likely to have a variety of Cbearing molecules in their CSEs (e.g. Gong et al. 2015), while oxygenrich Mtype stars, with C/O < 1, are more likely to contain a variety of Obearing molecules (e.g. Justtanont et al. 2012).
SO and SO_{2} are two such Obearing molecules. They are thought to exist in shells in the CSE, having been formed through the photodissociation of the parent molecule H_{2}S and subsequent reactions with O and OH (Cherchneff 2006; Willacy & Millar 1997). Observations of the red supergiant VY CMa by Adande et al. (2013) contradict this view, however; the modelling results show small, concentrated envelopes of SO and SO_{2} and indications that SO_{2} may itself be formed directly, rather than being a photodissociation product of H_{2}S or another molecule and are found in a hollow shell around the star. Similarly, when Decin et al. (2010a) modelled SO_{2} emission around the AGB star IK Tau, they had difficulty reconciling their observations with a shell model.
The Yamamura et al. (1999) ISO/SWS detections of the 7.4 μmν_{3}SO_{2} band in a few AGB stars suggest that SO_{2} is formed in the warmest regions of the CSE. Analysis of these data by Yamamura et al. (1999) and Cami et al. (1999) indicates that the SO_{2} is mostly likely formed within a few stellar radii of the star at a temperature of ~600 K. Cami et al. (1999) also find that the excitation of SO_{2} to the ν_{3} band varies with pulsation period. Their simple models put the outer radius of SO_{2} within ~ 5R_{∗}.
In this paper we present new observations of circumstellar SO and SO_{2} from an APEX spectral survey of the Mtype AGB star R Dor. We combine these results with SO and SO_{2} detections from Herschel/HIFI, developing comprehensive models of the SO and SO_{2} distributions around R Dor using 17 SO lines and 98 SO_{2} lines, all spectrally resolved.
We also model the sparse detections of SO and SO_{2} emission towards the other Mtype AGB stars observed with Herschel/HIFI, supplemented with archival data where available.
2. Sample and observations
The stars included in this study come from the sample of Mtype AGB stars observed as part of the HIFISTARS guaranteed time key programme (Justtanont et al. 2012, and see Sect. 2.2 for details). The OH/IR stars are excluded, as is Mira, which has a complicated and asymmetric CSE induced by a white dwarf companion (see Ramstedt et al. 2014). That leaves a sample of five Mstars, four of which had SO and SO_{2} lines detected by HIFI. The remaining star, TX Cam, has previously been detected in SO at lower frequencies.
Some basic information about the five sources is given in Table 1.
Basic information about our five sources.
2.1. APEX data
We performed a spectral survey of R Dor in the ranges 213−321.5 GHz and 338.5−368.5 GHz (λ = 0.8−1.4 mm) using the Swedish Heterodyne Facility Instrument (SHeFI; Vassilev et al. 2008) on the Atacama Pathfinder Experiment telescope (APEX). The data were observed over several observing seasons between May 2011 and June 2015. The observations were carried out using beam switching with a standard beam throw of 3′. A detailed description of this survey will be presented by De Beck et al. (in prep.).
Data reduction was carried out using the Gildas/Class^{1} package. Scans with very unstable baselines were ignored and bad channels were blanked. After masking the regions with line emission, polynomial baselines of typically first degree were subtracted from the averaged spectra to obtain a 0 K baseline. Rms noise levels throughout the survey are around 2−10 mK at a velocity resolution of 1 km s^{1}. The spectra were then converted to main beam temperatures using efficiency correction factors of η_{mb} = 0.75 for ν< 270 GHz, η_{mb} = 0.74 for 270 <ν< 320 GHz, and η_{mb} = 0.73 for ν> 320 GHz. The halfpower beamwidths were calculated using the general formula (1)where ν is in GHz and θ is in arcseconds. The beamwidths across our frequency range are between 17–29″.
The detections of SO using APEX are listed in Table 2 and the SO_{2} detections are listed in Table C.1. There were also some detections of SO and SO_{2} isotopologues: ^{34}SO, and tentative detections of ^{34}SO_{2} and SO^{18}O. We model ^{34}SO, but are unable to perform a full radiative transfer analysis for the other isotopologues. See Table 3 for a list of isotopologue detections and for the full discussion, see Sect. 3.2.4.
In terms of other Sbearing molecules, there were no conclusive detections of either CS (out of three possible transitions in the range J_{up} = 5 to J_{up} = 7) or SiS (out of nine possible transitions in the range J_{up} = 12 to J_{up} = 20). There was a tentative detection of CS (6 → 5) but it is blended with ^{29}SiO(7 → 6, v = 3) line and hence allows no reliable conclusion on the detection of CS. (We note that there are several other detections of ^{29}SiO in the survey, but none of CS.) No other Sbearing molecules were detected in this survey.
SO observations towards R Dor using APEX, listed in order of descending energy of the upper level.
SO and SO_{2} isotopologue observations towards R Dor using APEX.
2.2. HIFI data
R Dor, IK Tau, R Cas, TX Cam and W Hya were observed as part of the HIFISTARS guaranteed time key programme, using the Herschel/HIFI instrument (de Graauw et al. 2010) to observe emission lines with high spectral resolution. The full results are presented in detail in Justtanont et al. (2012). Since those data were published, there have been updates to the main beam efficiencies (Mueller et al. 2014^{2}) and for this work we have rereduced the HIFI data to take this into account (using HIPE^{3} version 12.1, Ott 2010). We have also identified three additional SO_{2} lines that were not included in Justtanont et al. (2012). The detected SO and SO_{2} HIFI lines are listed in Table C.2. We note that no SO or SO_{2} lines were detected with HIFI in TX Cam.
Stellar properties and input from CO models.
2.3. Archival data
To supplement the HIFI data for IK Tau, R Cas, W Hya, and TX Cam, we have used observations found in the literature. These are listed in Table C.3. As the older data generally covers lowerenergy transitions than those observed by HIFI, we are better able to constrain our models over a larger energy range. This is particularly important for R Cas, IK Tau and TX Cam where the HIFI lines (or nondetections in the case of TX Cam) are clustered close together energetically.
3. Modelling
3.1. Modelling procedure
We perform detailed radiative transfer modelling of the molecular emission lines using an accelerated lambda iteration method code (ALI), which has been previously described and implemented by e.g. Maercker et al. (2008), Schöier et al. (2011), Danilovich et al. (2014). ALI is particularly useful in this work as it is able to take into account extensive descriptions of molecular properties – such as large numbers of energy levels and transitions – while still fully solving the statistical equilibrium equations and taking temperature and velocity profiles into account.
We assume a smoothly expanding spherical CSE produced by a constant massloss rate. The molecules are located in this CSE until they eventually become photodissociated. They are excited by collisions with H_{2} molecules and through radiation from the star, the dust, and the cosmic microwave background. ALI input parameters such as the kinetic temperature distribution, dust temperature, and dust optical depth, are taken from CO modelling and, where applicable, are listed in Table 4. For R Dor, R Cas, IK Tau, and TX Cam Maercker et al. (in prep.) performed detailed radiative transfer modelling of the CO and H_{2}O lines and we use their results in our modelling. We based our CO model of W Hya on the results of Khouri et al. (2014a), but generated a CO model using the same code as in Maercker et al. (in prep.) for consistency between the stars.
We calculated the best fit model for each star and molecule using a χ^{2} statistic, which we define as (2)where I is the integrated line intensity, σ is the uncertainty in the observations, and N is the number of lines being modelled. We also calculate a reduced χ^{2} value such that where p is the number of free parameters.
After testing both centrallypeaked and shelllike abundance distributions, we came to the conclusion that the best radial abundance distribution profiles for both SO and SO_{2} in R Dor and W Hya were Gaussian profiles of the form (3)where f_{p} is the peak abundance at the inner radius, and R_{e} is the efolding radius, the radius at which the abundance has dropped by a factor of 1 /e.
Fig. 1 SO energy level diagram with levels labelled using the N_{J} convention. Transitions detected with HIFI and APEX towards R Dor are indicated in light blue and green, respectively. 

Open with DEXTER 
In the cases of IK Tau and R Cas, we found that a shell model was a better fit to the observed SO lines. As such, we modelled IK Tau and R Cas assuming a Gaussian shell for the abundance distribution of the form (4)where f_{p} is the peak abundance, R_{p} is the radial distance of the peak of the distribution from the centre of the star, and R_{w}, is the width of the shell at the efolding radius. Using a shell distribution for both IK Tau and R Cas rather than a central Gaussian distribution significantly improved the χ^{2} fits of the models.
Similarly, we can firmly rule out a centrallypeaked model for TX Cam, as for such a model to fit the archival data we would expect conclusive detections in the HIFI data. As the undetected HIFI lines are of higher energy than the archival detections, a lower abundance in the inner regions of the CSE is expected, than in the outer regions, which points to a shelllike abundance distribution.
3.1.1. SO
For the radiative transfer analysis of SO we include 182 rotational energy levels, denoted N_{J}, up to N = 30 in the ground and first excited vibrational states. There are 907 radiative transitions. These include pure rotational transitions in the X^{3}Σ^{−}v = 0 and v = 1 states as well as the v = 1 → 0 rovibrational lines. There are 8629 collisional transitions including collisions between all rotational states within a vibrational state, as well as between vibrational states. The rotational energy levels, transition frequencies, and Avalues have been adapted directly from the CDMS (Müller et al. 2001, 2005). The infrared line list has been computed directly from the rotational levels with the bandhead frequency adjusted to give very good agreement with the line positions measured by Burkholder et al. (1987). The vibrationrotation line strengths have been computed in intermediate coupling and have been verified by comparison with the pure rotational line strengths in the CDMS tables. The vibrationrotation transition dipole moment has been taken to be 0.08843 Debye, which yields inverse lifetimes of A_{tot} = 3.6 s^{1} for the v = 1 → 0 band as computed by Peterson & Woods (1990). The collisional rate coefficients for pure rotational transitions were adapted from the HeSO rates computed by Lique et al. (2006) with massscaling to H_{2} as in the smaller data set in the LAMDA database (Schöier et al. 2005). Rates for transitions within v = 1 were assumed to be identical to those within v = 0. Crude collision rates for v = 1 → 0 were scaled in proportion to normalised radiative line strengths for electricdipoleallowed transitions, with the largest values of the order of 1 × 10^{11} cm^{3} s^{1}.
In Fig. 1 we include an energy level diagram for SO. Here we have indicated all the transitions of SO detected towards R Dor with HIFI and APEX. These cover most of the transitions also detected in IK Tau, R Cas, W Hya, and TX Cam.
For the purposes of modelling the ^{34}SO emission in R Dor, we used a simpler molecular description than that for ^{32}SO, including the rotational energy levels up to N = 30, corresponding to those included for ^{32}SO, but only in the ground vibrational state. When adopting the corresponding simpler molecular description for ^{32}SO in the case of R Dor specifically, we found that the final best fit model only shifted by a few percent between the detailed and simpler descriptions, justifying this approach for ^{34}SO. There was, however, some shift in final model for the other, especially higher massloss rate, stars when changing between the detailed and simpler molecular description for SO.
3.1.2. SO_{2}
Our radiative transfer analysis of SO_{2} includes 2600 energy levels, denoted J_{Ka,Kc}, across the ground vibrational state and the ν_{1} = 1 (8.7 μm), ν_{2} = 1 (19.3 μm) and ν_{3} = 1 (7.3 μm) vibrationally excited states. Levels with energies up to 4830 K and J = 38 were included. This gives 15243 radiative transitions, with spectroscopic data taken from the HITRAN database (Rothman et al. 2013), and 15244 collisional transitions. The collision rates in the literature for SO_{2} are inadequate for our purposes. Green (1995) calculated rate coefficients for HeSO_{2} collisions in the infiniteorder sudden approximation for the lowest 50 rotational levels (up to 100 K excitation energy and J ≤ 13 only). Cernicharo et al. (2011) published rates for H_{2}SO_{2} collisions for the lowest 31 rotational levels at low temperatures, 5 to 30 K. The rates for H_{2} impact were found to be approximately 10 times higher than corresponding rates for He impact. For the much larger number of states in our models, we adopted instead a set of crude collision rates in which the downward rate coefficient is proportional to the radiative line strength and normalised to a total collisional quenching rate of 2.0 × 10^{10} cm^{3} s^{1}, which is comparable to the highest collision rates found by Cernicharo et al. (2011). We tested the impact of the chosen collisional transition rates by multiplying the rates, in stages, by up to two orders of magnitude in both directions. We find that such drastic changes had only a very small and barely detectable effect on the resulting models. Hence we conclude that SO_{2} excitation is radiatively dominated with the choice of collisional transition rates playing only a minor role in the radiative transfer modelling.
In Fig. 2 we include an energy level diagram for SO_{2}. Here we have indicated all transitions of SO_{2} detected towards R Dor with HIFI and APEX. As can be seen, SO_{2} has many close energy levels. This leads to a multitude of overlapping transitions, especially in AGB winds with typical expansion velocities of 5−25 km s^{1}. The number of levels, transitions and overlaps presents some computational challenges, especially when it comes to fully taking overlapping lines into account or running exhaustive grids. To reduce running time to a manageable interval we restrict the overlaps so that only those within the sampled frequency range, between 200−1200 GHz, are included. This reduced the total number of overlaps by more than an order of magnitude (down to 441 lines participating in overlaps), hence decreasing running time and memory usage. This, however, neglects possible overlaps in pumping lines, which could have a significant effect on some of the lines included in the model. From what tests we were able to run we believe that the overall impact of these omitted lines is relatively minor.
Fig. 2 SO_{2} energy level diagram. Levels are labelled J_{Ka,Kc}. Transitions detected with HIFI and APEX towards R Dor are indicated in light blue and green, respectively. 

Open with DEXTER 
3.2. R Dor
Our modelling is based on the radiative transfer results obtained by Maercker et al. (in prep) for CO in R Dor. They find a massloss rate of Ṁ = 1.6 × 10^{7}M_{⊙} yr^{1} and an expansion velocity of υ_{∞} = 5.7 km s^{1}. They also find an expansion velocity profile following (5)where υ_{min} = 3 km s^{1} is taken to be the sound speed at R_{in} = 1.6 × 10^{14} cm, the dust condensation radius. β = 1.5 governs the acceleration of the gas, having the most significant impact in the inner regions, and hence on the excitation of the higherenergy lines. The other relevant stellar properties of R Dor are listed in Table 4.
3.2.1. SO results
To model the 17 SO lines detected towards R Dor with APEX and HIFI, we set up a grid sampling different SO abundances and efolding radii. We then ran a finer grid with steps of 0.1 × 10^{6} in abundance and 0.1 × 10^{15} cm in efolding radius to find the best possible fit to the observations. The results of our χ^{2} analysis can be seen in Fig. 3. Our resulting bestfit model, with , has a peak SO abundance relative to H_{2} of (6.7 ± 0.9) × 10^{6} and efolding radius R_{e} = (1.4 ± 0.2) × 10^{15} cm and is plotted against the observed lines with respect to the LSR velocity in Fig. 4. A plot illustrating the goodnessoffit for all the lines is given in Fig. 5. The abundance profile for SO is plotted in Fig. 8 along with SO_{2} and the CO and H_{2}O results from Maercker et al. (in prep.) for comparison.
One of the detected SO lines, (8_{9} → 7_{8}), overlaps with SO_{2}(16_{4,12} → 16_{3,13}) in its wing. We note that this is the only SO line which is significantly overpredicted by the model. Our code is unable to properly take heteromolecular overlaps such as this into account. We suspect that although the SO_{2} line is much fainter than the SO line (in fact it is difficult to see even in Fig. A.1 where the SO_{2} model is overplotted), their interaction likely affects the flux from SO(8_{9} → 7_{8}).
Fig. 3 SO χ^{2} plots for R Dor, W Hya, IK Tau and R Cas. The contours show the confidence intervals and the shading represents the value for the corresponding model, with the colourbar indicating multiples of the minimum value. The white cross indicates our bestfit model (see Table 6). For IK Tau, the slice for which R_{w} = 1.8R_{p} is shown. For R Cas, the slice for which R_{w} = 1.0R_{p} is shown. 

Open with DEXTER 
Fig. 4 SO models (blue lines) and observations (black histograms) for R Dor. 

Open with DEXTER 
Fig. 5 SO goodness of fit plots for R Dor, R Cas, IK Tau, and W Hya. New HIFI lines as well as archival data listed in Table C.3 are included. The green points in the R Dor plots represent the observations from the APEX spectral survey. Undetected HIFI lines are shown as cyan points with arrows, in this case representing lower limits because the vertical axis is the ratio of model integrated intensities to observed integrated intensities or the upper limits thereof. 

Open with DEXTER 
3.2.2. SO_{2} results
We detected 100 SO_{2} lines in R Dor with APEX and HIFI. We exclude the v_{2} = 1 (25_{4,22} → 26_{1,25}) line at 279.497 GHz from our analysis since it is most likely a maser^{4}. We also concluded that it was not computationally viable to model the line with the highest energy level in the ground state, (40_{4,36} → 40_{3,37}) at 341.403 GHz, as the number of additional levels and transitions required to fully account for this line represented a significant increase in computation time. (We would have required 3583 levels and 19889 radiative transitions.) The two excluded lines are plotted in Fig. 6.
Fig. 6 SO_{2} lines excluded from modelling for R Dor. See text for full explanation. 

Open with DEXTER 
This leaves us with 98 detected SO_{2} lines with which to constrain our model. Our best fit model has f_{p} = 5.0 × 10^{6}, R_{e} = 1.6 × 10^{15} cm and . Due to the significant computational time in running SO_{2} models, we are unable to provide a comprehensive error analysis as we do for the SO model, hence the lack of formal uncertainties on our results. The model lines are plotted with the observed lines in Fig. A.1 with goodness of fit shown in Fig. 7. Overlaps are discussed in detail in Sect. 3.2.3. Figure 8 shows our bestfit abundance profiles for SO_{2} and SO, along with the results for CO and H_{2}O from Maercker et al. (in prep.).
There is a lot of scatter in the goodnessoffit plots in Fig. 7. There is no trend in goodnessoffit with upper energy level or J, but the observed lines that are most strongly underpredicted by the model are those lines for which the upper energy level has quantum number K_{a} ≥ 6 (see lower right plot in Fig. 7). This corresponds to the lines further away from the “backbone” of K_{a} = 0,1 energy levels in the energy level diagram in Fig. 2. We suspect this could be partially due to our exclusion of overlaps for lines outside of the observed frequency range (see Sect. 3.1.2). When testing models with and without overlaps enabled, we note that some lines that do not participate in overlaps can still be strongly affected by the inclusion (or not) of overlaps in our model. For example SO_{2}(27_{7,21} → 27_{6,22}) at 657.885 GHz was one such line, with the model predicting weaker emission by a factor of a few when overlaps were omitted. Unfortunately, due to computational limitations, it is not feasible to properly include overlaps in a full radiative transfer analysis, as discussed in Sect. 3.1.2. It should also be noted that the R Dor data were taken over a long observational campaign (see Sect. 2), so any variability in SO_{2} line brightnesses with pulsation period may contribute to the scatter.
Fig. 7 SO_{2} goodness of fit plots for R Dor. Top: goodness of fit with upper energy level of the transition. HIFI lines are shown as blue points and APEX lines are shown as green crosses. Error bars are excluded to make the plot clearer to read. Lower left: goodness of fit with J. Lower right: goodness of fit with K_{a}, a clear downwards trend for K_{a} ≥ 6. 

Open with DEXTER 
3.2.3. Overlapping lines
Overlapping lines in R Dor.
Table 5 contains an inventory of known line overlaps for the presented lines. In our radiative transfer modelling, we are able to take into account overlaps which occur between two lines of the same molecule – i.e. two SO_{2} lines. (For computational purposes we only include SO_{2} overlaps in the range 200 GHz−1.2 THz. Note, however, that all possible homomolecular overlaps are taken into account for SO in all modelled stars.) However, if there is a line overlap between two lines generated by different molecules, we are unable to properly treat this, as our code only allows for the modelling of one molecular species at a time. In R Dor we observe two such heteromolecular overlaps. The first between the SO(8_{9} → 7_{8}) and SO_{2}(16_{4,12} → 16_{3,13}) lines, where the much weaker SO_{2} line appears in the wing of the bright SO line, and the second between SO_{2}(13_{2,12} → 12_{1,11}) and H^{13}CN(4 → 3), where the two lines coincide very closely so as to be indistinguishable. Based on our model, we expect approximately half the flux to be due to the H^{13}CN(4 → 3) transition, which would agree with the H^{13}CN(3 → 2) line also covered by the APEX survey. However, without modelling H^{13}CN, it is not possible to fully gauge the impact of this overlap on our model.
The remaining line overlaps for lines modelled in this paper are homomolecular.
Three of the line pairs that are treated as overlapping in the code consist of a bright primary line in the vibrational ground state and a very weak secondary line in the ν_{2} = 1 vibrationally excited state. As can be seen in Fig. A.1, these secondary lines are not detectable above the noise in our observations, but are taken into account in our modelling.
The SO_{2} (24_{7,17} → 26_{6,18}) line at 659.898 GHz overlaps with the SO_{2} (40_{1,39} → 40_{0,40}) at 659.886 GHz and we would expect the latter to have an effect on the former. However, the SO_{2} (40_{1,39} → 40_{0,40}) line falls outside of the range of energy levels we included in our model. As noted in Sect. 3.1.2, it was not feasible to include a larger number of higher energy levels, hence this particular overlap is not taken into account in our modelling.
Fig. 8 Abundance profiles for R Dor, W Hya, IK Tau and R Cas. The abundances for CO and H_{2}O are taken from Maercker et al. (in prep.), except for W Hya, for which they are taken from Khouri et al. (2014a,b). The dashed line for the SO_{2} results for IK Tau and R Cas indicates that they are tentative. 

Open with DEXTER 
3.2.4. Isotopologue results
Based on the analysis of 6 ^{34}SO lines, and assuming the same efolding radius as found for ^{32}SO, we find a ^{34}SO abundance of (3.1 ± 0.8) × 10^{7} in a best fit model that has . This gives a ^{32}SO/^{34}SO ratio of 21.6 ± 8.5. The best fit model is shown in Fig. 9. The goodness of fit plot showing the ratio between the model and observed integrated intensities is shown in Fig. 5.
Modelling ^{34}SO_{2} in the same detailed manner as we have modelled ^{32}SO_{2} is impractical given the computational time required, the complexity of the molecular data file, and the low number of detected lines. However, all of the ^{32}SO_{2} lines we modelled are optically thin, so we can approximate the ^{32}SO_{2}/^{34}SO_{2} ratio by comparing the intensity ratios of two lines of the same transition. The best ^{34}SO_{2} transition for this purpose is 20_{0,20} → 19_{1,19}. Comparing the integrated intensities for this transition, we find a ^{32}SO_{2}/^{34}SO_{2} ratio of 21.6 ± 12.1, in good agreement with the result from ^{34}SO modelling.
The solar system value of ^{32}S/^{34}S is 22.5 (Cameron 1973) and Kahane et al. (1988) found a value of 20.2 for the carbon star CW Leo using SiS isotopologues, both in agreement with our results.
While a detailed model of ^{34}SO_{2} would be extremely time consuming, a detailed model of SO^{18}O would not be computationally feasible. Due to the asymmetry of the two oxygen atoms, SO^{18}O has approximately double the number of energy levels and transitions as SO_{2}, when looking at the same energy range, meaning that an SO^{18}O molecular data file would have to be approximately twice the size of our already very large SO_{2} file to probe a similar range of energies. The more complex energy level structure also means it is not possible to directly compare lines between SO^{18}O and SO_{2}, even when the transitions have the same quantum numbers. For ^{34}SO_{2} and SO^{18}O we present the (tentative) detections in Fig. A.2.
Fig. 9 ^{34}SO model (blue lines) and observations (black histograms) for R Dor. 

Open with DEXTER 
3.3. Other M stars
We model SO and SO_{2} line emission for the remaining stars using HIFI observations, as listed in Table C.2, and archival observations with different groundbased instruments, as listed in Table C.3. Several of these older observations probe energy levels significantly lower than the HIFI observations, allowing us to better constrain the size of the emitting molecular envelope. This is particularly important for IK Tau, where only the three N = 13 → 12 SO lines were detected with HIFI, as these are emitted from a similar region of the CSE.
The stellar parameters used in our SO and SO_{2} models, taken from CO model results, are listed in Table 4.
3.3.1. W Hya
In the case of W Hya we find an SO model that fits the data well using the Gaussian abundance distribution given in Eq. (3). We found f_{p} = (5.0 ± 1.0) × 10^{6} and R_{e} = (1.5 ± 0.5) × 10^{15} cm, with . This result is qualitatively similar to that of R Dor. As with R Dor, this suggests that SO in the CSE of W Hya is formed close to the star and is not found in a shell around the star as might be expected if it were a photodissociation product of another molecule such as H_{2}S. The HIFI observations and model line plots for SO are shown in Fig. 10. The corresponding χ^{2} plot is shown in Fig. 3.
Fig. 10 Models (blue lines) and observations (black histograms) for SO towards W Hya. 

Open with DEXTER 
The HIFI observations and model line plots for SO_{2} towards W Hya are shown in Fig. 11. The main difficulty we had in fitting an SO_{2} model was finding a model which fit the two highestenergy lines. As can be seen in Table C.2, the SO_{2}(37_{1,37} → 36_{0,36}) and SO_{2}(36_{1,35} → 35_{2,34}) lines are only ~ 3 K apart in upper energy level. Also we note that the lowerenergy line is almost a factor of 3 brighter than the higherenergy line. Our model invariably predicts a smaller difference in intensity with the higherenergy line being the brighter. The same is true for R Dor, however, in R Dor the detected lines reflect this (although the model fit is not perfect). This phenomenon is probably due in part to the noise in our observations but could also reflect a problem with our molecular description of SO_{2}. In this case, the most likely cause is the cutoff in included energy levels at J = 38. The variation in these lines cannot be due to variations in brightness due to stellar pulsations as both lines were observed simultaneously (and, indeed, all the SO_{2} lines in W Hya were observed within two days). In any case, the apparently outlying line of (37_{1,37} → 36_{0,36}) strongly contributes to the poorly fitting model we find for SO_{2} in W Hya. We are able to find a better fit by excluding this line, but do not have a strong basis for doing so, hence we leave it in.
Our best fit model for SO_{2} has f_{p} = 5.0 × 10^{6}, based on a small grid with steps of 0.5 × 10^{6}, and R_{e} = 3.0 × 10^{15} cm, based on a small grid with steps of 0.5 × 10^{15} cm. This model has . We also test an SO_{2} model using the parameters we found for SO. That model is not a significantly worse fit with almost the same .
The abundance distributions for SO and SO_{2}, along with the CO and H_{2}O abundance distributions from Khouri et al. (2014a,b) for comparison, are shown in Fig. 8.
Fig. 11 Models (blue lines) and observations (black histograms) for SO_{2} towards W Hya. 

Open with DEXTER 
3.3.2. IK Tau
When we try to fit the SO IK Tau observations with a centrally peaked Gaussian distribution, we cannot constrain the efolding radius with the available data. The χ^{2} analyses of centrallypeaked Gaussian models point towards very large efolding radii, significantly larger (by more than half an order of magnitude) than the halfabundance radius Maercker et al. (in prep.) found for the corresponding CO envelope. Since it is highly unlikely that the SO envelope is more extensive than that of CO, we conclude that a centrallypeaked Gaussian distribution is unlikely for SO in IK Tau. Instead, we run a threeparameter grid across f_{p}, R_{p}, and R_{w} (see Eq. (4)) to find the best model. We find f_{p} = (1.0 ± 0.2) × 10^{6}, R_{p} = (1.3 ± 0.2) × 10^{16} cm, and R_{w} = 1.8R_{p} (which we gridded in steps of 0.2R_{p}), with and the resultant lines are shown in Fig. 12.
The χ^{2} plot for SO in IK Tau is shown in Fig. 3. IK Tau has a significantly larger χ^{2} value for the best fit model (compared with R Dor and W Hya) because of some noisy observations. This is also seen in the goodness of fit plot in Fig. 5. In comparison, R Dor and W Hya have brighter and more uniform line observations, making it easier to find a good model fit.
Decin et al. (2010a) perform a radiative transfer analysis of IK Tau in a way that is similar to our method. They find an SO abundance distribution that is similar to our shelllike distribution, but with an increased abundance in the inner region. They find an abundance at 200R_{∗} (which corresponds to about 3 × 10^{15} cm) of ~ 2 × 10^{7} using two lines to fit the model. This did not change significantly in the follow up in Decin et al. (2010b) which included one of the HIFI lines as well. Our model results give a corresponding abundance about a factor of 2 higher at the same radius but using a different shape for the abundance distribution. We also use 10 lines with a broader range of energy levels to constrain the model.
In the case of SO_{2} in IK Tau we are unable to include overlaps as we do for R Dor and W Hya due to the larger expansion velocity of the circumstellar gas around IK Tau. The larger expansion velocity means there are a larger number of overlaps (since the lines are about three times wider than for R Dor) which quickly become computationally infeasible to fully account for.
We could not find a consistent model for IK Tau that matched all the available observed SO_{2} lines. In particular, there was a very large scatter in goodnessoffit for the lines with upper energy levels of 136 K or less (which is all of the lines other than the one HIFI observation). There was no way to simultaneously fit all these observed lines well. A centrallypeaked Gaussian model matches the data reasonably well – particularly the HIFI line, which according to the best shell model should have been a nondetection – and much better than the shell model. A Gaussian model with efolding radius located at the peak of the SO distribution is a better fit than a model with the SO distribution parameters, but we find that decreasing the efolding radius to R_{e} = 1 × 10^{16} cm gives a better fit again. We cannot constrain the efolding radius better than by a factor of ~ 2, however, because of the large scatter in the lowerenergy lines. The model we present in this paper, plotted in Fig. 13, has a peak SO_{2} abundance f_{p} = 2 × 10^{6}, and R_{e} = 1 × 10^{16} cm. This model has , the high value reflecting the poor overall fit. The large scatter in the IK Tau SO_{2} lines could be due to variability in line brightness with pulsation period. The data we used were observed at different times corresponding to different phases of pulsation. For example, the brightest lines (17_{1,17} → 16_{0,16}) and (13_{2,12} → 12_{1,11}), were observed less than two weeks apart close to maximum brightness in 2006. The most underpredicted line, (14_{3,11} → 14_{2,12}), was observed four months later when the star was approaching minimum brightness. On the other hand, the most wellfit lines – those with J = 5,4,3 as can be seen in Fig. 13 – were variously taken close to minimum and maximum brightness, so perhaps it is the higher J lines which are most strongly affected. Future monitoring of these lines observationally would allow us to confirm whether the effect on the higherJ lines is really due to variability over a pulsation period.
Fig. 12 SO models (blue lines) and observations (black histograms) for IK Tau. 

Open with DEXTER 
Fig. 13 SO_{2} model (blue line) and observations (black histograms) of IK Tau. For details on the archival observations, see Table C.3. 

Open with DEXTER 
The abundance distributions for SO and SO_{2} in IK Tau, along with the CO and H_{2}O abundance distributions from Maercker et al. (in prep.) for comparison, are shown in Fig. 8. In general, we do not consider our SO_{2} results for IK Tau conclusive. A more rigorous model which is properly able to take overlaps into consideration and which perhaps includes more lines in the intermediate to high energy range (with upper energy level > 136 K) is recommended.
Decin et al. (2010a) have similar issues modelling the SO_{2} in IK Tau, especially with the (17_{1,17} → 16_{0,16}) and (13_{2,12} → 12_{1,11}) lines which we also strongly underpredict, as can be seen in Fig. 13. When they exclude these two lines, Decin et al. (2010a) find a high inner abundance of SO_{2}, in general agreement with our results. The poor fit of our model could be a result of unusual structure in the CSE of IK Tau or could be a result of not being able to properly consider overlaps in the SO_{2} model. The higher wind velocity would also lead to more overlapping lines overall – including in regions we have not observed – which could have an effect on the overall energy distribution between all molecular energy levels.
Kim et al. (2010) use a combination of MonteCarlo radiative transfer modelling, to find CSE properties, and LTE formulations, to determine SO and SO_{2} abundance for IK Tau. For SO they find fractional abundances in the range 3 to 8 × 10^{7} and for SO_{2} their fractional abundances were in the range 4 × 10^{6} to 1 × 10^{5}. Their results are not accompanied by clear abundance distributions, making them difficult to compare with our results. Nevertheless, their SO result is very close to our peak abundance for SO, while their SO_{2} result is much higher than we found.
3.3.3. R Cas
As with IK Tau, we find that a model with a centrally peaked Gaussian distribution of SO does not match the observed data. We again run a threeparameter grid to find the best shellmodel fit to the data and find f_{p} = (6.0 ± 1.2) × 10^{6}, R_{p} = (3.2 ± 0.3) × 10^{15} cm, and R_{w} = 1.0R_{p} cm (gridded in steps of 0.2R_{p}), with . The resultant model lines are shown in Fig. 14 with the observations. The χ^{2} plot for SO in R Cas is shown in Fig. 3 and the goodness of fit plot is included in Fig. 5.
As there are only 2 SO_{2} lines observed towards R Cas, we are only able to find an approximate model for SO_{2}. As with IK Tau, a shelllike model based on the R Cas SO results does not fit the SO_{2} observations. Our best model has f_{p} = 7 × 10^{6} (best within steps of 1 × 10^{6}) and R_{e} = 6 × 10^{15} cm (best within steps of 1 × 10^{15} cm). In Fig. 15 we plot the HIFI detection with our model. An abundance plot for SO, SO_{2}, CO and H_{2}O towards R Cas is shown in Fig. 8.
We note that the HIFI detection in Fig. 15 has a central narrow peak, much narrower than the gas expansion velocity, This skews the overall integrated line intensity somewhat. Interestingly, IK Tau has a similar narrow peak in the same transition line (see Fig. 13), also at approximately the stellar velocity. R Dor does not have such a peak and W Hya may have one which is significantly less bright with respect to the rest of the emission line. The cause of this feature is unclear.
Fig. 14 SO models (blue lines) and observations (black histograms) for R Cas. 

Open with DEXTER 
Fig. 15 SO_{2} model (blue line) and observation (black histogram) for R Cas. 

Open with DEXTER 
3.3.4. TX Cam
In the case of SO towards TX Cam, we do not have sufficient constraints to run a full grid and perform a χ^{2} analysis as we did for the other stars. Instead we aim to fit the archival lines and find the largest R_{w} allowed by the HIFI nondetections. The best model with these assumptions has f_{p} = 1.7 × 10^{6}, R_{p} = 1.4 × 10^{16} cm and R_{w} = 1.6R_{p}. We plot the detected lines with our model in Fig. 16. The goodness of fit, represented by the ratio between model integrated intensities and observed integrated intensities, is plotted in Fig. 5. We stress that the dearth of observational results leaves our model poorly constrained and this is just one possible model that fits the available data. We can, however, rule out a centrally peaked model, as in that case we would expect the HIFI lines to have been detected, given the constraints on the fit from the archival lines.
There were no SO_{2} lines detected towards TX Cam.
Fig. 16 SO models (blue lines) and observations (black histograms) for TX Cam. We only plot the archival SO detections, not the nondetections from HIFI. 

Open with DEXTER 
SO and SO_{2} model results.
4. Discussion
4.1. SO distribution
Our results for circumstellar SO are summarised in Table 6. In Fig. 17 we plot the circumstellar SO abundance profiles of the stars we modelled. We also show the radial range probed by the available observational data for each line with a thicker line. These ranges were found by considering the brightness distributions for each emission line and the radii at which these fall to half of their maximum values. It is interesting to note that for R Cas, IK Tau, and TX Cam, the three stars with shelllike SO distributions, the location of the peak is found progressively further out with increasing massloss rate. The two low massloss rate stars, however, both seem to have centrally peaked SO distributions, which could be interpreted as shells with peaks close to the star, especially if the peaks are near or within our inner radii. Looking at the three stars with shelllike SO distributions, there also seems to be a trend of decreasing SO abundance with increasing massloss rate (or with the radius of peak SO abundance).
Fig. 17 SO abundance distributions for all stars modelled. The vertical lines represent the dust condensation radii, where our models stop. The thicker sections of the curves represent the area probed by our observations and for TX Cam the thick dashed line is the area probed by the upper limits imposed by the HIFI nondetections. 

Open with DEXTER 
In Fig. 18 we plot the peak positions against the wind density, Ṁ/υ_{∞}, and fit a power law to the three higher massloss rate stars (R Cas, TX Cam and IK Tau). The results for these stars are well fit by a power law (6)with α_{R} = 1.15 ± 0.24. We extend the power law to predict the peak positions for R Dor and W Hya, were they also to fit this trend. This predicts the peak in SO abundance for R Dor to lie at 1.0 × 10^{15} cm, which is close to the R_{e} we found, and for W Hya the predicted SO peak lies at 4.4 × 10^{14} cm, about three times smaller than the efolding radius. Running models with shelllike distributions at the predicted peaks, we found they could not provide as good a fit for either R Dor or W Hya as the starcentred Gaussian models. In general, the best shell models had at least twice the χ^{2} values of the best starcentred Gaussian models. Given the region probed by our observations as shown in Fig. 17, it is not surprising that changing the inner abundance would have an effect on the model fit.
We performed a similar fit for the peak abundance values against density for the three highest massloss rate stars. The results for these stars are well fit by a power law (7)with α_{f} = −1.29 ± 0.17, Fig. 18. Doing a similar extrapolation to predict the peak abundance values for R Dor and W Hya based on the power law, we find fractional abundance predictions of 2.2 × 10^{5} for R Dor and 5.8 × 10^{5} for W Hya. Both of these are higher than the values we find from our modelling and in the case of W Hya this represents more sulphur than should be available, i.e. it exceeds the solar and ISM abundances (see below). Because the SO abundance cannot increase with decreased massloss rate indefinitely, there must be a maximum SO abundance set by the abundance of sulphur.
Fig. 18 Trends in SO peak radius (top) and fractional abundance (bottom) against the circumstellar density measure, Ṁ/υ_{∞}. The green lines are the trends fitted to the three higher massloss rate stars (R Cas, TX Cam and IK Tau), while the yellow stars are the predicted locations of the lower massloss rate stars (R Dor and W Hya) based on the trend. In the fractional abundance plot, the red line shows the hard limit for SO abundance based on solar S abundance and the blue points in line with the yellow stars represent the real abundance values for R Dor and W Hya. 

Open with DEXTER 
The shelllike distributions of SO for the three stars with the highest massloss rates means that circumstellar chemistry, most likely related to photodissociation, must play an important role, in these cases. It is therefore interesting to compare with results on photodissociation for other species. H_{2}O is particularly interesting here, since, as will be discussed below, SO (and also SO_{2}) may owe its origin to the presence of circumstellar OH, which in turn is a photodissociation product of H_{2}O. Netzer & Knapp (1987) predict a peak OH radius that scales with both massloss rate and expansion velocity. Their formulation has been used by Maercker et al. (2008, 2009), Schöier et al. (2011), Danilovich et al. (2014) and others to define the efolding radius of H_{2}O, since OH is a photodissociation product of H_{2}O and peaks in abundance where H_{2}O drops off. In Fig. 19 we plot our SO peak abundance radii or efolding radii (as relevant) against the H_{2}O efolding radii of the same stars found by Maercker et al. (in prep.) and Khouri et al. (2014b, for W Hya). We find a strong correlation, which is close to being 1:1. The chemistry of SO will be discussed in Sect. 4.3.
Fig. 19 Peak abundance radii of SO (for R Cas, TX Cam, and IK Tau) and efolding radii of SO (for R Dor and W Hya), plotted against the efolding radii of H_{2}O found by Maercker et al. (in prep.) and Khouri et al. (2014b, for W Hya). The solid green line is the best fit to the data and the dashed red line traces a 1:1 relationship. 

Open with DEXTER 
4.2. SO_{2} distribution
Our results for circumstellar SO_{2} are summarised in Table 6. For R Dor and W Hya we find circumstellar SO_{2} distributions of similar size and abundance to the circumstellar SO distributions. Our results for the higher massloss rate stars, however, are less clear. For both IK Tau and R Cas, the best models are centrallypeaked Gaussian distributions of SO_{2}, rather than the shelllike models we found for SO. For IK Tau, we were unable to constrain the efolding radius better than by a factor of 2 and for R Cas we only had two observations, but in both cases shell models similar to the corresponding SO models are ruled out. Our SO_{2} results for IK Tau and R Cas suggest that SO_{2} is formed in the inner regions irrespective of the massloss rate. It appears that SO_{2} is formed more favourably than SO in the inner regions. There does not appear to be a strong correlation between efolding radius and massloss rates or H_{2}O efolding radii, as we found for SO, but this might become clearer if we add more SO_{2} observations to our models, especially for R Cas and W Hya. We emphasise that our SO_{2} results for the higher massloss rate stars, R Cas and IK Tau, are particularly uncertain.
4.3. Comparisons with chemical models
There exist two studies of the abundances of SO and SO_{2} in the extended atmospheres of AGB stars (Cherchneff 2006; Gobrecht et al. 2016). In both cases the effects of shockinduced chemistry, due to pulsational motion, are included. Cherchneff (2006) has chosen TX Cam as the representative star, and the modelling covers the region from 1 to 5 stellar radii, R_{∗}, which means that the outer reach of their model is approximately an order of magnitude smaller than our inner radii. For their model star with C/O = 0.75, they find an SO abundance at 5R_{∗} of 3 × 10^{7}, and an SO_{2} abundance several orders of magnitude lower than this. Gobrecht et al. (2016) use IK Tau as the example, and extend the calculations to 10 stellar radii (the outer radius of their model then approximately meets the inner radius used by us). Their model focuses on the shock chemistry in this region, much of which varies with the pulsation phase. Their abundances of SO and SO_{2} are much lower than we observe, at ~ 10^{8} and ~ 2 × 10^{9}, respectively. This means that the predicted SO and SO_{2} abundances close to the star, at least for these higher massloss rate stars, are substantially lower than we derive for our sample stars.
The shelllike SO distributions for the higher massloss rate stars suggest a circumstellar origin. Willacy & Millar (1997) describe circumstellar chemical models of four Mtype AGB stars, including R Dor, TX Cam, and IK Tau. Their models differ from ours in terms of CSE parameters, for example taking the inner radius to be 2 × 10^{15} cm, about an order of magnitude larger than our inner radii. They assume that all the sulphur is carried by H_{2}S that is eventually photodissociated. SO is subsequently formed through the following reactions which are favoured depending on the availability of OH and SH, respectively, and with Eq. 8 dominating at the lower gas temperatures in the CSE. Following this, SO can be destroyed through (10)and hence form SO_{2}. Unfortunately, they only visualise their model results for TX Cam, the least wellconstrained star of our sample. Nevertheless, the location of the peak of the SO distribution that we find for TX Cam agrees quite well with their predicted peak location. Our peak abundance is about 50% higher than theirs, but this must be considered to be within the errors. Their SO_{2} distribution for TX Cam peaks at roughly the same radius as SO, but with a peak abundance about an order of magnitude lower than that of SO (as we do not have any SO_{2} detections for TX Cam, we cannot compare this directly). The molecular column densities they list for R Dor, IK Tau, and TX Cam are, in general, a few orders of magnitude lower than those predicted by our models.
In conclusion, our SO results for the outer CSE are reasonably consistent with the results of Willacy & Millar (1997) for the higher massloss rate objects, although they do not predict a peak abundance that decreases with massloss rate. Thus, an origin through OH is likely, a result that is further strengthened by the correlation between SO and H_{2}O sizes that we found. However, we note that neither Cherchneff (2006) nor Gobrecht et al. (2016) predict high abundances of H_{2}S in the upper atmosphere. For the lower massloss rate objects, where the SO abundance is high close to the star, the models of Cherchneff (2006) and Gobrecht et al. (2016) fail by more than two orders of magnitude to reproduce our estimated abundances. They even fail to reproduce the inner SO abundances for the higher massloss rate stars.
In the case of SO_{2} we find no evidence for a photoinduced circumstellar origin along the Willacy & Millar (1997) model for any of our objects. Once again, we caution that the results for R Cas and IK Tau are uncertain. The SO_{2} abundances that we estimate for R Dor and W Hya are an order of magnitude higher than those predicted by Cherchneff (2006) and Gobrecht et al. (2016).
4.4. Sulphur chemistry: can we account for all the sulphur?
AGB stars and their progenitors do not produce S via nucleosynthesis. As such, the quantity of sulphur available to form molecules in the CSE of an AGB star is fixed and not dependent on the stage of evolution or mass of the star in question. Rudolph et al. (2006) find an S/H ratio in the ISM of ~ 10^{5} in the solar neighbourhood and Lodders (2003) indicate a solar S/H abundance of 1.5 × 10^{5}. All stars in our sample are at distances < 400 pc and hence can be assumed to trace a similar S/H abundance. In this work we refer to the fractional molecular abundance with respect to H_{2}. Hence, assuming all hydrogen is in the form of H_{2} in the CSE, we will take the S/H_{2} ratio to be ~ 2 to 3 × 10^{5}, which represents the maximum total amount of sulphur that should be found in an AGB star.
For R Dor and W Hya we find combined SO and SO_{2} abundances of ~ 1.2 × 10^{5} and ~ 1.0 × 10^{5}, respectively. Hence, in these cases most of the sulphur is locked up in SO and SO_{2} within the inner regions of the CSE and within the errors. This result is consistent with the nondetections (or lowlevel emission) for CS and SiS in the APEX spectral scan of R Dor, and no reported detections of these species towards W Hya. In the case of R Cas, the combined SO and SO_{2} abundances in the midCSE is ~ 1.4 × 10^{5}, suggesting that these two species carry all the sulphur, but here the uncertainty on the SO_{2} abundance is substantial. For the high massloss rate object IK Tau, the combined SO and SO_{2} abundance is well below that of sulphur.
In general, higher massloss rate stars also show definitive detections of other Sbearing molecules. For example, SiS was detected in several carbon and Mtype stars by Schöier et al. (2007) and Danilovich et al. (2015). Schöier et al. (2007) reported circumstellar SiS abundances of 4 × 10^{7}, 4 × 10^{7}, and 1 × 10^{7} for R Cas, TX Cam, and IK Tau, respectively, suggesting that SiS is less abundant than SO and SO_{2} by up to an order of magnitude, at least in the outer CSE. It should be noted that Schöier et al. (2007) assume a Gaussian distribution of SiS centred on the star, but find their model fit greatly improved when they include a highabundance inner component, which could represent the SiS reservoir before depletion through dust condensation. Decin et al. (2010a), who model IK Tau in detail, also find evidence of depletion of SiS.
Another Sbearing molecule, CS, has mainly been detected in carbon stars rather than Mtype stars. CS has been detected and modelled in IK Tau by Kim et al. (2010) and Decin et al. (2010a), detected in TX Cam and IK Tau by Bujarrabal et al. (1994, with nondetections in R Cas and W Hya) and Lindqvist et al. (1988, with a nondetection in R Cas). Derived CS abundances for Mtype stars have generally been low, in the range ~ 10^{8} to ~ 5 × 10^{8} (see Bujarrabal et al. 1994; Decin et al. 2010a, for examples).
H_{2}S is considered as a parent species of sulphur in the chemical modelling of Willacy & Millar (1997). However, H_{2}S has not been widely detected in AGB stars other than in OH/IR stars. For example, in the HIFISTARS project H_{2}S was only detected in AFGL 5379 (Justtanont et al. 2012), despite being in the observed range for all stars except TX Cam, while Justtanont et al. (2015) detected H_{2}S in all OH/IR stars observed with SPIRE and some observed with PACS. In a study of 25 stars, Ukita & Morris (1983) detected H_{2}S only in OH231.8+4.2 aka the Rotten Egg Nebula. Omont et al. (1993) detected H_{2}S in several high massloss rate stars, including several OH/IR stars. Of the stars we modelled, Ukita & Morris (1983) did not detect H_{2}S in W Hya, R Cas, and TX Cam, but it was detected in IK Tau by Omont et al. (1993) and De Beck et al. (in prep.). This suggests that H_{2}S may require high densities to form, or may be able to survive longer in the CSEs of high massloss rate stars, or that the excitation conditions are such that the emission is only bright enough in the very high massloss rate stars to be detectable. We note that Gobrecht et al. (2016) predict a fairly rapid decline of H_{2}S inside of the dust condensation radius (which is where our models start) even for IK Tau, which is a relatively high massloss rate object.
To check what the lack of detections predicts in terms of H_{2}S abundances, we run radiative transfer models for R Dor and IK Tau to find upper limits for the H_{2}S abundances based on the nondetections in HIFI, the nondetection in APEX for R Dor, and the SMA detection in IK Tau by De Beck et al. (in prep.). We used the orthoH_{2}S molecular data file available on LAMDA^{5} (Schöier et al. 2005) which includes the lowest 45 rotational energy levels, 139 radiative transitions with frequencies taken from JPL^{6} and 990 collisional transitions taken from Dubernet et al. (2009) for temperatures from 5−1500 K. For IK Tau, using the detection from De Beck et al. (in prep.) and the HIFI upper limit to also constrain the envelope size, we find a small envelope with R_{e} ≃ 4 × 10^{14} cm and f_{p} ≃ 4 × 10^{6}. This is consistent with a rapid destruction of H_{2}S. For R Dor, using both nondetections and assuming the R_{e} we find for SO, we find an upper limit on the abundance of f_{p} ≲ 2.5 × 10^{7}. If we instead use the H_{2}S envelope size found for IK Tau, the abundance upper limit increases slightly to f_{p} ≲ 6 × 10^{7}. In any case, these results limit the possibility that H_{2}S is a significant Scarrier in the inner CSE, certainly for the low massloss rate objects.
None of the molecules discussed thus far have been found in sufficient quantities towards higher massloss rate AGB stars to account for the full amount of expected sulphur. It is possible that the remaining sulphur is locked up in dust or left as atomic S or locked up in molecules that are difficult to detect for various reasons, such as the spectral region they are most likely to emit in, as is the case with HS. Both Cherchneff (2006) and Willacy & Millar (1997) predict a rapid decline of HS with radius as it is consumed by various chemical processes (although we note that the two studies make predictions for different regions around the star). The only detection of HS in the literature is through rovibrational lines identified by Yamamura et al. (2000) towards R And (an Stype AGB star). They estimate a molecular abundance of HS/H ~ 1 × 10^{7}, which is well below the sulphur limit. There have been no other detections of circumstellar HS, although it has been detected in the ISM (see e.g. Neufeld et al. 2015).
To fully study the issue of sulphur in the CSEs of AGB stars of different massloss rates, a more thorough investigation including more molecular species – such as SiS, CS, and H_{2}S in addition to SO and SO_{2} – across a larger sample of stars is needed.
5. Conclusions
We present new APEX observations of a very large number of SO and SO_{2} lines towards the low massloss rate Mtype AGB star R Dor. Combining these data with higherfrequency observations from Herschel/HIFI, we compute comprehensive radiative transfer models to determine the molecular abundances and distributions of the two molecules. For R Dor we find a Gaussian abundance distribution centred on the star, with a peak SO fractional abundance of (6.7 ± 0.9) × 10^{6} and efolding radius of (1.4 ± 0.2) × 10^{15} cm, and an SO_{2} fractional abundance of 5.0 × 10^{6} and efolding radius of 1.6 × 10^{15} cm. Our ^{34}SO model assumes the same efolding radius as for ^{32}SO and we find an abundance of (3.1 ± 0.8) × 10^{7}. This gives an ^{32}SO/^{34}SO ratio of 21.6 ± 8.5, which is in agreement with previous results from other nearby stars.
We also model SO in four other Mtype AGB stars that were observed as part of HIFISTARS: IK Tau, TX Cam, W Hya, and R Cas. For TX Cam for we are only able to provide an upper limit model since there are no SO lines detected with HIFI. Of these four stars only W Hya has a similar SO distribution to R Dor. The other three stars, all of which have higher massloss rates, are best fit with shelllike abundance distributions. We find that the radial position of the peak of the distributions increases with massloss rate, while the peak abundances decrease. The location of the peaks of the SO distributions correlates with the photodissociation of H_{2}O into OH (itself partly dependent on massloss rate), suggesting that the production of SO depends on the availability of OH to participate in the formation process.
We are only able to model SO_{2} in an additional three stars, IK Tau, W Hya, and R Cas, owing to the dearth of detections towards TX Cam. For W Hya we find an SO_{2} distribution similar to SO in abundance and envelope size. We have some difficulty fitting an SO_{2} model to observations for IK Tau and ultimately find an uncertain model which differs in shape from the SO distribution. For R Cas the SO_{2} model is also very uncertain because there are only two detected lines.
Overall, the circumstellar SO and SO_{2} abundances are much higher than predicted by chemical models of the extended stellar atmosphere. These two species may also account for all the available sulphur in the lower massloss rate stars. The Sbearing parent molecule appears not to be H_{2}S. The SO_{2} models for the higher massloss rate stars are less conclusive, but suggest an origin close to the star for this species. This is not consistent with present chemical models. The combined circumstellar SO and SO_{2} abundances are significantly lower than that of sulphur for these higher massloss rate objects.
To better constrain the behaviour of sulphur we need more observations of SO and SO_{2}, as well as other Sbearing species. Observations of a larger sample of stars will also allow us to confirm the trends we see in the SO abundance distributions.
The main evidence for this supposition is that it is in a vibrationally excited state, and that ΔK_{a,c} = 3 for this transition. Although this is an allowed transition, it is a very unlikely one under normal circumstances and, if included, our (nonmasering) model predicts almost no emission from this transition.
The Leiden Atomic and Molecular Database, found at http://home.strw.leidenuniv.nl/~moldata/
Acknowledgments
T.D. and K.J. acknowledge funding from the Swedish National Space Board. HO acknowledges financial support from the Swedish Research Council. This publication is based on data acquired with the Atacama Pathfinder Experiment (APEX). APEX is a collaboration between the MaxPlanckInstitut für Radioastronomie, the European Southern Observatory, and the Onsala Space Observatory. HIFI has been designed and built by a consortium of institutes and university departments from across Europe, Canada and the United States under the leadership of SRON Netherlands Institute for Space Research, Groningen, The Netherlands and with major contributions from Germany, France and the US. Consortium members are: Canada: CSA, U.Waterloo; France: CESR, LAB, LERMA, IRAM; Germany: KOSMA, MPIfR, MPS; Ireland, NUI Maynooth; Italy: ASI, IFSIINAF, Osservatorio Astrofisico di ArcetriINAF; Netherlands: SRON, TUD; Poland: CAMK, CBK; Spain: Observatorio Astronómico Nacional (IGN), Centro de Astrobiología (CSICINTA). Sweden: Chalmers University of Technology – MC2, RSS & GARD; Onsala Space Observatory; Swedish National Space Board, Stockholm University – Stockholm Observatory; Switzerland: ETH Zurich, FHNW; USA: Caltech, JPL, NHSC.
References
 Adande, G. R., Edwards, J. L., & Ziurys, L. M. 2013, ApJ, 778, 22 [NASA ADS] [CrossRef] [Google Scholar]
 Bujarrabal, V., Fuente, A., & Omont, A. 1994, A&A, 285, 247 [NASA ADS] [Google Scholar]
 Burkholder, J. B., Lovejoy, E. R., Hammer, P. D., Howard, C. J., & Mizushima, M. 1987, J. Mol. Spectr., 124, 379 [NASA ADS] [CrossRef] [Google Scholar]
 Cameron, A. G. W. 1973, Space Sci. Rev., 15, 121 [NASA ADS] [CrossRef] [Google Scholar]
 Cami, J., Yamamura, I., de Jong, T., et al. 1999, in The Universe as Seen by ISO, eds. P. Cox, & M. Kessler, ESA SP, 427, 281 [Google Scholar]
 Cernicharo, J., Spielfiedel, A., Balança, C., et al. 2011, A&A, 531, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cherchneff, I. 2006, A&A, 456, 1001 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Danilovich, T., Bergman, P., Justtanont, K., et al. 2014, A&A, 569, A76 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Danilovich, T., Teyssier, D., Justtanont, K., et al. 2015, A&A, 581, A60 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 de Graauw, T., Helmich, F. P., Phillips, T. G., et al. 2010, A&A, 518, L6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Decin, L., De Beck, E., Brünken, S., et al. 2010a, A&A, 516, A69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Decin, L., Justtanont, K., De Beck, E., et al. 2010b, A&A, 521, L4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dubernet, M.L., Daniel, F., Grosjean, A., & Lin, C. Y. 2009, A&A, 497, 911 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gobrecht, D., Cherchneff, I., Sarangi, A., Plane, J. M. C., & Bromley, S. T. 2016, A&A, 585, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gong, Y., Henkel, C., Spezzano, S., et al. 2015, A&A, 574, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Green, S. 1995, ApJS, 100, 213 [NASA ADS] [CrossRef] [Google Scholar]
 Guilloteau, S., Lucas, R., Omont, A., & NguyenQRieu. 1986, A&A, 165, L1 [NASA ADS] [Google Scholar]
 Justtanont, K., Khouri, T., Maercker, M., et al. 2012, A&A, 537, A144 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Justtanont, K., Barlow, M. J., Blommaert, J., et al. 2015, A&A, 578, A115 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kahane, C., GomezGonzalez, J., Cernicharo, J., & Guelin, M. 1988, A&A, 190, 167 [NASA ADS] [Google Scholar]
 Khouri, T., de Koter, A., Decin, L., et al. 2014a, A&A, 561, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Khouri, T., de Koter, A., Decin, L., et al. 2014b, A&A, 570, A67 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kim, H., Wyrowski, F., Menten, K. M., & Decin, L. 2010, A&A, 516, A68 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lindqvist, M., Nyman, L.A., Olofsson, H., & Winnberg, A. 1988, A&A, 205, L15 [NASA ADS] [Google Scholar]
 Lique, F., Spielfiedel, A., & Cernicharo, J. 2006, A&A, 451, 1125 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lodders, K. 2003, ApJ, 591, 1220 [NASA ADS] [CrossRef] [Google Scholar]
 Maercker, M., Schöier, F. L., Olofsson, H., Bergman, P., & Ramstedt, S. 2008, A&A, 479, 779 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Maercker, M., Schöier, F. L., Olofsson, H., et al. 2009, A&A, 494, 243 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Müller, H. S. P., Thorwirth, S., Roth, D. A., & Winnewisser, G. 2001, A&A, 370, L49 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Müller, H. S. P., Schlöder, F., Stutzki, J., & Winnewisser, G. 2005, J. Mol. Struct., 742, 215 [NASA ADS] [CrossRef] [Google Scholar]
 Netzer, N., & Knapp, G. R. 1987, ApJ, 323, 734 [NASA ADS] [CrossRef] [Google Scholar]
 Neufeld, D. A., Godard, B., Gerin, M., et al. 2015, A&A, 577, A49 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Olofsson, H., Lindqvist, M., Nyman, L.A., & Winnberg, A. 1998, A&A, 329, 1059 [NASA ADS] [Google Scholar]
 Omont, A., Lucas, R., Morris, M., & Guilloteau, S. 1993, A&A, 267, 490 [NASA ADS] [Google Scholar]
 Ott, S. 2010, in Astronomical Data Analysis Software and Systems XIX, eds. Y. Mizumoto, K.I. Morita, & M. Ohishi, ASP Conf. Ser., 434, 139 [Google Scholar]
 Peterson, K. A., & Woods, R. C. 1990, J. Chem. Phys., 93, 1876 [NASA ADS] [CrossRef] [Google Scholar]
 Ramstedt, S., Mohamed, S., Vlemmings, W. H. T., et al. 2014, A&A, 570, L14 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rothman, L. S., Gordon, I. E., Babikov, Y., et al. 2013, J. Quant. Spec. Radiat. Transf., 130, 4 [NASA ADS] [CrossRef] [Google Scholar]
 Rudolph, A. L., Fich, M., Bell, G. R., et al. 2006, ApJS, 162, 346 [NASA ADS] [CrossRef] [Google Scholar]
 Sahai, R., & Wannier, P. G. 1992, ApJ, 394, 320 [NASA ADS] [CrossRef] [Google Scholar]
 Schöier, F. L., van der Tak, F. F. S., van Dishoeck, E. F., & Black, J. H. 2005, A&A, 432, 369 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schöier, F. L., Bast, J., Olofsson, H., & Lindqvist, M. 2007, A&A, 473, 871 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schöier, F. L., Maercker, M., Justtanont, K., et al. 2011, A&A, 530, A83 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ukita, N., & Morris, M. 1983, A&A, 121, 15 [NASA ADS] [Google Scholar]
 Vassilev, V., Meledin, D., Lapkin, I., et al. 2008, A&A, 490, 1157 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Vlemmings, W. H. T., Humphreys, E. M. L., & FrancoHernández, R. 2011, ApJ, 728, 149 [NASA ADS] [CrossRef] [Google Scholar]
 Willacy, K., & Millar, T. J. 1997, A&A, 324, 237 [NASA ADS] [Google Scholar]
 Yamamura, I., de Jong, T., Onaka, T., Cami, J., & Waters, L. B. F. M. 1999, A&A, 341, L9 [NASA ADS] [Google Scholar]
 Yamamura, I., Kawaguchi, K., & Ridgway, S. T. 2000, ApJ, 528, L33 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: R Dor plots
Our best fit model lines for SO_{2} in R Dor are plotted along with the corresponding observations in Fig. A.1. For more details see Sect. 3.2.2.
The tentative detections of the isotopologues ^{34}SO_{2} and SO^{18}O from the APEX survey towards R Dor are plotted in Fig. A.2.
Fig. A.1 SO_{2} model (blue line) and observations (black histograms) of R Dor. In the case of overlapping lines, the top line listed is always the line centred at υ_{LSR} = 5.7 km s^{1}. 

Open with DEXTER 
Fig. A.1 continued. 

Open with DEXTER 
Fig. A.1 continued. 

Open with DEXTER 
Fig. A.1 continued. 

Open with DEXTER 
Fig. A.2 Detections of the isotopologues ^{34}SO_{2} and SO^{18}O towards R Dor. 

Open with DEXTER 
Appendix B: HIFI OBSIDs
The observation IDs for the HIFI observations used in this work are given in Table B.1, including the nondetections used to constrain our TX Cam SO model.
HIFI OBSIDs for observations used in this work.
Appendix C: Additional tables
SO_{2} observations towards R Dor using APEX, listed in order of descending energy of the upper level.
SO_{2} and SO observations using HIFI, listed by molecule in order of descending energy of the upper level.
Archival observations of SO and SO_{2} towards IK Tau, TX Cam, R Cas, and W Hya.
All Tables
SO observations towards R Dor using APEX, listed in order of descending energy of the upper level.
SO_{2} observations towards R Dor using APEX, listed in order of descending energy of the upper level.
SO_{2} and SO observations using HIFI, listed by molecule in order of descending energy of the upper level.
Archival observations of SO and SO_{2} towards IK Tau, TX Cam, R Cas, and W Hya.
All Figures
Fig. 1 SO energy level diagram with levels labelled using the N_{J} convention. Transitions detected with HIFI and APEX towards R Dor are indicated in light blue and green, respectively. 

Open with DEXTER  
In the text 
Fig. 2 SO_{2} energy level diagram. Levels are labelled J_{Ka,Kc}. Transitions detected with HIFI and APEX towards R Dor are indicated in light blue and green, respectively. 

Open with DEXTER  
In the text 
Fig. 3 SO χ^{2} plots for R Dor, W Hya, IK Tau and R Cas. The contours show the confidence intervals and the shading represents the value for the corresponding model, with the colourbar indicating multiples of the minimum value. The white cross indicates our bestfit model (see Table 6). For IK Tau, the slice for which R_{w} = 1.8R_{p} is shown. For R Cas, the slice for which R_{w} = 1.0R_{p} is shown. 

Open with DEXTER  
In the text 
Fig. 4 SO models (blue lines) and observations (black histograms) for R Dor. 

Open with DEXTER  
In the text 
Fig. 5 SO goodness of fit plots for R Dor, R Cas, IK Tau, and W Hya. New HIFI lines as well as archival data listed in Table C.3 are included. The green points in the R Dor plots represent the observations from the APEX spectral survey. Undetected HIFI lines are shown as cyan points with arrows, in this case representing lower limits because the vertical axis is the ratio of model integrated intensities to observed integrated intensities or the upper limits thereof. 

Open with DEXTER  
In the text 
Fig. 6 SO_{2} lines excluded from modelling for R Dor. See text for full explanation. 

Open with DEXTER  
In the text 
Fig. 7 SO_{2} goodness of fit plots for R Dor. Top: goodness of fit with upper energy level of the transition. HIFI lines are shown as blue points and APEX lines are shown as green crosses. Error bars are excluded to make the plot clearer to read. Lower left: goodness of fit with J. Lower right: goodness of fit with K_{a}, a clear downwards trend for K_{a} ≥ 6. 

Open with DEXTER  
In the text 
Fig. 8 Abundance profiles for R Dor, W Hya, IK Tau and R Cas. The abundances for CO and H_{2}O are taken from Maercker et al. (in prep.), except for W Hya, for which they are taken from Khouri et al. (2014a,b). The dashed line for the SO_{2} results for IK Tau and R Cas indicates that they are tentative. 

Open with DEXTER  
In the text 
Fig. 9 ^{34}SO model (blue lines) and observations (black histograms) for R Dor. 

Open with DEXTER  
In the text 
Fig. 10 Models (blue lines) and observations (black histograms) for SO towards W Hya. 

Open with DEXTER  
In the text 
Fig. 11 Models (blue lines) and observations (black histograms) for SO_{2} towards W Hya. 

Open with DEXTER  
In the text 
Fig. 12 SO models (blue lines) and observations (black histograms) for IK Tau. 

Open with DEXTER  
In the text 
Fig. 13 SO_{2} model (blue line) and observations (black histograms) of IK Tau. For details on the archival observations, see Table C.3. 

Open with DEXTER  
In the text 
Fig. 14 SO models (blue lines) and observations (black histograms) for R Cas. 

Open with DEXTER  
In the text 
Fig. 15 SO_{2} model (blue line) and observation (black histogram) for R Cas. 

Open with DEXTER  
In the text 
Fig. 16 SO models (blue lines) and observations (black histograms) for TX Cam. We only plot the archival SO detections, not the nondetections from HIFI. 

Open with DEXTER  
In the text 
Fig. 17 SO abundance distributions for all stars modelled. The vertical lines represent the dust condensation radii, where our models stop. The thicker sections of the curves represent the area probed by our observations and for TX Cam the thick dashed line is the area probed by the upper limits imposed by the HIFI nondetections. 

Open with DEXTER  
In the text 
Fig. 18 Trends in SO peak radius (top) and fractional abundance (bottom) against the circumstellar density measure, Ṁ/υ_{∞}. The green lines are the trends fitted to the three higher massloss rate stars (R Cas, TX Cam and IK Tau), while the yellow stars are the predicted locations of the lower massloss rate stars (R Dor and W Hya) based on the trend. In the fractional abundance plot, the red line shows the hard limit for SO abundance based on solar S abundance and the blue points in line with the yellow stars represent the real abundance values for R Dor and W Hya. 

Open with DEXTER  
In the text 
Fig. 19 Peak abundance radii of SO (for R Cas, TX Cam, and IK Tau) and efolding radii of SO (for R Dor and W Hya), plotted against the efolding radii of H_{2}O found by Maercker et al. (in prep.) and Khouri et al. (2014b, for W Hya). The solid green line is the best fit to the data and the dashed red line traces a 1:1 relationship. 

Open with DEXTER  
In the text 
Fig. A.1 SO_{2} model (blue line) and observations (black histograms) of R Dor. In the case of overlapping lines, the top line listed is always the line centred at υ_{LSR} = 5.7 km s^{1}. 

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

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

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

Open with DEXTER  
In the text 
Fig. A.2 Detections of the isotopologues ^{34}SO_{2} and SO^{18}O towards R Dor. 

Open with DEXTER  
In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.