c ESO 2013 Astronomy & Astrophysics manuscript no. 8680maer June 16, 2013 Circumstellar water vapour in M-type AGB stars: Radiative transfer models, abundances and predictions for HIFI M. Maercker1 , F. L. Schöier1,2 , H. Olofsson1,2 , P. Bergman2,3 and S. Ramstedt1 1 2 arXiv:0801.0971v3 [astro-ph] 14 Apr 2010 3 Stockholm Observatory, AlbaNova University Center, SE-106 91 Stockholm, Sweden e-mail: maercker@astro.su.se Onsala Space Observatory, SE-439 92 Onsala, Sweden European Southern Observatory, Casilla 19001, Santiago 19, Chile received 15 September 2007; accepted 12 December 2007 ABSTRACT Context. Surprisingly high amounts of H2 O have recently been reported in the circumstellar envelope around the M-type asymptotic giant branch star W Hya. This has lead to the speculation that evaporation of icy cometary or planetary bodies might be an effective ongoing mechanism in such systems. However, substantial uncertainties remain, as the required radiative transfer modelling is difficult due to high optical depths, sub-thermal excitation and the sensitivity to the combined radiation field from the central star and dust grains. Aims. By performing a detailed radiative transfer analysis, we determine fractional abundances of circumstellar H2 O in the envelopes around six M-type asymptotic giant branch stars. The models are also used to predict H2 O spectral line emission for the upcoming Herschel/HIFI mission. Methods. We use Infrared Space Observatory Long Wavelength Spectrometer spectra to constrain the circumstellar fractional abundance distribution of ortho-H2 O, using a non-local thermal equilibrium, and non-local, radiative transfer code based on the accelerated lambda iteration formalism. The mass-loss rates and kinetic temperature structures for the sample stars are determined through radiative transfer modelling of CO line emission based on the Monte-Carlo method. The density and temperature profiles of the circumstellar dust grains are determined through spectral energy distribution modelling using the publicly available code Dusty. Results. The determined ortho-H2 O abundances lie between 2 × 10−4 and 1.5 × 10−3 relative to H2 , with the exception of WX Psc, which has a much lower estimated ortho-H2 O abundance of only 2 × 10−6 , possibly indicating H2 O adsorption onto dust grains or recent mass-loss-rate modulations. The estimated abundances are uncertain by, at best, a factor of a few. Conclusions. The high water abundance found for the majority of the sources suggests that either the ‘normal’ chemical processes are very effective in producing H2 O, or else non-local thermal equilibrium atmospheric chemistry, grain surface reactions, or a release of H2 O (e.g. from icy bodies like Kuiper belt objects) play a role. However, more detailed information on the physical structure and the velocity field of the region where the water vapour lines are formed is required to improve abundance estimates. We provide predictions for ortho-H2 O lines in the spectral window of Herschel/HIFI. These spectrally resolved lines cover a wide range of excitation conditions and will provide valuable additional information on the physical and chemical properties of the inner stellar wind where H2 O is abundant. Key words. Stars: abundances - Stars: AGB and post-AGB - Stars: evolution - Stars: mass-loss 1. Introduction Asymptotic giant branch stars (AGB stars) are chemically characterised according to their relative abundances of carbon and oxygen. Carbon stars have a photospheric C/O ratio > 1, whereas M-type (‘oxygen rich’) AGB stars have C/O < 1 (e.g., Russel 1934; Beck et al. 1992). It is believed that M-type AGB stars evolve to carbon stars by the dredging up of nucleosynthesised carbon during the thermally pulsing AGB (TP-AGB) (e.g., Iben 1975; Boothroyd & Sackmann 1988). The evolution along the AGB is dominated by the mass loss, starting with a low massloss rate ∼ 10−8 M⊙ yr−1 and ending in an intense wind with rates up to 10−4 M⊙ yr−1 . The photosphere and expanding circumstellar envelope (CSE) are effective producers of a variety of molecular species and microscopic dust particles, providing up to 80% of the dust in galaxies (Whittet 1992). When the star reaches the end of the AGB, most of the initial mass has been lost, leaving behind only a C/O core that ionises the surrounding material and, Send offprint requests to: M. Maercker for a short while, creates a planetary nebula (Habing 1996). The lifetime on the AGB is estimated to be a few ∼ 106 yr (Marigo & Girardi 2007). To date, about 70 circumstellar molecular species have been detected, most of them through transitions at radio wavelengths (Olofsson 2007). H2 O is of special importance in the case of Mtype AGB stars, as it is possibly the dominant coolant in the inner part of the stellar wind (where it is abundant) with a large number of far-infrared (FIR) transitions (Truong-Bach et al. 1999). H2 O line emission is also a potentially important probe of the physical and chemical properties of the inner regions of CSEs. H2 O is the main reservoir of circumstellar O, and whereas CO is abundant in all AGB stars, chemical models predict H2 O to be one of the most abundant species only in O-rich stars (e.g., Cherchneff 2006). However, H2 O has also been detected in the nearby Crich AGB star IRC+10216 (Melnick et al. 2001; Hasegawa et al. 2006), but here the estimated abundance is comparatively low, and the origin of H2 O uncertain. 2 M. Maercker et al.: Circumstellar water vapour in oxygen-rich AGB stars The atmosphere of the Earth is mostly opaque at H2 O line wavelengths and it is therefore difficult to observe H2 O from the ground, in particular in its ground vibrational state. However, Menten et al. (2006) managed to observe two vibrationally excited lines in the ν2 = 1 bending mode (523 − 616 at 336.2 GHz and 661 − 752 at 293.6 GHz) towards the red supergiant VY CMa. H2 O masers from the ground and excited vibrational states can also be detected from the ground (Menten & Melnick 1989). Observations from the Infrared Space Observatory (ISO) are not hindered by the atmosphere, and the Short- and LongWavelength Spectrometer (SWS and LWS) data show a rich H2 O spectrum in M-type AGB stars (e.g., Truong-Bach et al. 1999; Barlow et al. 1996; Neufeld et al. 1996; Justtanont et al. 2005). The ortho-H2 O ground state line at 557 GHz was also observed by the SWAS (Melnick et al. 2000) and Odin satellites (Nordh et al. 2003). In contrast to the CO radio lines, the H2 O line emission originates from closer to the star, providing information on the warm, high-density innermost layers of the CSE (TruongBach et al. 1999). However, the high optical depths and the subthermal excitation of H2 O in CSEs make it a challenge to accurately model the observed H2 O spectral lines (e.g., Justtanont et al. 2005). We present here the results of radiative transfer modelling of the ortho-H2O line emission in six M-type AGB stars with mass-loss rates between 1 × 10−7 M⊙ yr−1 and 4 × 10−5 M⊙ yr−1 . Detailed radiative transfer models are created using the accelerated lambda iteration (ALI) method (Rybicki & Hummer 1991, 1992). The model results are fit to LWS ISO observations between 43 and 197 µm. The mass-loss rate and kinetic temperature structure are determined through CO radio line modelling, and the dust temperature and dust density structure are determined using Dusty (Ivezić et al. 1999). The results show the ability of ALI to model circumstellar H2 O lines. We also make predictions for H2 O lines in the frequency range of the upcoming Herschel/HIFI mission. These observations will include information on the line profiles and can be used to set firmer constraints on the physical structure of the inner CSE and the H2 O abundance. Previous attempts to determine the water vapour content in stars included in our sample have been made for W Hya (Barlow et al. 1996; Zubko & Elitzur 2000; Justtanont et al. 2005) and R Cas (Truong-Bach et al. 1999). González-Alfonso & Cernicharo (1999) used the large-velocitygradient (LVG) method to model the 183 GHz maser line in O-rich evolved stars, to study the line intensity dependence on source properties and physical conditions. Ryde & Eriksson (2002) model the 2.6 −3.6 µm spectrum for R Dor observed with the SWS on ISO, showing the dominance of water vapour in the spectrum. Determining the para-H2 O abundance would be a straightforward task. However, considering the uncertainty in the abundance estimates, this would not give any firm constraints on the ortho-to-para ratio or the total H2 O abundance. In Sect. 2 the ISO observations are briefly described. Section 3 contains a description of the dust and CO emission line modelling. In Sect. 4 we describe the H2 O line model, and in Sects. 5 and 6 we present and discuss the results, respectively. Appendix A contains a description of the accelerated lambda technique, adopted in our radiative transfer model, in more detail and benchmark tests are presented. All abundances mentioned refer to the fractional abundance compared to H2 , unless stated otherwise. Table 1. ISO observations for the 6 M-type AGB stars. Source TX Cam R Cas R Dor W Hya WX Psc IK Tau AOT LWS01 LWS01 LWS01 LWS01 LWS01 LWS01 Observer ID MBARLOW MBARLOW KERIKSSO MBARLOW MBARLOW MBARLOW Date 1997-10-10 1997-06-06 1997-07-01 1997-08-02 1997-06-15 1997-09-02 Int time (sec) 3428 2204 1330 1778 2796 3430 2. Observations We have used archived ISO LWS spectra (43 − 197 µm) for six M-type AGB stars (Table 1) to identify and measure the intensity of ortho-H2O spectral lines in the ground vibrational state. For all observations the LWS01 AOT was used (Clegg et al. 1996). The spectra were sampled at 1/4 of a resolution element, one element being 0.3 µm in second order (λ ≤ 93 µm) and 0.6 µm in first order (λ ≥ 80 µm). The analysis of the data was done using the ISO Spectral Analysis Package (ISAP). The data were cleaned from cosmic rays by excluding deviant points that were only present in individual scans. All scans for each detector were median averaged. In order to get a straight (continuumsubtracted) spectrum, a baseline (typically of order 2) was fit to the data and subtracted. Finally, the ten detectors were combined to form a complete spectrum. A section of the merged spectrum for R Dor with labelled ortho- and para-H2 O transitions is shown in Fig. 1. Table 2 gives the measured ortho-H2O line fluxes. The observed fluxes are determined by fitting a Gaussian profile to the observed ortho-H2O lines using available tools within ISAP. The error in the baseline fit dominates over the error in the Gaussian fit, and is ∼ 30% on average (in terms of the line intensity). The fit is more difficult in the short-wavelength part of the spectrum due to the high noise levels. The measured values lie within ± 20% of those of previous studies of H2 O emission using ISO spectra in W Hya (Barlow et al. 1996) and R Cas (Truong-Bach et al. 1999). 3. Circumstellar dust and CO line modelling 3.1. Dust emission modelling using Dusty The dust optical depth and dust temperature profile were determined using the continuum radiative transfer code Dusty (Ivezić et al. 1999). The observational constraints are in the form of spectral energy distributions (SEDs) and consist of 2MASS and IRAS fluxes. A ‘standard model’ containing amorphous silicate dust was used, with optical constants adopted from Justtanont & Tielens (1992). The standard model assumes the mass loss to be spherically symmetric with a constant rate and expansion velocity, resulting in a density profile with ρ ∝ r−2 . Prompt dust formation is assumed at the condensation temperature T c , which is treated as a free parameter. Additionally, the dust grains are assumed to be spherical with a radius of 0.1 µm. The dust optical depth and dust condensation temperature were adjusted in order to find the best fit to the SED using a χ2 -statistic. A more detailed description of the modelling procedure is given in Schöier et al. (2002). The luminosities of the Miras are derived using a periodluminosity (PL) relationship based on bolometric magnitudes (Feast et al. 1989). The distances are obtained from the apparent bolometric magnitude obtained in the SED modelling. For semi-regular objects, Knapp et al. (2003) give a PL relationship based on K-band magnitudes. Since this only uses information in M. Maercker et al.: Circumstellar water vapour in oxygen-rich AGB stars 3 Fig. 1. The ISO LWS spectrum of R Dor between 95 and 185 µm, including the position of ortho-H2O, para-H2 O lines, and CO lines (transitions J = 15 − 14 to J = 27 − 26) . Table 2. Measured and modelled ISO line fluxes of ortho-H2O [in units of 10−20 W cm−2 ]. Lines without fluxes were either too weak to detect or otherwise unusable, due to blending with other lines or unacceptable line fits. The dominant error in the observed fluxes is due to the baseline subtraction and is ∼ 30%. Trans JK− K+ − JK− K+ 221 − 212 212 − 101 303 − 212 330 − 321 514 − 505 423 − 414 432 − 423 734 − 643 414 − 303 221 − 110 505 − 414 616 − 505 423 − 312 321 − 212 707 − 616 330 − 221 532 − 423 523 − 414 λ [µm] 180.49 179.53 174.63 136.50 134.94 132.41 121.72 116.78 113.54 108.07 99.49 82.03 78.74 75.38 71.95 66.44 47.97 45.11 TX Cam F obs F mod — 10 10 — — — — — 9 11 10 — 14 23 — 18 — — — 7 6 — — — — — 8 11 10 — 13 23 — 22 — — R Cas F obs F mod — 8 7 3 2 4 4 — 14 16 16 — 25 28 — 35 — — — 9 8 3 2 3 4 — 10 14 17 — 17 28 — 25 — — R Dor F obs F mod 6 22 17 8 6 7 9 10 — 31 52 58 — 43 — 59 109 155 one point of the SED (the K-band magnitude), we instead used the revised Hipparcos distances for W Hya and R Dor (Knapp et al. 2003), and determined the luminosity through the SED modelling. 6 17 20 8 6 9 9 8 — 31 39 40 — 63 — 65 89 102 W Hya F obs F mod — 13 12 3 3 4 6 9 17 18 20 — — 35 24 26 52 70 7 8 4 4 4 5 7 11 15 16 — — 31 23 31 46 64 WX Psc F obs F mod — — — 1 — 2 1 — — — — — — — — 9 — — — — — 1 — 1 1 — — — — — — — — 12 — — IK Tau F obs F mod — 13 15 6 5 5 7 — 20 21 50 41 57 50 37 48 105 161 — 10 8 7 5 6 10 — 15 19 19 26 33 44 39 54 132 161 3.2. Modelling of CO lines The CO line modelling is based on the Monte-Carlo method presented in detail in Schöier & Olofsson (2001) and Olofsson et al. (2002). The radiative transfer model includes the lowest 40 rotational levels in the ground and first vibrational states. The collision rates are taken from Flower (2001) for CO-H2 , and they 4 M. Maercker et al.: Circumstellar water vapour in oxygen-rich AGB stars are extrapolated to higher rotational temperatures as was done in Schöier et al. (2005). The collisional rates between CO and ortho- and para-H2 were weighted assuming an ortho- to paraH2 ratio of 3. A CO fractional abundance (CO/H2 ) of 2 × 10−4 is adopted, a value typically used for O-rich AGB stars. The model includes cooling terms from the adiabatic expansion of the envelope and radiative line cooling through CO and H2 . The dominant heating term is given by dust-gas collisions. The density distribution and temperature structure of the dust grains are taken from the dust radiative transfer modelling (see Sect. 3.1), while the parameters for grain size, grain mass density and dust-to-gas ratio are combined into one parameter h. The h-parameter is, in principal, a free parameter in the CO modelling, and is chosen to be consistent with previous models (Olofsson et al. 2002: W Hya and R Dor; González Delgado et al. 2003: R Cas; Ramstedt et al. 2007: TX Cam, WX Psc, and IK Tau). The outer radius is set by the photodissociation of CO (Stanek et al. 1995, based on results by Mamon et al. 1988). The expansion velocity and micro-turbulent velocity (0.5 km s−1 doppler width), with local thermal (Gaussian) line broadening, are also set in the model. This leaves the mass-loss rate as the remaining free parameter that affects the model line fluxes. The best fit model is then determined by fitting the model line fluxes to the observed fluxes and minimising χ2 , where χ2 = X (Iobs,i − Imod,i )2 σ2i (1) and Iobs and Imod are the observed and modelled line fluxes, respectively, σi the uncertainty in observation i, and the sum goes over all lines included. The line profile is not taken into consideration, but we note that in the case of W Hya, Justtanont et al. (2005) had to reduce the CO photodissociation radius and increase the mass-loss rate in order to fit the line profile. However, there is sufficient uncertainty in the underlying assumptions of the circumstellar model that we prefer to treat all objects in the same way and avoid individual adjustments of, e.g., the CO envelope size. The uncertainty in the mass-loss rate is estimated to be ± 50% within the adopted model (Schöier & Olofsson 2001). Previous mass-loss-rate estimates by González Delgado et al. (2003) and Olofsson et al. (2002) differ somewhat from the results in Table 4. The difference is less than 10% for the lowmass-loss-rate and intermediate-mass-loss-rate objects, but more for the high-mass-loss-rate objects WX Psc and IK Tau (± 50%). This is due to the combination of new adopted distances, inclusion of new data, updated collisional rate coefficients and the inclusion of thermal radiation from dust grains in our CO models. The observed CO line fluxes, used in the model fitting, are presented in Table 3. IK Tau, WX Psc and TX Cam were modelled by Ramstedt et al. (2007). The parameters derived from the dust and CO modelling are presented in Table 4. 4. Circumstellar H2 O line modelling (ALI) We have developed a detailed radiative transfer code based on the accelerated lambda iteration method, in order to accurately model circumstellar H2 O line emission (Bergman, internal report). The Monte Carlo method is not well suited for such problems due to its slow convergence at high optical depths. The ALI code was previously used in Justtanont et al. (2005), but a number of modifications have subsequently been made to be able to treat the inclusion of dust emission and vibrationally excited energy levels. The accelerated lambda technique and its numerical Table 3. Observed CO line fluxes in main beam brightness scale. Observations were taken at Onsala Space Observatory, JCMT, and SEST (see references for details). Source TX Cam R Cas R Dor W Hya WX Psc IK Tau a J=1→0 [K km s−1 ] 20.1e 8.4a 5.0a 0.83a 44.4e 46.2e J =2→1 [K km s−1 ] 60.6e 32.1b 41.6a 18.7c 80.6e 95.1e J =3→2 [K km s−1 ] 71.1e 100.2b 70.4a 40.2d 70.1e 124.5e J=4→3 [K km s−1 ] 148.7e 89.6b 44.6c 83.8e 129.6e Kerschbaum & Olofsson (1999), b Delgado et al. (2003), c JCMT archive, d Justtanont et al. (2005), e Ramstedt et al. (2007) implementation is described in more detail in appendix A, including results of benchmarking tests against the standard problems presented in van Zadelhoff et al. (2002). 4.1. Dependence on numerical parameters ALI calculates the mean intensity at each radial point in the CSE, including the contribution from all other points in the envelope. This is done by paying special attention to a proper sampling of the central source. Typically, 50 radial points (distributed in a logarithmic way) and 32 angular points were used. The dependence on the radial and angular sampling was tested for the low mass-loss-rate object R Dor and intermediate mass-loss-rate object TX Cam by increasing the two parameters by a factor of two. The differences in the modelled line fluxes were small in both cases, less than 2%. After each iteration the new level populations are calculated, using a user-defined fraction of the old level populations (here set to 0.8) to prevent the models from diverging too easily. The models are said to converge when the average change in the level populations between iterations is less than a defined limit. This limit was also tested by running models with different limits for R Dor and IK Tau. For an average change between iterations in the level populations of less than 10−3 the models changed by less than 0.1% This limit was subsequently used for all models. 4.2. Modelling of H2 O lines The radiative transfer model includes the 45 lowest levels in the ground and first excited vibrational states (the bending mode ν2 = 1, 6.3 µm or ≈ 2300 K above the ground state) of orthoH2 O, while excitation through the stretching modes can be neglected. For the high mass-loss-rate objects IK Tau and WX Psc, only the ground vibrational state was included in order to improve the convergence. The effect of excluding the ν2 state on the resulting model line fluxes was tested and found to be insignificant (≤ 4%) for these high-mass-loss-rate objects (see Sect. 5.3). The collisional rates in the ground state are taken from the rates of H2 O with He, corrected by a factor of 1.4 to account for collisions with H2 (Green et al. 1993). Phillips et al. (1996) suggest that this may lead to rather different rates compared to those obtained when including H2 in the calculation. However, they only calculate collisional rates between H2 O and H2 up to a temperature of 140 K, while typical temperatures in the CSEs of AGB stars can be up to 1000 K. The importance of collisions within the ν2 state and between the ground and ν2 state was also tested (see Sect. 5.3). Rotational collision rates within the first excited M. Maercker et al.: Circumstellar water vapour in oxygen-rich AGB stars 5 5. Results 5.1. Modelled lines and ortho-H2 O abundances Fig. 2. Modelled ISO ortho-H2O lines for TX Cam. The dotdashed line marks the stellar vLSR , the solid line shows vLSR ± vexp . Self-absorption and weak P-Cygni profiles are present in the model line profiles. state are taken to be the same as for the ground state, while collisions between the ground and excited ν2 state are based on the ground state rotational collision rate coefficients scaled by a factor of 0.01. The molecular data for H2 O was taken from the Leiden Atomic and Molecular Database1 (LAMDA) (Schöier et al. 2005). As with the CO models, the density distribution and temperature structure of the dust grains are taken from the dust radiative transfer modelling (see Sect. 3.1). From the CO line modelling the mass-loss rate, expansion velocity, and kinetic temperature profile are included (see Sect. 3.2 and Table 4). The inner radius of the CSE is set to the inner radius of the dust envelope, i.e. the dust condensation radius. Setting a smaller radius does not affect the line fluxes. For all models, a Gaussian distribution of circumstellar H2 O is assumed. The maximum radius is defined by the e-folding radius re (i.e., where the abundance has decreased to 37%) and is set to 3 × re . By varying re , the size of the circumstellar H2 O envelope is treated as a free parameter. An estimate of re is given by the results of theoretical models of H2 O photodissociation (Netzer & Knapp 1987), where the radius at which the abundance of OH (the photodissociation product of H2 O) is highest is given by rNK = 5.4 × 1016 Ṁ 10−5 M⊙ yr−1 !0.7  vexp −0.4 cm. km s−1 (2) The remaining parameter to fit is the ortho-H2O abundance f o−H2 O relative to H2 . The model lines are fit to the line fluxes of the observed spectra in the same way as for the CO emission line models. The best fit model was determined, as for the CO modelling, by minimising χ2 as defined in Eq. 1, with re and fo−H2 O as free parameters. 1 available at http://www.strw.leidenuniv.nl/∼moldata For each object a grid, of models with varying ortho-H2O abundance and H2 O envelope size was calculated. The resulting χ2 maps are shown in Fig. 3. Between 4 and 16 lines are matched for each object. The fluxes of the best-fit modelled ortho-H2O lines are shown in Table 2, together with the observed fluxes. For the low-mass-loss-rate objects the abundance and re can be reasonably well determined, while setting constraints for re becomes increasingly difficult with increasing mass-loss-rate. In the high-mass-loss-rate objects, the emitting region is excitation limited. Consequently, the emission is not very sensitive to the envelope size. As a result, it is difficult to constrain the size of the of the H2 O envelope for these objects. For the low-mass-lossrate objects, the excitation can be high throughout the envelope, thus making the emission sensitive to the envelope size, i.e. the emission can be said to be photodissociation limited. The grids are terminated in envelope size when an increase in re does not lead to a significant change in line fluxes anymore, or when the 1σ contour closes. For IK Tau and WX Psc a lower limit could not be determined, as smaller envelope sizes require a higher abundance, for which the models could not be made to converge. The χ2 maps are truncated at ≈ 0.1 × rNK for these objects. In the cases where an outer radius could be determined, the radius given by Eq. 2 lies within the 1 σ contour of the lowest χ2 model radius. Table 4 gives the radius from the theoretical model (rNK ), the radius of the lowest χ2 model (re ), and the corresponding abundance estimates ( f NK and f e , respectively). From here on, the ‘best-fit’ models refer to the results obtained using rNK from Eq. 2. The resulting abundances of circumstellar ortho-H2O, from fitting the model lines to the observed ISO lines, lie between 2 × 10−6 and 1.5 × 10−3 relative to H2 (Table 4). Four of the modelled ISO lines of TX Cam are shown in Fig. 2 as an example of an intermediate-mass-loss-rate object. Self absorption is present in all the lines, whereas weak P-Cygni profiles are seen in the higher-frequency lines, where the continuum emission from thermal dust grains is stronger. Since the line shapes are not spectrally resolved in the ISO spectra, only the fluxes were fit. HIFI will be able to resolve the line profiles, and these can then be taken into consideration when fitting the models (Sect. 5.5). Figure 4 shows the kinetic temperature profile compared to the excitation temperature for four lines for W Hya and IK Tau (upper panels), and the tangential optical depth at line centre for the same lines as a function of distance from the central star (lower panels). It should be pointed out that the optical depth in the line wings usually is considerably higher. Two of the lines are observed by ISO, while three of the lines will be observable by HIFI. In general, the lines are subthermally excited, and the JK− K+ − JK− K+ = 532 − 441 transition at 620.9 GHz shows maser action. The lines are more thermalised in the high-massloss-rate object IK Tau. The tangential optical depth shows the regions where the lines are created, and this indicates that the ISO observations do indeed probe the inner regions of the H2 O envelope. W Hya has a smaller envelope due to the lower massloss-rate and consequently the lines are formed closer to the star. 5.2. Constraining the H2 O abundance distribution Several test models were calculated in order to see how the results depend on the ortho-H2 O abundance, mass-loss-rate, and 6 M. Maercker et al.: Circumstellar water vapour in oxygen-rich AGB stars Fig. 3. χ2 maps in ortho-H2O abundance and radius of the H2 O envelope for all stars in the sample. The contours give the 1, 2 and 3σ limits. The stars mark the best model based on the theoretical radius given by Eq. 2, the dots mark the models with the lowest χ2 values. The χ2red values in the figures refer to the best model fits. Table 4. Results of the dust emission, and CO and H2 O line modelling (see text for details). rNK and re are the theoretical model radius and the e-folding radius, respectively. The corresponding H2 O abundances for both radii are given under f N K and f e , respectively. The χ2red values refer to the best H2 O model fits using rNK . Likewise, diff gives the average difference of the modelled H2 O line fluxes using rNK to the observations in %. a b Source Type TX Cam R Cas R Dor W Hya WX Psc IK Tau M M SRb SRa M M La [L⊙ ] 11900 8725 6500 5400 14600 10400 Db [pc] 440 172 59 78 720 300 τdust 10 µm 0.4 0.05 0.03 0.07 3.0 1.2 Ṁ [10−6 M⊙ yr−1 ] 7.0 0.9 0.2 0.1 40 10 vexp [kms−1 ] 18.5 10.5 6.0 7.2 19.3 19.0 rNK [1015 cm] 13.1 3.9 1.7 1.1 43.0 17.0 re [1015 cm] 20.0 3.0 1.1 2.2 - f NK [10−4 ] 3.0 3.0 2.0 15.0 0.02 3.5 fe [10−4 ] 3.0 3.5 3.0 15.0 - χ2red 0.6 0.3 0.7 0.6 2.7 1.2 diff % -6.8 -4.8 -1.9 -10.4 -7.3 -8.5 From period-luminosity relation (Feast et al. 1989), except for R Dor and W Hya, where the Hipparcos distances and the dust modelling are used. Determined in the dust modelling using L, except for R Dor and W Hya, where the Hipparcos distances are used (Knapp et al. 2003). size of the circumstellar envelope. This was done by varying one of the parameters, while all other parameters were held constant. Table 5 shows the results in the cases of R Dor and TX Cam, as examples of low and intermediate-mass-loss-rate objects. Increasing and decreasing the abundance by a factor of two changes the line fluxes on average by ±20% and ±10% for the low- and high-mass-loss-rate objects, respectively (i.e. less than the observational error in the ISO lines). In general, the error in the H2 O abundances is estimated to be of the order ±50% within the adopted model. Varying the mass-loss rate by ±50% changes the line fluxes by ≈ ± 30%. The change in mass-loss rate affects the density distribution, and hence the efficiency of the collisional excitation, making the lines more sensitive for changes in the mass-loss rate than for changes in the abundance. Increasing and decreasing re by a factor of two changes the line fluxes by less than ±20%. The emitting region is excitation limited, instead of limited by the size of the H2 O envelope, making it difficult to provide good constraints on re . The position of OH masers for W Hya (0.8 − 1.4 × 1015 cm from the central star, see Szymczak et al. (1998)) and for R Cas (1.3 − 2.6 × 1015 cm from the central star, see Chapman et al. (1994)), indicate that the H2 O envelope might actually be smaller than given by Eq. 2. We conclude that the saturation of the lines, and hence the insensitivity of the integrated intensity to the abundance, is the major source of uncertainty in the determined abundances. Taking into account the uncertainty in the mass-loss rate (Ramstedt et al. 2007), the line fluxes, and in the model-dependent assumptions (such as the molecular distribution and velocity field), the determined abundances are accurate to within a factor of approximately 5. 5.3. Dependences on ‘fixed’ parameters Additional test models were calculated to examine the dependence of the line fluxes on the ‘fixed’ parameters (see Table 5). As before, only one parameter was varied, while all other parameters were held fixed. The effect of excitation to the first M. Maercker et al.: Circumstellar water vapour in oxygen-rich AGB stars 7 Fig. 4. Kinetic temperature profile and excitation temperatures for four modelled lines (indicated by their JK− K+ − JK− K+ quantum numbers) in W Hya and IK Tau (upper panels), and the tangential optical depths (at line centre) for the same lines (lower panels). The lines include a high transition line observed with ISO (523 − 414 at 45 µm), a low lying transition observed with ISO (303 − 212 at 174.6 µm), and the 557 GHz line (110 − 101 ) observed by Odin and SWAS (the latter two can also be observed by HIFI). A line, exhibiting maser action, which can be observed by HIFI is shown (532 − 441 at 621 GHz); the dashed line indicates the region where the level populations are inverted, and shows the absolute value of the optical depth in the lower panel. For radii larger than 3 × 1015 cm, numerical noise dominates the 532 − 441 transition in IK Tau and the model is consequently terminated at this radius. excited vibrational state was tested for the four stars for which models including this state converged. Excluding excitation to this state leads to a stronger decrease of the line fluxes for the low-mass-loss rate objects than for the intermediate and highmass-loss rate objects. Hence the error introduced by not including the first excited vibrational state for IK Tau and WX Psc is expected to be well below the uncertainty due to the observed line fluxes. For the low-mass-loss-rate objects, the inclusion of vibrational excitation results in significantly lower abundances compared to models that do not include the excited vibrational state (Justtanont et al. 2005). The variation of the line fluxes when excluding the dust radiation field and the central star was also examined in test models for the low-mass-loss-rate object R Dor and intermediate-mass-loss-rate object TX Cam (Table 5). It must be kept in mind that the changes in line fluxes all lie within (or close to) the observational uncertainties. This is especially true for TX Cam, where the change between different models is at most 12% on average. Including the rotational collision rates within the ν2 = 1 state and between the ground state and first excited state, as described in Sect. 4.2, does not change the fluxes significantly for any of the stars, showing that the dust radiation field and radiation from the central star are more important for the excitation of the ν2 = 1 state than collisions. The role of collisions within the vibrational ground state is tested by setting all collisional coefficients to zero, decreasing the observed line fluxes, on average, by 17% in the case of R Dor and 58% for TX Cam. As the lines are generally subthermally excited, the error in the collision rates in the ground state is likely to affect the modelled line fluxes (see Sect. 4.2). Since the ISO H2 O emission lines probe the inner parts of the CSE, the expansion velocity in the region where the lines are created may not have reached its terminal value as determined in the CO line modelling (see Sect. 3.2). [Note that changing the expansion velocity, while keeping the mass-loss-rate constant, effectively changes the density profile.] Decreasing the expansion velocity by 20% only changes the line fluxes by +6% for both TX Cam and R Dor, showing that the uncertainty in the expansion velocity does not contribute significantly to the error in the estimated abundances. 5.4. H2 O line cooling In addition to CO and H2 cooling by line emission, H2 O could also be an important coolant in the inner part of CSEs where it is 8 M. Maercker et al.: Circumstellar water vapour in oxygen-rich AGB stars Table 5. Parameter dependence of the model ISO line fluxes for the low-mass-loss-rate object R Dor and the intermediate-massloss-rate object TX Cam. Column 2 gives the change of the various parameters relative to their best-fit value (Table 4). Columns 3 and 4 give the average difference of all model line fluxes compared to the best-fit model. param. fe Rmod Ṁ vexp T kin collisions ν2 dust star change +100% –50% +100% –50% +50% –50% –20% ±20% excl. excl. excl. excl. R Dor +20% –20% +20% –20% +30% –30% +6% ±15% –17% –35% –11% –19% TX Cam +10% –10% +20% –20% +30% –30% +6% ±15% –58% –4% –12% –3% abundant. From the excitation analysis, the cooling (or heating) arising from line emission can be estimated using XX Λ = n H2 (nl γlu − nu γul )hνul , (3) l u>l where nu and nl are the level populations in the upper and lower levels participating in the transition (with rest frequency νul ) and γul and γlu are the collisional rate coefficients (van der Tak et al. 2007). The cooling rate Λ in erg s−1 cm−3 is defined as positive for net cooling. In Fig. 5 the line cooling from o-H2 O is compared with that from CO and H2 for the W Hya and IK Tau models used in the present analysis. In the region where H2 O is abundant its cooling rate will be more than an order of magnitude higher than that from CO. Due to the constraints set by the CO-line intensity ratios, the kinetic temperature profile determined through the CO line modelling is most accurate in the region in which the CO lines are emitted. Any major cooling due to H2 O must therefore be partly compensated by additional heating, by adjusting the h-parameter so that the kinetic temperature structure is consistent with the observed CO rotational line intensities. To include molecular cooling from water self-consistently, the energy balance equations for CO and H2 O must be solved simultaneously. For a complete description of the energy balance, other processes (such as a complicated heating term) would have to be included as well, which is beyond the scope of this paper. A crude estimate of how H2 O line cooling and other processes might affect the modelled line intensities is done by changing the temperature profile resulting from the CO modelling by ± 20%. This changes the line intensities on average by less than ±15%, without changing the line ratios. 5.5. Predictions for HIFI HIFI (Heterodyne instrument for the far infrared) is a heterodyne receiver onboard the Herschel space telescope, scheduled for launch in 2008. The instrument is a high-resolution (up to 107 , 0.03 − 300 km s−1) heterodyne instrument designed to provide a wide and continuous frequency coverage. The frequency range covered includes a large amount of molecular lines (de Graauw et al. 2005). We have modelled 19 ortho-H2O lines covering a range of excitation and optical depths within the range 557 GHz to 1885 GHz. Fig. 5. The emissivity Λ in erg s−1 cm−3 , scaled by (r/1014cm)4 , of the line cooling rate from o-H2 O, CO and H2 in the circumstellar envelopes around W Hya and IK Tau. There is a net cooling (Λ > 0) throughout the envelopes for these molecules, except very close to the star in the case of CO in W Hya (Λ < 0 indicated by a dashed line). The sensitivity limit of HIFI for a 5 σ line detection after 1 hour of integration at a resolution of R = 100000 is between 2.2 and 22.0 × 10−30 W cm−2 Hz−1 (de Graauw et al. 2005), depending on the frequency band. The modelled HIFI lines generally have fluxes of 10−27 W cm−2 Hz−1 or higher, making them observable without difficulty. Typically, between 1 and 10 hours of integration time is required to obtain a sufficient signal-to-noise ratio, even for weaker lines in the first excited vibrational state. Figure 7 shows a selection of the modelled lines obtained from the best-fit models for TX Cam, W Hya and IK Tau as given in Table 4 (based on the theoretical H2 O envelope radius given in Eq. 2). The dot-dashed line marks the stellar vLS R and the solid line marks ±vexp . Maser action is present in the 532 −441 and 634 −541 transitions at 621 GHz and 1158 GHz, respectively, in all three sources. The line profiles for models with lower mass-loss rates generally appear to suffer more from self-absorption and show significant emission outside the range vLS R ±vexp . The CSEs of these sources are relatively small and the optical depths at the outer parts of the CSEs are around, or larger than, unity (see Fig. 4). This leads to effective self-absorption on the blue-shifted side, but also to a contribution of the local line profiles to the line wings. This is particularly pronounced in the lines where the optical depths at line centre can reach up to ∼100 (Fig. 4) for these kinds of objects, such as 110 − 101 and 303 − 212 (Fig. 7). The broadening of the line profiles is further enhanced by the fact that the thermal motion, vth , is of the same order or larger when compared with the microturbulent motion, vturb , in the emitting region. The result is that optically thin lines can also have significant line broadening (Fig. 6). These effects are less significant in the high-mass-loss-rate objects, as the optical depths in the outer parts of these CSEs, which contribute most to the observed emission, are around, or below, one. Also, for the high-mass-loss-rate objects vth << vturb in the outer parts of the envelope. The sensitivity of the HIFI lines to changes in various parameters has also been tested. A change in abundance by ±50% changes the HIFI line fluxes on average by ≈ ±30% compared to the best-fit model. This compares to the 10 − 20% change in the ISO lines, indicating that the lines in the range of HIFI are somewhat more sensitive to the H2 O abundance (see below). Increasing and decreasing re by a factor of two change the line fluxes on average by less than ±16%. The line fluxes for the HIFI lines are more sensitive for a -20% change of the expansion velocity than the ISO lines (+20% on average; however, optically thin lines are more sensitive, see below). Resolved line profiles M. Maercker et al.: Circumstellar water vapour in oxygen-rich AGB stars Table 6. Model maximum intensities for 3 optically thin lines within the range of HIFI (in units of 10−30 W cm−2 Hz−1 ). Trans JK− K+ − JK− K+ 854 − 761 854 − 927 845 − 752 ν [GHz] 1169 1596 1885 R Dor TX Cam 10 3 320 4 1 40 Table 7. Parameter dependence of the model fluxes of three optically thin lines observable by HIFI for the low-mass-loss-rate object R Dor. Column 2 gives the change of the various parameters relative to their best-fit value (Table 4). The remaining columns give the change of the model fluxes compared to the best-fit fluxes for the respective lines. param. fe Rmod Ṁ vexp T kin change ±50% +100% –50% +50% –50% –20% +20% –20% 854 − 761 ±30% +4% –2% +77% –25% +28% +21% +9% 854 − 927 ±6% +5% –2% +10% –7% +5% +1% –1% 9 Table 8. Same as Table 7, but for the intermediate-mass-loss-rate object TX Cam. param. fe Rmod Ṁ vexp T kin change ±50% +100% –50% +50% –50% –20% +20% +20% 854 − 761 ±30% +20% –14% +70% –44% +48% +13% +13% 854 − 927 ±7% +3% –2% +14% –11% +8% +2% +2% 845 − 752 ±30% +6% –6% +71% –53% +40% +22% +22% 845 − 752 ±55% +2% –4% +153% –45% +55% +55% +20% will be crucial for getting information on the dynamical structure of the inner CSE. Of the 19 modelled lines, 5 are optically thin throughout the CSE. Two of these are maser lines (the 532 −441 and 634 −541 transitions at 621 GHz and 1158 GHz, respectively). The remaining three lines are shown in Table 6 with their predicted maximum intensities. The 5 σ line detection limit after 1 hour of integration time at R = 100000 is ≈ 2 × 10−29 W cm−2 Hz−1 , making it difficult to detect some of these lines. The sensitivity of these lines to the ortho-H2O abundance, the outer radius, mass-loss rate, kinetic temperature profile, and expansion velocity was tested for R Dor and TX Cam. The resulting change in line fluxes is presented in Tables 7 and 8 for R Dor and TX Cam, respectively. Except for the 854 − 927 transition, the optically thin lines are much more sensitive to the abundance, mass-loss rate and expansion velocity than the ISO lines are. The P-Cygni profile causes the insensitivity of the integrated flux in the 854 − 927 transition. The high sensitivity to the expansion velocity is due to the change in the density profile as a result of keeping the mass-loss rate constant while decreasing the expansion velocity. They are only slightly more sensitive to a change in the kinetic temperature profile compared to the ISO lines, and not at all sensitive to the outer radius. For R Dor, transitions 321 − 312 , 212 − 101 , and 221 − 212 at 1163 GHz, 1670 GHz, and 1661 GHz, respectively, are more sensitive to the change of the outer radius (with up to a 30% change in line flux), but the emission region is essentially excitation limited. An example of how the 845 − 752 line profile changes between the different models for R Dor (top) and TX Cam (bottom) is shown in Fig. 6. The left panels show models in which the parameters were reduced as indicated in Tables 7 and 8, and the right panels show the models with increased parameters. The line profiles change mainly in intensity. However, for TX Cam the optical depth in the 845 − 752 transition is just below one in the best fit model. The line therefore changes from optically thick to optically thin for the different models, and as the optical depth decreases, the line becomes increasingly flat-topped. The Fig. 6. Line profiles of the 845 − 752 transition for the various models indicated in Tables 7 and 8 for R Dor (top) and TX Cam (bottom). The left panels show the models with decreased parameters; the right panels show models with increased parameters.One parameter is varied, while the others are held fixed, thus leading to a change in the density profile when decreasing the expansion velocity. detailed abundance distribution is likely to have the largest effect on the shape of the line profile. Detailed examination of the predicted HIFI line profiles based on the spectrally resolved 557 GHz lines of R Cas, R Dor and W Hya, obtained with the Odin satellite, will be done in a forthcoming paper (Maercker et al., in prep.). Since this line is excited throughout the CSE, it will in particular help put better constraints on the size of the H2 O emitting region. 5.6. ISO CO-lines The ISO LWS spectrum also contains high-lying transitions of the CO molecule (Fig. 1). These CO lines probe the same region as the H2 O lines, while the radio lines probe the regions further out. Modelling these high excitation CO lines can be used as a consistency check for the adopted circumstellar model used in the ALI models. Due to weak lines and blending, in general only between one and three lines were found in the ISO spectra. W Hya, R Dor, and R Cas show the CO J = 16 → 15 and 10 M. Maercker et al.: Circumstellar water vapour in oxygen-rich AGB stars Fig. 7. Eight ortho-H2O lines that can be observed with Herschel/HIFI TX Cam (top two rows), W Hya (middle two rows) and for IK Tau (bottom two rows) models based on the theoretical radius from Eq. 2. The dot-dashed line marks the stellar vLSR , the solid line shows vLSR ± vexp . A number of line profile shapes are predicted, and the HIFI observations will help to verify or adjust the models. M. Maercker et al.: Circumstellar water vapour in oxygen-rich AGB stars J = 22 → 21 lines. R Cas also shows the CO J = 17 → 16 line. In the remaining objects only the CO J = 16 → 15 line could be found. For all objects except W Hya, the best-fit model fits the observations within ±20%. The modelled CO lines for W Hya are ∼ 80% weaker than the observations, in line with the higher mass-loss rate used by Justtanont et al. (2005) (see Sect. 3.2). 6. Discussion 6.1. Comparison with previous results Abundance estimates for two of the sources treated here have been made previously. For W Hya, Justtanont et al. (2005) find an ortho-H2O abundance of 1.0 ± 0.4 ×10−3 relative to H2 . In their paper they use an ALI code to determine the H2 O abundance by modelling of ISO SWS and LWS lines, and the Odin 557 GHz line. This agrees well with the abundance determined here (1.5 ± 0.75 × 10−3 ). Our higher abundance is due to the lower mass-loss rate used (see Sect. 3.2). Barlow et al. (1996) employ a large velocity gradient (LVG) model and determine a total H2 O/H2 of 8 × 10−4 in the inner part of the CSE (r ≤ 6 × 1014 cm) and 3 × 10−4 in the outer part (r > 6 × 1014 cm). Although the ortho-to-para ratio in their article is not very well constrained, between 1 and 3, their results are consistent with the abundances derived here. Zubko & Elitzur (2000) determine an abundance for ortho-H2O of 5.7 × 10−5 , one order of magnitude lower than the results presented here, and a mass-loss rate one order of magnitude higher (2.3 × 10−6 M⊙ yr−1 ). Their results are therefore consistent with ours, given the higher mass-loss rate. For R Cas, Truong-Bach et al. (1999) determine an ortho-H2O abundance of only 7.9 × 10−6 , two orders of magnitude lower than our result, using ISO LWS spectra. The mass-loss-rate used in their paper is somewhat lower than the one determined from the CO modelling here, further increasing the discrepancy in the abundance estimates. They do not model the dust emission explicitly, but instead use a power law to reproduce the far-infrared continuum spectrum. They also estimate cooling by CO emission by adopting a three level model, compared to the more complex scheme included here. On the other hand, they include a dust grain drift velocity (relative to the gas) that depends on the distance to the star. Including a velocity field in our models would probably result in somewhat lower abundances. Although the models are different, this is far from enough to explain the different results. 6.2. Evaporating planetary systems? Assuming an ortho-to-para ratio of 1 gives an estimate of the total H2 O abundance of ∼ 6 × 10−4 in the majority of our stars. [Justtanont et al. (2005) and Zubko & Elitzur (2000) found an ortho-to-para H2 O ratio of 1 and 0.8 for W Hya, respectively. Truong-Bach et al. (1999) derive a ratio of 3 for R Cas, in accordance with what is expected in LTE.] The estimated W Hya H2 O abundance is significantly higher than this, but, as outlined in Justtanont et al. (2005), the W Hya mass-loss rate could be significantly higher than the value we use, and consequently the H2 O abundance could be significantly lower. The cosmic abundance of carbon and oxygen (Anders & Grevesse 1989) limits the abundance of H2 O relative to H2 to not more than 1 × 10−3 , assuming that all the carbon goes into CO, and all the remaining O goes into H2 O. The abundance of H2 O in thermal equilibrium (TE) is approximately 1-3× 10−4 relative to H2 (Mamon et al. 1987; Nejad & Miller 1988; Willacy & Millar 1997; Duari et al. 1999). Taken at face value, our results indicate that the 11 abundances of H2 O observed in M-type AGB stars cannot be explained by TE chemistry. The uncertainty in our results is, however, at least a factor of a few, and therefore does not rule out a TE abundance. NLTE chemistry can affect the H2 O abundance, e.g., due to different densities and temperatures in shocked regions close to the star (Duari et al. 1999; Cherchneff 2006), and photodissociation further out in the envelope will strongly affect the abundance (Mamon et al. 1987; Nejad & Miller 1988; Willacy & Millar 1997). Other processes suggested to explain the relatively high water abundances are the evaporation of icy bodies (Justtanont et al. 2005) or grain surface reactions, such as Fischer-Tropsch catalysis on the surfaces of small metallic grains (Willacy 2004). The former process was first suggested to explain the surprising detection of H2 O in C-rich CSEs (Melnick et al. 2001; Hasegawa et al. 2006). Future observations with Herschel will make it possible to discriminate between the different formation scenarios, as suggested by e.g. González-Alfonso et al. (2007) in the case of Crich CSEs. A high abundance of H2 O in the inner regions of the CSE will be observable through a significant amount of midexcitation H2 O transitions. These are already observed in the ISO spectra for our objects, indicating that at least some of the observed H2 O is formed in the inner envelope. Both FischerTropsch catalysis and the evaporation of cometary bodies result in similar shell-like H2 O distributions. The evaporation hypothesis, however, implies the presence of HDO, the detection of which would rule out the formation of H2 O from material ejected by the star. 6.3. The very low abundance in WX Psc WX Psc has an estimated abundance that is two orders of magnitude lower than for the other objects. A possible explanation for this is depletion of H2 O on dust grains. For a sample of carbon and M-type stars, Schöier & Olofsson (2006) show that the abundance of SiO decreases in the CSEs with increasing Ṁ/vexp (a measure of the density), and explain this by depletion of SiO onto the dust grains. The freezing of water ice onto dust grains is also the explanation for the lack of H2 O lines in the spectra of a sample of OH/IR stars (Justtanont et al. 2006). Another possible explanation for the low abundance in WX Psc could be that the ISO water vapour observations are probing a present day massloss rate, whereas the CO transitions probe an earlier mass-loss epoch (Ramstedt et al. 2007). If the present day mass loss is much lower, this would require an increase in the H2 O abundance, thereby putting WX Psc more in line with the other objects. Mass-loss-rate modulations in WX Psc with a low present day mass-loss-rate have been suggested previously by Kemper et al. (2003) and are further discussed by Ramstedt et al. (2007). 7. Conclusions We have demonstrated the ability of our accelerated lambda iteration (ALI) code to model circumstellar H2 O emission lines in M-type AGB stars at optical depths > 1000, and ortho-H2O abundances are estimated. The models include the ground and first excited vibrational states, collisions within the ground vibrational state, the first excited vibrational state, and between these states, and excitation by the stellar and dust radiation fields. The sample of stars modelled spans one order of magnitude in mass-loss rate. The ortho-H2O abundances determined, in the range (0.02–15)×10−4, are estimated to be accurate to within 50%, within the adopted circumstellar model, while the absolute error (including model-dependent assumptions) is within a 12 M. Maercker et al.: Circumstellar water vapour in oxygen-rich AGB stars factor of 5. The estimated H2 O abundances are high in comparison to what can be produced by chemical models, but we explain that the estimates are rather uncertain. The outer radius of the H2 O envelope is better determined for the low-mass-loss-rate objects and corresponds to within 1σ to the radius given by the theoretical estimate in Eq. 2. Including excitation of the ground state through the vibrationally excited state affects the low-mass-loss-rate objects, increasing the model line fluxes by ≈ 40%. This results in a lower abundance for these objects compared to models that do not include the excited vibrational state (Justtanont et al. 2005). For higher mass-loss-rate objects the vibrational excitation is less important. Several formation processes for the observed H2 O abundances are possible. LTE or NLTE models may explain the formation of H2 O in the innermost layers of the CSE, whereas Fischer-Tropsch catalysis or the evaporation of Kuiper-belt like objects are possible sources of H2 O in the outer parts of the envelope. Future observations with the Herschel telescope will help to differentiate between these scenarios. The HIFI lines are generally less saturated and will therefor give a better estimate of the circumstellar H2 O abundance. In particular, optically thin lines, such as the 845 − 752 transition, are sensitive to the H2 O abundance, and will help constrain this parameter. The ground state line at 557 GHz, on the other hand, is excited throughout the envelope and will help in setting constraints on the size of the H2 O envelope. For WX Psc, we suggest that sticking of H2 O onto dust grains may explain the observed low amount of H2 O. Alternatively, the present day mass-loss rate is low. Acknowledgements. The authors would like to thank the anonymous referee for useful comments that helped improve the quality of the manuscript. The authors acknowledge the financial support from the Swedish Research Council and the Swedish National Space Board. References Anders, E., & Grevesse, N. 1989, Geochem. Cosmochim. Acta, 53, 197 Bains, I., Cohen, R.J., Louridas, A., et al. 2003, MNRAS, 342, 8 Barlow, M.J., Nguyen-Q-Rieu, Truong-Bach, et al. 1996, A&A, 315, L241 Beck, H.K.B., Gail, H.-P., Henkel, R., & Sedlmayr, E. 1992, A&A, 265, 626 Boothroyd, A.I., & Sackmann, I.J. 1988, ApJ, 328, 671 Chapman, J.M., Sivagnanam, P., Cohen, R.J., & Le Squeren, A.M. 1994, MNRAS, 268, 475 Cherchneff, I. 2006, A&A, 456, 1001 Clegg, P.E., Ade, P.A.R., Armand, C., et al. 1996, A&A, 315, L38 Colomer, F., Reid, M.J., Menten, K.M., & Bujarrabal, V. 2000, A&A, 355, 979 Deguchi, S. & Nguyen-Q-Rieu 1990, ApJ, 360, L27 Duari, D, Cherchneff, I., & Willacy, K. 1999, A&A, 341, L47 Feast, M.W., Glass, I.S., Whitelock, P.A., Catchpole, R.M. 1989, MNRAS, 241, 375 Flower, D.R. 2001, J. Phys. B.: At. Mol. Opt. Phys., 34, 2731 González-Alfonso, E., & Cernicharo, J. 1999, ApJ, 525, 845 González-Alfonso, E., Neufeld, D.A., & Melnick, G.J. 2007, A&A, accepted González Delgado, D., Olofsson, H., Kerschbaum, F., et al. 2003, A&A, 411, 123 de Graauw, Th., Caux, E., Guesten, R., et al. 2005, AAS, 37, 1219 Green, S., Maluendes, S., & McLean, A.D. 1993, ApJS, 85, 181 Habing, H.J. 1996, A&A Rv, 7, 97 Hasegawa, T.I., Kwok, A., Koning, N., et al. 2006, ApJ, 637, 791 Iben, I. 1975, ApJ, 196, 525 Ivezić, Z̆, Nenkova, M. & Elitzur, M. 1999, User Manual for DUSTY, (Univ. Kentucky Internal Rep.) Justtanont, K., & Tielens, A.G.G.M. 1992, ApJ, 389, 400 Justtanont, K., Bergman, P., Larsson, B., et al. 2005, A&A, 439, 627 Justtanont, K., Olofsson, G., Dijkstra, C., & Meyer, W. 2006, A&A, 450, 1051 Kemper, F., Stark, R., Justtanont, K., et al. 2003, A&A, 407, 609 Kerschbaum, F. & Olofsson, H. 1999, A&A Suppl. Ser., 138, 299 Knapp, G.R., Pourbaix, D., Platais, I., & Jorissen, A. 2003, A&A, 403, 993 Lane, A.P., Johnston, K.J., Bowers, P.F., Spencer, J.H., & Diamond, P.J. 1987, ApJ, 323, 756 Mamon, G.A., Glassgold, A.E., & Omont, A. 1987, ApJ, 323, 306 Mamon, G.A., Glassgold, A.E., & Huggins, P.J. 1988, ApJ, 328, 797 Marigo, P., & Girardi, L. 2007, A&A, 469, 239 Melnick, G.J., Stauffer, J.R., Ashby, M.L.N., et al. 2000, ApJ, 539, L77 Melnick, G.J., Neufeld, D.A., Saavik Ford, K.E., et al. 2001, Nature, 412, 160 Melnick, G.J., & Bergin, E.A. 2005, AdSpR, 36, 1027 Menten, K.M., & Melnick, G.J. 1989, ApJ, 341, L91 Menten, K.M., Philipp, S.D., Güsten, R., et al. 2006, A&A, 454, L107 Nejad, L.A.M., & Millar, T.J. 1988, MNRAS, 230, 79 Netzer, N., & Knapp, G.R. 1987, ApJ, 323, 734 Neufeld, D.A., Chen, W., Melnick, G.J., et al. 1996, A&A, 315, L237 Ng, K.C. 1974, J. Chem. Phys., 61, 2680 Nordh, H.L., von Schéle, F., Frisk, U., et al. 2003, A&A, 402, L21 Olofsson, H., González Delgado, D., Kerschbaum, F., & Schöier, F.L. 2002, A&A, 391, 1053 Olofsson H. 2007, in ”Science with ALMA: a new era for astrophysics”, in press Olson, G.L., Auer, L.H., & Buchler, J.R. 1986, J. Quant. Spectrosc. Radiat. Transfer, 35, 431 Phillips, T.R., Maluendes, S., & Green, S. 1996, ApJS, 107, 467 Ramstedt, S., Schöier, F.L., Olofsson, H., & Lundgren, A.A. 2007, submitted. Russel, H.N. 1934, ApJ, 79, 317 Ryde, N., & Eriksson, K. 2002, A&A, 386, 874 Rybicki, G.B., & Hummer, D.G. 1991, A&A, 245, 171 Rybicki, G.B., & Hummer, D.G. 1992, A&A, 262, 209 Schöier, F.L., & Olofsson, H. 2001, A&A, 368, 969 Schöier, F.L., & Olofsson, H. 2006, A&A, 454, 247 Schöier, F.L., Ryde, N., & Olofsson, H. 2002, A&A, 391, 557 Schöier, F.L., van der Tak, F.F.S., van Dishoeck, E.F., & Black, J.H. 2005, A&A, 432, 369 Stanek, K.Z., Knapp, G.R., Young, K., & Phillips, T. 1995, Ap&SS, 100, 169 Szymczak, M., Cohen, R.J., & Richards, A.M.S. 1998, MNRAS, 297, 1151 Truong-Bach, Sylvester, R.J., Barlow, M.J., et al. 1999, A&A, 345, 925 van der Tak, F.F.S., Black, J.H., Schöier, F.L., Jansen, D.J., & van Dishoeck, E.F. 2007, A&A, 468, 627 Vlemmings, W.H.T., van Langevelde, H.J., Diamond, P.J., Habing, H.J., & Schilizzi, R.T. 2003, A&A, 407, 213 Whittet, D.C.B. 1992, Dust in the Galactic Environment, IoP Publishing, Bristol Willacy, K., & Millar, T.J. 1997, A&A, 324, 237 Willacy, K. 2004, ApJ, 600, L87 van Zadelhoff, G.-J., Dullemond, C.P., van der Tak, F.F.S., et al. 2002, A&A , 395, 373 Zubko, V. & Elitzur, M. 2000, ApJ, 544, L137 M. Maercker et al.: Circumstellar water vapour in oxygen-rich AGB stars 13 Appendix A: Accelerated lambda iteration (ALI) A.1. General considerations The radiative transfer models used in this work are created using the accelerated lambda iteration (ALI) method (Rybicki & Hummer 1991, 1992). As for most NLTE codes, the ALI technique makes an initial guess of the level populations ni and solves the statistical equilibrium (SE) equations (for transition l → l′ ) X [nl All′ − (nl′ Bl′ l − nl Bll′ ) J¯ll′ ] − l′ <l X [nl′ Al′ l − (nl Bll′ − nl′ Bl′ l ) J¯ll′ ] + l′ >l X (nlCll′ − nl′ Cl′ l ) = 0, (A.1) l′ where All′ is the Einstein A-coefficient, Bll′ and Bl′ l are the Einstein B-coefficients for stimulated emission and absorption respectively, and Cll′ are the collisional rates. J¯ll′ is the mean integrated intensity Z Z 1 J¯ll′ = dΩ dνφll′ (µ, ν)I(µ, ν), (A.2) 4π where I(µ, ν) is the specific intensity at frequencey ν along direction µ, and φll′ (µ, ν) is the normalized line shape profile (Doppler profile) used as a weight function. In the spherically symmetric case it is enough to solve the SE equations by determining J¯ll′ at every radial point in the circumstellar shell, considering both the local and non-local contributions from all directions µ. This is done by introducing the lambda operator Λµν , which relates to the specific intensity I(µ, ν) as I(µ, ν) = Λµν S tot (µ, ν) + Lbg , (A.3) where S tot (µ, ν) is the total source function (including dust and all overlapping lines) and Lbg is the contribution from the background to I(µ, ν). The lambda operator can be seen as a N × N matrix with N radial points. The diagonal elements of the matrix describe the local contribution to the specific intensity, while the off-diagonal elements describe the contribution from all other parts of the cloud. If τ(µ, ν) is the local optical depth in direction µ, the diagonal elements become Λr,r (µ, ν) = 1 − e−τ(µ,ν) . (A.4) For high optical depths (τ(µ, ν) ≫ 1) Λr,r → 1, whereas in the optical thin case (τ(µ, ν) → 0) Λr,r → 0. Instead of applying the lambda iteration directly, the accelerated lambda iteration technique uses an approximate lambda operator Λ∗µν to take care of the high optical depths separately † I(µ, ν) = Λ∗µν S tot (µ, ν) + (Λµν − Λ∗µν )S tot (µ, ν). (A.5) S tot † (µ, ν) is the source function from the previous iteration and can be used, since convergence guarantees that S tot † → Stot . It is convenient to simply choose the diagonal elements of the lambda operator as the approximate lambda operator. This makes the matrix inversion easily possible, and the SE equations can be solved at every radial point individually, based on the approximate lambda operator only. Using the accelerated lambda iteration technique improves the rate of convergence compared to normal lambda iteration, particularly in cases of high optical depths. To further improve convergence, Ng acceleration (Ng 1974) is used, with a good description of the method given by Olson et al. (1986). Fig. A.1. Results from benchmarking our ALI code against the results found in van Zadelhoff et al. (2002) for their problems 2a and 2b, consisting of an inside-out collapsing spherical cloud. The figure shows the radial distributions of the relative number density of HCO+ molecules in the rotational levels J = 1 and J = 4 as obtained by our ALI code after convergence. Also shown are the deviations (in percent) of our results when compared to the mean of the results found by all other radiative transfer codes, as presented in van Zadelhoff et al. (2002). A.2. Benchmarking When developing a new radiative transfer code it is important to perform detailed tests as to its accuracy. A number of such test problems are available for download at www.strw.leidenuniv.nl/∼radtrans. Our ALI code has been benchmarked against the solutions to these problems as found by seven different radiative transfer codes (presented in van Zadelhoff et al. 2002). We find that our ALI code reproduces the results for all these test problems with an accuracy that matches all the other codes included in van Zadelhoff et al. (2002). As an example, Fig. A.1 shows the solution obtained from our ALI code for problems 2a and 2b of van Zadelhoff et al. (2002). These problems consist of an inside-out collapsing spherical cloud and contain radial gradients in most physical parameters such as velocity, density, temperature and microturbulent velocity. The excitation analysis is performed for HCO+ , and the final level populations (in relation to the total number density of HCO+ molecules, ntot ) for the J = 1 and J = 4 rotational levels are shown in Fig. A.1. Also shown are the deviations of our ALI code (in percent) from the mean of the results obtained by the other radiative transfer codes presented in van Zadelhoff et al. (2002). The spread in the results, around the mean, between various codes in van Zadelhoff et al. (2002) is for the J = 1 level about ± 1% and ± 2% for problems 2a and 2b, respectively. For the J = 4 level a larger spread of about ± 5% and ± 15% was obtained for problems 2a and 2b, respectively.