© ESO, 2010

1. Introduction

Stars of low and intermediate masses (1 M_⊙  ≤ M _⋆  ≤ 9 M_⊙) expel a large part of their outer layers in a dust-driven wind during the later stages of their evolution, when ascending the asymptotic giant branch (AGB). The expulsion of stellar material creates a cool circumstellar envelope (CSE) containing dust grains and molecular gas phase species. The rates at which mass is lost (Ṁ) vary from  ~10^-8 M_⊙  yr^-1 up to  ~10^-4 M_⊙  yr^-1. A violent ejection of the circumstellar envelope ends the AGB phase and initiates further evolution into a protoplanetary nebula (PPN) and later into a planetary nebula (PN), when the ejected material is ionised (Habing & Olofsson 2003).

High-mass stars (M _⋆   >9  M_⊙) do not ascend the AGB, but evolve into supergiants and hypergiants. Similar to AGB stars, they are typified by circumstellar envelopes of dusty and molecular material created by strong winds. The mechanisms behind the dense outflows of both types of objects are likely to be different (Josselin & Plez 2007), since the central stars have obtained other physical and chemical properties after going through somewhat different phases of evolution.

We want to investigate the physical and chemical properties of these CSEs by analysing molecular-line observations. The most important species for such investigations is carbon monoxide, CO. More specifically, observations of lines resulting from transitions from states with angular momentum quantum number J to J−1 can be used to derive density, velocity and kinetic temperature in the envelope. Previous surveys have focussed on the analysis of only a few low-excitation rotational transitions of CO (J = 1−0, 2−1, 3−2), or have used simple analytical expressions to derive the mass-loss rate from the parameters of these transitions (Knapp et al. 1982; Knapp & Morris 1985; Olofsson et al. 1993; Loup et al. 1993; Neri et al. 1998; Ramstedt et al. 2008). Since we have access to an extensive data set covering both low and high-J molecular lines, we want to provide mass-loss estimators based on multiple CO lines, including those of high excitation levels. These analytical expressions will link stellar and CO-line parameters to the (variable) mass-loss rates of the central stars.

The presence of both ¹²CO and ¹³CO line transitions in the data set allows us to estimate the ¹²C/¹³C isotope ratio for many stars in the sample. This abundance ratio provides a measure for the chemical evolution and hence of the evolutionary state of the stars. In this paper we give first order estimates for ¹²C/¹³C, derived from line intensity ratios.

The sample of evolved stars and the CO-line observations are presented in Sect. 2 and in the (online) appendix. We discuss some of the sample stars individually based on the observations. Section 3 deals with the radiative transfer analysis and the parameter study that will lead to the construction of analytical expressions to estimate mass-loss rates. The constructed formalism is then compared to those presented in the literature. In Sect. 4 and Appendix C, we discuss the determination of the basic stellar parameters needed to derive mass-loss rates. The results of our study in terms of Ṁ, ¹²C/¹³C and wind driving efficiency are presented in Sect. 5. Conclusions are given in Sect. 6.

2. Observations

2.1. The sample

The data for our sample of 69 galactic objects have been assembled over many years and mainly consists of targets originally selected as potential candidates to be included in guaranteed time key programs (GTKP) of the Herschel-HIFIand Herschel-PACSinstruments. Different evolutionary and chemical types are covered, but there is a rather strong bias towards oxygen-rich (C/O  <  1) AGB stars. Several categories of AGB stars, differing in pulsational and mass-loss properties, are represented in the sample, i.e. Mira-variables, semi-regular (SR) variables and OH/IR-type stars. SRs and Miras have low to intermediate mass-loss rates (10^-8−10^-5 M_⊙  yr^-1), while OH/IR stars have very dusty and optically opaque envelopes, formed by intense mass loss that can reach up to a few 10^-4 M_⊙  yr^-1. Typical for these objects are the strong OH masering at 1612 MHz and a very large infrared excess. Miras have very regular pulsations, with a fixed period of a few 100 days (average  ~400 days) and large amplitudes (≥2.5 mag in V-band). Semi-regulars, on the other hand, have smaller amplitudes and can exhibit irregularities in their pulsations. Their pulsation periods can range from 20 up to 2000 days (GCVS Samus et al. 2004). OH/IR stars have pulsation periods ranging from a few 100 days up to more than 1000 days with an average of  ~1000 days and they could be considered the evolutionary successors of Miras (Vassiliadis & Wood 1993). Their thick envelopes, obscuring the central star in the optical, are produced by the so-called superwind phases (Iben & Renzini 1983; Vassiliadis & Wood 1993), which occur during the late parts of the quiescent hydrogen burning phase, just before the central star goes through a thermal pulse (TP). During these phases, the pulsation period increases with increasing mass loss (up to a few 10^-5 M_⊙  yr^-1) and luminosity. The enhanced mass loss leads to a very optically thick and dusty CSE, characterising an OH/IR star. When a superwind phase is halted, the pulsation period decreases again, the expelled material diffuses outwards, and a less intense mass-loss process is initiated (Vassiliadis & Wood 1993). The combination of these factors can lead to the central star being again visible in the optical, and being again observed and classified as a Mira or SR variable. This scenario is more plausible for the less massive stars (M _⋆   ≤ 2.5 M_⊙), as the higher-mass stars (M _⋆   > 2.5 M_⊙) experience relatively modest variations in their high mass-loss rates, luminosities, and pulsation periods due to the TPs (see Figs. 3 to 9 in Vassiliadis & Wood 1993). Low-mass stars experience only very short periods of enhanced mass loss. Moreover, during these periods the absolute value of the mass-loss rate is significantly lower than for their massive counterparts. Therefore, the chance to observe them in their superwind phases, i.e. as OH/IR, is far smaller than for the higher-mass stars.

Supergiants and hypergiants and some post-AGB objects are also included in the sample. The latter fit in this study as the progeny of AGB stars. Studying their extended envelopes will lead to an improved knowledge of the envelopes of AGB-type stars and the late stages of the AGB evolution. Two young stellar objects (YSOs) – AFGL 5502 and the Gomez Nebula – were observed together with the sample of evolved stars. The obtained data on these YSOs were not yet published and are presented in this paper. We will not further discuss these objects, since the focus of this paper is on evolved stars.

2.2. The observations

The data set presented in this paper consists of observations of multiple rotationally excited lines of both ¹²CO and ¹³CO in 69 stars. The bulk of the presented lines were first published by Kemper et al. (2003). Transitions from J = 1−0 up to J = 7−6 for ¹²CO and from J = 2−1 up to J = 6−5 for ¹³CO are covered. Since the high-J lines have formation regions deeper within the CSEs, this multitude of rotational lines can sample a much larger part of the CSE than only the low-excitation transitions. Also, because of slightly different molecular properties, the ¹³CO-data probe regions somewhat different from those sampled with the ¹²CO-rotational lines. Moreover, since the abundance of ¹³CO is lower than that of ¹²CO, the circumstellar layers are optically thin(ner) for the rotational lines of the former, implying that ¹³CO-rotational lines provide good diagnostics for the envelope’s density structure. In this respect, we should be able to obtain a more complete view on the structure and properties of the envelopes and the mass-loss history of the objects, than was possible before with single-dish data only sampling low-excitation ¹²CO-lines.

All previously unpublished data in our sample were obtained with APEX¹ (Atacama Pathfinder EXperiment) and JCMT² (James Clerk Maxwell Telescope). Other data were obtained via private communication or retrieved from the literature.

APEX is a 12 m single-dish telescope with a frequency range from 210 up to 1500 GHz located at Llano Chajnantor, Chile. CO-line data were obtained with three heterodyne SIS-receivers mounted on the Nasmyth-A focus: APEX-2A, FLASH-I and FLASH-II. Instrument specifics are listed in Table 1, together with the CO-lines observable in the respective frequency ranges. All observations were performed in beam-switching mode.

Reduction of the APEX data was done with CLASS, part of the GILDAS-package. The process mainly consisted of combining different scans of the same molecular line towards one object into a single spectrum and subsequently removing baselines and spikes in the obtained spectra. Correcting the data for telescope efficiencies was carried out by the APEX online calibrator.

JCMT data were obtained with heterodyne receivers A, B, HARP, RxW (bands C and D), and E in beam-switching mode. Data reduction was performed with the reduction package SPECX. Spectra of standard stars taken during the observing runs were compared to the standard spectra available on the web site of JCMT. In case of large deviations in the measured line intensities of these standard stars (>10%), the scientific-programme data obtained around the time of the measurement of these standard stars were corrected for these discrepancies. In case of multiple available standard spectra with other magnitudes of deviation, the correction factor was determined via linear interpolation in time. As for APEX data, most corrections and conversions of the data are performed online. The correction for the main-beam efficiency η_MB had to be done manually in the course of the reduction process using the values given in Table 1. η_MB allows transforming antenna temperatures, T_A, into main-beam-brightness temperatures, T_MB. The latter is the equivalent of the brightness temperature of measurements performed with a perfect antenna outside the earth’s atmosphere. The main advantage of a T_MB-based temperature scale is that the data are no longer dependent of any of the instrumental properties, apart from the beam width. All data in this paper are presented on the T_MB-scale.

The uncertainties on the JCMT data, so-called absolute errors, are fairly well established for all receivers (Kemper et al. 2003). For APEX this is not so clear. No standard spectra were available to check the performance of the instruments at the time of the execution of the scientific programme and no standard values for the uncertainties on the data were found in the literature. Ramstedt et al. (2006) mention an uncertainty of 20% on the absolute intensity scale for the J = 3−2 transition, but give no reference for this number.

2.3. Results

Tables A.1 and A.3 in Appendix A (available online) contain main-beam-brightness temperatures at the centre of the line profile (T_MB,c), velocity integrated main-beam intensities (I_MB), and expansion velocities (v_∞) for all CO-data of the sample stars. Among the high-quality CO-data for over fifty stars, there are ¹³CO data for 29 targets.

Some of the observed lines towards AGB or RSG stars in the sample were not detected, e.g. ¹²CO (3−2) for α Sco (=IRAS 16262–2619), and some others were too noisy or heavily contaminated by interstellar lines, e.g. V1360 Aql (=IRAS 18432–0149). In these cases, the data are presented in the appendix, but are not used in the further analysis of the sample. The respective targets are listed in Table 5 together with Post-AGB objects and YSOs.

Characterisations of the sample targets in terms of chemical type (oxygen-rich, carbon-rich or S-type) and pulsational and/or evolutionary type (e.g. Mira, SR, OH/IR) are listed in Table 6.

2.3.1. Line profile diagnostics

Because of its high stability against photodissociation and its large molecular abundance, carbon monoxide exists throughout the largest part of the CSE. Rotational transitions of CO are easily excited in the cool envelopes and can therefore trace many properties of all different layers, e.g. density and temperature. These properties are linked to the optical thickness of the envelope layers and the mass-loss rate of the central star.

In case of optically thick, unresolved envelopes, the observed line profiles have parabolic shapes. Changing to the optically thin and/or spatially resolved case, the molecular lines are more flat-topped or even display two-horned profiles (Knapp & Morris 1985). Considering this simple rule, the shape of the line profiles is a useful first diagnostic for the (density) properties and the geometric extent of a star’s CSE relative to the telescope beam. The data of, e.g., IK Tau shown in Fig. 1 clearly reflect that the envelope layers are optically thinner for ¹³CO than for ¹²CO, since the line profiles of the latter are consistently more parabolic in shape. The mass loss of CW Leo is very strong (~2 × 10^-5 M_⊙  yr^-1; Mauron & Huggins 2000), causing optically thick CO emission lines. These are, however, flattened due to the large angular size of the envelope (~200 arcsec; Mauron & Huggins 2000), causing resolution by the telescope beams with half power beam widths between 7 and 20 arcsec (see Table 1).

A first analysis of the data consists in fitting a so-called soft parabola (Olofsson et al. 1993) to the observed line profiles. Every rotational line is specified by several line parameters: the main-beam-brightness temperature at the line centre, T_MB,c, the velocity at the line centre, i.e. the velocity of the star with respect to the Local Standard of Rest, v_LSR, and the half width of the line profile, i.e. the expansion velocity of the CSE in the outermost regions where the studied molecule is present, v_∞. These parameters are obtained by fitting the data with the soft parabola line profile function, given in Eq. (1), where β is a measure for the shape of the profile function: $\begin{array}{ccc}\mathit{T}\mathrm{\left(}\mathit{v}\mathrm{\right)}& \mathrm{=}& {\mathit{T}}_{\mathrm{MB}\mathit{,}\mathrm{c}}{\left[\mathrm{1}\mathrm{-}{\left(\frac{\mathit{v}\mathrm{-}{\mathit{v}}_{\mathrm{LSR}}}{{\mathit{v}}_{\mathrm{\infty }}}\right)}^{\mathrm{2}}\right]}^{\mathit{\beta }\mathit{/}\mathrm{2}}\mathit{,}\mathrm{|}\mathit{v}\mathrm{-}{\mathit{v}}_{\mathrm{LSR}}\mathrm{|}\mathrm{\le }{\mathit{v}}_{\mathrm{\infty }}\\ & \mathrm{=}& \mathrm{0}\mathit{,} \mathrm{|}\mathit{v}\mathrm{-}{\mathit{v}}_{\mathrm{LSR}}\mathrm{|}\mathit{>}{\mathit{v}}_{\mathrm{\infty }}\mathit{.}\end{array}$(1)When β = 2, the line profile has a parabolic shape, representative for optically thick lines, observed towards spatially unresolved CSEs. Smaller, positive values for β lead to more flat topped profiles. Negative β values can be used to fit two-horned profiles observed towards optically thin envelopes. The best fit to the line profile was determined through minimising the total absolute difference between the data and the soft-parabola fit with the IDL-routine AMOEBA.

