A&A 516, A69 (2010)

Circumstellar molecular composition of the oxygen-rich AGB star IK Tauri

II. In-depth non-LTE chemical abundance analysis

L. Decin^1,2, - E. De Beck¹ - S. Brünken^3,4 - H. S. P. Müller^4,5 - K. M. Menten⁵ - H. Kim^5,6 - K. Willacy⁷ - A. de Koter^2,8 - F. Wyrowski⁵

1 - Department of Physics and Astronomy, Institute of Astronomy, K.U.Leuven, Celestijnenlaan 200D, 3001 Leuven, Belgium
2 - Sterrenkundig Instituut Anton Pannekoek, University of Amsterdam, PO Box 9429, 1090 CE Amsterdam, The Netherlands
3 - Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA
4 - I. Physikalisches Institut, Universität zu Köln, Zülpicher Street 77, 50937 Köln, Germany
5 - Max-Planck Institut für Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany
6 - MPI für Gravitationsphysik, Calinstr. 38, 30167 Hannover, Germany
7 - Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA 91109, USA
8 - Astronomical Institute, Utrecht University, Princetonplein 5, 3584 CC Utrecht, The Netherlands

Received 26 January 2010 / Accepted 30 March 2010

Abstract
Context. The interstellar medium is enriched primarily by matter ejected from evolved low and intermediate mass stars. The outflow from these stars creates a circumstellar envelope in which a rich gas-phase chemistry takes place. Complex shock-induced non-equilibrium chemistry takes place in the inner wind envelope, dust-gas reactions and ion-molecule reactions alter the abundances in the intermediate wind zone, and the penetration of cosmic rays and ultraviolet photons dissociates the molecules in the outer wind region.
Aims. Little observational information exists on the circumstellar molecular abundance stratifications of many molecules. Furthermore, our knowledge of oxygen-rich envelopes is not as profound as for the carbon-rich counterparts. The aim of this paper is therefore to study the circumstellar chemical abundance pattern of 11 molecules and isotopologs (¹²CO, ¹³CO, SiS, ²⁸SiO, ²⁹SiO, ³⁰SiO, HCN, CN, CS, SO, SO₂) in the oxygen-rich evolved star IK Tau.
Methods. We have performed an in-depth analysis of a large number of molecular emission lines excited in the circumstellar envelope around IK Tau. The analysis is done based on a non-local thermodynamic equilibrium (non-LTE) radiative transfer analysis, which calculates the temperature and velocity structure in a self-consistent way. The chemical abundance pattern is coupled to theoretical outer wind model predictions including photodestruction and cosmic ray ionization. Not only the integrated line intensities, but also the line shapes are used as diagnostic tool to study the envelope structure.
Results. The deduced wind acceleration is much slower than predicted from classical theories. SiO and SiS are depleted in the envelope, possibly due to the adsorption onto dust grains. For HCN and CS a clear difference with respect to inner wind non-equilibrium predictions is found, either indicating uncertainties in the inner wind theoretical modeling or the possibility that HCN and CS (or the radical CN) participate in the dust formation. The low signal-to-noise profiles of SO and CN prohibit an accurate abundance determination; the modeling of high-excitation SO₂ lines is cumbersome, possibly related to line misidentifications or problems with the collisional rates. The SiO isotopic ratios (²⁹SiO/²⁸SiO and ³⁰SiO/²⁸SiO) point toward an enhancement in ²⁸SiO compared to results of classical stellar evolution codes. Predictions for H₂O emission lines in the spectral range of the Herschel/HIFI mission are performed.

Key words: astrochemistry - molecular processes - radiative transfer - submillimeter: stars - stars: AGB and post-AGB - stars: mass-loss

1 Introduction

Asymptotic giant branch (AGB) stars are well known to release significant amounts of gas and dust in the interstellar medium via (copious) mass loss. This mass loss dominates the evolution of the star and ultimately, when the stellar envelope is exhausted, causes the star to evolve off the AGB into the post-AGB phase. The outflow from these evolved stars creates an envelope, which fosters gas-phase chemistry. The chemical complexity in circumstellar envelopes (CSEs) is thought to be dominated by the elemental carbon to oxygen ratio: oxygen-rich M-stars have a C/O ratio less than unity, carbon-rich C-stars have C/O > 1, and for S-stars C/O is 1.

Many papers have focused on the CSEs of carbon-rich stars in which a rich chemistry takes place. This is reflected by the detection of over 60 different chemical compounds, including unusual carbon chain radicals in the CSE of IRC +10216, the prototype of carbon stars (e.g. Cernicharo et al. 2000). In contrast, only 10-12 compounds have been identified in the chemically most interesting oxygen-rich evolved stars, like IK Tau and VY CMa (e.g. Ziurys et al. 2007). The first observations of carbon-bearing molecules (other than CO) in oxygen-rich AGBs were somewhat unexpected (e.g. Jewell et al. 1986; Deguchi & Goldsmith 1985). Nowadays, the formation of carbon molecules is thought to be the result of shock-induced non-equilibrium chemistry in the inner circumstellar envelope (e.g. Duari et al. 1999) and/or a complex chemistry in the outer envelope triggered by the penetration of cosmic rays and ultra-violet radiation (e.g. Willacy & Millar 1997). Recently, a new interstellar molecule, PO (), has been detected toward the envelope of the oxygen-rich supergiant VY CMa (Tenenbaum et al. 2007). Phosphorus monoxide is the first interstellar molecule detected that contains a P-O bond, a moiety essential in biochemical compounds. It is also the first new species identified in an oxygen-rich, as opposed to a carbon-rich, circumstellar envelope. These results suggest that oxygen-rich shells may be as chemically diverse as their carbon counterparts.

Circumstellar molecules have been extensively observed both in the form of surveys of a single molecular species and in the form of searches for various molecular species in a limited number of carefully selected sources. The aim of these studies was to derive (i) the mass-loss rate (from CO rotational lines) or (ii) molecular abundances. For this latter purpose, several methods exist, each with varying degrees of complexity. (1.) For example, Bujarrabal et al. (1994) and Olofsson et al. (1998) showed that simple molecular line intensity ratios, if properly chosen, may be used to study the chemical behavior in CSEs. The use of line intensity ratios has the advantage of requiring no assumptions about a circumstellar model, but it also limits the type of conclusions that can be drawn. (2.) Several authors have derived new constraints on chemical and circumstellar models based on the simplifying assumption of unresolved optically thin emission thermalized at one excitation temperature (e.g. Bujarrabal et al. 1994; Lindqvist et al. 1988; Kim et al. 2010; Omont et al. 1993). (3.) Later on, observations were (re)-analyzed based on a non-LTE (non-local thermodynamic equilibrium) radiative transfer model (e.g. Teyssier et al. 2006; Bieging et al. 2000; Schöier et al. 2007a). In this study we will go one step further and abandon or improve few of the assumptions still made in many non-LTE analyses.

1.

Quite often, the temperature structure - being the most important factor determining the molecular line excitation - is approximated with a powerlaw (e.g. Teyssier et al. 2006; Bieging et al. 2000). Effects of different heating and cooling mechanisms are hence not properly taken into account. For instance, in the outermost parts of the envelope the temperature profile deviates from a power law distribution once the influence of photoelectric heating by the external interstellar radiation field becomes important (e.g. Crosas & Menten 1997; Justtanont et al. 1994; Decin et al. 2006).

2.

The shell is often assumed to expand at a constant velocity (e.g. Bieging et al. 2000; Schöier et al. 2007a). However, for molecular lines primarily formed in the wind acceleration zone, the effect of a non-constant velocity structure on the derived molecular abundance may be significant.

3.

The fractional abundances are estimated to follow an exponential or Gaussian distribution, assuming that the molecules are formed in the inner envelope and are photodissociated or absorbed onto dust grains further out (e.g. Bieging et al. 2000; González Delgado et al. 2003; Schöier et al. 2007a). The effect of extra formation and/or depletion processes in the envelope can hence not be taken into account.

4.

Often, a maximum of two molecules (CO and one other) is analyzed at once (e.g. González Delgado et al. 2003; Schöier et al. 2007a).

5.

Integrated line intensities are often used as a criterion to analyze the circumstellar chemical structure. However, line shapes provide us with strong diagnostic constraints to pinpoint the wind acceleration, which in turn has an influence on the deduced fractional abundances.

In this paper, we will study the circumstellar chemical abundance fractions of eleven different molecules and isotopologs in the oxygen-rich AGB star IK Tau based on the non-LTE radiative transfer code GASTRoNOoM (Decin et al. 2006,2007), which computes the temperature and velocity structure in the envelope in a self-consistent way. Chemical abundance stratifications are coupled to theoretical non-equilibrium (non-TE) predictions in the outer envelope by Willacy & Millar (1997) and compared to the shock-induced non-TE inner wind predictions by Duari et al. (1999) and Cherchneff (2006). The star IK Tau has been chosen for study because of the wealth of observations which are available for this target and because its envelope is thought to be (roughly) spherically symmetric (Lane et al. 1987; Hale et al. 1997; Kim et al. 2010; Marvel 2005).

The star IK Tau, also known as NML Tau, was discovered in 1965 by Neugebauer et al. (1965). It is an extremely red Mira-type variable with spectral type ranging from M8.1 to M11.2 and a period around 470 days (Wing & Lockwood 1973). From dust shell motions detected at 11 m with the ISI interferometer, Hale et al. (1997) deduced a distance of 265 pc. This agrees well with the results of Olofsson et al. (1998), who computed a distance of 250 pc from integrated visual, near-infrared and IRAS data using a period-luminosity relation. Estimated mass-loss rates range from 3.8  10^-6 (Neri et al. 1998) to 3  10^-5 /yr (González Delgado et al. 2003). IK Tau's proximity and relatively high mass-loss rate (for a Mira) facilitates the observation of molecular emission lines.

