EDP Sciences
Free Access
Issue
A&A
Volume 591, July 2016
Article Number A44
Number of page(s) 18
Section Stellar structure and evolution
DOI https://doi.org/10.1051/0004-6361/201628310
Published online 09 June 2016

© ESO, 2016

1. Introduction

Asymptotic giant branch (AGB) stars play an important role in the chemical evolution of galaxies and the Universe (Busso et al. 1999; Schneider et al. 2014). Elements created inside the star get returned to interstellar space in the form of molecules and dust grains in a wind from the stellar surface. This mass loss is driven by radiation pressure on dust grains formed in the upper atmosphere, which drag the gas along through collisions (Bowen 1988; Hoefner et al. 1998; Simis et al. 2001; Höfner 2008). The lifetime on the AGB is determined through the mass loss.

In order to determine stellar elemental yields to interstellar space from AGB stars, and to understand their role in the chemical evolution of the Universe, it is important to determine the molecular content of the circumstellar envelope (CSE), the chemical reaction paths, and the characteristics of the mass-loss process. Observations of molecular emission lines have been used to determine these properties. In particular, observations of CO emission lines have proven to be a useful tool to describe the mass loss (e.g., Schöier et al. 2002; González Delgado et al. 2003; Ramstedt et al. 2008; De Beck et al. 2010; Khouri et al. 2014a). Observations of other molecules have led to important insights in the chemistry and structure of AGB CSEs (e.g., González Delgado et al. 2003; Schöier et al. 2007; Maercker et al. 2008, 2009; Decin et al. 2010a; Cernicharo et al. 2010, 2011; De Beck et al. 2012; Justtanont et al. 2012; Khouri et al. 2014b; Lombaert et al. 2016).

A molecule of particular importance for M-type AGB stars is H2O. It is one of the most abundant molecules in the winds of these stars, rivalled only by CO and H2 (Cherchneff 2006; Maercker et al. 2008, 2009; Decin et al. 2010a,c; Khouri et al. 2014b). As an oxygen-bearing molecule it plays an important role in the oxygen chemistry and the formation of other molecules. Due to its many far infrared transitions it is a dominant coolant in the inner CSE (Truong-Bach et al. 1999). The H2O lines are emitted in the warm, inner envelope, and observations of emission lines also probe the acceleration of the wind (Decin et al. 2010a,c). As such, observations of H2O serve as an excellent diagnostic for the chemistry and physics of AGB CSEs. At the same time, a correct treatment of H2O and its effects on the physical structure of the CSE is essential in all radiative-transfer calculations and mass-loss descriptions.

Observations of H2O lines are difficult to obtain from ground-based telescopes due to the absorption in Earth’s atmosphere (although it is possible to observe masers and vibrationally excited lines in the ν2 = 1 bending mode; Menten et al. 2006). The Infrared Space Observatory (ISO) observed a large number of H2O lines in the wavelength range between 43 μm and 197 μm towards AGB stars (e.g., Barlow et al. 1996; Truong-Bach et al. 1999; Maercker et al. 2008). Radiative-transfer models generally require relatively high amounts of H2O to fit the observations of H2O emission lines towards the CSEs of M-type AGB stars. However, the ISO observations are spectrally unresolved, increasing the uncertainty in the determined H2O abundances, and entirely eliminating any information on the velocity field of the H2O line emitting region. The SWAS (Melnick et al. 2000) and Odin (Nordh et al. 2003) satellites provided spectrally resolved observations of the ground-state ortho-H2O (110−101) transition at 557 GHz towards M-type AGB stars. Modelling the line-width of this transition further constrained the determined H2O abundances, and demonstrated that the H2O line emitting region extends into the wind acceleration zone (Justtanont et al. 2005; Hasegawa et al. 2006; Maercker et al. 2009).

The Heterodyne Instrument for the Far-infrared (HIFI) aboard the Herschel Space Observatory (Herschel) for the first time provided high-sensitivity and spectrally resolved observations of a large number of H2O lines towards AGB CSEs. The observations used here were done as part of the HIFISTARS guaranteed time key program (P.I.: V. Bujarrabal) to study the CSEs around evolved stars of all chemical types, mainly targeting the CO and H2O lines (see Sect. 2.2). In particular, HIFISTARS contained a sample of M-type AGB stars for which, in addition to CO and H2O, a total of seven different molecules were detected (Justtanont et al. 2012). The HIFI observations allow us to determine the abundance and distribution of H2O in the CSEs of the observed sources, as well as to constrain the velocity profile, through radiative-transfer modelling of the observed CO and H2O lines.

In this paper we use advanced radiative-transfer models in order to determine the abundance of H2O in the CSEs of M-type AGB stars, constrain the velocity profile of the stellar wind, as well as explore the significance of H2O line cooling in the energy balance of the inner CSE. The basic modelling method has been used previously for models of H2O lines (Maercker et al. 2008, 2009), and has been upgraded to also include velocity profiles and the effect of H2O line cooling on the energy balance (Schöier et al. 2011; Danilovich et al. 2014). In addition to the HIFI observations, we use spectrally unresolved lines of CO and H2O observed with the Photoconductor Array Camera and Spectrometer (PACS) aboard Herschel to further constrain the models. All abundances are relative to H2. In Sect. 2 we present the observations and basic information on the modelled M-type AGB stars. In Sect. 3 the radiative-transfer models are presented. We present and discuss the results in Sects. 4 and 5, respectively, and summarise our conclusions in Sect. 6.

2. Observations

We model rotational lines from the main isotopologues of CO and H2O from both ground-based and space telescopes. The observations have been published previously, and are summarised below.

2.1. Ground-based observations of CO

We included ground-based observations in the molecular line radiative modelling of CO. The observations were taken at Onsala Space Observatory, the JCMT, and SEST (for details, see Kerschbaum & Olofsson 1999; González Delgado et al. 2003; Justtanont et al. 2005; Ramstedt et al. 2008; Maercker et al. 2008, and Table A.1). Transitions from CO(1−0) to CO(4−3) were included for all sources except R Dor, which lacks CO(4−3). The uncertainty in the integrated intensities is assumed to be 20%.