Figs. A.1 and A.2 show all ¹²CO and ¹³CO observations and the soft-parabola fits to these data. The β-values per line profile are listed in Tables A.1 and A.3. The parameter β provides information on the optical thickness of the envelope (and hence on Ṁ ) and the resolving power of the instrument. Its value will be used to improve the accuracy of the mass-loss determination (see Sect. 3.4).

2.3.2. Discussion on individual targets

If strong deviations from the soft-parabola fit are found, this fitting procedure could reveal significant deviations from the assumptions of, e.g., spherical symmetry or constant Ṁ, since detached shells, variability in the mass loss, vast clumps in the envelope, or jets can affect the symmetry and the overall shape of the CO lines. Some of the targets in the sample and their deviations from the soft-parabola fits are discussed individually in this section.

It is obvious from Figs. A.1 and A.2 that indeed not all profiles in the data set can be fitted with the simple type of line profile function presented in Eq. (1). R Cas (Mira) and V Hya (SRa) have line profiles exhibiting superimposed peaks – see Fig. 2 – that could be linked to bipolar outflows (Olofsson et al. 1993; Sahai et al. 2003). Olofsson et al. (1993) suggest that these bipolar outflows could be related to the possible presence of a binary companion in case of both targets. Sahai et al. (2003) suggest that V Hya could already be transitioning to the Post-AGB phase. The high-velocity bipolar outflow would be linked to this transition process.

U Ant and S Sct are discussed by Olofsson et al. (1993) as targets with convincing evidence in the CO data for highly episodic mass loss, possibly due to a thermal pulse. The shapes of their line profiles – see Fig. 2 – indeed imply the presence of a detached shell. The inner parts of the outflow have terminal velocities of respectively 7 km s^-1 and 8 km s^-1, while the outflow velocities of the outer parts of the envelopes reach 16 km s^-1 and 20 km s^-1. The carbon-rich semi-regular R Scl also has a detached shell (Olofsson et al. 1996), albeit not directly visible in the CO line profiles.

The CO line profiles of the SRb variable EP Aqr (Fig. 3) reveal the composite nature of its circumstellar environments. Except for the low-S/N observation of ¹³CO(4–3), all measured line profiles towards EP Aqr exhibit a composite profile with (a) a broad, low-T_MB profile and (b) a very narrow, high-T_MB profile centrally superposed on this broad profile. Only very little is known about the origin of these types of line profiles. The broad plateau emission has been suggested to originate from episodic mass loss, bipolar outflows, or circumstellar disks (Knapp et al. 1998; Kerschbaum & Olofsson 1999). Winters et al. (2007) mention that no obvious departure from circular symmetry can be seen, but that there is evidence for a ringed structure in the ¹²CO(2–1) map implying variation of Ṁ in time.

R Hya is well-known for its declining period and mass-loss rate (Decin et al. 2008a), which Wood & Zarro (1981) attributed to a possible recent thermal pulse. The detached shell that has been detected at 60 μm can be explained by a slowing down of the wind. Decin et al. (2008a) have modelled the wind of R Hya in detail, checking models produced with the radiative transfer code GASTRoNOoM against data of both rotational and vibration-rotation transitions of CO. WX Psc is a second object for which Ṁ could be considered to be under influence of TPs. This extreme oxygen-rich star could currently be going through a superwind phase (Decin et al. 2007), causing the large infrared excess. For both targets, however, these presumed deviations from constancy in mass-loss rate are not directly visible in the rotational CO line profiles.

The triply peaked low-J transition profiles in both ¹²CO and ¹³CO towards the supergiant VY CMa – see Fig. 2 – are most probably a superposition of an optically thin component (red-wing and blue-wing peaks) and an optically thick one (middle peak), as shown by Ziurys et al. (2007). For the higher-J transitions (J = 4−3 and up), the extra parabolic feature is no longer clearly present. This, in combination with the asymmetry of the profile, led to the assumption of a variable mass-loss history as described by Decin et al. (2006). NML Cyg and α Ori both have spiked line profiles, somewhat similar to those of VY CMa. In the case of NML Cyg, the lines could be contaminated with interstellar CO detection. The line profiles of IRC +10420 and μ Cep exhibit asymmetries, with a blue-wing feature superimposed on a more or less parabolic profile. The extra feature is far more pronounced for μ Cep and is also further out to the blue part of the line profile.

For some of the more evolved objects in our sample, i.e. Post-AGB stars, we also see deviations from the soft-parabola shape. The ¹²CO(3−2) line profile towards GSC 08608-00509 exhibits a sharp red-wing peak. For the Cotton Candy Nebula, we see a similar sharp feature in the blue wing. Since we have only few observations for these targets, it is hard to draw well-founded conclusions about the properties or geometry of their extended envelopes.

All three line profiles towards the Gomez Nebula are two-horned, supporting the hypothesis of a CO emitting disk in keplerian rotation, as mentioned by Bujarrabal et al. (2008, their Fig. 3). The CO line profiles towards the other young stellar object, AFGL 5502, bear evidence of asymmetry, since the higher-excitation lines have superimposed features in their blue wings.

We note here that T Cet was not included in the sample of red supergiants, contrary to what is often found in literature (e.g. Sloan & Price 1998), but in the sample of AGB stars. The low luminosity and the dust features in the ISO spectrum make classification of T Cet as AGB star more plausible (Speck et al. 2000; Verhoelst et al. 2009).

3. Radiative transfer analysis

Physical information on the envelope structure of evolved stars can be extracted from the observed line fluxes using a non-local thermodynamic equilibrium (non-LTE) radiative transfer code (see Sect. 3.1 for a short description). The availability of a large ¹²CO and ¹³CO line flux data base, as presented in Sect. 2.3, offers interesting possibilities for statistical analyses. However, this large database also forces us to construct simplified analytical expressions relating the mass-loss rate and the observed integrated CO-line intensities to deduce the mass-loss history of each target. We therefore derived new mass-loss rate formulae based on a grid computation. Similar formalisms that have been presented in the literature are summarised in Sect. 3.2. The grid construction and the resulting analytical expressions are discussed in Sects. 3.3 and 3.4. A discussion on the quality of the formalism can be found in Sect. 3.5.

3.1. Radiative transfer model

The structure of circumstellar envelopes is thought to be quite complex, due to the possible presence of shocks, bipolar outflows, or episodic mass-loss events. The complex heating and cooling processes due to excitation of molecules, and the unknown grain-size distributions contribute to the complexity of the envelopes. The often taken approach to model these envelope structures is therefore to choose a simplified approach with reasonable assumptions. We use the non-LTE radiative code GASTRoNOoM, which was presented and discussed by Decin et al. (2006, 2007, 2008a) and benchmarked according to the method described by van Zadelhoff et al. (2002), considering a simple 2-level molecule (their cases 1a and 1b). The CSE is assumed to be spherically symmetric and formed by a constant mass-loss rate. The code 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 Justtanont et al. (1994). Adopting a value for the gas mass-loss rate, the dust-to-gas ratio is adjusted until the required terminal velocity is obtained (Decin et al. 2006). The statistical equilibrium equations and the radiative transfer equation are solved in the co-moving frame using the Approximated Newton-Raphson operator as developed by Schönberg & Hempe (1986). Finally, a ray tracing program uses the computed non-LTE level populations to calculate the line profiles. For a full description of the GASTRoNOoM-code, we refer to Decin et al. (2006).

As described by Decin et al. (2006, 2007, 2008a), a detailed radiative transfer modelling of one target often results in the computation of many (>10⁵) models. Some prior knowledge on the input parameters helps to restrict the extensive parameter space. The main input parameters for this research are the stellar temperature T _⋆ , the stellar radius R _⋆ , the distance d, the CO fractional abundance f_H,CO = CO/H, the ¹²C/¹³C isotope abundance ratio, the terminal velocity v_∞, the mass-loss rate Ṁ, and the dust condensation radius R_inner. When each of these eight parameters would only be scanned over five different values, almost 400 000 models should be computed. Performing an in-depth line profile analysis for each target presented in this study is therefore too time-consuming and beyond the scope of this paper. The mass-loss rates of some twenty targets of the current sample will be individually modelled in a forthcoming paper (Lombaert et al., in prep.).

Before describing the grid computation in Sect. 3.3, we give a short overview of mass-loss rate formulae presented in the literature. This overview is for sure not complete, but provides a few often applied formulae, relevant to this study. These formalisms are also compared with the one we present in this paper.

3.2. Comparison with other studies

Several theoretical mass-loss rate formulae based on observed line intensities or fluxes were already presented in the literature. Knapp & Morris (1985) derived a formula to estimate Ṁ in terms of the CO(J = 1−0)³ line intensity. Their formula relates the mass-loss rate to an assumed CO abundance and a number of easily determined observables: the main beam temperature T_MB,c at the line centre (in K), the terminal velocity v_∞ (in km s^-1), and distance d (in pc). A distinction was made between optically thick and optically thin envelopes. The model for the CO emission assumed a spherically symmetric envelope expanding at constant velocity. The kinetic temperature profile throughout the envelope was assumed in all cases to be that derived by Kwan & Linke (1982) for CW Leo. The statistical equilibrium equations were solved for the lowest 15 rotational levels. Olofsson et al. (1993) transformed the formula of Knapp & Morris (1985) for an optically thick envelope to the slightly more general form $\mathit{Ṁ}\mathrm{KM}\mathrm{=}\mathrm{5.7}\mathrm{\times }{\mathrm{10}}^{-20}\frac{{\mathit{T}}_{\mathrm{MB}\mathit{,}\mathrm{c}}{\mathit{v}}_{\mathrm{\infty }}^{\mathrm{2}}{\mathit{d}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{b}}^{\mathrm{2}}}{\mathit{s}\mathrm{\left(}\mathit{J}\mathrm{\right)}{\mathit{f}}_{\mathrm{H}\mathit{,}\mathrm{CO}}^{\mathrm{0.85}}}{\mathit{M}}_{\mathrm{\odot }}\mathit{/}\mathrm{yr}\mathit{,}$(2)where θ_b is the beamwidth of the instrument (in arcsec) used for the observation of the line under study. The factor s(J) depends on the rotational transition J → J−1. For the J = 1−0 transition s(J) = 1, and for the J = 2−1 transition s(J) = 0.5 (see the discussion in Olofsson et al. 1993). In this equation – valid in the optically thick case – the CO envelope size was fixed at 3 × 10¹⁷ cm, appropriate for high-Ṁ targets. Schöier & Olofsson (2001) compared mass-loss rates derived using their Monte Carlo non-LTE radiative transfer model to those estimated using the formulae of Knapp & Morris (1985). They found that the formula underestimated the mass-loss rates when compared to those obtained from their radiative transfer analysis on average by a factor of about four, and the discrepancies were found to increase for lower mass-loss rate objects (Schöier & Olofsson 2001, see their Fig. 8). In the case of optically thin envelopes, the exact envelope size matters in estimating Ṁ, so Eq. (2) is indeed of lesser value. A correction for the envelope size (see Neri et al. 1998, and references therein) would decrease the discrepancy found by Schöier & Olofsson (2001). However, considering the simplicity of Eq. (2), the agreement is still remarkably good.

Loup et al. (1993) calculated the envelope size R_CO and the mass-loss rate Ṁ_L by numerically solving the system of Eqs. (3) and (4). $\begin{array}{ccc}\mathit{Ṁ}\mathrm{L}& \mathrm{=}& \frac{{\mathit{v}}_{\mathrm{\infty }}^{\mathrm{2}}{\mathit{d}}^{\mathrm{2}}{\mathit{T}}_{\mathrm{MB}\mathit{,}\mathrm{c}}}{\mathrm{2}\mathrm{\times }{\mathrm{10}}^{\mathrm{3}}\mathrm{(}\mathrm{2}{\mathit{f}}_{\mathrm{H}\mathit{,}\mathrm{CO}}{\mathrm{)}}^{\mathrm{0.85}}}\mathrm{\times }{\mathit{F}}^{-1}\mathrm{\left(}{\mathit{R}}_{\mathrm{CO}}\mathrm{\right)}{\mathit{M}}_{\mathrm{\odot }}{\mathrm{yr}}^{-1}\\ \mathit{Ṁ}\mathrm{L}& \mathrm{=}& {\mathrm{10}}^{-6}{\left(\frac{{\mathit{R}}_{\mathrm{CO}}}{\mathrm{7.3}\mathrm{\times }{\mathrm{10}}^{\mathrm{16}}}\right)}^{\mathrm{1}\mathit{/}\mathrm{0.58}}{\left(\frac{{\mathit{v}}_{\mathrm{\infty }}}{\mathrm{10}}\right)}^{\mathrm{1}\mathit{/}\mathrm{1.45}}\\ & & \mathrm{\times }{\left(\frac{\mathrm{4}\mathrm{\times }{\mathrm{10}}^{-4}}{\mathrm{2}{\mathit{f}}_{\mathrm{H}\mathit{,}\mathrm{CO}}}\right)}^{\mathrm{1}\mathit{/}\mathrm{1.16}}{\mathit{M}}_{\mathrm{\odot }}{\mathrm{yr}}^{-1}\mathit{.}\end{array}$They assumed f_H,CO of 2.5 × 10^-4 for oxygen-rich stars and 5 × 10^-4 for both carbon-rich and S-type stars. The function F(R_CO) introduces a correction factor for the envelope size and is normalised to unity for R_CO = 3 × 10¹⁷ cm.