In Sect. 2 we present the molecular line observational data used in this paper. Section 3 describes the background of the excitation analysis: the radiative transfer model used, the molecular line data and the theoretical ideas on molecular abundance stratification in the envelope. Section 4 describes the results: we first focus on the velocity structure in the envelope with special attention to the acceleration zone, after which the derived stellar parameters are discussed. Thereafter, the abundance structure for each molecule is derived and compared to the theoretical inner and outer wind predictions and observational results found in the literature. The time variability and SiO isotopic ratios are discussed in Sect. 5 and water line predictions are performed in Sect. 5.3. We end with some conclusions in Sect. 6.

Table 1:   Overview of the molecular line transitions used in this research, with indication of the frequency, the upper energy level, the telescope, the main beam half power beam width (HPBW) and main beam efficienciy ( ).

2 Observational data

Part of the observations were obtained from our own observing programs scheduled at the JCMT, APEX and IRAM. These observations and the data reduction are described in Sect. 2.1. Other data are extracted from the literature and summarized briefly in Sect. 2.2. An overview is given in Table 1.

2.1 Observations and data reduction

The ¹²CO(2-1), ¹²CO(3-2), ¹²CO(4-3) and the ¹³CO(2-1) observations were extracted from the JCMT archive. Additional data with the APEX 12 m telescope were obtained for the ¹²CO(3-2), ¹²CO(4-3), ¹²CO(7-6), and ¹³CO(3-2) molecular line transitions. With the 30 m telescope of the Institut de Radio Astronomie Millimetrique (IRAM) molecular line observations were performed in December 2006. During this observing run, data on the CO(2-1), SiS(8-7), SiS(12-11), SO₂( 14_3,11-14_2,12), SO₂( 4_3,1-4_2,2), SO₂( 3_3,1-3_2,2), SO₂( 5_3,3-5_2,4), and HCN(3-2) line transitions were obtained.

The JCMT and APEX observations were carried out in a position-switching mode. The IRAM observations were done in the wobbler-switching mode with a throw of 60''. The frequency resolution for the JCMT-data equals 0.0305 MHz, for the APEX data it is 0.1221 MHz. The resolution was 1.25 MHz for the 3 mm and 2 mm IRAM observations and 1 MH or/and 4 MHz at 1.3 mm and 1 mm, resulting in a resolution slightly higher than 1 km s^-1 for these observations.

The JCMT data reduction was performed with the SPLAT devoted routines of STARLINK, the APEX and IRAM data were reduced with CLASS. A polynomial of first order was fitted to an emission-free region of the spectral baseline and subtracted. To increase the signal-to-noise ratio, the data were rebinned to a resolution of 1 km s^-1 so that we have at least 40 independent resolution elements per line profile. The antenna temperature,  , was converted to the main-beam temperature ( ), using a main-beam efficiency,  as specified in Table 1. The absolute uncertainties are 20%.

2.2 Literature data

To have better constraints on the chemical abundance pattern in the wind region around IK Tau, additional data were taken from the literature (see Table 1).

High-quality observations were performed by Kim (2006) with the APEX telescope in Chile during observing periods in November 2005, April 2006 and August 2006 (see also Kim et al. 2010, hereafter referred to as Paper I). In total, 34 transitions from 12 molecular species, including a few maser lines, were detected toward IK Tau.

Schöier et al. (2007a) published the observations of four SiS lines: the (5-4) and (6-5) rotational line transitions were obtained with the Onsala Space Observatory (OSO) telescope, the (12-11) and (19-18) rotational line observations were performed with the JCMT telescope.

The SiO thermal radio line emission from a large sample of M-type AGB stars, including IK Tau, was studied by González Delgado et al. (2003). The SiO(2-1) line transition was obtained with the OSO telescope, the SiO(5-4) and (6-5) transitions with the Swedish ESO Submillimeterwave Telescope (SEST).

The ¹²CO line data were obtained by Teyssier et al. (2006) with the IRAM and JCMT telescope. The CO(1-0) line was observed by Ramstedt et al. (2008) in December 2003 with the 20 m OSO telescope. Olofsson et al. (1998) reported on the detection of the CS(2-1) line at 98 GHz with the OSO telescope with an integrated intensity of 0.5 K km s^-1.

3 Excitation analysis

3.1 Radiative transfer model

The observed molecular line transitions provide information on the thermodynamic and chemical structure in the envelope around IK Tau. The line profiles were modeled with our non-LTE radiative transfer code GASTRoNOoM (Decin et al. 2006). The code (1) calculates the kinetic temperature and velocity structure in the shell by solving the equations of motion of gas and dust and the energy balance simultaneously; then (2) solves the radiative transfer equation in the co-moving frame with the Approximate Newton-Raphson operator as developed by Schönberg & Hempe (1986) and computes the non-LTE level-populations; and finally (3) determines the observable line profile by ray-tracing. For a full description of the code we refer to Decin et al. (2006).

The main assumption of the code is a spherically symmetric wind. The mass-loss rate is allowed to vary with radial distance from the star. The local line width is assumed to be described by a Gaussian and is made up of a microturbulent component with a Doppler width of 1.5 km s^-1 and a thermal component, which is calculated from the derived kinetic temperature structure.

Two major updates have been made since the original publication in Decin et al. (2006).

• The code now iterates on steps (1) and (2) to obtain the kinetic temperature structure in a self-consistent manner from solving the energy balance equation, where the CO and H₂O line cooling (or heating) are directly obtained from the excitation analysis, i.e.

  (1)

  
  
  where n_l and n_u are the level populations in the lower and upper levels participating in the transition at rest frequency , and and are the CO-H₂ collisional rate coefficients. The cooling rate  (in erg s^-1 cm^-3) is defined as positive for net cooling. For IK Tau the water line cooling dominates the CO line cooling by more than one order of magnitude in the inner wind region; for regions beyond 10¹⁶ cm the CO line cooling dominates over H₂O cooling with the adiabatic cooling the dominant coolant agent.

• While the original version of the code approximates the stellar atmosphere with a blackbody at the stellar effective temperature, an additional option is now implemented to use the MARCS theoretical model atmospheres and theoretical spectra (Decin & Eriksson 2007; Gustafsson et al. 2008) to estimate the stellar flux. Molecular species in the CSE, less abundant than CO and with larger dipole moments are primarily excited by infrared radiation from the central star (with the possible exception of HCN, Jura 1983, and H₂O, see Sect. 5.3). For CO, the infrared radiation competes with rotational excitation by collisions and by trapped rotational line photons to determine the populations of the rotational levels (Knapp & Morris 1985). For those minor species, the blackbody approximation of the stellar flux may lead to inaccurate absolute intensity predictions on the order of 5 to 20%.

Decin et al. (2006,2007) and Crosas & Menten (1997) demonstrated that not only the integrated line intensities, but particularly the line shapes should be taken into account for a proper determination of the envelope structure. This is particularly true because (i) the line shapes contain valuable information on the line forming region, e.g. a Gaussian line profile points toward line formation partially in the inner wind where the stellar wind has not yet reached its full terminal velocity (Bujarrabal & Alcolea 1991), and (ii) while the absolute uncertainty of the intensity ranges between 20 to 50%, the relative accuracy (of the line shapes) is on the order of a few percent. In contrast to most other studies, we therefore will not use a classical -analysis using the integrated line intensity, but will perform the model selection and the assessment of goodness-of-fit using the log-likelihood function on the full line profiles as described in Decin et al. (2007).

3.2 Molecular line data

In this paper, line transitions of CO, SiO, SiS, CS, CN, HCN, SO, and SO₂ will be modeled, and H₂O line profile predictions for the Herschel/HIFI instrument will be performed. The molecular line data used in this paper are described in Appendix A.

3.3 Molecular abundance stratification

Theoretical chemical calculations clearly show that the fractional abundances (relative to H₂) vary throughout the envelope. Chemical processes responsible for the molecular content are dependent on the position in the envelope (see Fig. 1). In the stellar photosphere and at the inner boundary of the envelope, the high gas density and temperature ensure thermal equilibrium (TE). The TE is suppressed very close to the photosphere because of the action of pulsation-driven shocks propagating outwards. Furthermore the regions of strong shock activity correspond to the locus of grain formation and wind acceleration. This region is referred to as the inner envelope (or inner wind) which extends over a few stellar radii. At larger radii (5 to 100 ) the newly formed dust grains interact with the cooler gas. Depletion or formation of certain molecular/atomic species may result from this interaction and these layers are referred to as the intermediate envelope. This is also the region where parent molecules, injected in the envelope, may begin to break down, and daughter molecules are formed. At still larger radii (>100 ), the so-called outer envelope is penetrated by ultraviolet interstellar photons and cosmic rays resulting in a chemistry governed by photochemical processes.

Because our modeling results will be compared to chemical abundance predictions in the outer envelope by Willacy & Millar (1997) and in the inner envelope by Duari et al. (1999) and Cherchneff (2006), we first briefly describe these studies in Sects. 3.3.1 and 3.3.2 respectively. In Sect. 3.4 we discuss how we have implemented this knowledge in the modeling of the molecular line transitions.

3.3.1 Chemical stratification in the outer envelope

The chemistry in the outer envelope of IK Tau has been modeled by Willacy & Millar (1997). This chemical kinetic model aims at deriving the abundance stratification in the outer envelope (between 2  10¹⁵ and 2  10¹⁸ cm). The chemistry is driven by a combination of cosmic-ray ionization and ultraviolet radiation and starts from nine parent species injected into the envelope (see Table 2). The CSE of IK Tau was assumed to be spherically symmetric with a constant mass-loss rate and a constant expansion velocity of 19 km s^-1. The temperature was described by a power law

with r₀ = 2 10¹⁵ cm and T(r₀) taken to be 100 K for  > 5  10^-6 /yr and 300 K otherwise. For IK Tau T(r₀) was assumed to be 100 K. The derived fractional abundances for the molecules studied in this paper are represented in Fig. 2.

Table 2:   The fractional abundance (relative to H₂) taken for the parent species by Willacy & Millar (1997).