2.2. HIFI observations of CO and H2O

The HIFI observations covered a number of expected strong CO and H2O lines in the frequency range 556.6–1843.5 GHz. For CO the CO(6−5), CO(10−9), and CO(16−15) were observed. The observed integrated intensities for CO and the observed H2O lines are presented in Tables A.1 to A.3. These observations were already presented in Justtanont et al. (2012), but were reprocessed using the latest main beam efficiencies. For some of the bands this led to an increase in the integrated fluxes of up to 20%. We assume an average uncertainty in the integrated intensity for all lines of 20%.

2.3. PACS observations of CO and H2O

The PACS observations were taken as part of the MESS guaranteed-time key project (Groenewegen et al. 2011), which covered a sample of AGB stars of all chemical types. The observations and the detailed data reduction will be presented in a forthcoming paper dedicated to the M-type AGB stars of the MESS sample (Decin et al., in prep.). The PACS spectra cover the entire frequency range from 1580 up to 5450 GHz. For this study, we only include non-blended CO and H2O emission lines. For a detailed description on how line blends are identified, we refer to Lombaert et al. (2016). The observed integrated intensities are presented in Tables A.1 to A.3. We assume an average uncertainty in the integrated intensity for all PACS lines of 20%.

2.4. Sources

The HIFISTARS program included nine M-type AGB stars. Of these we include only fours stars: R Dor, R Cas, TX Cam, and IK Tau. R Dor is a SRb variable, and the remaining sources are Mira-type variables. The HIFI observations of IK Tau were modelled by Decin et al. (2010c) using a different radiative-transfer code. We include IK Tau here for a comparison of the modelling approaches. The HIFI observations of W Hya were modelled in detail by Khouri et al. (2014a,b) using the same code as Decin et al. (2010c). Mira has a very complex circumstellar environment due to binary interaction (Ramstedt et al. 2014), and the remaining sources are OH/IR stars with high mass-loss rates. These sources are therefore not included in this study.

The mass-loss rates and H2O abundances of the remaining sources have been modelled previously (Maercker et al. 2008, 2009). The lowest and highest mass-loss rates of the sample were estimated to be 2 × 10-7M yr-1 (R Dor) and 1 × 10-5M yr-1 (IK Tau), respectively, based on ground-based CO observations only. Ortho-H2O abundances were estimated to be between 2.0 × 10-4 (R Dor) and 3.5 × 10-4 (IK Tau) relative to H2, based on ISO observations. For the Mira-variables R Cas, TX Cam and IK Tau the luminosities were derived from a period-luminosity relation (Feast et al. 1989) and the distances were derived in the dust modelling, while for the SRb variable R Dor we assumed a distance of 59 pc and derived the luminosity in the dust modelling. The same procedure was used in Maercker et al. (2008).

3. Radiative-transfer modelling

The radiative-transfer modelling is done by combining results from models of the dust radiation field, CO radiative-transfer models, and H2O radiative-transfer models. All models assume a spherically symmetric, homogenous wind produced by a constant mass-loss rate. The procedure allows us to determine the dust temperature profile, the mass-loss rate, the kinetic temperature profile, the gas expansion velocity profile, and the fractional H2O abundance (relative to H2) distribution throughout the CSE. Since the aim of the paper is to compare the derived H2O abundances between the sources, derive gas-velocity profiles, and investigate the significance of H2O line cooling, we treat the sources as homogenously as possibly. Therefore basic properties, such as the dust grain sizes, dust composition, and the initial wind velocity, are assumed to be the same for all sources. Although detailed modelling of these parameters for individual sources may result in different values, the uncertainties in these individual models would complicate the comparative study that is the goal of this work. This basic modelling strategy has been thoroughly tested and used previously (e.g., Schöier & Olofsson 2001; Olofsson et al. 2002; Ramstedt et al. 2008; Maercker et al. 2008, 2009; Schöier et al. 2011; Danilovich et al. 2014).

3.1. The basic CO model

The code used for the modelling of the CO lines is described in detail in Schöier & Olofsson (2001). It employs the Monte Carlo method and is a non-local, non-LTE radiative-transfer code that solves the energy balance and derives a self-consistent kinetic temperature profile. We include radiation from dust by calculating the spectral energy distribution using DUSTY (Ivezić et al. 1999), constrained by 2MASS and IRAS fluxes. We used amorphous silicate dust with optical constants from Justtanont & Tielens (1992), assuming spherical grains with a radius of 0.05 μm and a density of 3 g cm-3. The dust model calculates the dust condensation radius, the dust optical depth, τd, at 10 μm and gives a dust temperature profile, both of which are included in the CO model. Our dust models are largely based on the results presented in Maercker et al. (2008), with possible slight adjustments in the luminosity and/or dust optical depth. The radius of the CO envelope is determined by models of photodissociation of CO (Mamon et al. 1988). The code and modelling method has been used extensively to model the emission lines from various molecules, including CO (e.g., Olofsson et al. 2002), SiO (Ramstedt et al. 2009), and HCN (Schöier et al. 2013).

Models of CO lines up to the J = 4−3 transition were used to determine the basic properties of the CSE in the modelling of H2O lines observed by ISO and Odin (Maercker et al. 2008, 2009). The emission of these lines is dominated by the outer parts of the envelope, where the wind has already reached its terminal velocity. HIFI provides new observations of higher-energy CO transitions. These are emitted in the inner CSE, and it is now possible to constrain the velocity profile. The models describe a velocity law for the gas and the dust separately. The gas-velocity profile follows the form (1)where υg,i is the gas expansion velocity at the inner radius (assumed to be 3 km s-1 in all models), υg, is the gas terminal velocity, Ri is the inner radius (set to the dust condensation radius), and β describes how fast the wind accelerates. The dust expansion velocity is connected to the gas expansion through the dust drift velocity υdr(r): (2)where υdr(r) follows the same shape as Eq. (1), with υdr,i the drift velocity at the inner radius. The terminal drift velocity υdr, is calculated by (3)where L is the stellar luminosity, υg, is the terminal gas expansion velocity, Q is the scattering efficiency assumed to be 3%, is the mass-loss rate, and c is the speed of light. The terminal dust velocity is then given by υd, = υdr, + υg,. Note that there are no constraints on the dust velocity profile, and we assume Eqs. (2) and (3) to give reasonable values. However, the exact shape of the drift velocity profile strongly affects the heating in the inner envelope (see Sect. 3.3).