Ramstedt et al. (2008) derived a formula relating the mass-loss rate to the integrated intensity I_CO,J of the ¹²CO (J → J−1) lines, the terminal velocity v_∞, the abundance f_H,CO, the distance d to the object, and the telescope’s beam size θ_b. Using the Monte Carlo radiative transfer code described by Schöier & Olofsson (2001), the estimator was derived from a grid of 60 models varying in Ṁ (1, 3, 10, 30, and 100  × 10^-7 M_⊙  yr^-1), v_∞ (5, 10, 15, and 20 km s^-1), and f_H,CO (1, 3, and 10  × 10^-4). The distance was taken to be 1000 pc. For each model, the velocity integrated intensities in CO (J = 1−0, 2−1, 3−2, 4−3) were calculated for a 20 m telescope (with corresponding beam sizes of 33, 16.5, 11, and 8.5). They use the so-called h-parameter to describe the free parameters of the dust (Schöier & Olofsson 2001): $\mathit{h}\mathrm{=}\left(\frac{\mathrm{\Psi }}{\mathrm{0.01}}\right)\left(\frac{\mathrm{2.0} \mathrm{g}\mathrm{c}{\mathrm{m}}^{-3}}{{\mathit{\rho }}_{\mathrm{d}}}\right)\left(\frac{\mathrm{0.05}\mathit{\mu }\mathrm{m}}{{\mathit{a}}_{\mathrm{d}}}\right)\mathit{,}$(5)where Ψ is the dust-to-gas mass ratio, and ρ_d and a_d are the average density and the radius of an individual dust grain, respectively. The value of h was assumed to be 0.2 for Ṁ in the range up to 3 × 10^-7 M_⊙  yr^-1, 0.5 for 1−3 × 10^-6 M_⊙ yr^-1, and 1.5 for 10^-5 M_⊙  yr^-1. The stellar temperature T _⋆ , stellar radius R _⋆  or inner dust condensation radius R_inner were not specified. An estimate for the mass-loss rate was determined from a minimisation procedure applied to Eq. (6): $\mathit{Ṁ}\mathrm{R}\mathrm{=}{\mathit{s}}_{\mathit{J}}\mathrm{(}{\mathit{I}}_{\mathrm{CO}\mathit{,J}}{\mathit{\theta }}_{\mathrm{b}}^{\mathrm{2}}{\mathit{d}}^{\mathrm{2}}\mathit{/}{\mathrm{10}}^{\mathrm{6}}{\mathrm{)}}^{{\mathit{a}}_{\mathit{J}}}{\mathit{v}}_{\mathrm{\infty }}^{{\mathit{b}}_{\mathit{J}}}{\mathit{f}}_{\mathrm{H}\mathit{,}\mathrm{CO}}^{\mathrm{-}{\mathit{c}}_{\mathit{J}}}{\mathit{M}}_{\mathrm{\odot }}{\mathrm{yr}}^{-1}\mathit{,}$(6)with I_CO,J in K km s^-1, v_∞ in km s^-1, d in pc, and θ_b in arcsec. A comparison between their input mass-loss rates and the values found from Eq. (6) was not shown, nor were estimates of the uncertainties on the derived Ṁ-values provided. The comparison between mass-loss rates derived using Eq. (6) and mass-loss rates derived using their Monte Carlo code for a sample of C-, M- and S-stars shows quite a good agreement. Close inspection, however, reveals that for low-J transitions (J = 1−0, 2−1) Eq. (6) tends to underestimate the mass-loss rate for Ṁ ≲ 10^-6 M_⊙  yr^-1 and to overestimate it for Ṁ ≳ 10^-5 M_⊙  yr^-1 (see Fig. A.1 in Ramstedt et al. 2008).

Since (1) our data set holds observations up to the high excitation J = 6−5 and J = 7−6 rotational line transitions of CO, and we want (2) to include the shape of the line profiles, represented by parameter β, in the estimates of Ṁ, and (3) to extend the computations to higher mass-loss rate values (Ṁ  > 1 × 10^-5 M_⊙ yr^-1), we elaborated on the studies discussed above by performing a grid computation covering more than 15 000 models.

3.3. Description of the grid

To study the influence of different stellar and envelope parameters – like stellar temperature, stellar radius (or stellar luminosity) and dust condensation radius – on the CO line intensities, the input parameters were varied to cover the relevant parameter space for evolved low-mass targets (see Table 2). The values for the f_H,CO were chosen to be 1 × 10^-4 for O-rich stars (Kahane & Jura 1994), 3 × 10^-4 for S-type stars (Ramstedt et al. 2006), and 5 × 10^-4 for C-rich stars (Zuckerman & Dyck 1986). All computations were performed for a 15 m class telescope, with beam sizes of 10′′, 20′′ and 30′′ as grid input values. This grid spans 15 309 star models, for which each rotational line intensity was computed for three values of the beamwidth.

As was done for the observational line profiles, the theoretical line profiles were fitted with a soft parabola. In the case of the theoretical line profiles, however, an extra offset should be added to the equation, resulting in $\begin{array}{ccc}{\mathit{T}}^{\mathit{g}}\mathrm{\left(}\mathit{v}\mathrm{\right)}& \mathrm{=}& {\mathit{T}}_{\mathrm{MB}\mathit{,}\mathrm{c}}^{\mathit{g}}{\left[\mathrm{1}\mathrm{-}{\left(\frac{{\mathit{v}}^{\mathit{g}}\mathrm{-}{\mathit{v}}_{\mathrm{LSR}}^{\mathit{g}}}{{\mathit{v}}_{\mathrm{\infty }}^{\mathit{g}}}\right)}^{\mathrm{2}}\right]}^{\mathit{\beta }\mathit{/}\mathrm{2}}\mathrm{+}{\mathit{T}}_{\mathrm{0}}^{\mathit{g}}\\ & & \mathrm{for}\mathrm{|}{\mathit{v}}^{\mathit{g}}\mathrm{-}{\mathit{v}}_{\mathrm{LSR}}^{\mathit{g}}\mathrm{|}\mathrm{\le }{\mathit{v}}_{\mathrm{\infty }}^{\mathit{g}}\\ & \mathrm{=}& {\mathit{T}}_{\mathrm{0}}^{\mathit{g}} \mathrm{for}\mathrm{|}{\mathit{v}}^{\mathit{g}}\mathrm{-}{\mathit{v}}_{\mathrm{LSR}}^{\mathit{g}}\mathrm{|}\mathit{>}{\mathit{v}}_{\mathrm{\infty }}^{\mathit{g}}\end{array}$(7)where T is the line intensity⁴ in main-beam temperature. The superscript “g” designates that we are dealing with theoretical profiles that were computed using the GASTRoNOoM-code. The reason for the offset is that, in case of optically thin envelopes, one still detects the almost unattenuated stellar continuum (see Fig. 4). This offset parameter is not present in Eq. (1) that was used to fit the data. During the data-reduction process, the offset is removed from the spectra when baselines are subtracted, since the observer does not know if and how much stellar light penetrates the envelope. When confronting observational data with theoretical integrated line intensities, one therefore has to compare to I, where I is the integrated line intensity ^∫T(v)dv in K km s^-1. The quantity will from now on be referred to as I^g,cc, with the superscript “cc” denoting the continuum correction. Optically thin lines allow one to trace the wind acceleration zone of which the emission is visible as a superimposed small peak around the line centre (see Fig. 4).

The soft-parabola fit is not ideal to reproduce optically thin resolved emission lines (see upper panels of Fig. 4). The β-value was determined in a way that the soft parabola should yield a good representation of the line depression. The integrated intensity of the fit was not required to match the integrated intensity of the GASTRoNOoM-predictions. This method provides information about the degree to which the line profiles are resolved by the telescope beam.

The values in Table 3 show that the dependence of the integrated intensity on the distance is strongest. For the unsaturated regime, the integrated intensity is quite sensitive on the terminal velocity and mass-loss rate, while this dependence is, as expected, much lower in the saturated regime. The intensity is moderately dependent on β, the stellar temperature, stellar radius, and dust condensation radius.

3.4. Grid results

To estimate the minimisation function to the integrated intensities, the dependence of the integrated intensity I^g on the different input parameters was inspected, and a power-law relation between the integrated intensity and all parameters was assumed. To estimate the outcome variable I^g,cc we therefore have minimised the function $\begin{array}{ccc}{\mathit{I}}^{\mathrm{D}}& \mathrm{=}& \frac{\mathrm{1}}{{\mathit{s}}_{\mathrm{1}\mathit{,J}}}\mathit{Ṁ}{\mathit{s}}_{\mathrm{2}\mathit{,J}}{\mathit{\theta }}_{\mathrm{b}}^{{\mathit{s}}_{\mathrm{3}\mathit{,J}}}{\mathit{d}}^{{\mathit{s}}_{\mathrm{4}\mathit{,J}}}{\mathit{v}}_{\mathrm{\infty }}^{{\mathit{s}}_{\mathrm{5}\mathit{,J}}}{\mathit{f}}_{\mathrm{H}\mathit{,}\mathrm{CO}}^{{\mathit{s}}_{\mathrm{6}\mathit{,J}}}{\mathit{T}}_{\mathit{\star }}^{{\mathit{s}}_{\mathrm{7}\mathit{,J}}}{\mathit{R}}_{\mathit{\star }}^{{\mathit{s}}_{\mathrm{8}\mathit{,J}}}{\mathit{R}}_{\mathrm{inner}}^{{\mathit{s}}_{\mathrm{9}\mathit{,J}}}{\mathit{\beta }}^{{\mathit{s}}_{\mathrm{10}\mathit{,J}}}\mathit{.}\end{array}$(8)The minimisation was carried out for each ¹²CO line transition, using the generalized reduced gradient method (Lasdon et al. 1978). The correlation between the GASTRoNOoM output variable I^g,cc and the integrated intensity I^D derived from Eq. (8) is not bad, but a clear saturation effect is visible for all transitions when I^g,cc exceeds a certain boundary value (see Fig. 5). These boundary values are higher for higher J. A higher mass-loss rate not only leads to an increase of the cooling of the circumstellar gas, but also to a decrease in drift velocity between the dust and the gas (Decin et al. 2006). This leads to a weaker dependence of the integrated intensity of the CO lines on Ṁ, visible as saturation. Not accounting for this saturation effect typically leads to underestimating Ṁ by a factor of 3−10 for Ṁ ≥ 3 × 10^-5 M_⊙  yr^-1.

Moreover, it is visible from Fig. 5 that the boundary values for saturation are distance dependent, with lower boundary values for larger distances. This dependence was removed by minimising Eq. (8) again, but with both the left-hand and the right-hand side multiplied with a factor d². The coefficients s_i,J for J = 1−0, 2−1, 3−2, 4−3, 6−5, 7−6, 10−9, 16−15 and the boundary values of I_MB × d² for saturation are listed in Table 3.

Replacing I^D in Eq. (8) with the observed integrated line intensity I_MB yields an estimate for the mass-loss rate, Ṁ_D, provided that the other parameters are known: $\mathit{Ṁ}\mathrm{D}\mathrm{=}{\left[\frac{{\mathit{s}}_{\mathrm{1}\mathit{,J}}\mathrm{\times }{\mathit{I}}_{\mathrm{MB}}}{{\mathit{\theta }}_{\mathrm{b}}^{{\mathit{s}}_{\mathrm{3}\mathit{,J}}}{\mathit{d}}^{{\mathit{s}}_{\mathrm{4}\mathit{,J}}}{\mathit{v}}_{\mathrm{\infty }}^{{\mathit{s}}_{\mathrm{5}\mathit{,J}}}{\mathit{f}}_{\mathrm{H}\mathit{,}\mathrm{CO}}^{{\mathit{s}}_{\mathrm{6}\mathit{,J}}}{\mathit{T}}_{\mathit{\star }}^{{\mathit{s}}_{\mathrm{7}\mathit{,J}}}{\mathit{R}}_{\mathit{\star }}^{{\mathit{s}}_{\mathrm{8}\mathit{,J}}}{\mathit{R}}_{\mathrm{inner}}^{{\mathit{s}}_{\mathrm{9}\mathit{,J}}}{\mathit{\beta }}^{{\mathit{s}}_{\mathrm{10}\mathit{,J}}}}\right]}^{\frac{\mathrm{1}}{{\mathit{s}}_{\mathrm{2}\mathit{,J}}}}\mathit{.}$(9)In Fig. 6 we show a comparison between the mass-loss rate Ṁ_input used as input for the grid models, and the mass-loss rate Ṁ_D, calculated from the CO line intensities that result from the models, using Eq. (9). In Table 3 and Fig. 6 we listed the mean values δ, standard deviations σ, and maximum absolute values μ of the quantity log (Ṁ_D/Ṁ_input). All listed δ-values are positive, meaning that on average we slightly overestimate the input-mass-loss rate. This is however limited to overestimates of a few percents only, except for the cases of 1−0 (15%, δ = 0.06), and 2−1 in the saturated regime (17 %, δ=0.07). These deviations are still well within the standard deviation σ and therefore not significant. In the panels of Fig. 6 we also mention the quantities π(3) and π(10). These are the percentages of models in the grid that result in a value of Ṁ_D that deviates by more than a factor 3, resp. 10, from the input mass-loss rate. The σ-values indicate uncertainties of a factor  ~1.3 for the unsaturated regime, and a factor  ~1.9 for the saturated regime when comparing Ṁ_D to Ṁ_input. When applying Eq. (9) to data, the error bars on the measured parameters, e.g. T _⋆ , add to these uncertainties, increasing the typical uncertainty to a factor of three.

The spread on the estimates reflects the interplay between the different stellar parameters and the line intensities, and the large parameter space covered by the grid. In the saturated regime the spread is somewhat larger than in the unsaturated regime, with standard deviations σ ranging up to 0.16 and 0.31 in the unsaturated and saturated regimes, respectively. These values are still rather low, and are also reflected in the values of π(3) and π(10) that do not exceed 10% and 0.45% in the saturated regime, clearly showing the very small number of outliers present in our results. We point out that we found no obvious relations between the (few) grid points showing the larger values of  |log (Ṁ_D/Ṁ_input)| . Therefore, we do not attach larger intrinsic uncertainties to specific input parameters. The quantities σ,μ,δ,π indicate a high accuracy of the estimator for the complete range of mass-loss rates included in the grid (Table 2).