The model of Willacy & Millar (1997) succeeded in reproducing the observed values of certain species, but failed for some other molecular abundances: the calculated abundance of HCN was too low and the injected abundance of the parent species SiS was about 10 times higher than observed. Duari et al. (1999) noted that the input molecular abundances of some parent species are sometimes questionably high because there exists no observational or theoretical evidence for the formation of these species in the inner and intermediate envelopes of O-rich Miras (see Sect. 3.3.2). More importantly, Duari et al. (1999) showed that HCN should form in the inner envelope or extended stellar atmosphere due to non-equilibrium shock chemistry and may be a parent species injected to the outer envelope. Recent observational studies also indicate that HCN must be formed in the inner envelope (Bieging et al. 2000; Marvel 2005). These results are in contrast to the modeling efforts of Willacy & Millar (1997), where HCN was not yet considered as a parent species.

3.3.2 Chemical stratification in the inner envelope

Carbon-bearing molecules have been identified in the envelopes of many oxygen-rich AGB stars (e.g., Bujarrabal et al. 1994) and it was first thought that the observed carbon species were produced in the outer wind of O-rich stars via photochemical processes. However, Duari et al. (1999) showed that shock-induced non-equilibrium chemistry models predict the formation of large amounts of a few carbon species, like HCN, CS and CO₂, in the inner envelope of IK Tau: these molecules are hence formed in the post-shocked layers and are then ejected in the outer wind as ``parent'' species. For some parent species, the non-equilibrium chemistry does not significantly alter the initial photospheric TE abundances. But other species, abundant in the TE photosphere, are quickly destroyed in the outflow by the non-equiblibrium chemistry generated by shocks (e.g., OH, SiS and HS). Again other species (like SO) appear to be absent in the inner regions of the wind, and are thought to be produced by ion-molecule reactions in the photo-dissociation regions of the outer wind.

Cherchneff (2006) continued the study of shock-induced non-equilibrium chemistry in the inner wind of AGB stars. She demonstrated that whatever the enrichment in carbon of the star (i.e. the C/O ratio), the atomic and molecular content after the passage of the first shock in the gas layers just above the stellar photosphere is very much the same, and in many cases totally different from what would be expected from thermodynamic equilibrium (TE) calculations. For the oxygen-rich envelope around TX Cam - which is almost a stellar twin of IK Tau, but with slightly lower luminosity -- Cherchneff (2006) found that while e.g. the TE abundance of HCN (CS) is predicted to be 1.9  10^-11 (2.5  10^-11), the non-TE fractional abundances at 2.5  are predicted to be 9  10^-5 (1.85  10^-5). The fractional abundances derived by Cherchneff (2006) differ from the abundances of the injected parent molecules in the study of Willacy & Millar (1997) (see Table 2 and Fig. 2): Willacy & Millar (1997) did not consider CS and HCN as parent molecules, and the (injected) abundance of SiS (3.5  10^-6) is much higher than the abundance stratification derived by Cherchneff (2006) (see their Fig. 8).

3.4 Modeling strategy

3.4.1 Envelope structure as traced by the CO lines

The physical properties of the circumstellar gas, like the temperature, velocity and density structure, are determined from the radiative transfer modeling of the multi-transitional (sub)millimetre CO line observations. Because higher-J lines are formed at higher temperature, different transitions offer the possibility to trace different regions in the envelope. The highest CO energy level traced is the CO   level at 154.8 K. The available rotational CO lines will hence be good tracers for the region beyond 100 , but they do not put strong constraints on the temperature in the inner CSE. The upcoming Herschel/HIFI mission will be crucial in the study of the temperature structure in this inner wind region.

An extensive grid has been calculated with parameters ranging from 2000 to 3000 K for the stellar temperature  , from 1  10¹³ to 6  10¹³ cm for the stellar radius , an inner (dust condensation) radius between 2 and 30 , distance between 200 and 300 pc, and a constant mass-loss rate between 1  10^-6 and 5  10^-5 /yr. As briefly explained in Sect. 3.1, a log-likelihood method (Decin et al. 2007) is used to find the best-fit model and derive a 95% confidence interval for the model parameters. The results will be presented in Sect. 4.

3.4.2 Abundance stratification through the envelope

From the descriptions of theoretical abundance estimates in the inner and outer envelope (Sect. 3.3.1, 3.3.2) it is clear that there is still some debate about the abundance structure in the envelope. SiS and HCN were already given as an example, but other molecules as e.g. CS and SO also pose a problem. This is illustrated in Fig. 2, where one notices for a few molecules a significant difference between the theoretically predicted fractional abundance in the inner envelope by Duari et al. (1999) and Cherchneff (2006) and the abundance of the parent molecules injected in the outer envelope by Willacy & Millar (1997). One of the big questions still existing concerns the modifications of the molecular abundances in the intermediate wind region due to gas-grain reactions. Currently, no theoretical efforts have been made to model this region in terms of molecular ``leftovers'' after the dust formation has occurred. For O-rich envelopes, it is thought that CO, HCN and CS are quite stable and travel the entire envelope unaltered until they reach the photo-dissocation region of the outer wind, because these molecules do not participate in the formation of dust grains like silicates and corundum (Duari et al. 1999). In contrast, SiO is a candidate molecule for depletion in the intermediate wind region due to the formation of SiO₂ (via a reaction with OH) whose condensation product, silica, is tentatively identified in post-AGB stars (Molster et al. 2002) and is claimed to be the carrier of the 13 m feature in low mass-loss rate AGB stars (Speck et al. 2000, but other studies argue that this feature is due to spinel). The theoretical modeling of Duari et al. (1999) and Cherchneff (2006) predict an SiS abundance 2 to 3 orders of magnitude lower than the observed value, indicating that SiS is produced in the outer envelope of IK Tau. However, recent observational results by Decin et al. (2008b) argue for a formation process in the inner envelope.

From the above arguments it is clear that we should allow for some variation in modeling the abundance structure in the envelope. However, one should also realize that we sometimes only have two rotational transitions of one isotopolog at our disposal with a restricted range in excitation temperature. The highest upper level energy traced is the SiS(20-19) transition at 183 K; hence none of the studied transitions is sensitive to the abundance in the inner envelope (R  5 ). In order to use some prior knowledge on the (theoretical) photo-dissociation rate in the outer regions and to allow for a depletion or an extra formation process in the intermediate/outer envelope, we therefore opted to divide the envelope in different regimes (see also Table 3 and Fig. 1): (i) in the dust-free zone ( ) the abundance is constant ( ); (ii) between and the abundance can decrease/increase from to linearly on a log-log scale, where both  and  are free parameters; (iii) from onwards, the abundance stratification follows the (photodissociation) results of Willacy & Millar (1997) scaled to at  . In that way, three parameters ( , , and ) have to be estimated to determine the abundance stratification of a species.

Table 3:   Modeling assumptions of the abundance stratification.

Most studies use the photodissociation results of Mamon et al. (1988) to describe the CO spatial variation in the outer envelope. For other molecules, the abundance pattern is often assumed to be described by a simple Gaussian or expontential distribution (e.g. Bieging et al. 2000; González Delgado et al. 2003; Schöier et al. 2007a). The e-folding radius then describes the photodissociation by ambient UV photons penetrating the dusty envelope or depletion of the molecules from the gas into dust grains in the outflowing stellar wind. That way, however, all molecules are assumed to be created in the extended atmosphere or inner wind region. Moreover, a combination of depletion and photodissociation or extra depletion/formation processes in the intermediate/outer region cannot be captured, and one cannot use the results by Willacy & Millar (1997) describing extra formation of a few molecules by ion-ion reactions in the outer wind region. The methodology outlined above (Table 3) captures these flaws, and may serve to considerably strengthen our knowledge on the abundance stratification in the envelope.

As was already alluded to in the previous paragraph, the line profiles in this study are not sensitive to a change in abundance in the inner wind region ( ). To assess the abundance stratification in this region, one either needs high-resolution near-infrared (see, e.g., Decin et al. 2008a) or far-infrared spectroscopy (as will be provided by the Herschel/HIFI instrument). Nonetheless, the derived abundance stratifications will be compared to the theoretical inner wind predictions by Duari et al. (1999) and Cherchneff (2006), because this comparison may yield hints on the (un)reactivity of the molecules in the dust-forming region and on uncertainties in the inner wind predictions.

4 Results

With the log-likelihood method the parameters for the model yielding the best-fit to the CO line profiles are derived (see parameters listed in the second column in Table 4, ``model 1''). The CO lines, however, only trace the envelope beyond 100 . One therefore should use other molecules to put constraints on the structure in the inner wind region. HCN is the only molecule for which we have observational evidence that it is formed in the inner wind region: using interferometric data Marvel (2005) deduced a maximum size for the HCN distribution of 3.85 (in diameter), or a radius of 7.2  10<sup>15</sup> cm at 250 pc. They concluded that the deduced size indicates a shock origin for HCN close to the star and a radius for the HCN distribution limited by photodissociation. The HCN line profiles (Fig. 7) are clearly Gaussian, indicating a line formation (at least partly) in the inner wind region, where the wind has not yet reached its terminal velocity. That way, HCN observations yield important clues on the velocity structure in the inner wind region. With the stellar parameters given in the second column in Table 4 (``model 1''), we were unable to derive a HCN-abundance structure yielding a satisfactory fit to the line profiles. While the integrated intensities could be well predicted, the line profiles were flat-topped with a FWHM (full width at half maximum) that was too large. The only way to reconcile this problem was (1) concentrating the HCN abundance in the inner 2  10¹⁵ cm with [HCN/H] = 9  10^-6; or (2) allowing for a smoother velocity law. While the abundance in the former solution is within the predictions of Cherchneff (2006), the angular distance is much smaller than the 3.85  diameter deduced by Marvel (2005).

4.1 Velocity structure