The gas expansion velocity at the inner radius, υg,i, equals typical values of the sound-speed in the inner winds of AGB stars (e.g., Decin et al. 2006). The gas-velocity profile is constrained by fitting the widths of the CO lines at the zero-intensity level. Low-J lines constrain the terminal velocity, while higher-J lines are emitted from closer to the star and constrain the shape of the expansion profile. This is not affected very much by the adopted mass-loss rate, and the expansion velocity profile can be constrained by initially assuming reasonable parameters for the mass-loss rate. We base our initial values on previous models of CO lines for these sources assuming constant expansion velocities. Once the new expansion velocity profile is determined, we proceed to model the mass-loss rate and temperature profile of the CSE to include in the H2O modelling. The best-fit CO model is determined by calculating a grid of models varying the mass-loss rate and dust-to-gas ratio. To find the best fit to the data, a χ2 analysis comparing the observed line intensities with the model intensities is then performed by minimising (4)where Iobs,i is the integrated intensity of the observation i, Imod,i is the modelled intensity, σobs,i is the uncertainty in the observation. An error of 20% in the integrated intensities is assumed for all the CO lines.

3.2. The basic H2O model

We use an accelerated lambda iteration (ALI) method to calculate the radiative transfer of H2O lines. The model is the same as described in Maercker et al. (2008), but now includes a velocity profile in the same form as Eq. (1). We include the same dust-radiation field and velocity profiles as for CO. The mass-loss rate and kinetic temperature profiles are included from the results of the CO modelling. The remaining free parameters in the H2O models are the fractional H2O abundance at the inner radius and the e-folding radius of the H2O envelope (see Maercker et al. 2008). Together these describe the fractional abundance distribution of H2O. We model both ortho- and para-H2O. The best-fit H2O model is found by performing the same χ2 analysis as for CO, but now on a grid of models in H2O-abundance and radius of the H2O envelope. An error of 20% is assumed for the integrated intensities of the HIFI and PACS lines. A good estimate of the H2O radius is obtained from photodissociation models as the radius where the OH abundance peaks (as a photodissociation product of H2O; Netzer & Knapp 1987; Maercker et al. 2008, 2009). From here on, we will refer to the photodissociation radius of H2O as RNK, given by (5)

3.3. Heating and cooling in the energy balance: the problem of H2O line cooling

The CO model self-consistently solves the energy balance throughout the CSE by including all (known) heating and cooling terms (Schöier & Olofsson 2001). The effect of heating due to dust-gas collisions, photoelectric heating, and heat exchange between the dust and the gas is included. The total cooling includes cooling due to adiabatic expansion, and molecular line cooling from CO and H2. In order to determine the effect of H2O line cooling on the energy balance in the inner CSE, the contribution to the cooling by H2O lines is given as output in the H2O model. This can be included in the energy balance when calculating the CO model.

To fully include H2O line cooling, a careful iterative procedure is necessary, slowly increasing the contribution from H2O line cooling in each step. Immediately including full H2O line cooling typically leads to over-cooling in the inner envelope and results in negative temperatures. This iterative procedure was followed previously for the S-type stars χ Cyg and W Aql using HIFI data as constraints (Schöier et al. 2011; Danilovich et al. 2014, respectively).

In the case of M-type AGB stars, the abundance of H2O is high enough to make it one of the dominant sources of cooling, in particular in the inner envelope. Here we attempt to include H2O line cooling in an iterative procedure for the M-type stars R Dor and R Cas. We start from the best-fit CO model that does not include H2O line cooling. The kinetic temperature profile, velocity profile, and mass-loss rate are included as input in the H2O models. The abundance of H2O is then adjusted to fit the observed HIFI and PACS lines, and the H2O line cooling is calculated (see Sect. 5.1 for details). Only a fraction of the resulting H2O line cooling is then included in the CO model. This allows us to adjust the CO model in a way that counteracts extreme cooling and to ensure that the model still fits the observations by, e.g., increasing the mass-loss rate, increasing/decreasing the dust-to-gas ratio, and adjusting the drift-velocity profile. Once a new satisfactory CO model has been calculated, a new H2O model is calculated based on the new parameters, giving a new H2O abundance and hence cooling function. This procedure is repeated, increasing the fraction with which H2O line cooling is included with every step. Full H2O line cooling is included when 100% of the cooling function is included in the energy balance, good fits have been obtained to the observed lines for both CO and H2O, and the temperature profile has converged between iterations.

While this procedure was successful for the S-type stars, we come to the conclusion that it is not feasible for the high H2O abundances in M-type AGB stars. The iterative procedure makes the modelling very computationally expensive. It is therefore unfeasible to run a grid of models to properly probe parameter space, making the model results very unreliable. For the same reason the H2O e-folding radius is treated as a fixed parameter and is set to RNK. We nevertheless attempted to include 100% H2O line cooling for the low mass-loss-rate objects R Dor and R Cas and find that either a complete understanding of the heating and cooling processes in the CSE is lacking, or that the numerical implementation of radiative line cooling may have to be revised (see Sect. 5.1). Although a sufficient number of CO-lines nevertheless constrain the temperature profile of the CSE, a correct treatment of H2O line cooling is necessary for a full understanding of the physical processes in the CSEs of M-type AGB stars. Consequently the final CO and H2O models do not include H2O line cooling.

thumbnail Fig. 1

Best-fit CO models for R Dor. The blue histograms are the observations. The red lines are the model lines, the green lines are the model lines scaled to the same integrated intensities as the observations. The velocities are given with respect to the νLSR = 6.5 km s-1.

Open with DEXTER