3.5. Discussion on the derived mass-loss rate formulae

3.5.1. Varying the number of free parameters

As was done in other studies (e.g. Ramstedt et al. 2008), one can lower the number of free parameters in Eq. (8). This is illustrated in Fig. 7 for the J = 2−1 transition. We show the estimates in two cases, omitting the dependence of the integrated intensity on β and R_inner (left panel) and on β, R_inner, T _⋆ , and R _⋆  (right panel). Mainly for the mass-loss rates exceeding 10^-5 M_⊙  yr^-1, the estimates exhibit larger uncertainties. An overview of the values for the standard deviations σ, and maximum absolute values μ of log (Ṁ_D/Ṁ_input) is given in Table 4 for the different cases. Both μ and σ increase significantly when lowering the number of parameters. This implies that it is better to include a well-founded estimate for some parameter, like β or R_inner, than to omit it from Eq. (8).

3.5.2. Sensitivity analysis

To further assess the quality of our estimator, we focus on the influence of (1) input parameters for Eq. (9), such as effective temperature, luminosity, dust condensation radius, photospheric CO abundance, and distance, (2) the envelope’s outer radius, and (3) the gas kinetic temperature.

Input parameters:

We tested the influence of d, T _⋆ , L _⋆  (R _⋆ ), f_H,CO and R_inner on the Ṁ-estimates for J = 1−0, 2−1, 3−2, 4−3, 6−5, 7−6. We assumed a standard stellar model with Ṁ_input  = 10^-7 M_⊙  yr^-1, T _⋆   = 2000 K, R _⋆   = 3.14 × 10¹³ cm, R_inner  = 3 R _⋆ , f_H,CO  = 10^-4, d = 300 pc, and a fixed beamwidth θ_b = 10 ′′. The integrated intensities I_MB are those calculated with the GASTRoNOoM-code for this model. In Fig. 8 we show the influence of the different parameters on the Ṁ-estimates by varying d, f_H,CO, T _⋆ , L _⋆  and R_inner one by one while keeping all others constant at the model values. We point out that a slight gradient is present in the Ṁ-estimates for the standard model, i.e. for higher J we seem to estimate somewhat lower values of Ṁ. The deviations and the spread of the estimates are, however, well within a factor three from the input mass-loss rate. Since the latter is used as a criterion to decide upon constancy of Ṁ, we can state that the estimator can very well reproduce the input-Ṁ.

As expected, the distance has the strongest influence on the estimated Ṁ-values, with a similar increase or decrease in mass-loss rate for all transitions. Altering f_H,CO has a pure scaling effect on the estimates, with a lower Ṁ for a higher f_H,CO. Changing T _⋆, L _⋆, or R_inner leads to the largest differences for the higher-J transitions, since changing these parameters affects the layers closest to the star more strongly than it does the more outward layers. This effect is reflected in the larger absolute values of exponents s_7,J,s_8,J and s_9,J in Eq. (9), and by the fact that the mentioned gradient is slightly shallower or steeper when considering other values for T _⋆, L _⋆, or R_inner. This implies that the estimates based on the higher-J transitions suffer from larger uncertainties and that more value can be attached to the estimates based on the lower-J transitions.

Inspecting the different panels in Fig. 8, we find no indications for changes in the Ṁ-estimates with an opposite character for the different transitions, in the sense that, e.g., the estimates for some J would go lower, while those for other transitions go higher. For a constant mass-loss rate we will therefore only see a gradient (if any) and not a random distribution of the estimates. This implies that a faulty parameter assumption will only lead to a general shift of the estimates and is unlikely to change the interpretation of possible variability or constancy of the mass-loss rate.

Outer envelope radius:

the influence of the outer radius of the CSE, R_outer, can not be tested in the same way since it is not an input parameter in Eq. (9). This variable has therefore no direct influence on the estimates. It did however influence the derivation of the coefficients s_i,J, since the geometrical extent of the envelope affects the integrated intensities and the shapes of the emission lines. The transition most sensitive to changes in R_outer is J = 1−0, implying that it is a safer choice to use Ṁ-estimates based on the J = 2−1 transition.

Clumping of CSE material causes a larger R_outer than a smooth mass outflow because of shielding properties of the clumps. This implies that higher integrated intensities can be reached with the same mass-loss rate. Therefore, if clumping is present, we could be overestimating the mass-loss rate. This will again be most strongly reflected in the J = 1−0 estimates.

Gas kinetic temperature:

in Fig. fig.9 we compare gas kinetic temperatures T, which have been consistently computed with the GASTRoNOoM-code, and integrated intensities I^g,cc of transitions J = 1−0 up to J = 7−6 for two different values of Ṁ_input (10^-7 M_⊙  yr^-1 and 10^-5 M_⊙  yr^-1) and five different effects in cooling. All other model parameters are the same as those mentioned above for the standard model. The default cooling scheme, which was used to derive the formalism presented in this paper, includes among others cooling due to rotational excitation of H₂O. The cooling schemes we considered (1) ignore the H₂O-cooling, (2) include H₂O-cooling (default), (3) include H₂O-cooling ten times stronger than assumed by default, (4) assume a power law for the gas kinetic temperature $\mathit{T}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathit{T}}_{\mathit{\star }}{\mathrm{(}{\mathit{R}}_{\mathit{\star }}\mathit{/}{\mathit{r}}^{\mathrm{)}}}^{\mathit{ϵ}}$(10)with ϵ = 0.5, and (5) assume a power law with ϵ = 0.7.

We point out that, except for the cases in which a power law was used to describe T(r), the temperature can increase or decrease with increasing radius r, i.e., instead of cooling there can be heating of the envelope; see Fig. fig.9.

We focus on the H₂O rotational cooling since this part of the cooling is still subject to large uncertainties. In GASTRoNOoM, H₂O is modelled as a three-level system, which is a gross simplification of the actual structure of the molecule. The rotational rate coefficients used here are based on the results of Green et al. (1993), as are those of Neufeld & Kaufman (1993). Faure & Josselin (2008) recently calculated new rotational rates and compared them to those of Neufeld & Kaufman (1993), finding higher cooling rates for temperatures  ≲ 1500 K and lower cooling rates for temperatures  ≳ 1500 K. They also found that the rotational rate coefficients can be uncertain by factors of a few up to an order of magnitude, which is reflected in the uncertainties on the cooling rates.

The contribution of H₂O-cooling is most important in the inner wind regions, since farther out in the envelope CO-rotational cooling and especially cooling due to adiabatic expansion dominate. This implies that the gas kinetic temperature in the inner wind regions is more uncertain than in the rest of the envelope.

As is visible in Fig. fig.9(b), I^g,cc decreases for increasing H₂O-cooling. Again, this effect is more explicit for higher J. Overestimating the H₂O-cooling therefore leads to an underestimate of the mass-loss rate, considering the relation between the integrated intensities and the Ṁ-estimates.

Since the excitation regions of the higher-J transitions are partly situated in the inner wind regions, the integrated intensities of these transitions will be subject to larger uncertainties than those of the low-J transitions J = 1−0 and 2−1. This different degree of influence on the integrated intensities for the different J-values is very clear in Fig. fig.9(b), where the spread on the integrated intensities of the transitions J → J−1 increases for higher J. A second implication is that the mentioned gradient of (slightly) lower Ṁ for higher J-values is possibly linked to the uncertainties on the H₂O-cooling in the inner regions. If this is indeed the case, then the H₂O-cooling rates used in our models are overestimates of the actual rates.

Considering the current models, the uncertainties on the intensities of higher-J transitions (as discussed above), and taking into account that the intensities of the J = 1−0 transition can be affected by e.g. masering and clumping, we put forward that the J = 2−1 transition is likely the most reliable transition to use in estimating mass-loss rates.

3.5.3. Comparison to Ramstedt et al. (2008)

In benchmarking our formalism against mass-loss rate estimators presented in the literature, we focus on the one derived by Ramstedt et al. (2008). A direct comparison to the estimators of Knapp & Morris (1985) or Loup et al. (1993) is not possible, since these have the main-beam temperature at the line centre, T_MB,c, as a fundamental parameter, while we consider the integrated intensity I_MB and the shape of the lines.

In Fig. 10 we show the mass-loss rate estimates derived for transitions J = 1−0, 2−1, 3−2 and 4−3 via Eq. (6), using the same parameter grid as was used to derive the formalism presented in this paper. Ramstedt et al. (2008) mention that their estimator Ṁ_R is only valid for mass-loss rates between 10^-7 and 10^-5 M_⊙  yr^-1, so we marked those outside this interval in grey, for the sake of the clarity. Ramstedt et al. (2008) also note the requirement of unresolved envelopes to use their estimator to high accuracy, and we therefore eliminated the grid points that result in resolved envelopes by determining the CO-photodissociation radius, R_CO, and subsequently the outer radius, R_outer, of the envelope according to Stanek et al. (1995) and Schöier & Olofsson (2001): $\begin{array}{ccc}& & {\mathit{R}}_{\mathrm{CO}}\mathrm{=}\mathrm{5.4}\mathrm{\times }{\mathrm{10}}^{\mathrm{16}}{\left(\frac{\mathit{Ṁ}}{{\mathrm{10}}^{-6}}\right)}^{\mathrm{0.65}}{\left(\frac{\mathrm{15}}{{\mathit{v}}_{\mathrm{\infty }}}\mathrm{\times }\frac{{\mathit{f}}_{\mathrm{H}\mathit{,}\mathrm{CO}}\mathrm{\left(}{\mathit{R}}_{\mathit{\star }}\mathrm{\right)}}{\mathrm{4}\mathrm{\times }{\mathrm{10}}^{-4}}\right)}^{\mathrm{0.55}}\\ & & \mathrm{+}\mathrm{7.5}\mathrm{\times }{\mathrm{10}}^{\mathrm{15}}\left(\frac{{\mathit{v}}_{\mathrm{\infty }}}{\mathrm{15}}\right)\mathrm{cm}\\ & & \mathit{\alpha }\mathrm{=}\mathrm{2.79}\mathrm{\times }{\left(\frac{\mathit{Ṁ}}{{\mathrm{10}}^{-6}}\mathrm{\times }\frac{\mathrm{15}}{{\mathit{v}}_{\mathrm{\infty }}}\right)}^{\mathrm{0.09}}\\ & & {\mathit{f}}_{\mathrm{H}\mathit{,}\mathrm{CO}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathit{f}}_{\mathrm{H}\mathit{,}\mathrm{CO}}\mathrm{\left(}{\mathit{R}}_{\mathit{\star }}\mathrm{\right)}\mathrm{\times }\exp \left[\mathrm{-}\ln \mathrm{\left(}\mathrm{2}\mathrm{\right)}\mathrm{\times }{\left(\frac{\mathit{r}}{{\mathit{R}}_{\mathrm{CO}}}\right)}^{\mathit{\alpha }}\right]\mathrm{·}\end{array}$The CO-photodissociation radius and the outer radius of the envelope are the radii at which the CO abundance has dropped to, respectively, 50% and 1% of its photospheric value, f_H,CO(R _⋆ ). α describes the abundance profile throughout the envelope. We adopted the Ṁ_input-value to calculate R_CO and R_outer. When the angular size of the envelope with radius R_outer exceeded the beam size θ_b, we excluded the respective grid point.

The mean values δ, standard deviations σ, and maximum absolute values μ of log (Ṁ_R/Ṁ_input) are printed in each panel of Fig. 10, both for the full grid and for only those grid points inside the interval  [ 10^-7;10^-5 ]  M_⊙  yr^-1. Very high values of μ are reached when the complete grid is considered. Confining the statistics to the interval  [ 10^-7;10^-5 ]  M_⊙  yr^-1, we find that these values are much lower. For all transitions μ and σ are higher than the values mentioned in Table 3, derived for our estimator. There is also a significant shift to higher Ṁ-values, for all but the J = 1−0 transition, using the estimator of Ramstedt et al. (2008) in Eq. (6). This shift is indicated by δ-values significantly deviating from zero (δ > σ) and is expected due to the differences in the temperature structure adopted by us and by Ramstedt et al. (2008). The uncertainties on the cooling, especially by rotational excitation of H₂O, are still rather large, especially so in the inner wind region, as discussed in Sect. 3.5.2. Furthermore, our approach utilises a higher number of parameters and expands the validity range of the Ṁ-estimator by correcting for saturation of the CO lines, ensuring that we get small intrinsic uncertainties for our estimator.

4. Stellar parameters

4.1. Temperature, luminosity and distance

Since we want to obtain empirical relations between wind parameters (e.g. Ṁ, v_∞), molecular-line parameters (e.g. line intensity), and basic stellar parameters (e.g. T _⋆ , L _⋆ ), the latter need to be well determined. To draw meaningful conclusions from sample statistics, a highly homogeneous determination of basic stellar parameters is needed. Effective temperatures were derived from the dereddened V−K colour, and luminosities from period-luminosity relations. A description of the methods used to obtain these parameters is given in Appendix C. The results for the sample are given in Table 6, together with the evolutionary stage, chemical type and pulsational type of the objects.

4.2. Inner envelope radius