The expansion velocity of SiO, H₂O, and OH masers can be used to put further constraints on the velocity structure, and specifically on the acceleration in the inner wind region (see Fig. 3). It is clear that the velocity structure as derived from the parameters of the best-fit model only based on CO lines (model 1), is far too steep in the inner wind region as compared to the velocity indications of the maser lines. This problem can be solved by either increasing the dust condensation radius or by allowing for a smoother velocity profile. This latter can be simulated with the classical -law (Lamers & Cassinelli 1999) with (see Fig. 3)

with v₀ the velocity at the dust condensation radius.

Table 4:   Parameters of the models with best goodness-of-fit for IK Tau.

One should realize that several assumptions are inherent to the velocity structure derived from solving the momentum equation: (i) all dust species at all different grain sizes are assumed to be directly formed at the dust condensation radius  . However, theoretical results from e.g. Gail & Sedlmayr (1999) show that formation and growth of (silicate) dust grains typically occur between 1100 and 900 K, i.e. extending over a few stellar radii. (ii) The extinction efficiencies used in the GASTRoNOoM-code represent the Fe-rich silicate MgFeSiO₄. Thanks to their high absorption efficiencies at optical and near-infrared wavelength Fe-rich silicates like MgFeSiO₄ (and solid Fe) are efficient wind drivers (Woitke 2006). However, other oxides or pure silicates like Al₂O₃, SiO₂, Mg₂SiO₄ and MgSiO₃ have low absorption efficiencies at optical and near-infrared wavelengths, resulting in a negligible radiative pressure on all glassy condensates. If these latter molecules were the most abundant in the envelope of IK Tau, the wind acceleration would be much lower. No medium resolution infrared (from the Infrared Space Observatory - Short Wavelength Spectrometer or the Spitzer - InfraRed Spectrograph) data are, however, available for IK Tau, hence we were unable to study the circumstellar dust composition in detail. (iii) ``Complete momentum coupling'' is assumed. This means that the grain motion everywhere in the flow can be computed by equating the local radiative and collisional drag forces, implying that virtually all of the momentum gained by a grain through the absorption of radiation from the stellar photosphere is transferred via collisions to the atmospheric gas (MacGregor & Stencel 1992). For physical conditions typical of the circumstellar envelopes of oxygen-rich red giants, MacGregor & Stencel (1992) found that silicate grains with initial radii smaller than about 5  10<sup>-2</sup> m decouple from the ambient gas near the base of the outflow. (iv.) The momentum equation used in the GASTRoNOoM-code (see Eq. (2) in Decin et al. 2006) implicitly assumes that the mass outflow is steady in time and that the circumstellar dust is optically thin to the stellar radiation (Goldreich & Scoville 1976). Dust emission modeling by Ramstedt et al. (2008) suggests a circumstellar envelope which is slightly optically thick at 10 m ( = 1.2). These results suggest that the acceleration of the gaseous particles in the inner wind might be slower than deduced from solving the momentum equation (``model 1'' and ``model 2''), because not all dust species take part in the momentum transfer.

Solving the momentum equation and taking both the CO and HCN line profiles into account, the inner radius of the dusty envelope is shifted towards higher values,   10¹⁴ cm, and the mass-loss rate slightly increases (``model 2'' in Table 4). Using the same stellar parameters as in ``model 2'', but simulating a smoother velocity law still compatible with the maser line mapping ( = 1) decreases the mass-loss rate to 8  10^-6 /yr (model 3 in Table 4), due to the condition of mass conservation (model 3 in Table 4). The narrow Gaussian line profiles of the HCN lines (Sect. 4.3.2) give more support to model 3 than to model 2, which is why the thermodynamic structure as deduced from model 3 (Fig. 4) will be used to model the other molecular line transitions.

4.2 Stellar parameters derived from the CO and HCN lines

The stellar parameters for the best-fit model (model 3) are listed in Table 4. The outer radius of CO is computed using the results of Mamon et al. (1988). The derived 95% confidence intervals in Table 4 are statistical uncertainties, which should be interpreted in the light of the model assumptions of a spherically symmetric wind. As described in Decin et al. (2007), the log-likelihood function can also be used to compare different models with a different number of parameters. For IK Tau, we have assessed the likelihood preference of a model with constant mass-loss rate compared to a model with mass-loss rate variations. The preferences pointed towards the simpler model, i.e. with a constant mass-loss rate (of 8  10^-6 /yr). We also derived the dust-to-gas mass ratio from the amount of dust needed to drive a wind at a terminal velocity of 17.7 km s^-1 for a gas mass-loss rate of 8  10^-6 /yr with the deduced velocity profile. A dust-to-gas mass ratio of 1.9  10^-2 (or a dust mass-loss rate of 1.52  10^-7 /yr) is obtained for model 3, with an estimated uncertainty of a factor 5.

4.3 Fractional abundances

Using the thermodynamic envelope structure derived above (see Fig. 4), the abundance stratification of all molecules is derived. A comparison to the theoretical inner and outer wind predictions (as discussed in Sects. 3.3.1, 3.3.2) is given in Fig. 5. The studied molecular line transitions are not sensitive to the full envelope size, but have a limited formation region. The part in the envelope we can trace by combining the different available rotational line transitions is indicated with vertical dashes in Fig. 5 and tabulated in Table 5.

Table 5:   Molecular fractional abundance relative to H  = (H₂) (see Fig. 5).

A comparison between observed and predicted line profiles and a discussion of the deduced abundance stratification are given for each molecule separately in the following subsections. We will always first briefly describe the deduced abundances, then compare the results to the theoretical inner and outer wind predictions and finally compare to other results found in the literature (see Table 6). For the literature results, a difference is made between studies based on the assumptions of optically thin unresolved emission and a population distribution thermalized at an excitation temperature that is constant throughout the envelope, and those based on a full non-LTE radiative transfer calculation. One should also realize that most studies make use of integrated line intensities, and do not deal with a full line profile analysis as is done here. With the exception of Omont et al. (1993), the other literature studies listed in Table 6 assume the shell to expand at a constant velocity. As discussed in Sect. 4.1, the Gaussian line profiles of the HCN and a few of the SiO lines are the result of line formation partially in the inner wind, where the stellar wind has not yet reached its full terminal velocity. As a result, more emission is produced at velocities near the line center than would be the case for a uniform-velocity wind. Hence observational studies assuming a constant expansion velocity will be unable to predict the line profiles properly.

Table 6:   Comparison of the deduced fractional abundances to other observational studies and theoretical predictions.

4.3.1 CO

Results:

a comparison between the observed rotational ¹²CO and ¹³CO lines and theoretical predictions is shown in Fig. 6. The ¹²CO and ¹³CO lines are very well reproduced by the GASTRoNOoM-predictions, both in integrated intensities and in line shapes. Only the IRAM ¹²CO(1-0) and ¹²CO(2-1) lines are slightly over-predicted. It is, however, not the first time that the non-compatibility of the IRAM absolute flux level is reported (see, e.g. Decin et al. 2008a; Schöier et al. 2006).

Comparison to theoretical predictions:

the CO fractional abundance assumed in all observational studies (see also Table 6) is always lower than the deduced inner wind theoretical non-TE values of Duari et al. (1999) and Cherchneff (2006). The non-TE theoretical values are comparable to the TE-value of 6.95  10^-4 at 1  (Duari et al. 1999). Increasing the CO fractional abundance by a factor 5 would decrease the mass-loss rate by a factor 2.8 to reproduce the observed CO rotational line profiles (with the scaling laws deduced by De Beck et al. 2010).

Comparison to other observational studies:

most observational studies listed in Table 6 assume a fractional CO abundance of [CO/H] = 1-1.5  10^-4. Using different CO rotational lines the mass-loss rate is then derived. The deduced ¹²CO/¹³CO ratio ranges between 9 and 14.

Using CO rotational line transitions, other studies have also estimated the mass-loss rate (see Table 7). The results depend on the assumed or derived temperature distribution, the distance, the adopted [CO/H₂] abundance ratio, and the radiative transfer model or analytical approximation used. All (scaled) mass-loss rate values are in the narrow range between 6.5  10^-6 and 9  10^-6 /yr, the exception being the result of González Delgado et al. (2003), which is a factor 4 higher. We note that the work of these authors was not devoted to the study of CO, and it remains unclear what line intensities they used in their modeling.

Table 7:   Mass-loss rate values derived from ¹²CO rotational line transitions for IK Tau.

4.3.2 HCN

Results:

the narrow Gaussian line profiles clearly point toward (at least) an inner wind origin of HCN. As discussed above, we imposed a value R₂ = 500  to simulate the mapping results of Marvel (2005). The deduced abundance around 50  is [HCN/H] = 2-2.5  10^-7.

Comparison to theoretical predictions:

theoretical predictions by Duari et al. (1999) and Cherchneff (2006) and observational studies by, e.g., Bieging et al. (2000) and Marvel (2005) indicate that HCN forms in the inner wind region of M-type envelopes by shock-induced non-equilibrium chemical processes. This is contrary to the photochemical models of Willacy & Millar (1997), where HCN is only produced in the outer envelope by photochemical reactions involving NH₃ and CH₃, which they assumed to be parent species originating close to the stellar photosphere and injected in the outer region. This formation route, however, leads to line shapes which are clearly non-Gaussian, but are flat-topped (for optically thin unresolved emission) because the wind is already at its full velocity in this region.

The maximum size distribution of 3.85 derived by Marvel (2005) is suggested to be caused by photodissociation of HCN. Bieging et al. (2000) used a parametrized formula to describe the photodissociation radius for HCN as a function of gas mass-loss rate and wind velocity (see their Eq. (2)). This estimate leads to a HCN photodissociation radius of 1.2  10¹⁶ cm, which agrees well with the result of Marvel (2005) of 7.6  10¹⁵ cm (at 265 pc). With HCN as a parent species with an injected abundance of 1.5  10^-7, new chemical outer wind models were calculated withthe code described in Willacy & Millar (1997). The derived photodissocation radius is around 2  10¹⁶ cm, which is a factor 3 higher than the observed value of Marvel (2005).