Note that Decin et al. (2010b) do manage to include H2O line cooling in the energy balance for the M-type AGB star IK Tau. However, they do not include any explicit discussion on the issues of H2O line cooling, and it is not clear how they avoid the problems with the energy balance that we encounter in our modelling.

4. Results

4.1. CO models

Figure 1 shows the best-fit CO models compared to the observed ground-based and HIFI lines for R Dor. Figures A.1 to A.3 show the equivalent for R Cas, TX Cam, and IK Tau, respectively. Table A.1 gives the observed and modelled integrated intensities of all CO lines, including the PACS lines, together with the upper-level energies (Eup) and the ratios between models and observations for individual lines (Δ). Figure 2 shows the maps for the different CO model grids and Table 1 gives the best-fit parameters. The integrated intensities of the individual observed lines are all reproduced to within 25%. The average ratios between all observed and model lines (δ) are all close to 1. Although observed by HIFI, the CO(16−15) transition towards TX Cam is extremely over-estimated in the model by an order of magnitude, while the PACS line is well reproduced. The HIFI line is therefore excluded in the χ2 analysis. The ratio between modelled and observed line intensities as a function of Eup is shown in Fig. 3 (red crosses). Generally the lines scatter around a ratio of one, although there might be a slight correlation, with models preferentially under-predicting the intensities for low-Eup transitions.

For all models we use the photodissociation radius described by Mamon et al. (1988). This generally leads to overly resolved model line-profiles for the lower transitions, as evidenced by the double-peaked model line-profiles (specifically discussed for R Dor by Olofsson et al. 2002). This is a common problem in CO radiative-transfer models (e.g., Justtanont et al. 2005; Maercker et al. 2008; Danilovich et al. 2015). We attempted to fit the low-J line profiles by reducing the CO envelope size. However, this generally leads to an underestimation of the total integrated intensity compared to the high-J transitions. A possible way of correcting for this is by changing the dust parameters and the velocity field. However, the angular size of the envelope also depends on the assumed and uncertain distance to the source. Since this is a degenerate problem, we chose to fix the radius to the CO photodissociation radius and use the integrated intensities when determining the best-fit models.

In principle the terminal expansion velocity of the gas, υg,, should be well constrained by the width at zero power (in the following referred to as the width of the line) of the low-J CO transitions, while the widths of the higher-J CO transitions and H2O transitions constrain the shape of the velocity profile (β in Eq. (1)). We attempt to constrain the velocity law using the observed CO and H2O lines, adjusting β in the radiative-transfer models to fit the observed line widths. For R Cas we derive β = 2.5, and for the remaining three sources β = 1.5. For IK Tau the momentum equation gives a steeper profile with β = 0.6 (Decin et al. 2010c). Although comparatively shallow velocity profile seem to be preferred, the β is not very well constrained by the observed CO and H2O lines (Fig. 4), and we estimate the uncertainty in β to be ±1.

thumbnail Fig. 2

χ2-maps of the CO model grids for R Dor (top left), R Cas (top right), TX Cam (bottom left), and IK Tau (bottom right). The contours show the 1σ, 2σ, and 3σ levels. The cross marks the lowest model.

Open with DEXTER

The estimated mass-loss rates are consistent with previous estimates using the same radiative-transfer code and based on ground-based data alone (Maercker et al. 2008). Only for TX Cam do we derive a mass-loss rate a factor of 2 higher than the previous estimate still within the absolute uncertainty of independent mass-loss rate estimates of a factor of three (Ramstedt et al. 2008). The uncertainty in the mass-loss rate within the adopted model is 30%. For IK Tau it proved difficult to find an upper limit to the mass-loss rate, likely due to a saturation effect of the CO lines at high mass-loss rates, caused by increasing line optical depths and cooler envelopes (Ramstedt et al. 2008). See Sect. 5.5 for a comparison with similar codes.

Table 1

Stellar parameters and results of the SED and CO line-emission models.

The derived dust-to-gas ratios have considerable uncertainties. For R Dor and R Cas we derive values of a few times 10-4 with relatively strong constraints on the upper limit. For TX Cam we derive a value a factor 50 larger. For IK Tau the best-fit value is 10 times larger than for R Dor and R Cas, and agrees with chemical models at 6 R (Gobrecht et al. 2016), however no limits can be set. Ramstedt et al. (2008) derive values of a few times 10-3. The dust-to-gas ratio mainly affects the energy balance, where it is part of the heating of the gas through dust-gas collisions. This heating term depends degenerately also on the grain size and density. For this reason Schöier & Olofsson (2001) chose to combine these parameters in the so-called h-parameter. The derived d/g-ratios give h-parameters of 0.05, 0.03, 1.7, and 0.3 for R Dor, R Cas, TX Cam, and IK Tau, respectively. Ramstedt et al. (2008) derive values of 1.0 and 0.3 for TX Cam and IK Tau, respectively.

Finally, the kinetic temperature profiles are shown in Fig. 5. Note again that although the energy balance is solved in the CO models, the contribution from H2O line cooling is not included in any of the models (see Sect. 3.3). For TX Cam, R Cas, and R Dor the profiles can be described well by power laws of the form (6)where Ti is the temperature at the inner radius Ri. The exponents α are between −0.6 and −1.0. For the low mass-loss-rate M-type AGB star W Hya Khouri et al. (2014a) derive an exponent of −0.65 ± 0.05. In models of the chemical network in the CSE of AGB stars, Willacy & Millar (1997) use α = −0.6 for R Dor, TX Cam, and IK Tau.

Our profile for IK Tau deviates more from a simple power law in the inner and outermost parts of the CSE, but can generally be fit with α = −0.98. Solving the energy balance including H2O line cooling in the modelling of CO lines towards IK Tau, Decin et al. (2010b) derive a temperature profile consistent with an exponent of −0.6, significantly different from our value. Their derived mass-loss rate is 50% higher than our value, and they use a dust-to-gas ratio a factor four higher. These values lie within the uncertainties quoted here and in Decin et al. (2010b), and it is not clear why there is a difference in temperature profile and whether this is related to the H2O line cooling.