As discussed in Sect. 3, the inner radius of the circumstellar envelope, R_inner, is used in determining the mass-loss rate. In estimating Ṁ for the sample targets, R_inner is calculated assuming that the dust temperature T_d at this radius is equal to the condensation temperature, and that ${\mathit{T}}_{\mathrm{d}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathit{T}}_{\mathit{\star }}{\left(\frac{{\mathit{R}}_{\mathit{\star }}}{\mathrm{2}\mathit{r}}\right)}^{\mathrm{2}\mathit{/}\mathrm{(}\mathrm{4}\mathrm{+}\mathit{s}\mathrm{)}}$(14)with s ≈ 1 (Olofsson in Habing & Olofsson 2003). The condensation temperature is taken to be 1500 K for O-rich and S-type stars and 1200 K for C-rich stars.

5. Results and discussion

5.1. Mass-loss rate

For all targets in the sample for which we could fit at least one rotational line profile with a soft parabola (Eq. (1)) we listed estimates for Ṁ in Table 7. The estimates from all ¹²CO lines for the sample AGB stars are shown in Fig. 11. For the supergiants and hypergiants they are shown in Fig. 12. In case the soft-parabola fit procedure yielded a negative β, the exponent s_10,J in Eq. (9) was assumed to be zero.

The listed estimates are averages over all values obtained via Eq. (9) for the fitted line transitions, and are written as . In Table 7 we also listed the number of lines used and the spread of the estimates around the mean.

We did not attempt to fit the line profiles for U Ant and S Sct, since these are composed of emission from the inner part of the CSE and from the detached shell at larger outflow velocities, respectively. As discussed in Sect. 2.3.2 some other targets show clear deviations from the fitted soft-parabola profile. The estimates for these targets, e.g. EP Aqr, are therefore subject to larger uncertainties than the spreads listed in Table 7 and have to be interpreted with caution.

If the mass-loss-rate estimates of all available lines are within or close to the typical factor three uncertainty region, it is reasonable to state that the mass-loss rate of the target has been constant throughout the regions of the envelope sampled by these transitions. For most AGB stars in the sample this is indeed the case, but there are some exceptions indicating a possible variability in mass loss in the part of the envelope that is traced.

5.1.1. Variability

Criterion:

Since higher-J transitions trace deeper, more recently formed envelope layers, a criterion to decide upon variability should consider the behaviour of the Ṁ-estimates in view of the J-values. The presence of a picket-fence structure in the estimates, i.e. with Ṁ varying strongly between transitions could suggest variability of the mass-loss rate, but on rather short time scales of some hundreds of years. The presence of a strong overall gradient can also imply mass-loss-rate variability: e.g. decreasing Ṁ for increasing J can be indicative for a mass-loss rate that has decreased over time. When discussing the sensitivity of our results to the modelling of the cooling of the envelope in Sect. 3.5.2, we pointed out that our estimates might intrinsically hold a shallow gradient. The differences associated with these model uncertainties between the highest and lowest estimates (for a constant Ṁ_input) are, however, not expected to exceed a factor three for J = 1−0 up to J = 7−6. Therefore, if a much steeper gradient is found in the estimates for some target, with a total coverage exceeding an order of magnitude in Ṁ, we are therefore inclined to state that the mass-loss rate of the respective target has varied over time.

Variability in the sample:

the estimates for V669 Cas, shown in Fig. 11, exhibit the picket-fence type of trend and suggest that variability in the mass-loss rate is traced for this star.

The estimates for NV Aur and GLMP 1048 based on data of the J = 2−1, 3−2, 4−3, 6−5 transitions, the estimates for R Scl and IRC +20370, for J = 1−0, 2−1, 3−2, 4−3, 7−6, and the estimates for CW Leo, for J = 1−0, 2−1, 3−2, 4−3, 6−5 exhibit a spread of about an order of magnitude. This spread is within the uncertainties we expect around a Ṁ-value and therefore points to a constant mass-loss rate throughout the traced regions of the envelope. The estimates for RW LMi, LP And and χ Cyg with J = 1−0, 2−1, 3−2, 4−3, 6−5, respectively cover factors about ten, twenty and thirty in Ṁ. The estimates for RX Boo, for J = 1−0, 2−1, 3−2, 6−5, cover a factor about forty in Ṁ. In all these cases, we find a strong downward gradient, which most likely indicates that the more inward regions, i.e. the excitation regions of the higher-J lines, were produced by a lower mass-loss rate than the outer regions. A strong gradient is also found for both WX Psc and LL Peg for J = 1−0 up to J = 6−5. For both targets, however, the Ṁ-estimate based on the J = 7−6 transition is higher by about an order of magnitude than the one based on the J = 6−5 transition, indicating a recent increase of the mass-loss rate, as was proposed by Decin et al. (2007) for WX Psc. Qualitatively, this phenomenon is independent of the modelling of the H₂O-cooling, i.e. an increase would still be present if another cooling mechanism was considered in deriving the estimates. Obtaining a result corresponding to detailed modelling using an estimator as in Eq. (9) can be used as a benchmark for the quality of the estimator.

It is likely that the OH/IR stars and the more extreme Mira stars in the sample have already undergone mass-loss modulations, but to verify this more data are needed, especially of J = 1−0, 6−5 and 7−6.

Data quality:

the IRAM measurements⁵ in our sample are subject to somewhat larger uncertainties than our own APEX and JCMT data since some questions have arisen about the data reduction. Decin et al. (2008a) already reported large uncertainties on the IRAM data of R Hya and Skinner et al. (1999) mentioned significant deviations between different measurements of the same lines towards CW Leo. The IRAM data Skinner et al. gathered from the literature and archives seemed to suffer from large uncertainties linked to pointing and calibration. They also discussed inconsistencies in the literature concerning the reported absolute line intensities. Very often confusion or unclarity arises about e.g. conversion from antenna temperatures to main-beam-brightness temperatures. For these reasons, we considered leaving the IRAM measurements out of the data set. In case of LL Peg this leads to a decrease of the estimated mass-loss rate with a factor 5.3, for IRC +20370 with a factor 3.9. For IRC +10420 and R Hya omitting the IRAM data has no effect on the Ṁ-estimate. For the eight other objects with IRAM observations, the decrease is with a factor between 1.3 and 2.7. The Ṁ-estimates ignoring IRAM data are always lower.

Omitting the IRAM data has an effect on the gradients seen in the estimates and discussed earlier. In many cases, the CO-sampling is now restricted to the transitions higher than J = 1−0 and/or J = 2−1, making possible trends less significant. In case of CW Leo, RX Boo, IRC +20370, RW LMi and LP And the supposed trends in the estimates are no longer present. For LL Peg the presence of an OSO J = 1−0 measurement still gives the trend seen when the IRAM data are included.

For the oxygen-rich semi-regular variable RX Boo we find two estimates that could hint towards variability, i.e. those based on the IRAM J = 1−0 and JCMT J = 6−5 measurements. The latter is duplicated by a measurement carried out with CSO, and the estimate based on that observation is about an order of magnitude higher than the one derived from the JCMT. We are however inclined to attach more value to our JCMT observation, since the CSO calibration as described by Teyssier et al. (2006) holds larger uncertainties, such as rather large pointing errors. We conclude from this that there is some, however somewhat uncertain, downward trend with lower Ṁ-estimates for higher J, indicating the possibility of a decrease in mass-loss rate over time.

Comparison to the literature:

Teyssier et al. (2006) note that a change in the mass-loss rate is needed to reproduce the data for three of their sample targets: R Hya, RW LMi and χ Cyg. Decin et al. (2008a) constructed a mass-loss model for R Hya and showed that a detached shell is needed to explain the rotation-vibrational lines. The purely rotational lines could however be reproduced with a constant mass-loss rate. Our estimates do, indeed, not point to any significant variations in the mass-loss rate, as can be seen from Fig. 11. The Ṁ-estimates for RW LMi cover one order of magnitude and show a (shallow) gradient with lower Ṁ for higher J. This could point to a mass-loss rate that is presently lower than the one that produced the outer layers of the circumstellar envelope. Teyssier et al. specify a difference of only a factor 1.5 between the two mass-loss rates, so we can say that our estimates support this hypothesis of mass-loss-rate variability. In case of χ Cyg Teyssier et al. mention a minimum change of a factor 4 between the former and the present mass-loss rate, with the present one being lowest. This result is again reflected in our estimates.

The detached-shell structure of the carbon-rich semi-regular variable R Scl was described by Bergeat & Chevallier (2005) as the consequence of a former mass-loss rate of  ~5.5 × 10^-6 M_⊙ yr^-1 and a present mass-loss rate of  ~3.5 × 10^-7 M_⊙ yr^-1. The range of mass-loss rates covered by our estimates agrees with these values.

The mass-loss rates estimated from the individual CO emission lines for WX Psc – shown in Fig. 11 – are in very good agreement with the detailed modelling of the star by Decin et al. (2007). They present a profile Ṁ(r) exhibiting a high Ṁ-value, followed in time by a strong decrease and again an increase of the mass-loss rate. We do not reproduce the same large difference of three orders in magnitude as mentioned by Decin et al. (2007). A possible explanation is that the excitation regions of the individual observed transitions are not restricted to either the high-Ṁ or the low-Ṁ regions. The estimates then represent more than one of these regions and are a kind of “mean” value of Ṁ. However, the estimates we present in this paper for WX Psc reflect the detailed mass-loss modelling very well.

Winters et al. (2000) reported on possible time variability of the mass-loss rate of CW Leo. The present-day mass-loss rate is suggested to be lower than the mass-loss rate that produced the outer envelope layers by a factor of up to five. The decrease we see in our estimates, however, is by a factor  ~30. Also, the Ṁ-values presented in the literature are about an order of magnitude higher than the -value we derived. A possible explanation for these discrepancies is the H₂O-cooling uncertainty, leading to underestimates of Ṁ. Another explanation is that the J = 7−6 transition for the first time reveals a stronger variation in Ṁ than could be traced with single-dish observations of CO transitions up to J = 4−3 or 6−5. The model presented by Ramstedt et al. (2008) fitting lines J = 1−0 up to J = 4−3 shows clear deviations from the data. The discrepancies are systematic, i.e. the model underestimates the observed intensities more strongly for increasing J. This indicates the need for consideration of Ṁ-variability for this target.

Neri et al. (1998) suggested that LL Peg might have a bipolar detached shell and Mauron & Huggins (2006) reported on a spiral structure in the envelope. The latter is indicative for variability in the outflow from the central star.

Supergiants and hypergiants:

the estimates for all five RSGs in our sample suggest constant mass-loss rates; see Fig. 12. In case of the yellow hypergiant IRC +10420, Humphreys et al. (1997); Castro-Carrizo et al. (2007) and Dinh-V.-Trung et al. (2009), among others, have put forward the possibility of a complex mass-loss history, with several episodes of high and low mass-loss rates. This resulted in multiple arc-like structures and shells around the central star. For AFGL 2343, the second hypergiant in our sample, we only used the J = 3−2 line of ¹²CO to estimate the mass-loss rate, so no conclusions on variability in mass-loss rate can be drawn from this. Castro-Carrizo et al. (2007) mention variations in the mass-loss rate for both hypergiants on time scales of  ~1000 years. In this study, where we only consider CO emission lines measured with single-dish telescopes, we lack sensitivity for the rapid density fluctuations. This is due to the overlap of the line forming regions of the different rotational transitions.

5.1.2. AGB mass loss in view of pulsational type and chemistry type

The histograms in Fig. 13 show that the spread on for AGB stars is quite large, but there is a noticeable difference between the semi-regular variables on the one hand and the Miras and OH/IR targets on the other – see Fig. 13(a). The spread on Ṁ for the 16 Miras is rather extended, with values ranging from 9 × 10^-8 M_⊙  yr^-1 (R Leo) up to 6 × 10^-5 M_⊙  yr^-1 (NX Ser) and a mean value of 10^-5 M_⊙  yr^-1.

The mass-loss rates for the 9 OH/IR targets were all estimated in the range  M_⊙  yr^-1, with an average of 2 × 10^-5 M_⊙  yr^-1. The estimates for these targets seem low, considering that OH/IR stars have optically thick envelopes and can be heavily obscured. However, it has already been pointed out by e.g. Delfosse et al. (1997) that the mass-loss rates derived from CO lines can be significantly lower than those derived from the infrared flux in the case of OH/IR stars. The most likely explanation for the discrepancy is found in the recent (~1000 yr ago) onset of a superwind phase. The latter would not be traced by the CO lines available in this data set, except possibly for WX Psc.

The SRb targets have lower mass-loss rates, with an average of 5 × 10^-7 M_⊙  yr^-1 and the highest value reaching 2 × 10^-6 M_⊙ yr^-1 (R Scl). These results are in accordance with the idea that the SR variables are likely progenitors of Miras (Whitelock & Feast 2000; Yeşilyaprak & Aslan 2004).

When considering the histogram representing the AGB targets grouped per chemistry type (oxygen-rich, carbon-rich or S-type) in Fig. 13b, no clear distinction between these groups can be made. Though this could be due to the bias of our sample towards oxygen-rich stars, there are no clear indications for significant discrepancies in Ṁ between the different types. Ramstedt et al. (2006) also reported no noticeable differences in Ṁ-values for the different chemistries. Putting all AGB targets together for which Ṁ-estimates could be made – irrespective of pulsational or chemical type – in Fig. 13c, we get a minimum value of 4 × 10^-8 M_⊙  yr^-1 (SRb V438 Oph), a maximum value of 6 × 10^-5 M_⊙  yr^-1 (O-rich Mira NX Ser) and a mean of 10^-5 M_⊙ yr^-1.

5.1.3. AGB mass loss as a function of pulsation period