Compared to the theoretical predictions for TX Cam at 5  by Cherchneff (2006) or for IK Tau at 2.2  by Duari et al. (1999), our deduced abundance is a factor 10 to 40 lower, respectively. There are a few possibilities for the origin of this difference. (i) Contrary to what is thought (e.g. Duari et al. 1999), HCN may participate in the formation of dust grains in the inner envelope. (ii) The formation mechanism of HCN in the inner wind is directly linked to its radical CN by

The destruction route is the reverse reaction. Shocks trigger CN and further HCN formation in the gas (Cherchneff 2006). The formation processes of both molecules depend critically on the physical parameters of the shocked gas, specifically on the physics in the ``very fast chemistry zone'', which is the narrow region after the shock front itself. The modeling of this zone is still subject to many uncertainties (e.g. cooling rate, velocity, shock strength), yielding uncertainties on the theoretical fractional abundances of at least one order of magnitude.

Comparison to observational studies:

all observational deduced values agree within a factor 2, but are clearly lower than the inner wind non-TE theoretical predictions (see Table 6).

4.3.3 SiS

Results:

next to CO, SiS is the molecule with the most molecular line transitions at our disposal. From the SiS(5-4) up to SiS(20-19), excitation temperatures from 13 to 183 K are covered, and one can trace the envelope abundance spatial variations between 40 and 6000 . The lower lying transitions (SiS(5-4) and (6-4)) probe material at larger radial distances from the star; the SiS(19-18) and (20-19) line intensities are sensitive to the abundance in the inner wind region. Modeling the strength of the low and high-excitation lines gives evidence of a depletion of SiS: the abundance at 40  is estimated to be   10^-6 and decreases to   10^-9 at 120 . This is also evident from the extra sub-structure in the high signal-to-noise high-resolution line profiles of the higher excitation lines observed with JCMT and IRAM, where the small peak traces a high abundance component at low velocity close to the star and the broader plateau the lower abundance component further away in the envelope. SiS can be depleted due to the adsorption of SiS molecules onto dust grains before it is photodissociated much farther away. Using ISI 11 m interferometric data Hale et al. (1997) detected dust excess, with a first intensity peak around 0.1 , going out up to at least 0.7 (or 2.8  10¹⁵ cm  180 ) around IK Tau. SiS condensation onto dust species has been postulated already (e.g. Bieging & Tafalla 1993), but not all atomic Si or molecular SiS should condense onto dust species (e.g. Dominik et al. 1993). The different rotational lines give us the potential to constrain the compact, pre-condensation, SiS fractional abundance component quite well, but the uncertainty on the extended post-condensation component is higher.

Comparison to theoretical predictions:

the abundance at 40  is higher than the theoretical value predicted in the inner wind envelope by Duari et al. (1999) and Cherchneff (2006). This may either indicate that the destruction reaction which occurs at larger radii in the inner O-rich envelope via (Cherchneff 2006)

is not as efficient, or that an extra formation route is not yet taken into account in the theoretical modeling. In general, the uncertainties on chemical reaction rates involving sulfur are still very high (I. Cherchneff, priv. comm.). The high abundance around 40  is in accordance with the observational suggestion by Decin et al. (2008b) that SiS forms close to the star, whatever the C/O ratio of the target.

Comparison to observational studies:

similar results concerning the depletion of SiS in the intermediate wind region are obtained by Schöier et al. (2007a). From the simplified analyses assuming optically thin emission thermalized at one excitation temperature (see Table 6) it is not possible to derive this kind of abundance depletion pattern.

4.3.4 SiO

Results:

for the ²⁸SiO isotopolog, five transitions were observed, with excitation energies ranging between 6 K and 75 K. The line strengths and profile shapes of all ²⁸SiO, ²⁹SiO, and ³⁰SiO lines are well predicted, except for the ²⁸SiO(6-5) line as observed with the SEST by González Delgado et al. (2003) (see Fig. 10). Because the strength of both the ²⁸SiO(5-4) and ²⁸SiO(7-6) are well reproduced, and both lines share the line formation region with the ²⁸SiO(6-5) transition, an absolute calibration uncertainty can be the cause of this discrepancy, but time variability of the lines should also be considered (see Sect. 5.1).

Although not as pronounced as for SiS, the modeling of the different rotational transitions indicates an abundance decrease with a factor 40 around 180 . SiO is a parent molecule and a strong candidate to be depleted in the wind of O-rich envelopes: at larger radii in the inner envelope, OH alters the SiO abundance via (Cherchneff 2006)

SiO₂ may condense as silica. It may also participate in the formation of amorphous or crystalline silicates.

The derived isotopic ratios in the envelope are [²⁸SiO/²⁹SiO]  27 and [²⁸SiO/³⁰SiO] = 80, [²⁹SiO/³⁰SiO] = 3. They are discussed in Sect. 5.2.

Comparison to theoretical predictions:

the SiO abundance at the inner dust condensation radius is slightly below the theoretical predictions of Duari et al. (1999) and Cherchneff (2006), which is very reasonable taking the assumptions of both the theoretical chemical kinetic calculations and our modeling into account. It possibly points toward the condensation of SiO onto dust grains in the intermediate wind zone, before 70 , a region where our observed lines are not very sensitive to the exact abundance distribution. The observational study by Decin et al. (2008b) indicates that SiO is formed close to the star, in support of the theoretical predictions by Duari et al. (1999) and Cherchneff (2006). In the outer envelope, SiO is very stable and only photodissociated around a few thousand stellar radii (Willacy & Millar 1997).

Comparison to observational studies:

Lucas et al. (1992) mapped the ²⁸SiO(2-1) v=0 flux distribution, showing that it has a circular geometry. The half-peak intensity radius has a diameter of 2.2 +/- 0.1 (or a radius of 4.35  10¹⁵ cm at 265 pc, in our model being 190 ). The SiO(2-1) emission regions indicate that the final expansion velocity is not yet reached in the SiO emission region, suggesting that grain formation must still take place as far as 10¹⁶ cm from the star. This supports our results on the velocity structure based on the study of the HCN line profiles (Sect. 4.1).

Compared to other observational studies, the deduced abundance around 70  is quite high, while the outer wind abundance agrees with the result by González Delgado et al. (2003).

4.3.5 CS

Results:

for CS, we only have two lines at our disposal, the (6-5) and (7-6) rotational transitions, both with a low S/N-ratio. The fractional abundance at 300  is estimated to be   10^-8.

Comparison to theoretical predictions:

the derived abundance of   10^-8 is much higher than the TE-abundance of 2.5  10^-11 for an oxygen-rich star. Cherchneff (2006) predicts CS to be a parent molecule, with a non-TE abundance for TX Cam around 1.8  10^-5 at 2  and around 2  10^-6 at 5  away in the envelope, while Duari et al. (1999) predicts a value of 2.75  10^-7 at 2.2 . In their theoretical modeling of the carbon-rich AGB-star IRC+10216 Millar et al. (2001) also argue for the need of CS as a parent species to account for the vibrationally excited CS lines detected in IRC+10216.

As for HCN (Sect. 4.3.2), our derived abundances are a factor 50 lower compared to the predictions of Cherchneff (2006) and a factor 8 lower compared to Duari et al. (1999). The dominant formation pathways of both CS and HCN occur in the fast chemistry zone of the gas parcel excursion involving CN. Knowing that this zone is very difficult to model (see Sect. 4.3.2) and that the sulfur reaction rates are not well known (see Sect. 4.3.3), this difference is not so cumbersome. However, as suggested for HCN, it may also be the case that CS or the radical CN are involved in dust formation, altering its abundance in the intermediate wind region. The low-resolution of the two observed CS lines do not provide enough information to firmly prove this. In the outer envelope, CS is first formed from H₂S. Somewhat farther away, the reaction of atomic carbon with SO and HS forms CS, before it is photodissociated by UV radiation (Willacy & Millar 1997).

Comparison to observational studies:

while our deduced CS fractional abundance agrees with the values derived by Bujarrabal et al. (1994) and Kim et al. (2010), the observational result by Lindqvist et al. (1988) is higher by a factor 4.

4.3.6 CN

Results:

the CN lines display a peculiar profile, probably related to the hyperfine structure of the molecule. Although the signal-to-noise of the individual components is low, Kim et al. (2010) noted already that the strength of the different peaks do not agree with the optical thin ratio of the different hyperfine structure compoments and hint to hyperfine anomalies as already reported by Bachiller et al. (1997). Simulations with the GASTRoNoOM-code taking all the hyperfine components into account confirm this result. We therefore opted to simulate both CN lines with the strongest component only. I.e., for the N=3-2, J=5/2-3/2 line we used the F=7/2-5/2 component at 340 031.5494 MHz, for the N=3-2, J=7/2-5/2 line the F=9/2-7/2 component at 340 248.5440 MHz was used.

Due to low signal-to-noise ratio of the lines and the problems with the different hyperfine componets, the derived abundance fractions are loosely constrained. To illustrate this, two model predictions are shown in Fig. 14. For one model, the inner abundance ratio is taken to be 3  10^-8 and from 1000  onward, the abundance stratification follows the predictions by Willacy & Millar (1997) (dashed line in Fig. 5). For the other model, the inner abundance ratio, f₁, is lowered to 1  10^-10 yielding a peak fractional abundance around 2000  of   10^-6, which is a factor 8 higher than the peak fractional abundance derived by Willacy & Millar (1997) (see dotted line in Fig. 5).

Comparison to theoretical predictions:

the predicted inner wind abundance fractions give higher preference to the second model described in the previous paragraph: Duari et al. (1999) predicts a fractional abundance of CN around 2.4  10^-10 at 2.2  for IK Tau, while Cherchneff (2006) gives a value of 3  10^-10 at 5  for TX Cam. Further out in the wind, CN is produced by neutral-neutral reactions and the photodissociation of HCN, yielding peak fractional abundances around 3  10^-7 (Willacy & Millar 1997).