thumbnail Fig. 3

Ratio between observed and modelled integrated intensities vs. upper-level energy for the different transitions (CO: red crosses, ortho-H2O: blue dots, para-H2O: green circles).

Open with DEXTER

thumbnail Fig. 4

Velocity profiles derived for R Dor, R Cas, TX Cam and IK Tau (top to bottom) compared to the measured line-widths from lines emitted at different radii in the CSE (CO: red crosses, ortho-H2O: blue dots, para-H2O: green circles. The widths are measured at zero power. The radii are determined by the peak in the brightness profiles of each transition). The solid black lines give the profiles with the β given in Table 1, the dotted black lines give the profile for β = 0.6.

Open with DEXTER

thumbnail Fig. 5

Kinetic temperature profiles (solid lines) calculated in the energy balance of the CO modelling for R Dor (left panel, blue), R Cas (left panel, red), TX Cam (right panel, blue), and IK Tau (right panel, red). The dashed lines show power-law fits to the individual profiles (Eq. (6)). The numbers give the exponents for the respective fits.

Open with DEXTER

4.2. H2O models

Figure 6 shows the best-fit ortho-H2O and para-H2O models for the HIFI lines compared to the observed lines for R Dor. Figures A.4 to A.6 show the equivalent for R Cas, TX Cam, and IK Tau, respectively. Tables A.2 and A.3 list the same parameters as Table A.1, but for the HIFI and PACS lines of ortho-H2O and para-H2O, respectively. Figure 7 shows the -maps for ortho-H2O and para-H2O for all objects. We indicate both the best-fit radius in the grid, as well as the best-fit H2O abundance using the radius determined by Netzer & Knapp (1987), RNK. It is not possible to set upper limits to the H2O radii. The H2O line emitting region is limited by radiative excitation, and the H2O lines hence do not probe the outer regions very well. In all cases RNK is consistent with the best-fit radius when allowing R to vary freely, and we use this radius for the H2O envelope. All model results refer to models using RNK. Table 2 gives the best-fit values for the abundance, the average ratios between all integrated model and observed line intensities (δ), and the values. The dependence of the model and observed line ratios on the upper level energy Eup is shown in Fig. 3. As with the CO lines, the scatter lies around one, with a slight under-prediction by the model lines. However, no clear trend with Eup can be seen.

The derived ortho-H2O abundances lie between 10-510-4 relative to H2. For R Dor and R Cas this is significantly lower than previously derived values based only on ISO observations (Maercker et al. 2008). Including the spectrally resolved 557 GHz line observed with Odin required a reduction of the (constant) expansion velocity to fit the observed line-width, and led to a slight reduction in the ortho-H2O abundance for R Dor, albeit well within the uncertainties (Maercker et al. 2009). Including a velocity profile forces a significant reduction of the H2O abundance to reproduce the observed intensities. Observations of several para-H2O lines allow us to set constraints on the para-H2O abundances and derive ortho/para-H2O ratios (o/p-ratio). The uncertainties in the derived abundances lead to considerable uncertainties in the derived o/p-ratios. This is especially true for TX Cam, where the low number of observed para-H2O lines leads to a bad fit to the observations with the by-far highest value.

5. Discussion

5.1. The significance of H2O line cooling

The current implementation of H2O line cooling in the radiative-transfer codes implies that line cooling by H2O significantly affects the temperature in the inner envelope of M-type AGB stars (at radii Ri<r< a few Ri). The energy balance in this inner part of the CSE is very sensitive to the amount of H2O line cooling included in the model, quickly leading to extreme cooling. Models of the S-type stars χ Cyg and W Aql, that include H2O line cooling, had H2O abundances factors 2−35 lower than the abundances estimated for the M-type stars here (Schöier et al. 2011; Danilovich et al. 2014). For our high mass-loss-rate objects, H2O line cooling is too strong to be included in the energy balance. For R Dor and R Cas we nevertheless managed to include 100% H2O line cooling. While the iterative procedure does not allow us to realistically calculate a grid of models due to computational limitations, the resulting fit to the observed lines gives similar values to the models without H2O line cooling. The derived mass-loss rates and H2O abundances do not change significantly between models with and without cooling, since the observed CO lines effectively constrain the temperature profile. However, to balance the additional cooling in the inner envelope the d/g ratio needs to be increased by a factor of approximately two. Much more importantly, the only way we managed to include full H2O line cooling was to change the drift-velocity profile from starting low and increasing with radius to starting high and decreasing with radius. This changes the dust-velocity profile to a constant value throughout the entire CSE. Although the dust velocity profile and the drift velocity are not constrained observationally, theoretical models of dust-driven winds generally predict very low drift velocities between the dust and the gas for high mass-loss rates, and increasing drift-velocity with increasing distance (as assumed in this paper), or a constant drift velocity. A situation where the dust velocity is constant throughout the envelope therefore seems unphysical.

A constant dust velocity drastically increases the heating due to the dust-gas interaction in the inner CSE, and hence balances the cooling due to H2O. While this scenario is unrealistic, it indicates that some physics leading to heating in the inner envelope may be missing. Since the dust-velocity profile is not constrained by observations, changing the profile in an ad hoc manner becomes similar to adding an arbitrary heating term to counter-balance the cooling by H2O. We have no suggestion as to what this additional heating may be, or whether this indicates a fundamental problem in our understanding of the physics in the inner CSEs of AGB stars.

Another possible explanation for the problem of extreme H2O line cooling is the way the line cooling is calculated. In the current version of the radiative transfer, the contribution due to line cooling (for any molecule) is included following Sahai (1990): (7)where ΔEul is the energy difference of a transition from upper level u to lower level l, kB the Boltzmann constant, pl and pu the lower and upper level population densities, and clu and cul the collisional excitation and de-excitation rates from the lower to upper and upper to lower levels, respectively. If the total sum is positive, the rate of downward collisions is lower than the rate of upward collisions. The difference in energy is assumed to be radiated away, leading to net cooling. The contribution to the cooling is calculated in each radial bin. While this approach may be valid for lines that are optically thin throughout the CSE, the fact that no description of the optical depth is included may be problematic at the very high optical depths of the majority of the H2O-transitions. The net cooling assumes that the lost energy leaves the radial bin (and CSE). For very high optical depths, however, the radiated energy may be re-absorbed – either in a radial bin further out, or possibly still in the same radial bin. This would feed energy back into the envelope and reduce the net cooling. The current implementation may hence overestimate the amount of line cooling from H2O.