The -values for all AGB stars in the sample – except for the RV Tauri star R Sct – are plotted versus pulsation period in Fig. 14 with error bars indicating the spread as listed in Table 7. In case only one or two estimates could be obtained, we assumed a spread of a factor three, which is appropriate considering the overall spread in the sample. The red data points in the graph were taken from He & Chen (2001) and cover different types of OH/IR stars with a large spread on the periods of pulsation. Five objects in the sample of He & Chen (2001) are also present in our sample. Their pulsation periods and mass-loss rates are in accordance with our estimates, except in the case of V669 Cas, where He & Chen use P = 1995 days (Slootmaker et al. 1985), while we adopted the more recently determined P = 1525 days (Groenewegen et al. 1999). In constructing a linear relation between the pulsation period and the logarithm of the mass-loss rate, we only took into account the data points in our sample – i.e. we did not consider the data points of He & Chen (2001) – with periods shorter than 1000 days. The positions in the graph of the objects with periods longer than 850 days are represented by a constant Ṁ of 3.4 × 10^-5 M_⊙  yr^-1. The linearly increasing part of the fit (P ≲ 850 days) was constructed to represent the Miras, the short-period OH/IR stars and the AGB semi-regulars. The constant fitting value represents a levelling-off or saturation of the AGB mass-loss rate for stars with long periods, as proposed by, e.g., Vassiliadis & Wood (1993). The lack of data points for stars with longer periods complicates the determination of the cut-off period and the mean Ṁ-value for the long-P objects. Both were more or less arbitrarily chosen to yield a good fit. The fit in Fig. 14 is represented by $\begin{array}{ccc}\begin{array}{c}\\ & & & & \\ & & & & \end{array}& & \end{array}$with Ṁ in units of M_⊙  yr^-1 and P in days. Although He & Chen (2001) note that the quality of their derived mass-loss rates may be questionable, we find that their data points of both short and long-period targets are well represented by our fit. The targets that deviate strongest from the fit are semi-regular variables, e.g. W Hya and U Hya, and the Miras NX Ser (O-rich) and W Aql (S-type). We note that in constructing the fit we did not take into account the effect of differences in stellar mass that are likely to cause extra spread on the mass-loss rate for fixed pulsation periods, as shown by Vassiliadis & Wood (1993). We also point out the very good agreement between our fit and the results presented by Schöier & Olofsson (2001), who produced a similar diagram of Ṁ versus pulsation period for a sample of 68 optically bright carbon stars (their Fig. 12), and found a relation Ṁ ∝ P^2.3 for periods P between 100 and 800 days.

5.2. ¹²C/¹³C isotope abundance ratio

For 29 sample stars we have observations of both ¹²CO and ¹³CO rotational lines, allowing us to estimate the ¹²C/¹³C isotope abundance ratio. This ratio is a tracer of the star’s chemical evolutionary state and nucleosynthesis (Sopka et al. 1989; Schöier & Olofsson 2000; García-Hernández et al. 2007; Milam et al. 2009). ¹³C is produced during the interpulse phases when, after dredge-up, protons are partially mixed in the region between the hydrogen and helium burning shells. Reactions of the protons with ¹²C in this intershell region form ¹³C. In massive AGB stars (3 M_⊙   ≤  M _⋆    ≤  7 M_⊙), “hot bottom burning” (HBB) can be activated when the convective envelope penetrates the hydrogen-burning shell. During HBB, ¹²C/¹³C is lowered, since the net production of ¹³C is higher than that of ¹²C in the re-activated CN-cycle. For HBB-stars, a ratio ¹²C/¹³C around 10 is expected, but could possibly range down to values about 3 or even 1 (Ohnaka & Tsuji 1999). HBB is expected to be active in the most massive and luminous oxygen-rich AGB stars. In this respect, OH/IR stars – typified by high luminosities, large infrared excesses and high mass-loss rates – are the best candidates for having very low ¹²C/¹³C isotope ratios (García-Hernández et al. 2007, and references therein).

Observational ¹²C/¹³C-values for AGB stars range from as low as 1 (Ohnaka & Tsuji 1999) up to about 90 (Milam et al. 2009). Lebzelter et al. (2008, and references therein) mention some clustering of empirical values of ¹²C/¹³C between 50 and 70 for galactic C-rich stars, while their evolutionary models even predict values higher than 100. The solar value is about 89 (Anders & Grevesse 1989).

Line strength ratios of ¹²CO and ¹³CO transitions with identical J provide a first order estimate of the isotope abundance ratio in case both lines are optically thin (Schöier & Olofsson 2000). In case of optically thick lines, these estimates are only tentative, and lower than the actual isotopic ratio. Schöier & Olofsson (2000) mention that the ¹²C/¹³C values derived from their radiative transfer models are mostly higher than those estimated from intensity ratios. This effect is largest for their high Ṁ-objects, where the difference amounts to a factor 6.

We note here that the line intensity ratios for J = 2−1 are consistently larger than those for J = 3−2 for our sample stars. A linear least-squares fit for the targets with line intensity ratios available for both J = 2−1 and J = 3−2 shows that