If HCN is indeed photodissociated around 400 to 500  (see Sect. 4.3.2), the peak fractional abundance of CN is not expected to occur around 2000  (see the model predictions by Willacy & Millar 1997, in Fig. 5), but slightly beyond 500  because the photodissociation of HCN is the main formation route to CN in the outer envelope. Shifting the CN peak fractional abundance in the second model to 500  with an abundance value of 1.5  10^-7 also yields a good fit to the (noisy) data.

Comparison to observational studies:

this is the first time that the non-LTE CN fractional abundance for IK Tau has been derived, although the uncertainty on the derived abundance is large.

4.3.7 SO

Results:

as for of CN (Sect. 4.3.6) the noisy profiles prevent a good determination of the abundance stratification. A fractional abundance of   10^-7 is derived around 200 .

Comparison to theoretical predictions:

an abundance stratification compatible with both the inner wind predictions of Cherchneff (2006) and the outer wind model of Willacy & Millar (1997) can be derived, yielding a good representation of both SO lines observed with APEX. Willacy & Millar (1997) assumed no SO injection from the inner wind at large radii, but in-situ formation processes only. Assuming that SO is indeed injected to larger radii can increase the predicted peak fractional abundance computed by Willacy & Millar (1997), which was somewhat too low compared to the observed value listed in their Table 6.

Comparison to observational studies:

for the first time, the SO abundance fraction is derived by a non-LTE radiative transfer analysis, although the low S/N prevents an accurate abundance determination. The result agrees with the LTE analysis by Kim et al. (2010), but is a factor of a few lower than Omont et al. (1993) and Bujarrabal et al. (1994).

4.3.8 SO₂

As for SO₂ the main formation channel in both the inner and outer wind region is

In the outer wind, SO₂ rapidly photodissociates back to SO (Willacy & Millar 1997). The inner wind predictions for TX Cam only yield an abundance of   10^-12, but the detection of SO₂ at 7.3 m in several O-rich giants (Yamamura et al. 1999) would imply a formation site close to the star and an abundance in the range  10^-8-10^-7, slightly below the SO values. Cherchneff (2006) argues that a limited number of reactions in the involved SO₂ formation scheme and the lack of measured reaction rates may explain this discrepancy.

Results:

the availability of 10 different transitions gives some hope that we can shed light on the discussion about the inner wind fractional abundance. However, it turns out that we are unable to fit the high-excitation SO₂( 17_1,17-16_0,16) and SO₂( 13_2,12-12_1,11) lines observed with APEX (see below). The high-excitation SO₂( 14_4,10-14_3,11) and SO₂( 14_3,11-14_2,12) can be predicted quite well. These two lines are narrower than the other SO₂ lines in the sample, indicating that their full formation region is in the inner wind region where the wind has not yet reached its full expansion velocity.

Extensive modeling efforts were made to predict the high-excitation SO₂ observed with APEX. While the high-excitation lines involving the J = 14-14 levels are reasonably well predicted with the models proposed above, the  J =13-12 and J = 17-16 are far too weak. In a study of SO₂ in star forming regions, van der Tak et al. (2003) encountered a similar problem, which they solved by introducing a high temperature, high-abundance component. The increase in abundance could be a factor of 100-1000. Introducing an unrealistically compact, very high-abundance component with f₁ > 1  10^-4 up to 200  reproduces the APEX J = 13-12 and J = 17-16 within a factor 2, but the J = 14-14 line, involving similar excitation levels, is a factor 15 too strong. A possible cause of the discrepancy could be a misidentification of the observed lines. However, different line data bases always point towards an identification as SO₂ transitions. But also the collision rates from the LAMDA database may be problematic. In the LAMDA database, Schöier et al. (2005) extrapolated the collisional rates as computed by Green (1995). Green (1995) computed the collisional rates for temperatures in the range from 25 to 125 K including energy levels up to 62 cm^-1; Schöier et al. (2005) extrapolated this set of collisional rates to include energy levels up to 250 cm^-1 and for a range of temperatures from 10 to 375 K.

Neglecting the SO₂( 17_1,17-16_0,16) and SO₂(13_2,12-12_1,11) lines, the strength of the other eight SO₂ lines can only be explained with a high inner abundance ratio of 1  10^-6, clearly pointing toward an inner wind formation region (in accordance with Yamamura et al. 1999).

Comparison to observational studies:

as for SO, the SO₂ fractional abundance is derived for the first time using a full non-LTE radiative transfer analysis. The derived abundance value is somewhat lower than the results from classical studies, assuming optically thin emission and one excitation temperature, although we have to state clearly that the SO₂ modeling still poses many problems.

5 Discussion

The derived fractional abundances are already discussed in Sect. 4.3. In this section, we focus on the possible time variability of the emission lines and on the derived SiO isotopic ratios. In the last part, H₂O line profile predictions for the Herschel/HIFI mission are performed.

5.1 Time variability

The observed molecular emission lines could be time-variable. Unfortunately, no dedicated study has yet been performed to study the time variability of the molecular lines in IK Tau. Carlström et al. (1990) reported on a monitoring program of the SiS(J = 4-3, 5-4, and 6-5) emission from the Mira-type carbon-rich AGB star IRC +10216. It was found that the circumstellar J = 5-4 and J = 6-5 line emission toward IRC +10216 varies both in line intensity and in line shape. A clear correlation between the variations and the infrared flux (as measured with the K-band magnitude) is found for the J = 5-4 and J = 6-5 lines (Bieging & Tafalla 1993), but not for the J = 4-3 line. This indicates that at least the population of a few levels varies in phase with the stellar flux. A change in the pumping mechanism of the infrared vibrationally excited levels of a molecule will modify the excitation in the ground vibrational state. A change in dust emission may also alter the excitation of a molecule, because dust emission has the potential of affecting the level populations in the ground vibrational state, in particular through pumping via excited vibrational states.

Cernicharo et al. (2000) has looked for time-related intensity variations in a line survey at 2 mm of IRC+10216. Among the 2-mm lines, the most likely lines to be affected are (i) CS, HC₃N, SiO and SiS, four species whose IR lines are known to be optically thick, as well as (ii) the vibrationally excited lines of C₄H and HCN. During the 10-year long run, these lines were observed at several occasions. The ground-state mm lines were found to have stable shapes and intensities (within 20% which is consistent with the calibration uncertainty). Only the strong HCN, , J=2-1 line, which is known to be masering, showed a factor of 2 intensity variation with time.

Currently the effect of time-variability on circumstellar line emission is unknown for AGB stars in general.

5.2 SiO isotopic ratios

The SiO abundance isotopic ratios derived for IK Tau are ²⁸SiO/²⁹SiO = 27, ²⁸SiO/³⁰SiO = 80, and ²⁹SiO/³⁰SiO = 3, with an uncertainty of a factor of 2 due to the low signal-to-noise ratio of the ²⁹SiO and ³⁰SiO lines. The ²⁹SiO/³⁰SiO is similar to the simple ratio of integrated intensities corrected for a difference in transition strength and beam filling factor (the combined effect is a frequency factor of ), but the ²⁸SiO/²⁹SiO and ²⁸SiO/³⁰SiO ratios are a factor 7 and 10 larger, respectively, due to neglect of optical depth effects of the ²⁸SiO lines if the simplified intensity ratio is used. Compared to the solar isotopic ratios of (²⁸SiO/²⁹SiO) = 19.6, (²⁸SiO/³⁰SiO) = 29.8 and (²⁹SiO/³⁰SiO) = 1.52, IK Tau is underabundant in neutron-rich isotopologs or overabundant in ²⁸Si.

Due to the weakness of the ²⁹SiO and ³⁰SiO lines, only few results on the silicon isotopic ratios in the circumstellar envelopes around AGB stars are reported in literature. By fitting the SiO maser emissivity Jiang et al. (2003) estimated the ²⁸SiO/²⁹SiO isotopic ratio in two oxygen-rich Miras, R Cas and NV Aur, to be 29 and 32, respectively. Ukita & Kaifu (1988) measured the relative intensities of the ²⁹SiO/³⁰SiO lines ( J=2-1,v=0), which was 2.4 for the S-type Mira  Cyg (with C/O ratio slightly lower than 1), 1.5 for IK Tau, and 2.9 for the oxygen-rich Mira V1111 Oph. For the carbon-rich AGB star IRC+ 10216, Cernicharo et al. (2000) and He et al. (2008) derived respectively ²⁸SiO/²⁹SiO = 15.4  1.1 (17.2  1.1) and ²⁸SiO/³⁰SiO = 20.3  2.0 (24.7  1.8) and ²⁹SiO/³⁰SiO = 1.45  0.13 (1.46  0.11) from the ratios of integrated line intensities, which is close to the solar values. The results on IRC+10216 are, however, lower limits, because no correction for opacity effects has been done.

Lambert et al. (1987) analyzed high-resolution spectra of the SiO first overtone band around 4 m. They obtained estimates of the atmospheric ²⁸SiO/²⁹SiO abundance ratios for four red giants. For the M-type  Peg and the S-type star HR 1105, the ²⁸SiO/²⁹SiO ratio is close to the solar ratio. ²⁹SiO appears to be underabundant in the MS star o¹ Ori (²⁸SiO/²⁹SiO = 40) and the M-type star 10 Dra (²⁸SiO/²⁹SiO  53). The ³⁰SiO isotope appears to be underabundant by a factor of 2 in all four red giants.

Tsuji et al. (1994) reported on high spectral resolution observations of the 4 m SiO first overtone band in six late-type M giants and two M supergiants. The atmospheric ²⁸Si/²⁹Si and ²⁸Si/³⁰Si and ²⁹Si/³⁰Si ratios in the M giants are always slightly lower than the terrestrial values, i.e. more neutron-rich nuclei tend to be more abundant. This is opposite to the results of Lambert et al. (1987), and assuming that the isotopic ratios are not modified in the circumstellar envelope, the result of Tsuji et al. (1994) is also in contrast to the circumstellar isotopic ratios of M-type giants listed above.