thumbnail Fig. 6

Best-fit o-H2O and p-H2O models of the HIFI lines for R Dor. The blue histograms are the observations. The red lines are the model lines, the green lines are the model lines scaled to the same integrated intensities as the observations. The velocities are given with respect to the νLSR = 6.5 km s-1.

Open with DEXTER

thumbnail Fig. 7

χ2-maps of the o-H2O (left) and p-H2O (right) model grids for R Dor, R Cas, TX Cam, and IK Tau (top to bottom). The cross marks the best-fit model using the H2O radius as a free parameter, the circle the best-fit model using RNK.

Open with DEXTER

This problem could possibly be solved by simply scaling the contribution to the cooling in Eq. (7) by the escape probability (1−eτ) /τ. The implementation of this in the code is however not trivial. The optical depth depends on the direction and is affected by the surrounding density and velocity profiles. Since the code only works along the radial direction, an angle-averaged optical depth would have to be calculated at each radial point, which is not straight-forward. We are currently working on this problem, trying to address the issue of properly including radiative line cooling into the code.

Note that this does not only affect the cooling due to H2O. In principle CO can also reach high optical depths, and in carbon AGB stars cooling due to HCN lines affect the temperature structure in the inner CSE in a very similar way to H2O in M-type AGB stars. Including line cooling correctly in the radiative-transfer codes, and identifying any potentially missed heating (or cooling mechanisms) is essential for a complete understanding of the physical conditions in the inner CSEs of AGB stars. This in turn affects chemical and dynamical models of the inner winds, and our understanding of the evolution of the star on the AGB in general.

5.2. The velocity profiles

Although not very well constrained, we derive comparatively shallow velocity profiles for all sources, with β = 1.5 ± 1.0. For IK Tau, Decin et al. (2010c) derive a β = 1−1.8. The velocity profiles derived here and by Decin et al. (2010c) are consistent with measurements of maser lines in the CSE of IK Tau (Bains et al. 2003). The momentum equation for IK Tau on the other hand gives β = 0.6 (Decin et al. 2010c). Khouri et al. (2014a) adopt a β = 1.5 to describe the velocity profile in W Hya, and find evidence that higher values may fit the high-J transitions of CO better. Note that using the upper-state energy Eup as a measure of the distance in the CSE and the line width to constrain the velocity profile as done by Justtanont et al. (2012) can be misleading. For R Dor the measured line width vs. Eup for all molecules detected in the HIFISTARS observations implies that R Dor has a decelerating wind (Justtanont et al. 2012), showing the importance of using full radiative transfer in order to derive correct velocity profiles.

All measurements indicate that the wind follows a more shallow profile than generally predicted by dynamical models solving the coupled momentum equations of the dust and gas (e.g., Decin et al. 2006; Ramstedt et al. 2008; Decin et al. 2010c). Note however, that Ramstedt et al. (2008) derive a velocity profile for IK Tau that is consistent with the results by Bains et al. (2003), showing that it is possible to derive velocity profiles from dynamical models that are consistent with the observed, shallow profiles. At the same time this shows the large uncertainty in the dynamical models, likely due to a sensitive dependence on the input parameters, and the assumption of a fully momentum coupled wind. The choice of υg,i does not significantly affect our results. The observed lines are emitted from too far out in the wind to effectively constrain the wind velocity at the inner radius. Any reasonable changes in υg,i would not change the derived β values, but would mainly affect the energy balance in the inner CSE. However, compared to the uncertainty in the treatment of the line cooling, this effect is small.

5.3. H2O abundances

For all sources we derive total H2O abundances between (0.3−4.0) × 10-4 relative to H2 (at the inner radius of our model CSE, i.e. Ri). Previous modelling of ISO observations of o-H2O resulted in the same o-H2O abundances for R Dor and IK Tau, while the values for R Cas and TX Cam were factors 5 and 10 higher, respectively. For IK Tau, Decin et al. (2010c) derive a total H2O abundance of 6.6 × 10-5, uncertain by a factor of 2. Within our uncertainties we just manage to meet the value derived by Decin et al. (2010c). However, they use a mass-loss rate that is higher by 50% than our value, likely resulting in a lower H2O abundance. For the M-type AGB star W Hya, previously derived H2O abundances are 8 × 10-4 and 3 × 10-4 (at radii smaller and larger than 4.5 × 1014 cm, respectively; Barlow et al. 1996), and 9 × 10-4 (Khouri et al. 2014b). For non-LTE chemical models, the predicted H2O abundance in the inner CSE (<5 R) of M-type AGB stars is 4 × 10-4 relative to H2, while LTE models predict H2O abundances closer to 7 × 10-5 (using TX Cam as an example; Cherchneff 2006). Non-LTE chemical models of IK Tau predict H2O abundances of 2 × 10-4 at 6 R (Gobrecht et al. 2016), consistent with our value derived here. Although the H2O abundances derived here and in other publications based on radiative-transfer modelling cover a large range, values of a few times 10-4 appear to be an upper limit. Taken at face value, and assuming that H2O abundances remain largely unaffected by processes in the intermediate wind, the results here favour the formation of H2O under non-LTE conditions in the inner wind. However, within the uncertainties, the derived abundances seem also to be consistent with the formation of H2O under equilibrium conditions at 1 R. The derived ortho- to para-H2O ratio is very uncertain, but is generally consistent with a value of 3, implying formation of H2O under warm conditions.

5.4. Discrepancies between models and observations