$\frac{\mathit{I}{\mathrm{(}}^{\mathrm{12}}\mathrm{CO}\mathrm{(}\mathrm{3}\mathrm{-}\mathrm{2}\mathrm{)}\mathrm{)}}{\mathit{I}{\mathrm{(}}^{\mathrm{13}}\mathrm{CO}\mathrm{(}\mathrm{3}\mathrm{-}\mathrm{2}\mathrm{)}\mathrm{)}}\mathrm{=}\mathrm{4.3}\mathrm{+}\mathrm{0.48}\mathrm{\times }\frac{\mathit{I}{\mathrm{(}}^{\mathrm{12}}\mathrm{CO}\mathrm{(}\mathrm{2}\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{)}}{\mathit{I}{\mathrm{(}}^{\mathrm{13}}\mathrm{CO}\mathrm{(}\mathrm{2}\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{)}}\mathrm{·}$

The standard deviation from this linear fit is 2.1 units in I(¹²CO(3−2))/I(¹³CO(3−2)). The slope and offset in this relation clearly indicate that individual line intensity ratios will not lead to accurate values of the actual ¹²C/¹³C abundance ratio and can only give rough estimates.

We have calculated ¹²C/¹³C-values for our sample targets, correcting for line strength differences and – for the sake of simplicity – assuming the targets were unresolved by the beam, in the way this was done by Schöier & Olofsson (2000) and as is given by Eq. (15). $\mathrm{12}\mathrm{C}{\mathit{/}}^{\mathrm{13}}\mathrm{C}\mathrm{=}\frac{{\mathit{I}}_{\mathrm{MB}}{\mathrm{(}}^{\mathrm{12}}\mathrm{CO}\mathrm{(}\mathit{J}\mathrm{\to }\mathit{J}\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{)}}{{\mathit{I}}_{\mathrm{MB}}{\mathrm{(}}^{\mathrm{13}}\mathrm{CO}\mathrm{(}\mathit{J}\mathrm{\to }\mathit{J}\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{)}}\mathrm{\times }{\left(\frac{{\mathit{\nu }}_{\mathrm{12}\mathrm{CO}\mathrm{(}\mathit{J}\mathrm{\to }\mathit{J}\mathrm{-}\mathrm{1}\mathrm{)}}}{{\mathit{\nu }}_{\mathrm{13}\mathrm{CO}\mathrm{(}\mathit{J}\mathrm{\to }\mathit{J}\mathrm{-}\mathrm{1}\mathrm{)}}}\right)}^{-3}\mathrm{·}$(15)The frequency correction factor is equal to 0.87 for all values of J considered. If multiple line intensity ratios were available for one J, these were averaged. The line intensity ratios and the resulting values for ¹²C/¹³C are listed in Table 8. Figure 15 shows a histogram plot of the ¹²C/¹³C-values for the AGB stars in the sample. It is clear from this graph that we derive very low values for ¹²C/¹³C, but as mentioned before, these estimates are lower limits in case of large mass-loss rate objects due to optical depth effects, that are not taken into account using this approach. Since high-Ṁ objects are more easily observed, the sample is biased towards the higher mass-loss rates. The values for the 12 Miras range from 3.6 (carbon-rich LL Peg) up to 14.4 (oxygen-rich R Cas) and for the 3 OH/IR stars they range from 4.9 (OH 44.8-2.3) up to 7.2 (AFGL 5379). The spread on ¹²C/¹³C for the other stars in the sample is quite large and we do not see any clear differences in the first-order estimates of ¹²C/¹³C when considering the different types of stars in the sample.

Milam et al. (2009) found ¹²C/¹³C   ~ 13 for NML Cyg and ¹²C/¹³C   ~ 14 for IRC +10420, while for VY CMa they derived a much higher value between 25 and 46. We find similar ratios with a first order estimate for VY CMa of 14.8, and values of 8.7 and 8.1 for NML Cyg and IRC +10420, respectively. We note that Milam et al. (2009) deconvolved the line profiles for VY CMa and V Hya into separate velocity components, since the different structures probed could have discrepant chemical properties (see also Ziurys et al. 2007, for VY CMa). They analysed the multiple components for ¹²C/¹³C separately, resulting in three different values for VY CMa (46, 33 and 25) and 2 values for V Hya (71 and 72). Milam et al. (2009) claim a  ± 30% uncertainty on their ¹²C/¹³C-values.

5.3. Wind driving efficiency

Olofsson et al. (1993) and Oudmaijer et al. (2009) suggest that the acceleration of winds around C-rich stars is more efficient than that of winds around O-rich stars because of differences in dust absorption properties. Comparing samples of 40 S-type, 77 oxygen-rich and 61 carbon-rich stars, Ramstedt et al. (2006) also find indications for a higher terminal velocity in carbon-rich stars.

Figure 16 shows histograms of v_∞ in our sample for the three different chemical types of objects: oxygen-rich, carbon-rich and S-type stars. The mean v_∞-values for respectively 24 oxygen-rich and 13 carbon-rich AGB targets in our sample are 14.5 km s^-1 and 15.4 km s^-1. These values are  ~4–10 km s^-1 higher than those of Ramstedt et al. (2006). We note that the largest v_∞-values for our carbon-rich and S-type sample exceed those of Ramstedt et al. (2006) by about 10 km s^-1. These extreme values are however the only ones in both samples – V Hya in the C-rich and π¹ Gru in the S-type sample – so they are only weakly reflected in the samples’ mean values. Our results do not lead to the same clear conclusion that the wind driving efficiency is higher for carbon-rich than for oxygen-rich stars. They do, however, point in the same direction considering the sampling bias we have in favour of oxygen-rich objects and the occurrence of higher v_∞-values for carbon-rich than for oxygen-rich targets in our sample.

6. Conclusions

The data set presented in this paper covers a very large number of high-quality measurements of CO rotational line transitions J → J−1 with J ranging from 1 up to 7. Observations of both ¹²CO and ¹³CO are included. It is the first time that a CO data set allowing one for such a detailed sampling of circumstellar envelopes of evolved objects is presented in the literature.

A grid of radiative transfer models calculated with GASTRoNOoM covering AGB stellar and wind parameters was used to derive a formalism allowing one to make estimates of mass-loss rates. The Ṁ-estimator is valid for J up to 16, while until now, only estimators valid for J up to 4 have been presented in the literature. We also extended the validity range of the Ṁ-estimator by correcting for saturation of the CO line intensities and therefore accounting for optically thick winds. The intrinsic spread of the estimator ensures good Ṁ-estimates when all input parameters are well determined. We like to point out that, when only one transition is used to estimate a mass-loss rate, the J = 2−1 transition will give the most reliable results.

For many of the sample AGB stars, the estimates based on the different rotational lines result in constant Ṁ-values. The estimates for 7 targets do however hint at a variability of the mass-loss rate, by showing a large spread (more than an order of magnitude in Ṁ) and clear trends. Uncertainties linked to the implemented H₂O-cooling could produce the effect of seemingly decreasing Ṁ in time, but there are two convincing cases in which these uncertainties cannot account for the variations seen in the estimates. This is the case for LL Peg and WX Psc, and we have shown that the Ṁ-estimates for the latter target reflect the mass-loss history constructed through detailed radiative-transfer modelling by Decin et al. (2007). The fact that we (1) do not have full sampling of all rotational J-levels⁶ for all sample targets and that (2) did not perform full radiative-transfer modelling for the sample, of course limit our ability to prove conclusively that a star has experienced a constant mass loss or mass-loss modulations. This kind of detailed modelling will be presented in a forthcoming paper by Lombaert et al. (in preparation). The examples discussed in this paper show, however, that the use of the estimator allows one not only to establish good mass-loss-rate estimates in case of constant mass-loss targets, but also to trace strong variability in the mass-loss rate of AGB stars, given a sufficiently large number of CO line transitions. The fact that our estimates bear quite close resemblance to the mass loss derived through detailed radiative transfer modelling, gives confidence that our estimator is of good quality.

Further observations of mostly high-excitation rotational lines, such as J = 6−5 and J = 7−6 would make it possible to make a better distinction between constant or (possibly) variable mass loss. The J = 5−4 transition will be measured by HIFIfor many stars in our sample. These additional measurements will provide a more complete sampling of the circumstellar envelopes and the mass-loss history of these evolved stars.

Online material

Appendix A: CO data – figures and tables

All ¹²CO and ¹³CO data of the sample and the soft-parabola fits (Sect. 2.3.1) to these data are shown in Figs. A.1 and A.2. The most important line parameters – main-beam temperature at the line centre and integrated intensity – and the β parameter from the soft-parabola fitting procedure as described in Sect. 2.3.1 are given in Tables A.1 and A.3. References for, or remarks on, the data are listed in Tables A.2 and A.4 for the ¹²CO and ¹³CO data respectively.

Appendix B: Dereddening

Since the effective temperatures of late-type stars are well-determined by the intrinsic V−K colour (Sect. C.1), we need a way to derive this colour accurately. The magnitudes of both the V and K bands – at 0.550 μm and 2.2 μm, respectively – are reddened due to interstellar and circumstellar contributions. The data were corrected for interstellar extinction using the dust maps of Schlegel et al. (1998), which give the interstellar reddening at infinity, E_B−V,∞, for a given set of galactic latitudes b and longitudes l. This value was transformed into interstellar reddening at a given distance d using ${\mathit{E}}_{\mathit{B}\mathrm{-}\mathit{V}}\mathrm{=}{\mathit{E}}_{\mathit{B}\mathrm{-}\mathit{V,}\mathrm{\infty }}\left[\mathrm{1}\mathrm{-}\exp \left(\mathrm{-}\mathrm{10}\frac{\mathit{d}}{\mathrm{kpc}}\sin \mathrm{\left|}\mathit{b}\mathrm{\right|}\right)\right]\mathit{,}$(B.1)following Feast et al. (1990). The resulting value for E_B−V was combined with R_V = 3.1 and A_K = 0.112A_V (Schlegel et al. 1998) to give interstellar extinction coefficients for the V and K bands.

Circumstellar extinction was not accounted for, even though it is assumed to be significant for high mass-loss rate LPV targets. Knapp et al. (1998) use the mass-loss rate Ṁ to determine the circumstellar extinction in the K band. Since the aim of the paper is to establish empirical relations between stellar parameters and Ṁ, this method can not be applied to the data for the sake of consistency. Knapp et al. (2003) suggest a formalism to calculate the circumstellar extinction in the K band, A_K,CS, as a function of S₁₂/S_2.2, the ratio of flux densities at 12 μm and 2.2 μm, respectively. This method is only valid for low optical depths (Knapp et al. 2003). When applied to the sample, unrealistic values for A_K,CS as high as several tens (LP And) or hundreds (CW Leo) of magnitudes are indeed obtained for the objects with presumed high Ṁ values. To decide, however, whether or not this correction should be applied for a given target seems arbitrary. For these reasons we did not correct for circumstellar reddening.

Appendix C: Basic stellar parameters

In this appendix we describe the methods that we used to obtain the effective temperature, stellar luminosity and distance towards objects. The results are given in Table 6 in the paper.

C.1. Effective temperature

Since broad-band near-infrared colours are very good temperature criteria, being relatively independent of stellar gravities and abundances, they are often used to determine effective temperatures. The V−K colour is a good T_eff-indicator for cool, late-type stars, e.g. GKM giants and supergiants (di Benedetto 1993; Bessell et al. 1998; Bergeat et al. 2001; Levesque et al. 2005). For hotter A through K stars the V−I colour is more suitable (Bessell et al. 1998).

In all T_eff versus V−K relations the Johnson-Cousins-Glass photometric system is used. K magnitudes for the sample stars that were obtained in the systems of 2MASS (K_s) and SAAO (K_SAAO) were transformed via $\mathit{K}\mathrm{=}{\mathit{K}}_{\mathrm{s}}\mathrm{+}\mathrm{0.044}\mathrm{=}{\mathit{K}}_{\mathrm{SAAO}}\mathrm{-}\mathrm{0.0275}\mathit{,}$(C.1)following (Bessell 2005). Since the accuracy of the K_SAAO-measurements (~0.03 mag; Olivier et al. 2001) is almost a factor 10 better than that of K_s-data (~0.2 mag), preference was given to the SAAO-measurements when both were available. The errors on the magnitudes are expected to be small relative to the amplitude in the K band caused by pulsations of the stars. The error bars on the T_eff values determined here, were calculated assuming an error on (V−K)₀ of 1 mag, leading to rather conservative temperature estimates.

In all equations listed (V−K)₀ is the intrinsic V−K broad-band colour with (1) both the V and K magnitudes at maximum light, or (2) both at their mean value if maximum light was not available, and (3) corrected for interstellar reddening in the way described in Sect. B.

C.1.1. Long-period variables: Miras, semi-regulars and irregular variables

Bessell et al. (1998) defined a T_eff versus V−K calibration for KM giants (luminosity class III), given by $\begin{array}{ccc}{\mathit{T}}_{\mathrm{eff}}& \mathrm{=}& \mathrm{9037.597}\mathrm{-}\mathrm{3101.282}\mathrm{\times }\mathrm{(}\mathit{V}\mathrm{-}\mathit{K}{\mathrm{)}}_{\mathrm{0}}\\ & & \mathrm{+}\mathrm{717.7044}\mathrm{\times }\mathrm{(}\mathit{V}\mathrm{-}\mathit{K}{\mathrm{)}}_{\mathrm{0}}^{\mathrm{2}}\\ & & \mathrm{-}\mathrm{85.83809}\mathrm{\times }\mathrm{(}\mathit{V}\mathrm{-}\mathit{K}{\mathrm{)}}_{\mathrm{0}}^{\mathrm{3}}\\ & & \mathrm{+}\mathrm{5.021194}\mathrm{\times }\mathrm{(}\mathit{V}\mathrm{-}\mathit{K}{\mathrm{)}}_{\mathrm{0}}^{\mathrm{4}}\\ & & \mathrm{-}\mathrm{0.1137841}\mathrm{\times }\mathrm{(}\mathit{V}\mathrm{-}\mathit{K}{\mathrm{)}}_{\mathrm{0}}^{\mathrm{5}}\mathit{.}\end{array}$(C.2)This relation is valid for 1.5 mag ≤ (V−K)₀ ≤ 10 mag (see Fig. C.1). However, Bessell et al. (1998) mention that their static model atmospheres cannot adequately represent the M giants cooler than 3000 K, since most of these are long-period variables (LPVs) with shock waves traversing and greatly extending their atmospheres.

di Benedetto (1993) presented a temperature calibration in which T_eff is derived from angular-diameter measurements based on Michelson interferometry. In case of variable stars, measurements at maximum light were used. The effective temperature of GK stars of luminosity classes II-III-IV-V is described by $\begin{array}{ccc}\log {\mathit{T}}_{\mathrm{eff}}& \mathrm{=}& \mathit{a}\mathrm{-}\mathrm{(}\mathit{b}\mathrm{+}\mathrm{0.0118}\mathrm{)}\mathrm{(}\mathit{V}\mathrm{-}\mathit{K}{\mathrm{)}}_{\mathrm{0}}\\ & & \mathrm{+}\mathrm{0.0118}\mathrm{(}\mathit{V}\mathrm{-}\mathit{K}{\mathrm{)}}_{\mathrm{0}}^{\mathrm{2}}\end{array}$(C.3)for 1.4 mag ≤ (V−K)₀ ≤ 3.7 mag, for FGKM supergiants by $\log {\mathit{T}}_{\mathrm{eff}}\mathrm{=}\mathrm{(}\mathit{a}\mathrm{-}\mathrm{0.205}\mathrm{)}\mathrm{-}\mathrm{(}\mathit{b}\mathrm{-}\mathrm{0.086}\mathrm{)}\mathrm{(}\mathit{V}\mathrm{-}\mathit{K}{\mathrm{)}}_{\mathrm{0}}$(C.4)for (V−K)₀ ≥ 0.6 mag and for M-type giants of type IV, Miras and carbon stars by Eq. (C.4) for 3.7 mag ≤ (V−K)₀ ≤ 10 mag (see Fig. C.1). The constants a and b in Eqs. C.3 and C.4, and the validity ranges in terms of (V−K)₀ are listed in Table C.2.

For oxygen-rich LPVs, except for OH/IR stars, the effective temperature was derived via Eq. (C.4) (di Benedetto 1993). The uncertainty on the temperature would then be of the order of 1 % according to di Benedetto (1993).

Bergeat et al. (2001) focussed on carbon-rich objects and derived a relation for T_eff versus V−K given by Eqs. C.5, C.6 and C.7. The notation T_eff,d denotes that the relation is based on angular-diameter measurements, T_eff is the final temperature adopted for their sample stars. $\begin{array}{ccc}\log {\mathit{T}}_{\mathrm{eff}\mathit{,}\mathrm{d}}& \mathrm{=}& \mathrm{-}{\mathrm{0.079}}_{\mathrm{(}\mathrm{\pm }\mathrm{0.008}\mathrm{)}}\mathrm{(}\mathit{V}\mathrm{-}\mathit{K}{\mathrm{)}}_{\mathrm{0}}\mathrm{+}{\mathrm{3.91}}_{\mathrm{(}\mathrm{\pm }\mathrm{0.02}\mathrm{)}}\\ & & \mathrm{for}\mathrm{4.99}\mathrm{\le }\mathrm{(}\mathit{V}\mathrm{-}\mathit{K}{\mathrm{)}}_{\mathrm{0}}\mathrm{\le }\mathrm{7}\\ \log {\mathit{T}}_{\mathrm{eff}\mathit{,}\mathrm{d}}& \mathrm{=}& \mathrm{-}{\mathrm{0.061}}_{\mathrm{(}\mathrm{\pm }\mathrm{0.008}\mathrm{)}}\mathrm{(}\mathit{V}\mathrm{-}\mathit{K}{\mathrm{)}}_{\mathrm{0}}\mathrm{+}{\mathrm{3.79}}_{\mathrm{(}\mathrm{\pm }\mathrm{0.02}\mathrm{)}}\\ & & \mathrm{for}\mathrm{7}\mathrm{\le }\mathrm{(}\mathit{V}\mathrm{-}\mathit{K}{\mathrm{)}}_{\mathrm{0}}\mathrm{\le }\mathrm{8.31}\\ \log {\mathit{T}}_{\mathrm{eff}\mathit{,}\mathrm{d}}& \mathrm{=}& \mathrm{1.003}\log {\mathit{T}}_{\mathrm{eff}}\mathrm{-}\mathrm{0.009.}\end{array}$These T_eff values are estimated to be accurate to within an error margin of  ± 140 K (Bergeat et al. 2001). Because this T_eff-calibration is based on a sample of 54 carbon-rich stars, preference was given to this calibration, rather than to the one by di Benedetto (1993), in case of carbon stars. The validity range was adopted to be 4.99 ≤ (V−K)₀ ≤ 8.31 mag, consistent with the minimum and maximum in Table 5 of Bergeat et al. (2001).

It is clearly visible from Fig. C.1 that the carbon stars have lower predicted effective temperatures than the oxygen-rich stars. This is in agreement with the results of Marigo (2002), showing a clear dichotomy in T_eff between both chemical types of AGB stars. This dichotomy is caused by different molecular opacities in stellar atmospheres with differing C/O abundance ratios.

When using these T_eff estimators, one should be aware of the effect of circumstellar reddening on the V and K band data. This reddening is stronger for objects with higher values of Ṁ and depends on the chemistry of the central star (Knapp et al. 1998, 2003). As mentioned in Appendix B, the data were not corrected for circumstellar reddening.

For oxygen-rich targets with (V−K)₀ outside the valid ranges for the above mentioned relations, the temperature corresponding to the spectral type listed in Table 6 was used (Table C.1; Marigo et al. 2008). An error bar of 500 K was assumed for these values. For the carbon-rich targets with (V−K)₀ outside the valid ranges the effective temperatures found in literature were used. References are given in Table 6.

C.1.2. OH/IR stars

In the case of OH/IR stars the photosphere is not easily probed since the CSE is optically thick and consequently the infrared sources often do not have optical counterparts. Therefore, the V−K colour is no longer a meaningful diagnostic of T_eff. For these objects an effective temperature of 2750 ± 750 K was assumed.

C.1.3. Red supergiants

The T_eff-calibration of Bessell et al. (1998) is not valid for red supergiants (RSGs) because of their low log g-values. Levesque et al. (2005) derived a T_eff-scale for RSGs in the Milky Way given by Eq. (C.8). This relation between the broad-band colour V−K and the effective temperature is consistent with the results from their grid of marcs models and is valid for 2.9mag ≤ (V−K)₀ ≤ 8.0mag. Adopting this new T_eff-calibration for RSGs yields a better agreement between evolutionary tracks and observationally determined locations of RSGs in the HR-diagram (see Massey et al. 2008, Fig. 3). This T_eff-scale for RSGs closely follows the relation of Bessell et al. (1998) defined for KM giants (see Fig. C.1), but since Eq. (C.8) was derived especially for RSGs, preference was given to this relation over Eq. (C.2). It is given by $\begin{array}{ccc}{\mathit{T}}_{\mathrm{eff}}& \mathrm{=}& \mathrm{7741.9}\mathrm{-}\mathrm{1831.83}\mathrm{(}\mathit{V}\mathrm{-}\mathit{K}{\mathrm{)}}_{\mathrm{0}}\\ & & \mathrm{+}\mathrm{263.135}\mathrm{(}\mathit{V}\mathrm{-}\mathit{K}{\mathrm{)}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{-}\mathrm{13.1943}\mathrm{(}\mathit{V}\mathrm{-}\mathit{K}{\mathrm{)}}_{\mathrm{0}}^{\mathrm{3}}\mathit{.}\end{array}$(C.8)Levesque et al. (2005) expect that the effective temperatures for M-type supergiants were obtained with a precision of 50 K.

For RSGs with (V−K)₀ < 2.9 or (V−K)₀ > 8.0, the temperature was determined using the scale of T_eff versus spectral type – in case of galactic targets – presented by Levesque et al. (2005).

C.1.4. Yellow hypergiants

The adopted T_eff-calibration of di Benedetto (1993) for FGKM targets of luminosity class I seems in accordance with the results presented by de Jager (1998) and Oudmaijer et al. (2009) for the yellow hypergiants (YHGs) IRC +10420 and AFGL 2343. Since both objects are known to exhibit large T_eff-variations on time scales of a few decades (Oudmaijer et al. 2009), this calibration is somewhat uncertain.

C.1.5. Alternative methods

Haniff et al. (1995) derive effective temperatures from angular diameters and bolometric fluxes obtained through fits to optical and infrared photometry. The infrared flux method (IRFM), presented by Blackwell & Shallis (1977), consists of determining the ratio of the bolometric flux to the flux in a selected photometric band and comparing these results to predictions from model atmospheres. Since the goal of this paper is not to individually model all sample targets, we did not pursue these methods, nor did we perform full fits of the spectral energy distributions (SEDs). The latter would require detailed radiative transfer analysis for each star.

C.2. Luminosity

The luminosities of Miras, OH/IR stars, semi-regular pulsators of types a and b and red supergiants were determined through empirical period-bolometric magnitude relations presented in the literature (Feast et al. 1989; Groenewegen & Whitelock 1996; Whitelock et al. 1991; Barthès et al. 1999; Whitelock & Feast 2000; Yeşilyaprak & Aslan 2004). These relations are then used to calculate luminosities adopting the solar bolometric magnitude M_bol, ⊙   = 4.75 mag.

C.2.1. Miras and OH/IR stars

Feast et al. (1989) derived period-luminosity, P−L _⋆ , relations for oxygen-rich and carbon-rich Miras in the LMC with periods ranging from 115 up to 420 days. Groenewegen & Whitelock (1996) revised the P−L _⋆ relation for carbon-rich Miras and extended the period range up to periods of 520 days. Whitelock et al. (1991) developed a similar relation for OH/IR targets with periods of pulsation exceeding 1000 days. These relations – and their extrapolations to shorter or longer pulsation periods – are shown in Fig. C.2.

The relation derived by Feast et al. (1989) for oxygen-rich Miras was assumed to be valid also for galactic oxygen-rich Miras and was adjusted with a distance modulus of 18.5 mag. It is given by ${\mathit{M}}_{\mathrm{bol}\mathit{,}\mathrm{O}\mathrm{-}\mathrm{MIRA}\mathit{,}\mathrm{OH}\mathit{/}\mathrm{IR}}\mathrm{=}\mathrm{-}{\mathrm{3.00}}_{\mathrm{\pm }\mathrm{0.24}}\log \mathit{P}\mathrm{+}{\mathrm{2.85}}_{\mathrm{\pm }\mathrm{0.57}}\mathit{,}$(C.9)with P in days.

Whitelock et al. (1991) mention a rather large scatter of 0.39 mag on their derived P−L _⋆ relation for OH/IR variables with periods ranging up to 2000 days. The large range in predicted luminosities associated with this scatter has been made visible in the figure. Note that the extrapolated relation for oxygen-rich Miras derived by Feast et al. (1989) matches the higher luminosity edge derived for OH/IR targets very well. Whitelock et al. (1991) also point out that there are suggestions that the luminosities of OH/IR sources with P > 1000 days probably follow this extrapolated P−L _⋆ relation that was established for oxygen-rich Miras by Feast et al. (1989). Because of these arguments and the strong evolutionary connection between oxygen-rich Miras and OH/IR stars (see Sect. 2.1) we adopt the P−L _⋆ relation of Feast et al. (1989) for oxygen-rich Miras, given in Eq. (C.9), for the determination of luminosities for all oxygen-rich Miras and OH/IR targets, regardless of their period of pulsation. Relation C.9 was also used in determining the luminosity of S-type Miras in the sample.

Luminosities of carbon-rich Miras were determined via ${\mathit{M}}_{\mathrm{bol}\mathit{,}\mathrm{C}\mathrm{-}\mathrm{MIRA}}\mathrm{=}\mathrm{-}\mathrm{2.59}\log \mathit{P}\mathrm{+}{\mathrm{2.02}}_{\mathrm{(}\mathrm{\pm }\mathrm{0.26}\mathrm{)}}$(C.10)following Groenewegen & Whitelock (1996). This relation predicts luminosities substantially higher than does the relation derived by Feast et al. (1989). Preference was given to the former since the sample of Groenewegen & Whitelock (1996) is almost three times as large as the one used by Feast et al. (1989) (54 versus 20 objects). For longer-period carbon-rich Miras, this P−L _⋆ relation was extrapolated, similar to what was done in the oxygen-rich case.

C.2.2. Semi-regular pulsators: SRa & SRb

Many semi-regular type a (SRa) objects differ from Miras only by their smaller amplitudes in the V band (ΔV < 2.5 mag). Objects classified as semi-regular b (SRb) variables have rather poorly defined periodicity, show a superposition of multiple pulsation periods, or have alternating intervals of periodic and slow irregular changes.

Semi-regulars and Miras are linked in terms of evolution, since SRs of types a and b are thought to be the progenitors of Miras (Whitelock & Feast 2000; Yeşilyaprak & Aslan 2004). Therefore some type of P−L _⋆ relation, analogous to those for Miras is probably justified for these LPVs.

Yeşilyaprak & Aslan (2004) presented period-magnitude relations for all magnitudes in the Johnson system and the IRAS system for SRa and SRb stars. To determine the bolometric magnitude, and hence the luminosity, only the J and K period – magnitude relations were used $\begin{array}{ccc}\mathit{J}& \mathrm{=}& \mathrm{-}{\mathrm{1.03}}_{\mathrm{(}\mathrm{\pm }\mathrm{0.56}\mathrm{)}}\log \mathit{P}\mathrm{-}{\mathrm{3.47}}_{\mathrm{(}\mathrm{\pm }\mathrm{0.98}\mathrm{)}}\\ \mathit{K}& \mathrm{=}& \mathrm{-}{\mathrm{1.38}}_{\mathrm{(}\mathrm{\pm }\mathrm{0.58}\mathrm{)}}\log \mathit{P}\mathrm{-}{\mathrm{3.99}}_{\mathrm{(}\mathrm{\pm }\mathrm{1.00}\mathrm{)}}\mathit{.}\end{array}$Yeşilyaprak & Aslan mention a scatter of 0.63 mag and 0.85 mag for the J and K relations, respectively. The bolometric correction on the K magnitude was calculated via $\begin{array}{ccc}\mathit{B}{\mathit{C}}_{\mathit{K}}& \mathrm{=}& \mathrm{10.86}\mathrm{-}\mathrm{38.10}\mathrm{(}\mathit{J}\mathrm{-}\mathit{K}{\mathrm{)}}_{\mathrm{0}}\\ & & \mathrm{+}\mathrm{64.16}\mathrm{(}\mathit{J}\mathrm{-}\mathit{K}{\mathrm{)}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{-}\mathrm{50.72}\mathrm{(}\mathit{J}\mathrm{-}\mathit{K}{\mathrm{)}}_{\mathrm{0}}^{\mathrm{3}}\\ & & \mathrm{+}\mathrm{19.48}\mathrm{(}\mathit{J}\mathrm{-}\mathit{K}{\mathrm{)}}_{\mathrm{0}}^{\mathrm{4}}\mathrm{-}\mathrm{2.94}\mathrm{(}\mathit{J}\mathrm{-}\mathit{K}{\mathrm{)}}_{\mathrm{0}}^{\mathrm{5}}\end{array}$(C.13)taken from Whitelock et al. (2000), where $\begin{array}{ccc}{\mathit{M}}_{\mathrm{bol}}& \mathrm{=}& \mathit{K}\mathrm{+}\mathit{B}{\mathit{C}}_{\mathit{K}}\mathit{.}\end{array}$(C.14)Relations (C.11) and (C.12) were however only adopted in estimating the luminosities for SRa variables, and not for the SRb variables in the sample. Barthès et al. (1999) presented P−(J−K) relations for long-period variables, including a large number of SRb type stars. The sample of Barthès et al. (1999) is divided in four groups, and the SRb type stars are largely split up in two of those – Groups 2 and 3, corresponding to younger and older kinematic properties respectively. Their Group 2 is far less luminous than their Group 3 and preference is given to the relations for their Group 3 – Eqs. (C.15) and (C.16) – since these estimates are in best accordance with luminosities derived in the literature for the SRb objects in the sample. Bolometric corrections were again derived via Eq. (C.13). $\begin{array}{ccc}& & {\mathit{K}}_{\mathrm{SRb}}\mathrm{=}\mathrm{-}{\mathrm{2.11}}_{\mathrm{(}\mathrm{\pm }\mathrm{0.13}\mathrm{)}}\mathrm{(}\log \mathit{P}\mathrm{-}\mathrm{1.75}\mathrm{)}\mathrm{-}{\mathrm{6.51}}_{\mathrm{(}\mathrm{\pm }\mathrm{0.05}\mathrm{)}}\\ & & \mathrm{(}\mathit{J}\mathrm{-}\mathit{K}{\mathrm{)}}_{\mathrm{0}}\mathrm{=}{\mathrm{0.17}}_{\mathrm{(}\mathrm{\pm }\mathrm{0.04}\mathrm{)}}\log \mathit{P}\mathrm{+}{\mathrm{0.85}}_{\mathrm{(}\mathrm{\pm }\mathrm{0.08}\mathrm{)}}\mathit{.}\end{array}$The uncertainties on calculated luminosities only take into account the uncertainties on the coefficients in Eqs. C.9 through C.16. The error bars on the pulsation periods are assumed to be small compared to the errors on the coefficients in the above relations.

C.2.3. Red supergiants

The luminosity of the RSGs in the sample was determined as a function of period of pulsation and effective temperature. Equation (C.17) gives the relation between the K magnitude and the pulsation period for RSGs in the milky way galaxy as presented by Dambis (1993). This relation was completed with the bolometric correction on the K band magnitudes given by Levesque et al. (2005) (Eq. (C.18)). This method leads to luminosities ranging from  ~10⁴ L_⊙ up to  ~510⁵ L_⊙. It is given by $\begin{array}{ccc}& & \mathit{K}\mathrm{=}\mathrm{-}{\mathrm{10.10}}_{\mathrm{(}\mathrm{\pm }\mathrm{0.08}\mathrm{)}}\mathrm{-}\mathrm{3.44}\mathrm{(}\log \mathit{P}\mathrm{-}\mathrm{2.7}\mathrm{)}\\ & & \mathit{B}{\mathit{C}}_{\mathit{K}}\mathrm{=}\mathrm{5.574}\mathrm{-}\mathrm{0.7589}\mathrm{\times }\frac{{\mathit{T}}_{\mathrm{eff}}}{\mathrm{1000}\mathrm{K}}\mathrm{·}\end{array}$

C.2.4. Yellow hypergiants

No general period-luminosity relations can be established for these targets because of their very strong variability. The latter is reflected in the multiple motions from blue to red and vice versa in the HR diagram (de Jager 1998; Oudmaijer et al. 2009). For both AFGL 2343 and IRC +10420, the two YHGs in the sample, luminosities were taken from de Jager (1998), yielding values for log L _⋆ /L_⊙ of 5.30 and 5.80, respectively.

C.2.5. Alternative methods

More recently, Barthès et al. (1999) and Knapp et al. (2003) derived P−K relations for oxygen-rich LPVs, including Miras and semi-regulars of types a and b. Using these relations in determining the luminosities of Miras would introduce extra uncertainties, since bolometric corrections for K band magnitudes would have to be derived. This is the main reason why the above mentioned P−M_bolrelations of Feast et al. (1989) and Groenewegen & Whitelock (1996) were used for Miras.

C.3. Distance

Distances towards the sample stars have been determined in various ways. (1) OH-maser phase lag distances were used for most OH/IR stars, since these provide a very accurate distance estimate (van Langevelde et al. 1990). (2) Hipparcos-data were used when available with a relative error lower than 50 %. (3) values in the literature derived from SED fitting were used and assumed to have uncertainties of 50 %, unless stated otherwise in the original papers.

C.3.1. Alternative method

An alternative method to derive distances is to compare observed magnitudes with absolute magnitudes predicted through formalisms as presented in Sect. C.2. The relations that could be used for distance determination from magnitudes for Miras are $\begin{array}{ccc}{\mathit{K}}_{\mathrm{O}\mathrm{-}\mathrm{MIRA}}& \mathrm{=}& \mathrm{-}{\mathrm{3.47}}_{\mathrm{(}\mathrm{\pm }\mathrm{0.19}\mathrm{)}}\log \mathit{P}\mathrm{+}{\mathrm{0.98}}_{\mathrm{(}\mathrm{\pm }\mathrm{0.45}\mathrm{)}}\\ {\mathit{K}}_{\mathrm{C}\mathrm{-}\mathrm{MIRA}}& \mathrm{=}& \mathrm{-}\mathrm{3.56}\log \mathit{P}\mathrm{+}\mathrm{1.14.}\end{array}$Since the K magnitude is less sensitive to metallicity than the bolometric magnitude m_bol, it is a better distance indicator (Whitelock & Feast 2000, and references therein). The strong variability in the K band of most sample stars, however, gives rise to large uncertainties on absolute magnitudes. An additional constraint to the application of this method to determine distances is that the circumstellar material can redden the K band magnitudes significantly.

---

Acknowledgments

The authors wish to thank Sofie Dehaes and Sacha Hony for their contribution to the data reduction and Hans Olofsson and David Teyssier for their contributions to the data set presented in this paper. Carlos De Breuck from ESO is acknowledged for the helpful discussions on APEX data reduction. The authors thank Jan Cami for proposing and planning the observations of EP Aqr. We thank Joris Blommaert and Rens Waters for the many fruitful discussions. E.D. and L.D. acknowledge support from the Fund for Scientific Research of Flanders (FWO) under grant number G.0470.07. The computations for this research have been done with the VIC HPC Cluster of the KULeuven. We are grateful to the LUDIT HPC team for their support. This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center, funded by the National Aeronautics and Space Administration and the National Science Foundation. APEX data were taken for programmes O-77.F-0605-2006, E-77.D-4004-2006, E-77.D-0781-2006, E-078.D-0534A-2007, and E-081.D-0207A-2008. The JCMT data were obtained in the framework of the science programmes M00AN09, M00BN19, M01AN04, M01BN10, M02BN05, M03AN02, M03BN06, M04AN20, M06BN03, M08AN01, and M08BN11.

References

All Tables

All Figures