The silicon isotopic ratios are not obviously correlated with other stellar properties. In the literature it is conventional to express the silicon (and other element) isotope ratios in parts per thousand deviation from the solar silicon isotope ratio:

yielding values for IK Tau of -274 and -627, respectively. These values reflect both the initial isotopic composition of the star and possible effects due to nucleosynthesis.

For an AGB star, the initial Si isotopic composition in the envelope is altered by slow neutron capture reactions (s-process) in the He intershell and subsequent third dredge-up (TDU) events, increasing the ²⁹Si and ³⁰Si abundance fractions. A low-metallicity star is expected to have initial ²⁹Si/²⁸Si and ³⁰Si/²⁸Si ratios that are smaller than the solar ratios. The inferred isotopic shifts of the Si isotopes are smaller for an O-rich than for a C-rich AGB star because a C-rich star goes through more dredge-up events, increasing the ¹²C and s-processed material (Vollmer et al. 2008; Zinner et al. 2006). Using two different stellar evolution codes, Zinner et al. (2006) studied the change in silicon isotopic ratios in AGB stars: the shift in Si isotopic ratios and the increase of the ¹²C/¹³C in the envelope during third dredge-up are higher for higher stellar mass, lower metallicity, and lower mass-loss rate, but their predicted silicon isotopic shifts are always much higher than the observational values derived for IK Tau. The mimimum values plotted in their Fig. 6 correspond to the value and at C/O = 1, and are higher than -200. They find that no noticeable changes in the Si isotopes occur when the star is still O-rich. Consequently, the isotopic anomalies in silicon found for several M-type giants probably reflect those of the interstellar medium out of which stars were formed.

For the solar system material, the silicon isotopic ratios are thought to be understood by a mixture of the nuclear products by types I and II supernovae (Tsuji et al. 1994). The predicted silicon isotopic ratios by type II supernovae are around ²⁸Si/²⁹Si  15 and ²⁸Si/³⁰Si  35 according to Hoppe et al. (2009), while Hashimoto et al. (1989) arrive at lower values of 8.9 and 12.6 respectively. Type I supernovae produce mostly ²⁸Si with little ²⁹Si and ³⁰Si (Thielemann et al. 1986). Non-terrestrial silicon isotopic ratios can then be reasonably explained in the same way as for the solar system but by assuming a different contribution of types I and II supernovae.

The silicon isotopic shifts reported here for IK Tau are much lower than values deduced from presolar silicate grains (Mostefaoui & Hoppe 2004; Vollmer et al. 2008), of which the origin spans the range from red giant branch (RGB) and AGB stars up to supernovae. Looking to silicon isotopic ratios derived from presolar SiC grains (Fig. 2 in Zinner et al. 2006), the silicon isotope ratios of IK Tau correspond to the X-grains, which are thought to originate in type II supernovae (but major discrepancies between model predictions and observed isotopic ratios still exist; Nittler et al. 1995). The ¹²C/¹³C ratio inferred for IK Tau (=14) is at the lower limit of the values derived for X-type grains (see Fig. 1 in Zinner et al. 2006).

Hence, if the atmospheric (and circumstellar) silicon isotopic ratio is indeed not changed due to nucleosynthesis and subsequent dredge-ups, the above arguments seem to suggest that the interstellar medium out of which IK Tau was born has a mixture analogous to X-type grains of which supernovae type II are thought to be the main contributors. The measurement of other isotopic ratios can shed new light on this discussion.

5.3 H₂O line profile predictions

Line profile predictions are performed for a few water lines, which will be observed by Herschel/HIFI in the framework of the Guaranteed Time Key Programme HIFISTARS (P.I. V. Bujarrabal) (see Table 8). This key program focusses on the observations of CO, H₂O and HCN lines in a well-selected sample of evolved stars to gain deeper insight into the structure, thermodynamics, kinematics and chemistry of CSEs and into the mass-loss history of evolved stars. The inner wind abundance fraction is assumed to be [H₂O/H₂] = 3.5  10^-4 (Cherchneff 2006). The photodissociation radius is taken from the modeling of Willacy & Millar (1997), which is around 1600 . Applying the analytical formula from Groenewegen (1994) deduced from the results of Netzer & Knapp (1987), a photodissociation radius of 2.8  10¹⁶ cm or 1870  would be obtained. As standard set-up, the Barber H₂O line list (see Appendix A) is used, including 45 levels in both the ground state and first excited vibrational state (the bending mode at 6.3 m).

Table 8:   Line frequencies, upper energy levels and Einstein A-coefficients for the ortho-H₂O lines, which will be observed with Herschel/HIFI.

Description of the line profiles:

most H₂O lines displayed in Fig. 17 have a parabolic shape, characteristic for optically thick unresolved emission. A few lines suffer from self-absorption in the blue wing (e.g., the 8_5,4-9_2,7 line in the ground-state). Particularly for lines where the optical depths at the line centre can be up to 100, effective self-absorption on the blue-shifted side can be seen (e.g., the 1_1,0-1_0,1 and 3_0,3-2_1,2 lines in the ground-state). Lines in the first vibrational state allow one to trace the wind acceleration zone. These lines are considerably narrower than  (e.g., the 8_2,7-7_3,4 line in the -state).

Comparison to other line lists, including or omitting the and  vibrational state

(see Fig. 18): Using the Barber H₂O or the LAMDA linelist yields comparable results for the lines displayed here. Omitting the first vibrational state of the bending mode () yields a decrease in line flux of a few percent up to 60 percent at maximum, depending on the transition. This discrepancy will increase for lower mass-loss rate objects (Maercker et al. 2008). The inclusion of excitation to the first excited vibrational state of the asymmetric stretching mode () yields a change in line flux of 20% at maximum. This agrees with the recent results by González-Alfonso et al. (2007) and Maercker et al. (2009), who found that the inclusion of the -state is particularly important for low mass-loss rate objects. Because the Einstein A-coefficients of the symmetric stretching mode ( = 1) are an order of magnitude lower than for the asymmetric stretching, this state is not likely to affect the models.

Dependence on parameters:

the H₂O line fluxes which will be observed with HIFI are very sensitive to certain (stellar) parameters, and their observations will give us a handle to constrain the circumstellar structure to even higher accuracy. To illustrate this, some simulations are shown in Fig. 19.

• If the temperature structure is approximated with a power law (dotted line in Fig. 19), a differential change is seen for the predicted line fluxes. Observing a few water lines with different excitation levels will pin down the temperature structure in the CSE.

• A second simulation shows the effect of using a velocity structure which is computed from solving the momentum equation (dashed-dotted line in Fig. 19). Applying another velocity law results in another gas number density and slightly different temperature structure. Using the momentum equation, a steeper velocity gradient is obtained (see, e.g., Fig. 3), and the velocity reaches the terminal velocity at shorter distances from the star. This results in slightly broader line profiles and a flux enhancement in the blue wing, because that part of the CSE contributing to the line profile at a certain velocity v will be shifted somewhat inward, hence attaining a higher source function.

• Using a blackbody to represent the stellar radiation instead of a high-resolution theoretical spectrum calculated from a MARCS model atmosphere (see Sect. 3.1) only induces a change in the predicted line fluxes smaller than 2% (not shown in Fig. 19). For wavelengths beyond 200 m the stellar flux is always represented by a blackbody in the GASTRoNOoM-calculations (note that the ground-state of ortho-water is at 23.794 cm^-1 or around 420 m), the flux difference between the blackbody and the theoretical high-resolution spectrum around the  bending mode is shown in Fig. 20. The reason for this negligible difference is that the stellar radiation field is not important (in this case) for the H₂O excitation. Excluding the stellar radiation field only yields a reduction of the line emission by 2% at maximum.

• The dashed line in Fig. 19 shows the model predictions using the same stellar and envelope parameters as in Maercker et al. (2008):  = 2600 K,  = 3.53  10¹³ cm, D = 300 pc and  = 1  10^-5 /yr. This example shows how maser action in the 5_{3 2}-4_{4 1} transition at 621 GHz is very sensitive to the structural parameters.

6 Conclusions

We have for the first time performed a self-consistent, non-LTE radiative transfer analysis on 11 different molecules and isotopologs (¹²CO, ¹³CO, SiS, ²⁸SiO, ²⁹SiO, ³⁰SiO, HCN, CN, CS, SO, SO₂) excited in the circumstellar envelope around the oxygen-rich AGB star IK Tau. In contrast to previous studies, the temperature and velocity structure in the envelope are computed self-consistently, the circumstellar fractional abundances are linked to theoretical outer wind non-chemical equilibrium studies and the full line profiles are used as criteria to deduce the abundance structure. The Gaussian line profiles of HCN and SiO clearly point toward formation partially in the region where the wind has not yet reached its full velocity. With the HCN line profiles as a criterion, we can deduce that the wind acceleration is slower than deduced from classical theories (e.g. Goldreich & Scoville 1976). For a few molecules, a significantly different result is obtained compared to previous, more simplified, studies. SiO and SiS seem to be depleted in the intermediate wind region due to adsorption onto dust grains. The HCN and CS intermediate wind abundance around 50-300  is clearly below the inner wind theoretical predictions by Duari et al. (1999) and Cherchneff (2006), which may either indicate a problem in the theoretical shock-induced modeling or possibly that, contrary to what is thought, HCN and CS do participate in the dust formation, maybe via the radical CN through which both molecules are formed. The lack of high signal-to-noise data for CN and SO prevents us from accurately determining the circumstellar abundance stratification. It turned out to be impossible to model all the SO₂ line profiles, particularly a few of the high-excitation lines. This may be due to a misidentification of the lines or to problems with the collisional rates. The SiO isotopic fractions point toward high ²⁸SiO/²⁹SiO and ²⁸SiO/³⁰SiO ratios, which are currently not understood in the framework of nucleosynthesis altering the AGB isotopic fractions, but seem to reflect the chemical composition of the interstellar cloud out of which the star is born. Finally, in Sect. 5.3, we present H₂O line profile predictions for a few lines which will be observed with the Herschel/HIFI instrument (launched on May, 14 2009).