Although the 1D radiative-transfer models generally reproduce the overall line shapes (line-width and double-peaked, flat-topped, parabolic, or triangular lines), it is clear that deviations from these simple profiles are not reproduced. In addition to over-resolved low-J lines in CO, asymmetries between the red- and blue-shifted sides of the lines, and bumps and additional peaks (e.g. in R Dor and R Cas, CO(6−5) at 3 km s-1), are not reproduced. In the H2O lines the self-absorption on the blue side due to the high optical depths of the lines is generally well reproduced, however also with some deviations.

These differences may be caused by several effects. Varying mass-loss rates, velocity profiles, and photodissociation efficiencies will affect the spatial density distribution of the circumstellar wind, and the fractional abundance distribution of the molecules. This will change the optical depths of the lines and the excitation conditions. Such effects could, in principle, be included in the 1D radiative transfer. However, effects due to the three dimensional structure may have an equal, or even stronger, effect on the observed lines. These include asymmetric winds, disks and tori, and clumpiness of the circumstellar matter. Shocks in the inner wind and asymmetric illumination of the CSE due to convection cells on the stellar surface and atmospheric molecular absorption bands can also strongly affect the observed line shapes. Since the beginning of science observations with ALMA (Atacama Large Millimeter/submillimeter Array), an increasingly complex picture of the inner winds of AGB stars is being established, including the effects mentioned above. VLT-SPHERE observations in the optical of R Dor resolve the stellar disk, and show clear signs of an irregular brightness distribution from the stellar surface due to convection cells (Khouri et al., in prep.).

Table 2

Results of the H2O modelling.

5.5. Comparison to radiative transfer with GASTRoNOoM

Solving the molecular line radiative transfer for CSEs around AGB stars encounters several physical (e.g. energy balance, radiation fields, velocity profiles, molecular data) and numerical (e.g. modelling strategies, matrix inversion, optical depths) problems that have to be solved. As such, the number of advanced, non-LTE radiative-transfer codes in use is limited. This is unfortunate, since this makes benchmarking, and comparison and verification of model results difficult. One of the radiative-transfer codes that approaches the modelling of AGB CSEs in a similar way to our code is GASTRoNOoM (Decin et al. 2006, 2007, 2008).

The estimate of the mass-loss rate for IK Tau based on GASTRoNOoM lies within 50% of our results (Decin et al. 2010c). Khouri et al. (2014a) derive a mass-loss rate for W Hya of 1.5 × 10-7M yr-1, compared to 1.0 × 10-7M yr-1 in our previous modelling (Maercker et al. 2008). Based on a grid of models using GASTRoNOoM, De Beck et al. (2010) derive a formula to estimate the mass-loss rates from AGB stars based on the intensities of the CO lines. Using their formula results in mass-loss rates that are less than a factor two different from the values derived here. This is well within the limit of the absolute uncertainty expected from this modelling method (Ramstedt et al. 2008). Compared to what we find in our modelling, the derived H2O abundances using GASTRoNOoM appear to be lower by a factor of two for IK Tau. Comparing the results by Khouri et al. (2014b) with Maercker et al. (2008), they are a factor of five lower for W Hya. However, the absolute uncertainties are relatively large, and a more systematic comparison between the two radiative-transfer codes is necessary to identify whether the discrepancy in H2O abundances is due to numerical differences in the codes, or simply reflects the uncertainty and difficulty in modelling H2O line emission.

6. Conclusions

We have successfully modelled CO and H2O lines observed towards four M-type AGB stars, including spectrally resolved ground-based and HIFI observations for CO, spectrally resolved HIFI observations of o-H2O and p-H2O, and spectrally unresolved observations of CO and H2O with PACS. These are the first models of a larger sample of sources with spectrally resolved H2O line profiles. We derive H2O abundances, o/p-H2O ratios, the sizes of the H2O line emitting regions, velocity and temperature profiles, and attempt to include full H2O line cooling in the energy balance in the radiative transfer. Based on the modelling we arrive at following conclusions.

  • Our modelling indicates that the winds around these four M-type AGB stars have a comparatively shallow acceleration profile, in contrast to what is generally expected from dynamical models of dust-driven winds. However, we also note that the widths of the line profiles only moderately constrain the velocity profile.

  • Although we managed to include H2O line cooling in the energy balance for the low mass-loss rate objects R Dor and R Cas, this required an unphysical description of the dust expansion velocity in order to counter-balance the extreme cooling from H2O in the inner CSE. This implies that the physical processes that describe the heating and cooling in the inner envelope are not completely understood and/or that the formal implementation of line cooling is not correct.

  • We derive generally high H2O abundances, up to a few times 10-4 relative to H2. Within the uncertainties, all derived abundances are consistent with both NLTE and LTE abundances derived from chemical models, emphasizing the importance of improved radiative transfer models to reduce uncertainties on abundances derived from observations.

  • The derived o/p-ratios are uncertain, however consistent with a value of three, indicating formation of H2O in the warm, inner CSE.

Computational limitations currently do not allow us to self-consistently solve the physics and dynamics throughout the CSE in three dimensions. 1D radiative-transfer models as presented here are important for constraining the basic physical parameters of the envelopes (e.g., temperature, velocity, average mass-loss rates). These can be used as input for 3D models (e.g. LIME; Brinch & Hogerheijde 2010) to study the effect of structure in three dimensions. In order to arrive at a self-consistent model that fully describes the physical conditions throughout the CSE, however, it is important to solve the problem of radiative line cooling, in particular for H2O.

Acknowledgments

M.M. has received funding from the People Programme (Marie Curie Actions) of the EU’s FP7 (FP7/2007-2013) under REA grant agreement No. 623898.11. H.O. acknowledges financial support from the Swedish Research Council. M.M. would like to thank Fredrik Schöier for his contribution at the beginning of this project, his advice and support. Sadly he could not be here to see the results.

References

Appendix A: Observed and model line intensities, and model line profiles

Table A.1

Observed and modelled CO lines.

Table A.2

Observed and modelled ortho-H2O lines.

Table A.3

Observed and modelled para-H2O lines.

thumbnail Fig. A.1

Best-fit CO models for R Cas. The blue histograms are the observations. The red lines are the model lines, the green lines are the model lines scaled to the same integrated intensities as the observations. The velocities are given with respect to the νLSR = 24.5 km s-1.

Open with DEXTER

thumbnail Fig. A.2

Best-fit CO models for TX Cam.The blue histograms are the observations. The red lines are the model lines, the green lines are the model lines scaled to the same integrated intensities as the observations. The velocities are given with respect to the νLSR = 12.0 km s-1.

Open with DEXTER

thumbnail Fig. A.3

Best-fit CO models for IK Tau.The blue histograms are the observations. The red lines are the model lines, the green lines are the model lines scaled to the same integrated intensities as the observations. The velocities are given with respect to the νLSR = 34.0 km s-1.

Open with DEXTER

thumbnail Fig. A.4

Best-fit o-H2O and p-H2O models of the HIFI Lines for R Cas. The blue histograms are the observations. The red lines are the model lines, the green lines are the model lines scaled to the same integrated intensities as the observations. The velocities are given with respect to the νLSR = 24.5 km s-1.

Open with DEXTER

thumbnail Fig. A.5

Best-fit o-H2O and p-H2O models of the HIFI Lines for TX Cam. The blue histograms are the observations. The red lines are the model lines, the green lines are the model lines scaled to the same integrated intensities as the observations. The velocities are given with respect to the νLSR = 12.0 km s-1.

Open with DEXTER

thumbnail Fig. A.6

Best-fit o-H2O and p-H2O models of the HIFI Lines for IK Tau. The blue histograms are the observations. The red lines are the model lines, the green lines are the model lines scaled to the same integrated intensities as the observations. The velocities are given with respect to the νLSR = 34.0 km s-1.

Open with DEXTER

All Tables

Table 1

Stellar parameters and results of the SED and CO line-emission models.

Table 2

Results of the H2O modelling.

Table A.1

Observed and modelled CO lines.

Table A.2

Observed and modelled ortho-H2O lines.

Table A.3

Observed and modelled para-H2O lines.

All Figures

thumbnail Fig. 1

Best-fit CO models for R Dor. The blue histograms are the observations. The red lines are the model lines, the green lines are the model lines scaled to the same integrated intensities as the observations. The velocities are given with respect to the νLSR = 6.5 km s-1.

Open with DEXTER
In the text
thumbnail Fig. 2

χ2-maps of the CO model grids for R Dor (top left), R Cas (top right), TX Cam (bottom left), and IK Tau (bottom right). The contours show the 1σ, 2σ, and 3σ levels. The cross marks the lowest model.

Open with DEXTER
In the text
thumbnail Fig. 3

Ratio between observed and modelled integrated intensities vs. upper-level energy for the different transitions (CO: red crosses, ortho-H2O: blue dots, para-H2O: green circles).

Open with DEXTER
In the text
thumbnail Fig. 4

Velocity profiles derived for R Dor, R Cas, TX Cam and IK Tau (top to bottom) compared to the measured line-widths from lines emitted at different radii in the CSE (CO: red crosses, ortho-H2O: blue dots, para-H2O: green circles. The widths are measured at zero power. The radii are determined by the peak in the brightness profiles of each transition). The solid black lines give the profiles with the β given in Table 1, the dotted black lines give the profile for β = 0.6.

Open with DEXTER
In the text
thumbnail Fig. 5

Kinetic temperature profiles (solid lines) calculated in the energy balance of the CO modelling for R Dor (left panel, blue), R Cas (left panel, red), TX Cam (right panel, blue), and IK Tau (right panel, red). The dashed lines show power-law fits to the individual profiles (Eq. (6)). The numbers give the exponents for the respective fits.

Open with DEXTER
In the text
thumbnail Fig. 6

Best-fit o-H2O and p-H2O models of the HIFI lines for R Dor. The blue histograms are the observations. The red lines are the model lines, the green lines are the model lines scaled to the same integrated intensities as the observations. The velocities are given with respect to the νLSR = 6.5 km s-1.

Open with DEXTER
In the text
thumbnail Fig. 7

χ2-maps of the o-H2O (left) and p-H2O (right) model grids for R Dor, R Cas, TX Cam, and IK Tau (top to bottom). The cross marks the best-fit model using the H2O radius as a free parameter, the circle the best-fit model using RNK.

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

Best-fit CO models for R Cas. The blue histograms are the observations. The red lines are the model lines, the green lines are the model lines scaled to the same integrated intensities as the observations. The velocities are given with respect to the νLSR = 24.5 km s-1.

Open with DEXTER
In the text
thumbnail Fig. A.2

Best-fit CO models for TX Cam.The blue histograms are the observations. The red lines are the model lines, the green lines are the model lines scaled to the same integrated intensities as the observations. The velocities are given with respect to the νLSR = 12.0 km s-1.

Open with DEXTER
In the text
thumbnail Fig. A.3

Best-fit CO models for IK Tau.The blue histograms are the observations. The red lines are the model lines, the green lines are the model lines scaled to the same integrated intensities as the observations. The velocities are given with respect to the νLSR = 34.0 km s-1.

Open with DEXTER
In the text
thumbnail Fig. A.4

Best-fit o-H2O and p-H2O models of the HIFI Lines for R Cas. The blue histograms are the observations. The red lines are the model lines, the green lines are the model lines scaled to the same integrated intensities as the observations. The velocities are given with respect to the νLSR = 24.5 km s-1.

Open with DEXTER
In the text
thumbnail Fig. A.5

Best-fit o-H2O and p-H2O models of the HIFI Lines for TX Cam. The blue histograms are the observations. The red lines are the model lines, the green lines are the model lines scaled to the same integrated intensities as the observations. The velocities are given with respect to the νLSR = 12.0 km s-1.

Open with DEXTER
In the text
thumbnail Fig. A.6

Best-fit o-H2O and p-H2O models of the HIFI Lines for IK Tau. The blue histograms are the observations. The red lines are the model lines, the green lines are the model lines scaled to the same integrated intensities as the observations. The velocities are given with respect to the νLSR = 34.0 km s-1.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.