Appendix A: Molecular line data

For each of the treated molecules, we briefly describe the molecular line data used in this paper. Quite often, data from the Leiden Atomic and Molecular Database (LAMDA) are used (Schöier et al. 2005). When appropriate, transition probabilities are compared to the relevant data in this database.

CO - carbon monoxide.

For both ¹²CO and ¹³CO energy levels, transition frequencies and Einstein A coefficients were taken from Goorvitch & Chackerian (1994). Transitions in the ground and first vibrational state were included up to J = 40. The CO-H₂ collisional rate coefficients at kinetic temperatures from 10 to 4000 K are from Larsson et al. (2002). Figure A.1 shows a good match for the transition frequencies (better than 1%). The Einstein A coefficients for the rotational transitions in the ground-state and the vibra-rotational transitions correspond to better than 0.5%, but the Einstein A coefficients for the rotational transitions in the v=1 state may differ by up to a factor 3 for the high-lying rotational transitions, i.e. the ones with the largest J quantum number. For the CO lines of interest to this study, i.e. rotational transitions in the v=0 state with , the effect on the predicted line fluxes is small: calculating the temperature stratification self-consistently (see Sect. 3.1) using the Einstein A coefficients of the LAMDA-data base, the largest deviation (of 6%) occurs for the ¹²CO(1-0) line.

SiO - silicon monoxide.

The SiO linelist of Langhoff & Bauschlicher (1993) was used to extract the frequencies, energy levels and (vibra-)rotational radiative rates for ²⁸SiO, ²⁹SiO, and ³⁰SiO. Both ground and first vibrational state were included, with rotational quantum number J up to 40. The SiO-H₂ collisional rates in the ground state are taken from the LAMDA-database. For the rotational transitions in the first vibrational state, it is assumed that the collisional rates are equal to those in the ground state. The vibra-rotational collisional rates are assumed to be zero.

The LAMDA database only lists the frequencies and transition probabilities for the first 40 levels in the ground state. Comparison with the line list of Langhoff & Bauschlicher (1993) shows that both databases nicely agree for the rotational transitions in the ground state (see Fig. A.2).

SiS - silicon monosulfide.

Frequencies, transition probabilities and energy levels were taken from the Cologne Database for Molecular Spectroscopy (CDMS, Müller et al. 2005). Transitions in the ground and first vibrational state were included up to J = 40. The pure rotational transitions of silicon monosulfide (²⁸Si³²S) and its rare isotopic species have been observed in their ground state as well as vibrationally excited states by employing Fourier transform microwave (FTMW) spectroscopy by Müller et al. (2007). The LAMDA database only lists the lowest 40 levels in the ground vibrational state, for which the energy levels, transition frequencies and Einstein A coefficients were taken from the JPL catalog (Pickett et al. 1998). The values listed in the JPL and CDMS database correspond wel (see Fig. A.3). As for SiO, the collisional rates for the rotational transitions in the ground state are taken from LAMDA, the collisional rates in the first vibrational state are assumed to be equal to those in the ground state, while the vibra-rotational collisional rates are assumed to be zero.

CS - carbon monosulfide.

Because the LAMDA database only lists the 41 lowest levels in the ground state (up to J = 40, as extracted from the CDMS database), the level energies, frequencies and Einstein A coefficients for the first vibrational state up to rotational quantum number J=40 and the vibra-rotational transitions between v = 1 and v = 0 were extracted from CDMS. The collisional rates are taken from the LAMDA database and are treated as in the case of SiO and SiS.

CN - cyanogen, cyanide radical, .

Transition rates for rotational transitions in the ground vibrational state of the cyanide radical are extracted from the CDMS catalog. The CN energy levels are indicated by three rotational quantum numbers: N is the total rotational quantum numbers excluding electron and nuclear spin, J the total rotational angular momentum including electron spin, and F designating the spin quanta. All energy levels with are included, yielding 235 energy levels and 508 transitions. For lack of anything better, we have applied the H₂-CS collisional deexcitation rates from the LAMDA database (Black & van Dishoeck 1991; Hogerheijde & van der Tak 2000). The typical deexcitation rates are estimated not to be in error by more than factors 2-3 (Black & van Dishoeck 1991).

HCN - hydrogen cyanide.

A plethora of vibrational states of the HCN molecule are relevant for astronomical observations. The states are designated by ( ). The (100) mode is the CH stretching mode at 3311.5 cm^-1; the (010) mode is the doubly degenerate bending mode at 712.0 cm^-1; and the (001) mode is the CN stretching mode at 2096.8 cm^-1. The excitation analysis includes radiative excitation through the stretching mode at 3 m and the bending mode at 14 m. The stretching mode at 5 m includes transitions that are about 300 times weaker and is therefore not included in the analysis. In each of the vibrational levels, we include rotational levels up to J = 29. Hyperfine splitting of the rotational levels were included only in the J = 1 levels, where the splitting is larger than the local turbulent width. l-type doubling in the 14 m tranistions was included. Data used are from Schöier et al. (2007b).

SO - sulfur monoxide, .

While both the JPL and CDMS database (and hence also the LAMDA database for which the values were extracted from the JPL catalog) list the rotational transitions in the v = 0 and v = 1 vibrational state, the transition probabilities for the vibra-rotational transitions could not be found. We therefore have restricted the excitation analysis to the first 70 levels in the ground vibrational state as given by the LAMDA database. We note that the collisional rates as given by the LAMDA database are only listed in the range between 50 and 350 K, and no extrapolations to higher temperatures are provided. For temperatures in the envelope higher than 350 K, the rotational rates at T = 350 K were used.

SO₂ - sulfur dioxide.

Energy levels, frequencies, Einstein A coefficients and collisional rates are taken from the LAMDA database (where they were extracted from the JPL database). The first 198 levels in the ground state are included in the analysis. We note that the collisional rates are taken from Green (1995) and were calculated for temperatures in the range from 25 to 125 K including energy levels up to 62 cm^-1 for collisions with He. This set of collisional rate coefficients, multiplied by 1.4 to represent collisions with H₂, was extrapolated in the LAMDA database to include energy levels up to 250 cm^-1 and for a range of temperatures from 10 to 375 K.

H₂O - water.

The radiative transfer modeling includes the 45 lowest levels in the ground state and first vibrational state (i.e. the bending mode at 6.3 m). Level energies, frequencies and Einstein A coefficients are extracted from the high-accuracy computed water line list by Barber et al. (2006), currently the most complete water line list. For instance, while the LAMDA database lists 158 transitions between the 45 lowest levels in the ground state, the line list of Barber et al. (2006) contains 164 transitions. The Einstein A coefficients for the common transitions agree within 40% (see Fig. A.4).

The H₂O-H₂ collisional rates in the ground state are taken from the H₂O-He rates by Green et al. (1993), corrected by a factor 1.348 to account for collisions with H₂. Rotational collision rates within the first excited state are taken to be the same as for the ground state, while collisions between the ground  and first excited state  are based on the ground state rotational collision rate coefficients scaled by a factor 0.01 (Deguchi & Nguyen-Q-Rieu 1990). In their analysis, Deguchi & Nguyen-Q-Rieu (1990) show that the uncertainties in the vibrational collisional rates have no effect on the calculated H₂O lines.

In Sect. 5.3, we compare the H₂O line profile predictions using the Barber and LAMDA line list and the effect of excluding the bending mode and including the asymmetric stretching mode is discussed.

Acknowledgements

We thank I. Cherchneff for useful discussion on the circumstellar non-TE chemistry, and F. Schöier for providing us with an updated HCN linelist in the LAMDA database. L.D. acknowledges financial support from the Fund for Scientific Research - Flanders (FWO). E.D.B. acknowledges support from the FWO under grant number G.0470.07. H.S.P.M. is very grateful to the Bundesministerium für Bildung und Forschung (BMBF) for financial support aimed at maintaining the Cologne Database for Molecular Spectroscopy, CDMS. This support has been administered by the Deutsches Zentrum fur Luft- und Raumfahrt (DLR).The computations for this research have been done on the VIC HPC Cluster of the KULeuven. We are grateful to the LUDIT HPC team for their support.

References

Footnotes

...1,2

Postdoctoral Fellow of the Fund for Scientific Research, Flanders.

... JCMT

The James Clerk Maxwell Telescope (JCMT) is operated by The Joint Astronomy Centre on behalf of the Science and Technology Facilities Council of the United Kingdom, the Netherlands Organisation for Scientific Research, and the National Research Council of Canada.

... APEX

APEX, the Atacama Pathfinder Experiment, is a collaboration between the Max-Planck-Institut fur Radioastronomie, the European Southern Observatory, and the Onsala Space Observatory. Program IDs are 077.D-0781 and 077.D-4004.

... (IRAM)

IRAM is supported by INSU/CNRS (France), MPG (Germany), and IGN (Spain).

... HCN(3-2)

When the isotopic notation is not given, the molecular line transition is from the main isotopolog.

...

Estimates of these ranges are found by considering the place where I(p) p³, with I the intensity and p the impact parameter, is at half its maximum value.

All Tables

Table 1:   Overview of the molecular line transitions used in this research, with indication of the frequency, the upper energy level, the telescope, the main beam half power beam width (HPBW) and main beam efficienciy ( ).

Table 2:   The fractional abundance (relative to H₂) taken for the parent species by Willacy & Millar (1997).

Table 3:   Modeling assumptions of the abundance stratification.

Table 4:   Parameters of the models with best goodness-of-fit for IK Tau.

Table 5:   Molecular fractional abundance relative to H  = (H₂) (see Fig. 5).

Table 6:   Comparison of the deduced fractional abundances to other observational studies and theoretical predictions.

Table 7:   Mass-loss rate values derived from ¹²CO rotational line transitions for IK Tau.

Table 8:   Line frequencies, upper energy levels and Einstein A-coefficients for the ortho-H₂O lines, which will be observed with Herschel/HIFI.

All Figures

---

Copyright ESO 2010