text:
# The European Large Area ISO Survey V: a BeppoSAX hard X-ray survey of the S1 region

## 1 INTRODUCTION

Unified models of Active Galactic Nuclei (AGN) propose that all types of AGN are fundamentally the same although the presence of a dusty molecular torus prevents the observation of the broad line region for particular line of sights. In this picture Type 1 AGN (e.g. Seyfert 1/QSO) are those for which the nucleus can be directly observed, while in Type 2 AGN (e.g. Seyfert 2) it is obscured by the torus (e.g. Lawrence and Elvis, 1982 and Antonucci, 1993). Strong support for this model has been given by X-ray, near-IR and mid-IR observations showing that Type 2 AGN are characterized by strong absorption whilst Type 1 AGN are relatively unabsorbed (Turner et al. 1997, Alonso-Herrero, Ward, Kotilainen, 1997, Maiolino et al. 1998 and Clavel et al. 2000). In this framework the infrared radiation is interpreted in terms of thermal emission from hot dust grains heated by the high energy central source emission (optical to X-ray continuum, see e.g. Granato, Danese, Franceschini 1997).

Barcons et al. (1995 hereafter B95) investigated the 12$`\mu `$m to hard X-ray (HX, 2 to 10 keV) flux ratios of a sample of nearby AGN and galaxies using IRAS and HEAO1 A-2 observations. They found that Type 1 AGN have larger HX/12$`\mu `$m flux ratios than Type 2 AGN whilst normal galaxies have very weak HX emission. This implies high column densities (log(NH)$`>`$22 cm-2) in the Type 2 objects, supporting the basic unified model hypothesis. The study of B95 was necessarily restricted to nearby objects. The investigation of this flux ratio to deeper fluxes and higher redshifts will probe the properties of luminous sources such as QSOs and Broad Absorption Line QSOs (BALQSOs) and provide clues to the constancy of AGN activity in the universe.

The ASCA and BeppoSAX X-ray satellites have been successful in determining the HX properties of AGN and have resolved $``$30% of the HX background into discrete sources, mostly Type 1 and Type 2 AGN (e.g. Ueda et al. 1998, 1999, Cagnoni et al. 1998, Fiore et al. 1999 and Akiyama et al. 2000). The recent launch of the Chandra and XMM-Newton X-ray satellites have now made it possible to probe to fainter fluxes at HX energies. Recent Chandra observations (e.g. Mushotzky et al. 2000, Fiore et al. 2000a, Brandt et al. 2000, Giacconi et al. 2000 and Hornschemeier et al. 2000) have uncovered a number of HX emitting, but optically apparently normal galaxies. The HX slopes of the sources are sufficiently hard that, when combined with the HX emitting AGN sources, they can account for $`>`$75% of the HX background (e.g. Mushotzky et al. 2000, Giacconi et al. 2000 and Hornschemeier et al. 2001). Comparisons to deep sub-mm surveys have shown that these HX sources are almost always unassociated with sub-mm sources (e.g. Hornschemeier et al. 2000; Severgnini et al. 2000; Fabian et al. 2000, Barger et al. 2000), which instead appear to be star forming systems. Due to the expected association between HX and MIR emission a closer correlation should be found in the mid-IR band.

We present here BeppoSAX hard X-ray observations of the S1 region of the European Large Area ISO Survey (ELAIS) at 15$`\mu `$m. ELAIS (Oliver et al. 2000) was the largest program undertaken by ISO and, covering an area of $``$12 sq. deg, is the largest 15$`\mu `$m (Serjeant et al. 2000) and 90$`\mu `$m (Efstathiou et al. 2000) survey before SIRTF and FIRST come into play. In this paper we focus on observations in the Southern S1 region. In addition to the ISO observations this region has been covered in the U, B, R and I (La Franca et al. in prep., Heraudeau et al. in prep.) optical bands and the 1.4 GHz radio (Gruppioni et al. 1999) band. Selected objects have been observed spectroscopically in the optical (La Franca et al. in prep., Gruppioni et al. in prep.) and photometrically in the near-IR (Heraudeau et al. in prep.). The field of view of the BeppoSAX MECS (25 arcminutes of radius; Boella et al, 1997a,b) matches well with that of the ELAIS sub-fields, covering a considerably larger area than both Chandra and XMM-Newton (see Fig. 1). It is thus ideal for a shallow X-ray survey of ELAIS fields. Throughout this paper $`H_0`$=50 kms-1, $`\mathrm{\Lambda }`$=0 and $`q_0`$=0.5 are used.

## 2 OBSERVATIONS

The BeppoSAX observations were taken with both the LECS and MECS instruments. The LECS observations are not considered here due to their lower sensitivity and much poorer PSF below 2 keV.

The MECS observations (Boella et al. 1997b), from 1.3 to 10 keV, cover 5 of the 9 sub-fields in the S1 region of ELAIS (see Fig. 1). The total area covered is $``$1.7 sq.deg., although the sensitivity of MECS varies with off-axis distance with the most sensitive observations covering 0.7 sq.deg. The on-source exposure time for these 5 pointings was on average $``$36 ksec, corresponding to a flux limit of $``$10-13 $`\mathrm{erg}\mathrm{s}^1\mathrm{cm}^2`$ on-axis (see Table 1).

The BeppoSAX positional accuracy is dependent on many factors (see Fiore et al. 2000b for details) leading to a 95% error radius of one arcmin at off-axis distances $`<`$12 arcmin and 1.5 arcmin at greater distances. Slightly better accuracies have been achieved in those fields where we have been able to identify more than one source at another wavelength and perform a positional shift of the X-ray sources to match the more accurate optical positions (see Table 1).

## 3 DATA

### 3.1 ELAIS data

As discussed in the introduction, the characteristics of the ELAIS selection and the IR counts have been presented by Oliver et al. (2000) and Serjeant et al. . (2000). The survey has covered an area of $``$ 12 deg2 down to $``$1 mJy in the 15$`\mu `$m lw3 filter (referred to here as MIR, Cersarksky et al, 1996). In this paper we focus on observations in the S1 region centered at $`\alpha `$(2000) = 00h 34m 44.4s, $`\delta `$(2000) = -43 28 12′′, which covers an area of about $`2^{}\times 2^{}`$. About 100 sources of the Preliminary Analysis of ELAIS (Serjeant et al. 2000) down to R$``$20.5 mags have been spectroscopically identified during three spectroscopic identification campaigns in 1998 (La Franca et al. in prep., Gruppioni et al. in prep).

### 3.2 BeppoSAX data

The BeppoSAX observations were cleaned and linearized using the SDC on-line analysis. The MECS2 and MECS3 images were co-added together and binned by a factor of 2. Source detection was performed using XIMAGE with the interactive routine sosta. When measuring the counts of a source a box of 13 pixels (each pixel corresponds to 8 arcsecs) was used. This was found to maximize the signal to noise. As the PSF is much wider than the positional uncertainty, the effect on the flux from poor centroiding with weak sources is negligible. The net counts were automatically corrected for vignetting, sampling dead time and the point spread function. For each source the background was measured in large apertures at nearby orthogonal sky positions.

The conversion from counts to fluxes was determined using a power law of $`\alpha `$=0.7 ($`f_\nu \nu ^\alpha `$). We experimented with a range of column densities (20$`<`$log(NH)$`<`$23) but found small differences ($``$10%) in the conversion factors as compared to the uncertainty in the counts of the X-ray sources and the flux uncertainties of the MIR sources, which are of the order of 30% (Serjeant et al. 2000). The overall conclusions are not affected by these uncertainties.

We initially found 13 sources at the SNR$`>`$3 Poisson level in either the 1.3 to 10 keV (T) or 2 to 10 keV (HX) band although 2 of these sources could be unreliable because of possible source confusion (sources s2 and s4), and are included in a supplementary table (see Tables 2 and 3 and Fig. 1). Based on the source box size we would expect $`<`$1 spurious source due to background fluctuations. Considering the relatively small cosmic volume sampled, the SNRHX$`>`$3 number counts of sources are reasonably consistent with those found by Giommi, Perri and Fiore (2000), see Fig. 2. We have added to the supplementary list, Table 3, 2 sources just below the detection threshold, with SNRT=2.9, but which are close ($``$60 arcsec) to the position of an ELAIS 15$`\mu `$m QSO. These sources are used in all the further analysis.

The softness ratios of the sources were determined following the procedure of Fiore et al. (2000b) with the ratio defined as (S-H)/(S+H) where S=count rate in 1.3 to 4.5 keV band and H=count rate in 4.5 to 10 keV band, see Fig. 3 and Tables 2 and 3. The distribution of ratios is consistent with that found in the BeppoSAX HELLAS survey (i.e. $``$1$`<`$(S-H)/(S+H)$`<`$1, Fiore et al. 2000b).

### 3.3 Positional correlation between the BeppoSAX and the ISO sources

We have positionally correlated the HX sources to the MIR sources finding 6 matches (Table 4, Fig. 4). One further match is made with a 2.6$`\sigma `$ ISOCAM source. This source is included in the further analysis although it should be noted that this is not a significant ISOCAM detection. The separation between the X-ray and optical centroids are $`<`$ 60.5 arcsecs in all cases. The optical counterparts of the ISOCAM sources have magnitudes of R$`<`$20 mags. This optical magnitude limit is consistent with that found for $`f_{HX}>`$10-13 $`\mathrm{erg}\mathrm{s}^1\mathrm{cm}^2`$, the depth of this survey, in other HX surveys (e.g. Schmidt et al, 1998, Akiyama et al, 2000 and Hornschemeier et al, 2001).

The surface density of the ISOCAM sources is $``$150 sources/deg2, and therefore the chance of a mis-association with a BeppoSAX error box (60 arcsec radius) is 13%. In previous BeppoSAX and ASCA surveys (e.g. Akiyama et al. , 2000 and La Franca et al. 2000) sources have been matched to optical counterparts. In these studies the high optical source density often results in $`>`$1 possible optical counterpart within the large BeppoSAX and ASCA error boxes. The relatively low source density of MIR in this study results in less ambigious source matching although we are implictly assuming that an ISOCAM source within a BeppoSAX error box is associated with the HX source. This association is expected within the Unified model of AGN (Antonucci, 1993) and is found from both observational and theoretical tests of this model (e.g. B95, Granato, Danese and Franceschini, 1997, Krabbe, Boker and Maiolino, 2000 and Alonso-Herrero et al. , 2001). These X-ray and infrared observations are two orders of magnitude deeper than those of B95 (see Fig. 5).

### 3.4 Hard X-ray upper limits of ISO sources

We measured 3$`\sigma `$ X-ray upper limits for all the spectroscopically classified MIR sources, labeled in Fig 1. This was performed by measuring the counts at an MIR source position in the same fashion as for the detected sources. The measured upper limits depend on the off-axis distance and range from 10-13 $`\mathrm{erg}\mathrm{s}^1\mathrm{cm}^2`$ on-axis to 3-5x10-13 $`\mathrm{erg}\mathrm{s}^1\mathrm{cm}^2`$ off-axis.

## 4 DISCUSSION

In this discussion we comment on the HX-MIR source identifications and compare the observed HX-MIR properties of the BeppoSAX and ISOCAM objects to those of other samples previously reported in the literature (e.g. B95 and Elvis et al. 1994, hereafter E94).

### 4.1 The spectroscopic object identifications

As already discussed we found 7 MIR sources (including the 2.6$`\sigma `$ ISOCAM source) within the 16 error-boxes ($``$60 arcsec, 95% significance, see section 2) of the BeppoSAX detections. We have obtained optical spectroscopy for 6 of these sources (La Franca et al. in prep. and Gruppioni et al. in prep.), see Table 4. As the optical spectroscopic campaign targetted those sources with 17$`<R<`$20.5 mags, we did not obtain a spectroscopic identification for the brightest source. We find 4 normal broad emission line QSOs, with 0.4$`<z<`$2.6 and $`10^{44}<L_{HX}<10^{46}`$ erg s-1 and 1 Broad Absorption Line QSO (BALQSO) at $`z`$=2.2 and $`L_{HX}`$=6x1045 erg s-1. The 2.6$`\sigma `$ ISOCAM source is associated with a galaxy at $`z`$=0.325 and has $`L_{HX}`$$`10^{44}`$ erg s-1 if the association is real. Based on the distribution of object types found in other BeppoSAX and ASCA surveys (e.g. Akiyama et al, 2000 and La Franca et al. 2000) we would expect to find two Type 2 sources in our study. We find none which seems surprising but is not significant (1$`\sigma `$) given the small number of sources. One of the two low redshift QSOs ($`z<`$1) is radio detected (source 3, $`z`$=0.559, $`f_{1.4GHz}`$ = 1.5 mJy, Gruppioni et al. 1999). BALQSOs are rarely detected in the X-ray, particularly in the soft band, possibly due to large absorption of the X-ray emission (Gallagher et al. 1999). This is the first time a high redshift BALQSO has been detected in both the MIR and HX.

The unclassified source 8 is a bright optical source with an extended optical profile and is possibly detected at 90$`\mu `$m ($`f_{90\mu m}`$ = 81 mJy, Efsthathiou et al. 2000). Its soft HX ratio (section 4.2) is consistent with either an AGN or a thermal Bremsstrahlung emitting source. Its high HX/MIR flux ratio (section 4.3) suggests the former and it is probably a nearby Seyfert galaxy.

The detection of a normal galaxy seems surprising although many apparently normal galaxies have been detected in recent Chandra surveys (e.g. Mushotzky et al. 2000, Fiore et al. 2000a, Brandt et al. 2000, Giacconi et al. 2000 and Hornschemeier et al. 2000) although none have such a high X-ray luminosity and many tend to have harder X-ray spectral slopes. This object has an HX/MIR flux ratio, X-ray spectral slope and luminosity similar to QSOs although lacks the typical QSO broad optical emission lines. Other interesting possibilities are that it is a BL Lac object, although the lack of radio emission and a strong Calcium break suggest against this, or a dominant cD cluster galaxy with the X-ray emission coming from the cluster medium. The expected log(HX/MIR) ratio for this latter possibility is $``$4.7$`\pm `$0.4 (e.g. Edge and Stewart, 1991 and Bregman, McNamara and O’Connell, 1990). Due to the uncertain nature of this source these possibilities are not further pursued.

### 4.2 Softness ratios

We have compared the softness ratios (see Fiore et al. 2000b) to those expected for various HX SEDs, see Fig. 6. The softness ratios of the QSOs, although extremely uncertain, show a variety of soft and hard sources 0.5$`<`$(S-H)/(S+H)$`<`$-0.5, similar to those found in the BeppoSAX HELLAS survey (Fiore et al. 2000b).

### 4.3 The HX/MIR flux ratios

We have compared the observed HX/MIR flux ratios to those expected for a variety of object types. B95 statistically studied the correlation between HX and 12$`\mu `$m for the IRAS 12$`\mu `$m sample of AGN and emission line galaxies (Rush, Malkan, Spinoglio, 1993, hereafter RMS). The 12$`\mu `$m sample represents the properties of MIR sources in the local universe whilst the redshifts of these ISOCAM/BeppoSAX objects are substantially higher. Therefore in order to compare this study to that of B95 it is necessary to take into account the effects of redshift and the construction of spectral energy distributions (SEDs) are necessary. In the generation of SEDs we have decided to take an empirical approach due to the uncertain contribution of AGN and galactic activity at mid-IR wavelengths (e.g. see the detailed modelling of Cen A which is clearly an AGN although its infrared emission appears to be dominated by galactic processes, Alexander et al. , 1999). To do this we firstly evaluated and constructed HX and infrared SEDs of QSOs, Seyfert 2 and HII galaxies and then normalized the X-ray SED to reproduce the HX/MIR flux ratios of the local value at $`z`$0 as determined by E94 for QSOs, and by B95 for Seyfert 2 and HII galaxies. Only a few Seyfert 2s were detected at HX energies by B95 although the HX/MIR flux ratio distribution is very similar to the IRAS 60$`\mu `$m selected, HX detected 16 Seyfert 2 galaxies in Alexander (2001) (log(HX/MIR)=-7.0$`\pm `$0.7). The B95 HX/MIR flux ratio for the HII galaxies is an upper limit.

The HX SEDs were taken from Pompilio, La Franca and Matt (2000) for the cases of QSO and Seyfert 2; the QSO SED is a two part power law whilst the Seyfert 2 SED uses this same spectrum convolved with the Seyfert 2 absorbing column density distribution found by Maiolino et al. (1998). The Seyfert 2 galaxies in the Maiolino et al. study were not MIR selected. However, the column density distribution is very similar to that found for far-IR selected Seyfert 2 galaxies (i.e. log($`N_H`$)$``$22 to 25 cm-2, Alexander, 2001). The HII galaxy HX SED was determined assuming a Bremsstrahlung thermal spectrum with kT=5.8 keV, as found for the starburst galaxy NGC253 (Cappi et al. 1999).

The Seyfert 2 and HII galaxy infrared SEDs were determined using the Xu et al. (1998) empirical algorithm which takes the four IRAS band fluxes to predict an overall 2 to 120$`\mu `$m SED assuming three basic components: AGN, starburst and cirrus. The $`<`$12$`\mu `$m emission is predicted from the IRAS colours and produced using AGN, starburst and cirrus observational templates which include PAH emission and dust absorption features. The sample used to determine these SEDs was the RMS sample as classified by Alexander and Aussel (2000). These SEDs do not account for starlight, which can potentially be a large contributor for $`\lambda `$$`<`$7$`\mu `$m and would lead to underestimates of the HX/MIR ratios. Therefore we have only used these SEDs for z$`<`$1, which is sufficient for the non-QSO sources. For the QSOs we have used the empirically determined SED of E94 which includes the contribution from starlight and gives a good match to the low redshift optically selected PG QSOs (Schmidt and Green, 1983). As this SED is derived mostly from low redshift sources, we are assuming little SED evolution for z$`<`$3 (see section 4.3.1 for discussion).

The HX/MIR flux ratios of the BeppoSAX/ISO sample are shown in Fig. 7a. The errorbars associated with the SEDs take into account the statistical spread in the HX/MIR ratios measured locally by B95. In the case of Seyfert 2s and QSOs the effect of redshift is to increase the HX flux with respect to the MIR flux. The difference with redshift is not dramatic for either object type because the assumed HX column densities are quite low, although see section 4.3.2 for the case of a Compton thick source (log(NH)$`>`$24 cm-2). The HII galaxies are difficult to detect in the HX at any redshift but particularly at high redshift where the HX K-correction is positive. In the case of a higher temperature of the X-ray emitting gas than that assumed here, the HX/MIR ratio will stay constant to higher redshifts, dropping off significantly at the exponential cut-off energy (i.e. at $`z=(T_{HX}/2)1`$, where $`T_{HX}`$ is the Bremsstrahlung temperature in keV).

The flux ratio distribution of the HX/MIR upper limits are shown in Fig. 7b. The mean HX upper limit and log(HX/MIR) are 3.3x10-13 $`\mathrm{erg}\mathrm{s}^1\mathrm{cm}^2`$ and $``$5.1 respectively. To detect all the sources requires significantly deeper X-ray observations. HX observations of the Northern ELAIS fields will be carried out by Chandra and XMM-Newton to limiting fluxes of $``$3x10-15 and 10-14 $`\mathrm{erg}\mathrm{s}^1\mathrm{cm}^2`$ respectively. Assuming the mean flux and HX/MIR flux ratio of this BeppoSAX survey, the Chandra and XMM-Newton surveys should reach mean flux ratios of log(HX/MIR)$``$$``$7.1 and log(HX/MIR)$``$$``$6.6 although clearly some objects with lower ratios will also be detected. The Chandra survey (two ACIS-I 16 arcmin pointings) covers a small area but should detect all the AGN and some HII galaxies within the field of view whilst the XMM-Newton survey covers a similar area to this survey and should detect all the QSOs, most of the Seyfert 2 galaxies and a few HII galaxies.

#### 4.3.1 The QSOs

The majority of the HX-MIR associated sources are QSOs. The observed QSO HX/MIR flux ratios appear consistent with that of the median E94 QSO SED. The E94 UVSX sample includes 47 quasars selected to have at least 300 counts in the Einstein IPC and with V$`<`$17 mags to be observable by IUE; 29 members are radio quiet and 18 radio loud. The sample is consequently biased toward objects at low redshifts (although 7 objects have $`z>`$1, 4 $`z>`$2, 1 $`z>`$3), of moderate luminosity and with high X-ray to optical flux ratios.

In order to increase the statistics of our comparison we added to the MIR sample QSOs from Andreani, Franceschini and Granato (1999) with $`f_{12\mu m}>`$200 mJy, and created as a comparison an optically selected sample using 12$`\mu `$m detections or upper limits of the PG optically selected QSOs from Sanders et al. (1989). To increase the number of objects at high redshift we have taken those objects from Neugebauer et al. (1986) with z$`>`$2. We have determined which of these are BALQSOs using the list of Junkkarinen, Hewitt and Burbidge (1991). The data are shown in Fig. 8. The mean log(HX/MIR) ratios of the low redshift QSOs (z$`<`$0.6) are $``$5.7$`\pm `$0.7 and $``$5.6$`\pm `$0.3 for the MIR and optically selected objects respectively. These flux ratios are consistent with each other and with the E94 SED. This suggests that MIR selected QSOs of R$`<`$20 mags come from the same population as optically selected QSOs and that there is little SED evolution for z$`<`$3.

As a further and more definitive test we have compared the log(B/MIR) ratios of all the MIR detected QSOs in the S1 region to those of the PG QSO sample (Schmidt and Green, 1983 and Sanders et al, 1989), see Fig. 9. The advantage of this comparison is that we can use the full QSO set from both the ELAIS and PG surveys. The MIR QSOs are consistent with the flux ratios of the E94 SED and the PG QSOs. This further strengthens the evidence that MIR selected QSOs of R$`<`$20 mags come from the same population as optically selected QSOs.

#### 4.3.2 The BALQSO

Although the statistics are still poor, low redshift BALQSOs seem to be located in a different region of the HX/MIR plane to QSOs, see Fig. 8. This is consistent with that expected from an absorbed source. The low-redshift QSO, that lies close to the BALQSOs in the HX/MIR plane shown in Fig. 8, is a narrow line Seyfert 1 galaxy (IZw1), an unusually red source (e.g. see E94) which may be different to the general QSO population (e.g. see Brandt and Gallagher, 2000). Our BALQSO has a HX/MIR flux ratio similar to that found for normal QSOs, apparently in contradiction with that found for low redshift BALQSOs. This difference could be due to the negative HX K-correction effect of an absorbed HX spectrum, see Fig. 6. We can test whether the observed HX/MIR flux ratio is compatible with that expected for an absorbed HX source by constructing a possible BALQSO SED. We have chosen to use the Compton thick HX spectrum of the Type 2 QSO object IRAS09104+4109 (Franceschini et al. 2000) and the QSO infrared SED of E94, normalised to the HX/MIR ratio of IRAS09104+4109. The positive effect of redshift on the HX emission can clearly be seen (see Fig. 8) showing that a high redshift HX absorbed source can have a similar HX/MIR flux ratio to an unabsorbed source. Therefore, in contrast to low redshift BALQSOs, high redshift BALQSOs should be comparatively easy to detect in the HX band, allowing a determination of their true fraction in the high-z QSO population.

## 5 CONCLUSIONS

We have presented shallow BeppoSAX observations, reaching an on-axis ($``$0.7 sq. deg) 2-10 keV sensitivity of $``$10-13 $`\mathrm{erg}\mathrm{s}^1\mathrm{cm}^2`$, of the Southern S1 region of ELAIS, reaching a 15$`\mu `$m sensitivity of 1 mJy. This is the first HX analysis of an ISOCAM survey. We have constructed HX and infrared SEDs to determine the expected flux ratios for QSOs and BALQSOs up to $`z`$=3.2 and for Seyfert 2s and normal galaxies up to $`z`$=1.0. Our main findings are:

(i) we detect 13 sources with SNR$`>`$3 in the 1.3-10 keV or 2-10 keV X-ray bands and a further 2 sources with less reliable detections that have positions close ($``$60 arcsec) to QSOs. The number densities of the SNR$`{}_{HX}{}^{}>`$3 sources are consistent with the ASCA and BeppoSAX logN-logS function.

(ii) 6 of these sources have a reliable ISOCAM counterpart and one further source has a less reliable (2.6$`\sigma `$) ISOCAM counterpart. We have optical spectroscopic classifications for 6 of these sources finding 4 QSOs, 1 BALQSO at $`z`$=2.2 and 1 apparently normal galaxy (the 2.6$`\sigma `$ ISOCAM source). The unclassified object has X-ray and photometric properties consistent with that of a nearby Seyfert galaxy. The galaxy has properties suggesting either an unusual QSO or a galaxy cluster.

(iii) the QSOs cover a wide redshift range (0.4$`<`$z$`<`$2.6), and have HX/MIR flux ratios consistent with those found for nearby IRAS and optically selected PG QSOs and the QSO SED of Elvis et al (1994). By further comparing the B/MIR flux ratios of the MIR QSOs to those of the blue band selected PG sample, we suggest that MIR selected QSOs of R$`<`$20 mags come from the same population as optically selected QSOs.

(iv) the high redshift BALQSO has a HX/MIR ratio similar to that of QSOs, but different to that found for low redshift BALQSOs. This difference can be explained as the negative K-correction effect of an absorbed X-ray spectrum observed at high redshift. This suggests that high redshift BALQSOs should be comparatively easy to detect in the HX band, allowing the true fraction of BALQSOs in the high redshift QSO population to be determined.

We acknowledge the EC TMR network (FMRX-CT960068) for financial support and MURST Cofin-98-032 and ASI contracts for partial support. DMA additionally thanks the NSF CAREERS grant AST-998 3783 for post-doctoral support. This research has made use of SAXDAS (SAX Data Analysis System) linearized and cleaned event files (Rev.2.0) produced at the BeppoSAX Science Data Center and the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with NASA. We gratefully thank Belinda Wilkes, the referee, for an efficient and thorough reading of this manuscript. We further thank Cong Xu for providing the code to determine the mid-IR SEDs used in this paper and Sarah Gallagher, Ann Hornschemeier and Niel Brandt for valuable comments on earlier drafts of this paper.


format: text

text:
# Models of circumstellar molecular radio line emissionPresented in this paper is observational data collected using the Swedish-ESO submillimetre telescope, La Silla, Chile, the 20 m telescope at Onsala Space Observatory, Chalmers Tekniska Högskola, Sweden, and the NRAO 12 m telescope located at Kitt Peak, USA.

## 1 Introduction

The carbon star phenomenon, i.e., the presence of stars for which the abundance of carbon exceeds that of oxygen in the atmospheres, occurs during the final evolutionary stage of low to intermediate mass stars ($``$1$``$10 M; Wallerstein & Knapp 1998), and it is believed to be due to a dredge-up process driven by quasi-periodic He-shell flashes (Straniero et al. 1997). These stars, located on the asymptotic giant branch (AGB), lose copious amounts of matter ($``$10-8$``$10-4 M yr-1) in an intense stellar wind. In fact, it is the mass loss that determines the AGB lifetime and not the nuclear burning processes. Hence, a determination of the mass loss characteristics, i.e., the dependence on time, mass, chemistry, metallicity, etc., is of utmost importance for our understanding of late stellar evolution.

The mass loss creates a circumstellar envelope (CSE) of gas and dust around the star. The low temperature of the central star allows the formation of a wide variety of molecular species in its atmosphere, and the expanding gas is therefore mainly in molecular form. The chemistry in the CSE itself can be very rich, and it depends on the C/O-ratio, the thickness of the envelope, and the strength of the ambient ultraviolet radiation field. Presently, $``$60 different molecular species have been identified in CSEs around AGB-stars (Olofsson 1997).

The study of AGB-CSEs is important for the understanding of the late stages of stellar evolution as mentioned above. In addition, an understanding of the atmospheric and circumstellar chemistry is required to determine the elemental abundances in these winds, which contribute to the chemical evolution of the interstellar medium (Busso et al. 1999). Finally, since AGB stars are the progenitors of planetary nebulae (PNe), the CSEs are very likely important ingredients in the formation of PNe (Kwok 1993).

Observations of molecular millimetre-wave line emission have proven to be one of the best tools for studying the structure, kinematics, and chemistry of the CSEs, and CO and OH observations have proven particularly useful for the determination of accurate mass loss rates \[see Olofsson Olofsson (1996) and references therein\].

While the overall picture of mass loss on the AGB is fairly well understood, some of the underlying physical principles are not. The prevailing theory, which needs to be thoroughly tested observationally, is that the mass loss occurs in a two stage process. First the pulsations of the star deposit energy in the atmosphere, leading to a considerably increased scale height, and sufficient matter at low enough temperatures for efficient dust formation. Radiation pressure on the dust grains accelerate the dust wind to its terminal expansion velocity in the second stage. The gas is momentum coupled, through collisions, to the dust, and it will be effectively dragged along. The mechanisms that cause the mass loss to vary on a wide range of time scales still need to be pinpointed. Eventually, the possibility of mass ejected in clumps has to be addressed, although in this paper we will consider only a smooth spherically symmetric wind.

This study of an essentially complete sample of visually bright, relatively unobscured carbon stars will hopefully increase our knowledge of the modelling of circumstellar molecular radio line emission, the mass loss characteristics, the circumstellar chemistry, the elemental abundances in the winds, and possibly the evolutionary status of these objects. It will also provide us with the tools required to determine the properties of higher mass loss rate objects, for which the central star may be completely obscured.

In this paper the basic physical characteristics of the CSEs, e.g., the stellar mass loss rate, the gas kinetic temperature, and the gas expansion velocity, are determined using a radiative transfer analysis of the observed 12CO (hereafter CO) millimetre-wave line emission. The CO molecule is particularly well suited for this purpose, since it is difficult to photodissociate and easy to excite through collisions, and thus is a very good tracer of the molecular gas density and temperature. Furthermore, the molecular excitation is of a quasi-thermal nature which simplifies detailed modelling.

## 2 Observations

### 2.1 The sample

We have selected a sample of 68 bright N- and J-type carbon stars for which circumstellar CO emission was detected by Olofsson et al. (1993a). Their original sample consisted of carbon stars brighter than $`K`$$`=`$2 mag., in total 120 sources. The detection rate for all the carbon stars searched for circumstellar CO was as high as 72%. The original sample is thought to be close to complete to distances below 1 kpc \[it is estimated that about a third of the carbon stars within this distance are missing, presumably the ones with thick and dusty CSEs, i.e., the high mass loss objects (Olofsson et al. 1993a)\], and for distances below 500 pc all sources (a total of 41) were detected in CO (note that the distances used in this paper are generally lower than those adopted in Olofsson et al. 1993a). These facts taken together led to the conclusion that the great majority of all N-type carbon stars are losing mass at a rate higher than about 10-8 M yr-1. They do not, however, constitute a homogeneous group in terms of circumstellar properties. Instead they show a wide range of gas expansion velocities and mass loss rates, presumably due to differences in mass, evolutionary stage, etc.. In addition, the reasonably well studied, high mass loss rate carbon star, LP And (also known as IRC+40540), is included in our analysis, since it provides a good test case for the radiative transfer model.

The CO ($`J`$$`=`$$`1`$$``$$`0`$ and $`J`$$`=`$$`2`$$``$$`1`$) observations presented in Olofsson et al. Olofsson et al. (1993a) were obtained using the 20 m telescope at Onsala Space Observatory (OSO), Sweden, the 15 m Swedish-ESO sub-millimetre telescope (SEST) at La Silla, Chile, and the IRAM 30 m telescope at Pico Veleta, Spain, during the years 1986$``$1992.

### 2.2 New observations

We have supplemented the original data with new CO observations of some of the sample stars. In particular, in the $`J`$$`=`$$`3`$$``$$`2`$ line using the SEST (August 1992). Also some $`J`$$`=`$$`1`$$``$$`0`$ and $`J`$$`=`$$`2`$$``$$`1`$ data were obtained at SEST (October 1998). A few sources were observed in the $`J`$$`=`$$`1`$$``$$`0`$ and $`J`$$`=`$$`2`$$``$$`1`$ lines with the NRAO 12 m telescope at Kitt Peak, USA (June 1998), and some $`J`$$`=`$$`1`$$``$$`0`$ data were obtained using OSO (1998 to 1999).

A summary of telescope and receiver data \[the receiver temperature $`T_{\mathrm{rec}}`$ in single sideband mode (SSB), the full half power main beam width $`\theta _{\mathrm{mb}}`$, the main beam efficiency $`\eta _{\mathrm{mb}}`$, and the corrected main beam efficiency $`\eta _\mathrm{m}^{}`$\] at the observational frequencies are given in Table 1.

At SEST, a dual channel SIS mixer receiver was used to simultaneously observe at 115 GHz (the $`J`$$`=`$$`1`$$``$$`0`$ line) and 230 GHz (the $`J`$$`=`$$`2`$$``$$`1`$ line). At 345 GHz (the $`J`$$`=`$$`3`$$``$$`2`$ line), a single polarization SiS mixer receiver was used. As backends, two acousto-optical spectrometers (AOS) were used; one high resolution spectrometer (HRS) and one low resolution spectrometer (LRS). The wideband (1 GHz) LRS had 1440 channels separated by 0.7 MHz, whereas the narrow band (84 MHz) HRS used 2000 channels separated by 42 kHz.

For the NRAO 12 m observations we used a dual polarization SIS receiver. As backends two filterbanks, each with 256 channels and a channel separation of 1 MHz, were used. The spectra obtained at each of the polarizations were added to increase the signal-to-noise ratio. In addition, for strong sources, a millimetre autocorrelator (MAC) was used in order to obtain a higher velocity resolution. The MAC was used in a configuration where the usable bandwith was 300 MHz and the resolution was 98 kHz.

At OSO, a single polarization SIS receiver was used for the observations. As backends two filterbanks, with bandwidths of 512 MHz (MUL A) and 64 MHz (MUL B), were used. The MUL A used 512 channels separated by 1 MHz, and the MUL B filterbank used 256 channels with a separation of 250 kHz.

The SEST and OSO observations were made in a dual beamswitch mode, where the source is alternately placed in the signal and the reference beam, using a beam throw of about 12$`\mathrm{}`$ at SEST and 11$`\mathrm{}`$ at OSO. This method produces very flat baselines. The intensity scales are given in main beam brightness temperature, $`T_{\mathrm{mb}}`$. $`T_{\mathrm{mb}}`$$`=`$$`T_\mathrm{A}^{}/\eta _{\mathrm{mb}}`$, where $`T_\mathrm{A}^{}`$ is the antenna temperature corrected for atmospheric attenuation using the chopper wheel method, and $`\eta _{\mathrm{mb}}`$ is the main beam efficiency, see Table 1.

At the NRAO 12 m telescope the observations were carried out using a position switching mode, with the reference position located $`+`$10$`\mathrm{}`$ in azimuth. This is the preferred observation mode for spectral line observations at the 12 m telescope. The raw spectra, which are stored in the $`T_\mathrm{R}^{}`$ scale, were converted using $`T_{\mathrm{mb}}`$$`=`$$`T_\mathrm{R}^{}/\eta _\mathrm{m}^{}`$, where $`\eta _\mathrm{m}^{}`$ is the corrected main beam efficiency given in Table 1. The $`T_\mathrm{R}^{}`$ scale is related to $`T_\mathrm{A}^{}`$ through $`T_\mathrm{R}^{}`$$`=`$$`T_\mathrm{A}^{}`$/$`\eta _{\mathrm{fss}}`$, where $`\eta _{\mathrm{fss}}`$ is the forward scattering and spillover efficiency.

Regular pointing checks were made on strong SiO masers (OSO, SEST) and strong continuum sources (NRAO). The pointing was usually consistent within $``$3$`\mathrm{}`$ for SEST and OSO and $``$5$`\mathrm{}`$ at NRAO. The uncertainties in the absolute intensity scales at the various telescopes are estimated to be about $`\pm `$15$``$20%. However, due to the low efficiency of the SEST at 345 GHz, and the narrow beam, we consider these data to be particularly uncertain in terms of calibration.

The new observational results are summarized in Table 2 and the spectra are shown in Fig. 1. The line parameters, i.e., the main beam brightness temperature at the line centre ($`T_{\mathrm{mb}}`$), the line centre velocity ($`v_{}`$), and half the full line width ($`v_\mathrm{e}`$), are obtained by fitting the following line profile to the data (Olofsson et al. 1993a)

$$T(v)=T_{\mathrm{mb}}\left[1\left(\frac{vv_{}}{v_\mathrm{e}}\right)^2\right]^{\beta /2},$$
(1)
where $`\beta `$ is a parameter describing the shape of the line ($`\beta `$=2 represents a parabolic line shape, and $`\beta `$$`<`$0 means a profile with horns at the extreme velocities). The integrated intensity ($`I_{\mathrm{mb}}`$) is obtained by integrating the emission between $`v_{}\pm `$$`v_\mathrm{e}`$.

In addition, we have obtained publicly available data from the James Clerk Maxwell Telescope (JCMT) at Mauna Kea, Hawaii. The JCMT data are taken at face value. However, in the cases where there are more than one observation available, the derived line intensities are generally consistent within $`\pm `$20%. The good agreement with corresponding SEST intensities further support the reliability of the JCMT public data. Furthermore, interferometer observations of the CO($`J`$=$`1`$$``$$`0`$) brightness distribution around some of our carbon stars have been obtained by Neri et al. Neri et al. (1998) using the IRAM Plateau de Bure interferometer (PdBI). The data are publicly available and have been used in this paper.

## 3 Radiative transfer

### 3.1 The Monte Carlo method

In order to model the circumstellar molecular line emission and to derive the basic characteristics of the CSEs we have developed a non-LTE radiative transfer code based on the Monte Carlo method \[Bernes Bernes (1979); see also Schöier Schöier (2000), for details\]. An accurate treatment of the molecular excitation is needed for a correct description of the radiative transfer in an expanding circumstellar envelope, e.g., the Sobolev approximation has been shown to introduce significant errors in these circumstances (e.g., Schönberg 1985). The Monte Carlo method is well suited for the study of CSEs, since it is very flexible and close to the physics and yet simple in principle. For instance, one may include very complex geometries and velocity fields without being forced away from the physics of the problem. Crosas & Menten 1997 have recently used the Monte Carlo method to model the circumstellar CO radio line emission of the prominent carbon star IRC+10216 (from hereon called CW Leo).

The aim of the radiative transfer is to obtain the steady-state level populations, of the molecule under study, using the statistical equilibrium equations (SE). In the Monte Carlo method information on the radiation field is obtained by simulating the line photons using a number of model photons, each representing a large number of real photons from all transitions considered. These model photons, emitted both locally in the gas as well as injected from the boundaries of the CSE, are followed through the CSE and the number of absorptions are calculated and stored. Photons are spontaneously emitted in the gas with complete angular and frequency redistribution (CRD), i.e., the local emission is assumed to be isotropic and the scattering are assumed to be incoherent. The weight of a model photon is continuously modified, to take the absorptions and stimulated emissions into account, as it travels through the CSE. When all model photons are absorbed in, or have escaped, the CSE the SE are solved and the whole process is then repeated until some criterion for convergence is fulfilled. Once the molecular excitation, i.e., the level populations, is obtained the radiative transfer equation can be solved exactly.

The main drawback of the Monte Carlo method is its slow convergence ($``$$`\sqrt{N_{\mathrm{iter}}}`$). This is, however, usually outweighed by its great adaptability. One of our motivations for choosing the Monte Carlo method was to be able to treat varying mass loss rates, departures from spherical symmetry, and the possibility of a highly clumped medium. These complicating issues will be investigated in a future version of our Monte Carlo code.

### 3.2 The standard model

The numerical modelling of line emission from CSEs, where intricate interplays between physical and chemical processes take place, is a challenge. This applies, in particular, to the important inner regions of the CSE model, where very few observational constraints are available and our knowledge is very limited. In this analysis, we will neglect many of these complexities and assume a relatively simple, yet reasonably realistic, CSE model. By applying the model to a large sample we will derive fairly reliable mass loss rates, and we will also be able to pick out objects for which there appears to be clear deviations from the simple picture.

In what will be referred to as ”the standard model” we assume a spherically symmetric CSE that expands at a constant velocity. In our data we see no evidence for gas acceleration in the CO emitting region. The CSE is divided into a number of discrete shells, each with its own set of physical parameters describing the state of the molecular gas. The density structure (as a function of distance, $`r`$, from the central star) can then be derived from the conservation of mass

$$\rho _{\mathrm{H}_2}=n_{\mathrm{H}_2}m_{\mathrm{H}_2}=\frac{\dot{M}}{4\pi r^2v_\mathrm{e}},$$
(2)
where $`\dot{M}`$ is the hydrogen gas mass loss rate, and $`v_\mathrm{e}`$ is the gas expansion velocity taken from the line profile fits of Olofsson et al. (1993a) or Table 2. We assume that in the CSEs of interest here the hydrogen is in molecular form in the region probed by the CO emission (Glassgold & Huggins, 1983). The turbulent velocity is assumed to be equal to 0.5 km s-1 throughout the entire CSE.

The excitation of the CO molecules were calculated taking into account 30 rotational levels in each of the ground and first vibrational states. The radiative transition probabilities and energy levels are taken from Chandra et al. Chandra et al. (1996), and the rotational collisional rate coefficients (CO-H2) are based on the results in Flower & Launay Flower and Launay (1985). These are further extrapolated for $`J`$$`>`$11 and for temperatures higher than 250 K. Collisional transitions between vibrational levels are not important due to the low densities and their short radiative lifetimes.

The kinetic gas temperature is calculated after each iteration in a self-consistent way, using the derived level populations, by solving the energy balance equation (e.g., Goldreich & Scoville 1976),

$$\frac{dT}{dr}=(22\gamma )\frac{T}{r}+\frac{\gamma 1}{n_{\mathrm{H}_2}kv_\mathrm{e}}(HC),$$
(3)
where $`\gamma `$ is the adiabatic index \[assumed to be equal to 5/3, but the results are not very sensitive to its exact value, see also Groenewegen Groenewegen (1994) and Ryde et al. Ryde et al. (1999)\], $`k`$ is the Boltzmann constant, $`H`$ is the total heating rate per unit volume, and $`C`$ is the total cooling rate per unit volume. The first term on the right hand side is the cooling due to the adiabatic expansion of the gas. Additional cooling is provided by molecular line emission from CO and H2. The molecular cooling due to CO is calculated from the derived level populations using the expression of Crosas & Menten Crosas and Menten (1997). For the H2 cooling we use the approach by Groenewegen Groenewegen (1994). HCN could be an important coolant in the inner parts of the envelope (Cernicharo et al. 1996), but Groenewegen Groenewegen (1994) demonstrates that HCN cooling is only of minor importance. In the present version of the code HCN cooling is therefore not included (in a forthcoming paper, Schöier & Olofsson in prep., this issue will be addressed).

The mechanism responsible for the observed mass loss is probably (at least for the stars of interest here) radiation pressure acting on small dust grains, which in turn are coupled to the gas (e.g., Höfner 1999). The radiation pressure on the dust grains will give them a drift velocity, $`v_{\mathrm{dr}}`$, relative to the gas (Gilman 1972; Goldreich & Scoville 1976; Kwan & Hill 1977)

$$v_{\mathrm{dr}}=\left(\frac{Lv_\mathrm{e}Q}{\dot{M}c}\right)^{1/2},$$
(4)
where $`L`$ is the luminosity, $`Q`$ is the averaged momentum transfer efficiency, and $`c`$ is the speed of light \[see Sahai (1990) for a more elaborate treatment\]. As a result of the dust-gas drift, kinetic energy of the order of $`\frac{1}{2}m_{\mathrm{H}_2}v_{\mathrm{dr}}^2`$ will be transfered to the gas each time a particle collides with a dust grain. This is assumed to provide the dominating heating of the gas (Goldreich & Scoville 1976; Kwan & Hill 1977),

$$H_{\mathrm{dg}}=(n_\mathrm{d}\sigma _\mathrm{d}v_{\mathrm{dr}})\frac{1}{2}\rho _{\mathrm{H}_2}v_{\mathrm{dr}}^2,$$
(5)
where $`n_\mathrm{d}`$ is the number density of dust grains, and $`\sigma _\mathrm{d}`$ the geometrical cross section of a dust grain. Introducing the dust-to-gas mass ratio $`\mathrm{\Psi }`$ this equation can be written

$$H_{\mathrm{dg}}=\frac{3}{8}m_{\mathrm{H}_2}^2n_{\mathrm{H}_2}^2\frac{\mathrm{\Psi }}{a_\mathrm{d}\rho _\mathrm{d}}\frac{v_{\mathrm{dr}}^3}{1+\frac{v_{\mathrm{dr}}}{v_\mathrm{e}}},$$
(6)
where $`a_\mathrm{d}`$ and $`\rho _\mathrm{d}`$ are the average size and density of a dust grain, respectively.

When solving the energy balance equation free parameters describing the dust, i.e., via the dust-gas collisional heating, are introduced. These are highly uncertain, but affect the derived line intensities. Here we assume that the $`Q`$-parameter, i.e., the efficiency of momentum transfer, to be equal to 0.03 \[for details see Habing et al. (1994)\], and define a new parameter that contains the other dust parameters,

$$h=\left(\frac{\mathrm{\Psi }}{0.01}\right)\left(\frac{2.0\mathrm{g}\mathrm{cm}^3}{\rho _\mathrm{d}}\right)\left(\frac{0.05\mu \text{m}}{a_\mathrm{d}}\right),$$
(7)
where the normalized values are the ones we used to fit the CO line emission of CW~Leo using our model, i.e., $`h`$$`=`$1 for this object.

Additional heating is provided by the photoelectric effect (Huggins et al. 1988)

$$H_{\mathrm{pe}}=k_{\mathrm{pe}}n_{\mathrm{H}_2},$$
(8)
where $`k_{\mathrm{pe}}`$ is a constant that depends on the properties of the dust grains. Following Huggins et al. Huggins et al. (1988) we adopt $`k_{\mathrm{pe}}`$$`=`$1$`\times `$10-26 erg s-1. This heating mechanism is important in the cool, tenuous, outer parts of CSEs around high mass loss rate stars.

The spatial extent of the molecular envelope is generally an important input parameter, and the derived mass loss rate will depend on this. The radial distribution of CO in the CSE was estimated using the modelling presented in Mamon et al. Mamon et al. (1988). It includes photodissociation, taking into account the effects of dust-, self- and H2-shielding, and chemical exchange reactions. Mamon et al. Mamon et al. (1988) find that the radial abundance distribution of CO with respect to H2, $`f(r)`$, can be described by

$$f(r)=f_0\mathrm{exp}\left[\mathrm{ln}2\left(\frac{r}{r_\mathrm{p}}\right)^\alpha \right],$$
(9)
where $`f_0`$ is the initial (photospheric) abundance, $`r_\mathrm{p}`$ is the photodissociation radius (where the abundance has dropped to $`f_0/2`$), and $`\alpha `$ is a parameter describing the rate at which the abundance declines. Both $`r_\mathrm{p}`$ and $`\alpha `$ depend on the mass loss rate ($`\dot{M}`$), the expansion velocity ($`v_\mathrm{e}`$), and $`f_0`$. In our modelling we use analytical fits to these results from Stanek et al. Stanek et al. (1995)

$`r_\mathrm{p}`$ $`=`$ $`5.4\times 10^{16}\left({\displaystyle \frac{\dot{M}}{10^6}}\right)^{0.65}\left({\displaystyle \frac{v_\mathrm{e}}{15}}\right)^{0.55}\left({\displaystyle \frac{f_0}{8\times 10^4}}\right)^{0.55}`$ (10)
$`+\mathrm{\hspace{0.17em}7.5}\times 10^{15}\left({\displaystyle \frac{v_\mathrm{e}}{15}}\right)\mathrm{cm}`$
and Kwan & Webster Kwan and Webster (1993)

$$\alpha =2.79\left(\frac{\dot{M}}{10^6}\right)^{0.09}\left(\frac{15}{v_\mathrm{e}}\right)^{0.09},$$
(11)
where the units of $`\dot{M}`$ and $`v_\mathrm{e}`$ are $`\mathrm{M}_{\mathrm{}}\mathrm{yr}^1`$ and $`\mathrm{km}\mathrm{s}^1`$, respectively. In what follows we will assume $`f_0`$=1$`\times `$10-3. This is an average of the $`f_0`$:s estimated by Olofsson et al. Olofsson et al. (1993b) for this sample of optically bright carbon stars. The scatter in the estimated $`f_0`$:s is about 40%. We note here that a somewhat more sophisticated photodissociation model was developed by Doty & Leung Doty and Leung (1998) in which they include scattering by dust. For the CSE of CW~Leo they derive a CO envelope size which is 30% smaller than what is obtained using the Mamon et al. model. It is however not clear by how much the envelope size changes for an object with a significantly lower mass loss rate, and hence we use the results of Mamon et al. here.

In our models we include both a central source of radiation and the cosmic microwave background radiation at 2.7 K. The central radiation emanates from the star, and was estimated from a fit to its spectral energy distribution (SED), usually by assuming two blackbodies. This method is described in Kerschbaum Kerschbaum (1999). A fit to the SED gives the two blackbody temperatures, $`T_{}`$ and $`T_\mathrm{d}`$, and the relative luminosities of the two blackbodies, $`L_\mathrm{d}`$/$`L_{}`$. Any dust present around the star will absorb parts of the stellar light and re-emit it at longer wavelenghts. Thus, of the two blackbodies used, one represents the stellar contribution and one represents the dust. However, it should be noted that, e.g., the stellar temperature $`T_{}`$ derived in this manner is generally lower than the true temperature of the star (Kerschbaum 1999). In any case, the two blackbodies used give a good description of the radiation that the CSE is subjected to. The temperatures and luminosities used in the modelling are presented in Table 3. Only in the cases where $`L_\mathrm{d}`$/$`L_{}`$$`>`$0.1 the dust component was retained. The inner boundary of the CSE, $`r_\mathrm{i}`$, was set to reside outside the radius of the central blackbody(s), but never lower than 1$`\times `$1014 cm ($``$3 R), i.e., generally beyond both the sonic point and dust condensation radius.

The distances, presented in Table 3, were estimated using one of the following methods: the observed Hipparcos parallax, a period-luminosity relation (Groenewegen & Whitelock 1996), or an assumed bolometric luminosity. In the two former cases the luminosities were estimated using apparent bolometric magnitudes and the distances. In the rare cases when both the distance obtained from Hipparcos and the period-luminosity relation result in extremely high or low luminosities, a value of 4000 L was adopted as the luminosity and the distance was estimated using this value. For the irregular variables with no Hipparcos data we also assumed a luminosity of 4000 L. For a statistical study of a large sample of stars these distance estimates are adequate, although the distance estimate for an individual star has a large uncertainty. There is no apparent systematic difference between the distances derived from the period-luminosity relation and the Hipparcos parallaxes, Fig. 2. The distances in Olofsson et al. (1993a) were estimated from an adopted absolute $`K`$ magnitude. The distances presented here are systematically lower by, on the average, a factor of 1.4 than those used in Olofsson et al.

This was a summary of the assumptions used in ”the standard model”. In what follows, tests will be made to see how sensitive the model is to the various assumptions.

## 4 Model results

### 4.1 The fitting strategy

With the assumptions made in the standard model there remains two principal free parameters, the mass loss rate ($`\dot{M}`$) and the $`h`$-parameter. These two parameters were allowed to vary simultaneously, $`\dot{M}`$ in steps of $``$10% and $`h`$ in somewhat larger steps, until the model with the smallest deviations from the observed intensities was found. The quality of a particular model with respect to the observational constraints can be quantified using the chi-square statistic defined as

$$\chi ^2=\underset{i=1}{\overset{N}{}}\frac{[I_{\mathrm{mod},i}I_{\mathrm{obs},i})]^2}{\sigma _i^2},$$
(12)
where $`I`$ is the total integrated line intensity, $`\sigma _i`$ is the uncertainty in observation $`i`$, and the summation is done over all independent observations $`N`$. The errors in the observed intensities are generally dominated by the calibation uncertainty of $``$20%. In some cases the line profiles were used to discriminate between models. For each new $`\dot{M}`$ assigned, a consistent envelope size ($`r_\mathrm{p}`$) was calculated based upon the results of Mamon et al. Mamon et al. (1988).

It is found that the intensity ratios between the different lines are sensitive to various parameters, and hence can be used to constrain them (Sec. 4.3). Therefore, in order to accurately determine $`h`$, and ultimately $`\dot{M}`$, multi-transition observations are needed. Observed radial brightness profiles provide additional important information. For a limited number of the sources in our sample there are enough observational constraints that a reliable $`h`$-parameter may be estimated. These sources were first modelled, see Table 3.

### 4.2 The h-parameter

The derived $`h`$-parameters, as a function of the adopted luminosity of the central source, are shown in Fig. 3. A division between low luminosity objects and those with higher luminosities is evident. For luminosities lower than 6000 L the median value of the estimated $`h`$-parameter is 0.2 (based on 9 sources), while for the higher luminosities it is 0.5 (12 sources) (the J-type stars, i.e., stars with a high content of 13CO, were not included in these estimates since 13CO line cooling is not treated in the standard model). The large scatter for the high luminosity sources is unfortunate, since these are generally the sources where the derived mass loss rates are particularly sensitive to the adopted value of $`h`$. For the majority of objects, where a reliable estimate of the $`h`$-parameter is not possible, we have adopted a value of $`h`$$`=`$0.2 for the low luminosity objects (below 6000 L) and $`h`$$`=`$0.5 for the high luminosity objects (above 6000 L).

The $`h`$-parameter is related to the dust-to-gas mass ratio, $`\mathrm{\Psi }`$ (see Eq. 7). If we assume that the properties of a dust grain are equal to those used to fit the line emission from CW~Leo, and that they are the same for all stars, an $`h`$-parameter of 0.2 means that $`\mathrm{\Psi }`$$`=`$2$`\times `$10-3. Groenewegen et al. Groenewegen et al. (1998b) have modelled the dust emission towards a large sample of carbon stars and find a value of $`\mathrm{\Psi }`$ which is fairly constant at 2.5$`\times `$10-3 for stars with luminosities up to about 7900 L, and which increases drastically with luminosity up to 0.01 and higher. This is in good agreement with our results and lends further support to our distinction between low- and high-luminosity stars as regards the value of the $`h`$-parameter. Furthermore, Hiriart & Kwan Hiriart and Kwan (2000) derive, for 17 of our sample stars, on the average $`\mathrm{\Psi }`$$``$7$`\times `$10-4. The main reason for this discrepancy is the difference in the adopted dust parameters. In fact, it can be shown that for $`h`$$`=`$0.2 we get essentially the same heating term (cf. Eq. 6) in the energy balance equation as Hiriart & Kwan Hiriart and Kwan (2000).

### 4.3 The mass loss rates

We have managed to model reasonably well 61 of the 69 sample sources using our standard model and fitting strategy defined above. The reduced $`\chi ^2`$ for these models is estimated from

$$\chi _{\mathrm{red}}^2=\frac{\chi _{\mathrm{min}}^2}{Np},$$
(13)
where $`\chi _{\mathrm{min}}^2`$ is obtained from Eq. 12 for the best fit model, and $`p`$ is the number of adjustable parameters (one or two) in the model. Typically, $`\chi _{\mathrm{red}}^2`$$``$1$``$3 assuming a calibration uncertainty of 20%. The derived mass loss rates and $`h`$-parameters are presented in Table 3, where also the velocities by which the CSEs expand, $`v_\mathrm{e}`$, and their spatial extents, $`r_\mathrm{p}`$, are given.

The derived mass loss rates presented in Table 3 span more than three orders of magnitude, 6.5$`\times `$10-9 M yr-1 to 1.5$`\times `$10-5 M yr-1, over which the physical conditions of the CSEs vary considerably, posing a challenge to the model. The intensities, overall line shapes, and when available radial brightness distributions, of the circumstellar CO lines produced by the radiative transfer model generally agree well with those observed. This is illustrated here for three of our sample stars LP~And (Fig. 4), R~For (Fig. 5), and W~Ori (Fig. 6). These carbon stars span a large range in mass loss rate and serve to illustrate the various physical conditions in these CSEs.

A full, detailed, error analysis of the estimated mass loss rates is not possible due to the relatively large number of free parameters entering the model. Instead, we will vary some of the more important parameters in order to illustrate the sensitivity of the model, and to be able to get a rough estimate of the errors involved in the mass loss rate estimates. The results of these sensitivity tests are shown in Table 5 for our three example stars. In addition, in Fig. 7 we present chi-square contour plots for these stars, produced by varying the mass loss rate and the $`h`$-parameter. This will give an estimate of the accuracy in the determination of these two adjustable parameters when all other parameters are held fixed (the size and shape of the CO envelope is however allowed to change in accordance with Eqs 10 and 11). In Fig. 7 we indicate various confidence levels with the contour at $`\chi _{\mathrm{min}}^2`$$`+`$2.3 marking the 68% confidence limit (i.e., the ”1$`\sigma `$” level for two adjustable parameters). The information contained in the shape of the line profiles (to which the $`\chi ^2`$ defined above is not very sensitive) have occasionally been used to further constrain the most probable parameter space.

LP~And is a high mass loss rate Mira variable where the excitation of 12CO is dominated by collisions. Consequently, the line intensities are very sensitive to the temperature structure, i.e., the $`h`$-parameter. It is interesting to note that the assumed envelope size $`r_\mathrm{p}`$ does not significanly affect the derived mass loss rate in the high mass loss rate regime. The reason for this being that the density is too low to excite the CO molecules effectively (the lines are sub-thermally excited) in the cool outer parts of the envelope (cf. Fig. 4), i.e., the emission is excitation limited. In the high mass loss rate regime the line intensities are also insensitive to the adopted mass loss rate since an increase in $`\dot{M}`$ leads to more cooling ($`C`$$``$$`\dot{M}`$ while $`H`$$``$$`\dot{M}^{0.5}`$), which compensates for the increase of molecular density. This ”saturation”-effect of the line intensities for high mass loss rates has been noted before in other models where the cooling by CO is treated in a self-consistent manner (Sahai 1990; Kastner 1992). Additional sensitivity tests in the case of a high mass loss rate object (CW~Leo) can be found in Groenewegen (1994).

W~Ori, on the other hand, is a low mass loss rate object with low optical depths in the lowest rotational transitions. In fact, some of the low lying transitions are inverted over parts of the CSE, i.e., maser action is taking place. In particular, the $`J`$$`=`$$`1`$$``$$`0`$ is inverted over the entire CSE producing a line profile that is clearly double peaked, i.e., this is due to maser action and not due to resolution effects. The observations of this transition are complicated by possible contamination by interstellar lines, and the signal-to-noise ratio is not particularly good, erasing any possible signs of maser action. As expected the radiation emitted by the central star plays an important role in the excitation. The transitions exhibiting population inversion are particularly sensitive (both the intensity and line profile) to the physical conditions prevailing in the CSE. In this case where radiative excitation is important the model becomes sensitive to the choice of inner radius, $`r_\mathrm{i}`$, as well as tubulent velocity, $`v_\mathrm{t}`$. Here, also the size of the envelope is important when determining the mass loss rate, i.e., the emission is photodissociation limited (at least for the lowest transitions). Generally, the line intensities scale roughly linearly with $`\dot{M}`$.

R~For presents an intermediate case between LP~And and W~Ori. It is also a typical sample star in the sense that the majority of the stars have envelope parameters that resemble those of R~For. From Table 5 and Fig. 7 we see that R~For shares properties of both the low and the high mass loss rate star. An intermediate mass loss rate model is generally more sensitive to the temperature structure than the radiation field, thus resembling the high mass loss rate objects. However, as in the low mass loss rate regime the line intensities scale roughly linearly with $`\dot{M}`$. We also note that the derived line intensities, especially the $`J`$$`=`$$`1`$$``$$`0`$ transition, start to become sensitive to the outer radius of the CO envelope. The exact mass loss rate where the transition from photodissociation- to excitation-limited emission occurs is hard to pinpoint, since it depends on the envelope characteristics. We have made tests which show that for typical envelope parameters the transition lies at $``$5$`\times `$10-7 Myr-1 for the $`J`$$`=`$$`1`$$``$$`0`$ transition. Higher $`J`$-lines sample hotter and denser gas, located closer to the star, and are therefore less sensitive to the choice of the envelope size.

In summary, we believe that within the adopted circumstellar model, the estimated mass loss rates are accurate to about $`\pm `$50% (neglecting errors introduced by the uncertain CO abundance and the distance estimates) when good observational contraints are available, but one should keep in mind that the causes of the uncertainty varies with the mass loss rate.

Another rough way to estimate the errors involved is to compare with the results obtained from other self-consistent models. When comparing our mass loss rate estimate for CW Leo to those obtained from detailed radiative transfer models (Kastner 1992; Crosas & Menten 1997; Groenewegen et al. 1998a; Skinner et al. 1999) we find a very good agreement, within 20 %, when adjustments for differences in $`f_0`$ and distance have been made. Kastner Kastner (1992) also modelled the high mass loss rate object RW LMi (a.k.a. CIT 6) and obtained a mass loss rate in excellent agreement with our estimate (when corrected for the difference in distance). Recently, Hiriart & Kwan Hiriart and Kwan (2000) presented a method to combine a circumstellar dust model with a CO model, keeping physical quantities consistent between the models. They present mass loss rate estimates for 17 of our sample stars. Based upon their Eq. 17 and Table 6 we derive hydrogen mass loss rates for these stars that are, on the average, 1.3 times larger than those derived from our radiative transfer analysis (for an individual star, however, the discrepancy can be as large as a factor of 3).

We also compare our derived mass loss rates with those obtained by Olofsson et al. (1993a) using a formula based on the work by Knapp & Morris Knapp and Morris (1985),

$$\dot{M}_{\mathrm{KM}}=5.7\times 10^{20}\frac{T_{\mathrm{mb}}v_\mathrm{e}^2D^2\theta _{\mathrm{mb}}^2}{s(J)f_0^{0.85}}\mathrm{M}_{\mathrm{}}\mathrm{yr}^1,$$
(14)
which relates the mass loss rate $`\dot{M}_{\mathrm{KM}}`$ to easily determined observables ($`T_{\mathrm{mb}}`$ is given in Kelvins, $`v_\mathrm{e}`$ in km s-1, $`D`$ in pc, and $`\theta _{\mathrm{mb}}`$ in arcseconds). In Eq. 14 $`s(J)`$ is a factor that depends on the transition ($`J`$) used. For $`J`$$`=`$$`1`$$``$$`0`$ $`s(J)`$$`=`$1 and for $`J`$$`=`$$`2`$$``$$`1`$ $`s(J)`$$`=`$0.5 (see the discussion in Olofsson et al. 1993a). When estimating the mass loss rates from Eq. 14 we have used the observed $`J`$$`=`$$`1`$$``$$`0`$ and $`J`$$`=`$$`2`$$``$$`1`$ line emission, and averaged the results. This mass loss rate was compared to that obtained from the radiative transfer analysis, Fig. 8. It is found that Eq. 14 significantly underestimates the mass loss rates when compared to those obtained from the radiative transfer analysis by, on the average, a factor of about four. However, considering the simplicity of Eq. 14 the agreement is remarkably good. There is also a trend that the discrepancy increases with lower mass loss rate (indeed Olofsson et al. suspected that their simple approach led to a systematic underestimate of the mass loss rate that became worse the lower the mass loss rate). In this context it should be noted that Eq. 14 was derived assuming a fixed CO envelope size of 3$`\times `$1017 cm, which is appropriate for high mass loss rate objects. In the high mass loss rate regime the line intensities are not very sensitive to the exact envelope size since the emission is excitation limited. This is, however, generally not the case for the optically bright carbon stars in our sample, where the line intensities are sensitive to the size of the CO envelope set by photodissociation. A correction for the envelope size \[see Neri et al (1998) and references therein\] would decrease the discrepancy seen in Fig. 8 for the stars losing mass at lower rates. Taken at face values the discrepancies between the mass loss rates derived here and those derived by Olofsson et al. (1993a) are smaller, but this is due to the fact that on the average Olofsson et al. used larger distances. This also means that the trend in the dust-to-gas mass ratio, decreasing with the mass loss rate, that Olofsson et al. reported was due to underestimated gas mass loss rates \[see also Hiriart & Kwan (2000)\]. Indeed, in Sec. 4.2 we find the opposite behaviour, i.e., indications of a dust-to-gas mass ratio that increases with the mass loss rate.

### 4.4 The CO envelope sizes

The size of the CO envelope is an important parameter, and the derived mass loss rate will depend on it to an extent that depends on the mass loss rate (see discussion in Sec. 4.3). Neri et al. Neri et al. (1998) presented IRAM PdB interferometer (PDBI) CO($`J`$$`=`$$`1`$$``$$`0`$) maps of varying quality for 9 of our sample stars. From their CLEANed brightness maps we have obtained the averaged radial brightness distributions. In addition, we have obtained two single-dish CO($`J`$$`=`$$`1`$$``$$`0`$) maps using the OSO 20 m telescope. In Fig. 9 we compare the observed radial brightness distributions to those calculated from our model, which uses the Mamon et al. Mamon et al. (1988) results to estimate the CO envelope size, for the parameters presented in Tab. 3. We find that the estimated CO envelope sizes are generally consistent with the observations. In four cases (V384~Per, RW~LMi, T~Dra, and S~Cep) the calculated radial brightness distributions fail to reproduce the observations. At least in the case of V384~Per and RW~LMi this appears to be due to deviations from a simple $`r^2`$ density law, in the sense that two distinct components are clearly visible in the brightness maps.

We are not able to distinguish between the different CO envelope sizes obtained by Mamon et al. Mamon et al. (1988) and Doty & Leung Doty and Leung (1998) for the high mass loss rate object CW~Leo. This is simply because the emission in this case is excitation limited with the external cool and tenuous parts producing very little emission.

### 4.5 The kinetic temperature structure

As was demonstrated in Sec. 4.3 the derived mass loss rate can be very sensitive to the kinetic temperature structure of the gas. Adiabatic cooling alone results in $`T(r)`$$``$$`r^\delta `$, where $`\delta `$$`=`$$``$4/3 (assuming $`\gamma `$$`=`$5/3), but from Figs. 4$``$6 it is evident that other cooling processes as well as heating mechanisms are important. If CO line radiation is the dominant cooling mechanism and $`n`$$`>`$$`n_{\mathrm{crit}}`$ one would expect $`T(r)`$$``$$`r^1`$ with a somewhat shallower decline if $`n`$$`<`$$`n_{\mathrm{crit}}`$ (see Doty & Leung 1997).

Generally, it is not possible to assign a simple power law to the radial behaviour of the temperature throughout the entire CSE. However, in the case of LP~And we find $`\delta `$$``$$``$0.9 for $`r`$$``$5$`\times `$1016 cm. Further out in the cool, tenuous, outer layers of this CSE, heating due to the photoelectric effect is important and raises the temperature. For R~For and W~Ori we estimate, over the regions where the observed emissions emanate, i.e., from $``$1$`\times `$1015 to $``$5$`\times `$1016 cm, $`\delta `$ to be $``$$``$0.7 for both sources. For comparison, Doty & Leung Doty and Leung (1997) derive $`\delta `$$``$$``$1.0 for the high mass loss rate object CW~Leo, while we obtain $`\delta `$$``$$``$0.9.

From the fact that $`C`$$``$$`\dot{M}`$ and $`H`$$``$$`\dot{M}^{0.5}`$ one would expect high mass loss rate objects to have significantly cooler envelopes than low mass loss rate objects (e.g., Jura et al. 1988). In our modelling we do not see this trend. The reason being that in addition $`H`$$``$$`h`$$`(Lv_\mathrm{e})^{3/2}`$, where all of the parameters have a mass loss rate dependence in the sense that their values are generally smaller for the low mass loss rate objects.

### 4.6 13CO line cooling

Carbon stars generally have 12C/13C-ratios in the range 30$``$70 (Lambert et al. 1986; Schöier & Olofsson 2000). This means that for the large majority of all stars in our sample line cooling by 13CO is insignificant, even when optical depth effects are taken into account. However, there exists a small population of carbon stars, the J-type stars, which have 12C/13C-ratios of $``$3. For these stars 13CO line cooling could be important.

To test this we have included 13CO line cooling for one of our sample stars, Y~CVn, a known J-type star. We have used the 13CO model as presented in Schöier & Olofsson Schöier and Olofsson (2000) to calculate the the contribution from this isotopomer to the total amount of cooling, and used this extra cooling term as input to the standard model. We find that, after we have iterated between the 12CO and 13CO models a few times, that 13CO line cooling is somewhat less effective as a coolant than expected, based on the isotope ratio (probably due to slightly different excitation conditions between the two molecules). The effect on the derived line intensities is small; the largest effect is on the $`J`$$`=`$$`3`$$``$$`2`$ line emission which is lowered by $``$10%. The changes in line intensity may be compensated for by increasing the $`h`$-parameter.

Moreover, Ryde et al. Ryde et al. (1999) present a model of the high mass loss rate object IRAS+15194-5115 using the radiative transfer code outlined here. The 12C/13C-ratio for this object was estimated to be 5.5, i.e., it resembles the J-type stars in this respect. It was found that 13CO line cooling in this case was approximately one third of that of 12CO (due to optical depth effects the 13CO line cooling is in a relative sense more important). The decrease in kinetic temperature due to this extra cooling was compensated for by the fact that the $`h`$-parameter (heating) had to be somewhat increased to produce the observed line intensities. Consequently, the derived mass loss rate was not affected by introducing the 13CO line cooling.

## 5 Deviations from the standard model

For some stars, where there are enough observational constraints, it is possible to establish deviations from the standard model.

### 5.1 Detached CSEs

In our sample there exist a number of stars that exhibit a 60 $`\mu `$m excess, placing them in the regions VIa and VIb in the IRAS colour-colour diagram (see Olofsson et al. 1993a), suggesting that their CSEs are abnormal in some way. Previous CO observations of five of these 60 $`\mu `$m excess stars (see Table 6) have revealed that their CSEs cannot have been produced by a single smooth expanding wind, but rather suggested a drastic modulation of the mass loss on a time scale as short as a thousand years, which produced detached CSEs (dCSEs; Olofsson et al. 1996b). Presently, these stars are losing matter at a significantly lower rate, forming what is refered to as an attached circumstellar envelope (aCSE).

To be able to put constraints on the mass loss and its variation with time, high spatial resolution observations are needed. Such observations have been carried out using the IRAM PdB interferometer, with an angular resolution of $``$1′′, for two of the stars in Table 6, U~Cam and TT~Cyg. It was found that the dCSE around U~Cam is fairly young while that around TT~Cyg is significantly older. Both dCSEs show a remarkable overall spherical symmetry, and they are geometrically very thin. The CO line emission towards these two sources, supplemented by a detailed modelling using the radiative transfer code outlined in this paper, were presented in Lindqvist et al. Lindqvist et al. (1999) (U~Cam) and Olofsson et al. (1998, 2000) (TT~Cyg).

In this paper we have attempted to model the remaining three stars with known detached shells (R~Scl, S~Sct, and U~Ant), based upon the results obtained from the high angular resolution observations described above. The results for both the aCSEs and dCSEs around these stars are presented in Table 6 (where also the results for TT~Cyg and U~Cam are included). Due to the limited observational constraints available, the dCSEs around S~Sct and U~Ant were modelled using the same geometrical thickness and kinetic temperature structure as obtained for TT~Cyg (Olofsson et al. 2000). For R~Scl we used the same geometrical thickness as adopted by Olofsson et al. 1996b and assumed a kinetic temperature decrasing linearly with radius from a value of 50 K at the inner shell radius. The derived intensities from the model generally agree well with those observed (Table 7). However, we note that in the case of R~Scl we only get about half of the observed $`J`$$`=`$$`3`$$``$$`2`$ line emission from the model, but here the aCSE could produce a significant amount of the observed emission. We find that the estimated shell masses, $`M_\mathrm{s}`$, are relatively similar, but there is a trend of increasing shell mass with age (the R~Scl shell mass may be overestimated due to the difficulty in separating the aCSE and dCSE emission), which may indicate that gas is being swept up.

The present mass loss rate was modelled for U~Ant and S~Sct using the standard model. The results are presented in Table 6 and Table 7. For R~Scl the emission from the aCSE and dCSE cannot be separated. Again, the predicted line intensities agree, with one exception, well with those observed, considering the uncertainties of the observations. For S~Sct the model of the aCSE fails to reproduce the IRAM CO($`J`$$`=`$$`2`$$``$$`1`$) observations, in contrast to the dCSE model where the fit to the observations is excellent. An explanation to this inconsistency could be that there are pointing problems in these observations (the IRAM observations are more sensitive to the pointing due to its smaller beam), which would mainly affect the emission observed from the aCSE. We find that the present mass loss rates of the stars with the oldest shells are low, and that the present mass loss rate of U~Cam (which has a much younger dCSE) is significantly higher.

The central stars are semi-regular (SR) or irregular variables (Lb), and they have low present mass loss rates. They appear to have been subject to drastic modulations of their mass loss rates over relatively short periods of time (a few thousand years). We estimate that the CO dCSE emission is observable for $``$104 years (Bergman et al. 1993). If this mass ejection is a repeatable phenomenon, and all our sample stars go through this phase, we estimate that the time scale between mass ejections is about 105 years. Therefore, the mass ejections may possibly be linked with the helium-shell flashes (the process that leads to nuclear-processed matter being transported to the surface of the star) predicted to occur regularly in AGB stars with time intervals of the order of 105 years for the relevant masses (e.g., Blöcker 1995). We note that all carbon stars that have been found to have clearly detached gas shells are members of our sample, and that no M-type star with a detached CO shell has been detected.

### 5.2 V~Hya and TX~Psc

The sucessful modelling of the vast majority of our sample stars, assuming a spherically symmetric CSE, suggests that axi- or non-symmetric mass loss is not a common phenomenon among these optically bright carbon stars. However, in the case of V~Hya and TX~Psc, stars with complex radio line profiles, the possibility of bipolar outflows has been considered (Heske et al. 1989; Kahane et al. 1996). Due to the complexity of these outflows, no radiative transfer analysis was attempted here.

Kahane et al. Kahane et al. (1996) suggest that V~Hya, in addition to the high velocity bipolar outflow, has a slowly expanding spherical component. Using a 2D radiative transfer code, based upon the Sobolev approximation, they derive a total mass loss rate of $``$1.5$`\times `$10-6 Myr-1 for this object. The authors further suggest that V~Hya is a transition object, between carbon stars and PNe, having just developed a highly axi-symmetric mass loss. The slowly expanding spherical component is then interpreted as the fossil CSE created during the AGB stage. Knapp et al. Knapp et al. (1999), however, argue, based upon the spectral type, period, colours, and lack of ionizing radiation, that V~Hya is still on the AGB. Moreover, they estimate the total mass loss rate of this object to be $``$4$`\times `$10-5 Myr-1. This, in addition to the small dynamical age of the envelope, suggests that V~Hya has entered its ”superwind” phase.

The observed line profiles of the CO line emission towards TX~Psc are somewhat peculiar (Heske et al. 1989; Olofsson et al. 1993a). As a consequence, the expansion velocity derived from the observed radio line profiles range from 7.5 to 12.2 km s-1 (see Table 2 and Olofsson et al. 1993a) depending on the transition and the telescope used. Heske et al. Heske et al. (1989), mapped the circumstellar CO line emission around this object, and interpreted the observed asymmetries as being produced either by a bipolar outflow or a highly clumped wind. No detailed modelling exist for the CSE around TX~Psc, and its evolutionary status remains uncertain. It might be in a state similar to that of V~Hya.

### 5.3 DR~Ser

DR~Ser has the highest 60 $`\mu `$m excess of all the sample stars, which places it in the IRAS two-colour diagram among objects having detached envelopes (see Table 6). Indeed, Kerschbaum Kerschbaum (1999) suggests that this might be a detached shell source based upon his modelling of the SED of this object. However, the CO single dish observations reveal no signs of double-peaked profiles (confusion with interstellar lines complicates the interpretation), which are indicative of a detached CO shell. Moreover, DR~Ser has a relatively low 12C/13C-ratio of 6 (Abia & Isern 1997). Thus, it resembles the J-type stars.

Attempting to model this source with a standard model of its CSE fails to explain both the observed $`J`$$`=`$$`1`$$``$$`0`$ and $`J`$$`=`$$`2`$$``$$`1`$ line intensities in that they come out too weak, by more than a factor of two. This is due to the extremely cold envelope, resulting from the adopted $`h`$-parameter of 0.2 (due to the luminosity being lower than 6000 L), in combination with a relatively high mass loss rate $``$2$`\times `$10-6 Myr-1 required to maximize the line intensities. Further increasing the mass loss rate will make the envelope even cooler due to the increase in CO line cooling. If one instead assigns an $`h`$-parameter of $``$1 it is possible to obtain a reasonable fit to the data using a mass loss rate of $``$5$`\times `$10-6 Myr-1. This estimate must be regarded as highly uncertain.

## 6 Discussion

### 6.1 The mass loss rate distribution

Fig. 10 shows the distribution of the mass loss rates obtained from the radiative transfer model for all sources (also shown is the distribution for the stars within 500 pc of the Sun), as well as divided into variability groups (note that objects with known detached shells have not been included in this analysis). Since the statistics of the SRa stars are poor we have chosen to group them together with the SRb stars. The mass loss rate distribution for all stars is sharply peaked around the median mass loss rate of 2.8$`\times `$10-7 M yr-1 (the median mass loss rate for the stars within 500 pc of the Sun is 1.6$`\times `$10-7 M yr-1). It is interesting to note that this is roughly the rate at which the core mass is expected to grow due to nuclear burning (Schönberner 1983). Mira variables generally have larger mass loss rates than other variability types, while irregular variables and semiregulars appear to have very similar mass loss rate characteristics.

Based on the mass loss rate distribution for the sample of stars within 500 pc, which we believe to be close to complete and for which we have detected all sources, we are able to draw some general conclusions. The sharp decline at mass loss rates below $``$5$`\times `$10-8 Myr-1 is very likely real. Netzer & Elitzur Netzer and Elitzur (1993) estimate that a mass loss rate in excess of $``$10-7 M yr-1 is required to get a dust-driven wind. The exact limit is, however, sensitive to the adopted stellar and dust parameters. The drastically decreasing number of high mass loss rate objects is also real, although our selection criterion bias against these objects, and can be explained by the fact that carbon stars of the type discussed here, i.e., mainly low mass ones (Claussen et al. 1987), only for a limited time, or possibly never, reach high mass loss rates.

### 6.2 Mass loss and envelope kinematics

The $`v_\mathrm{e}`$-distribution for this sample of stars has already been shown and discussed by Olofsson et al. Olofsson et al. (1993a). However, a new comparison between the mass loss characteristics $`\dot{M}`$ and $`v_\mathrm{e}`$ is warranted considering the more reliable mass loss rates obtained in this paper, Fig. 11. We find a clear trend that $`v_\mathrm{e}`$ increases with $`\dot{M}`$, $`v_\mathrm{e}`$$``$$`\dot{M}^{0.40}`$ (with a correlation coefficient of 0.78). In comparison, Olofsson et al. Olofsson et al. (1993a) derived $`v_\mathrm{e}`$$``$$`\dot{M}^{0.53}`$, but their mass loss rates were calibrated using Eq. 14, which includes a $`v_\mathrm{e}^2`$-dependence. The correlation is much tighter than that obtained by Olofsson el al. Olofsson et al. (1993a), which is reassuring. Thus, the mass loss mechanism operates such that mass loss rate and expansion velocity increase together. However, the scatter appears larger than the uncertainties in the estimates, and so the mechanism also produces widely different mass loss rates for a given expansion velocity. At high mass loss rates the scatter is larger with a possible divison into objects with high mass loss rates but only moderately high velocities, and objects with moderately high mass loss rates but high velocities. We note also that for the detached shell sources the mass loss rates and expansion velocities for the two mass loss epochs follow the general trend.

Habing et al. Habing et al. (1994) studied the momentum transfer from the photons, via the dust, to the gas in the CSEs of AGB stars. They found that for mass loss rates close to the minimum mass loss rate for a dust-driven wind, $``$3$`\times `$10-8 M yr-1, $`v_\mathrm{e}`$ increases linearly with $`\dot{M}`$, while for mass loss rates above $``$10-6 M yr-1 the dependence weakens considerably. At 10-5 M yr-1 they found $`v_\mathrm{e}`$$``$$`\dot{M}^{0.04}`$. Habing et al. attribute the increase in $`v_\mathrm{e}`$ with mass loss rate (for low mass loss rates) to a higher efficiency in the coupling between gas and dust. This can explain the observed behaviour in Fig. 11 for mass loss rates below about 10-6 M yr-1, but for the higher mass loss rates it appears that a dependence of $`v_\mathrm{e}`$ (and $`\dot{M}`$) on luminosity gives the most reasonable explanation, see Sect. 6.3.

### 6.3 Dependence on stellar properties

In Fig. 12 we plot the circumstellar characteristics, $`\dot{M}`$ and $`v_\mathrm{e}`$, against the stellar characteristics luminosity ($`L`$), period ($`P`$), effective temperature ($`T_{\mathrm{eff}}`$), and the photospheric C/O-ratio. The effective temperatures and C/O-ratios used in this analysis are those presented in Olofsson et al. Olofsson et al. (1993b). We have looked for dependences assuming that the ordinate scales as the abscissa to the power of $`\xi `$, and a normal correlation coefficient $`r`$ was calculated as an estimate of the quality of the fit, see Fig. 12. The uncertainties in the estimated quantities are of the order: $`\pm `$50% ($`\dot{M}`$), $`\pm `$2 km s-1 ($`v_\mathrm{e}`$), $`\pm `$10% ($`P`$; but some periods may be poorly determined), a factor of 2 ($`L`$), and $`\pm `$200 K ($`T_{\mathrm{eff}}`$).

Clearly, both the mass loss rate and the expansion velocity increase with the pulsation period of the star. There is a weaker trend with the luminosity of the star. The latter dependence is not completely independent of the former, since some of the luminosities are estimated from a $`PL`$ relation. The existence of a $`PL`$ relation is usually attributed to a distribution in mass (Jones et al. 1994), i.e., the higher the mass the longer the period and the higher the luminosity. When looked at in detail, at least the apparent mass loss rate dependence on period may be attributed to a change from semiregular pulsation at short periods to regular pulsations at longer periods. We may therefore infer that for the mass lass rate it is not clear whether it is the regularity of the pulsation or the luminosity that causes the increase with period. For the expansion velocity the increase with period may be a combined effect of increasing mass loss rate with period (for low mass loss rates) and an increase in luminosity for the longer periods, e.g., Habing et al. Habing et al. (1994) derive $`v_\mathrm{e}`$$``$$`L^{0.35}`$ in their dust-driven wind model. In addition, there is evidence of a weak trend in the sense that higher mass loss rate objects have lower effective temperatures. These correlations are all consistent with a dust-driven wind, where the pulsation may play an important role. There appears to be no correlation between the mass loss rate and the C/O-ratio, which is surprising considering that the dust-to-gas mass ratio in a C-rich CSE should be sensitively dependent on this.

The present mass loss rates for the detached shell sources are typical for their periods but somewhat low for their luminosities. However, during the formation of the dCSEs the mass loss rate must have been atypically high for the present periods and luminosities. We also find that the expansion velocities of the dCSEs lie at the very high end of expansion velocities found for other stars with the same luminosity, indicating that these stars had higher luminosities during the shell ejection.

Two other stars, SZ~Car and WZ~Cas, stand out in these plots. The properties of SZ~Car resembles the stars with known dCSEs, but the CO line profiles give no indication of a detached shell. WZ~Cas has by far the lowest estimated mass loss rate in the sample, as well as the lowest C/O-ratio, 1.01, which classifies it as an SC-star. It is also a Li-rich J-star (Abia & Isern 1997).

### 6.4 Enrichment of the ISM

Carbon stars on the AGB are important in returning processed gas to the interstellar medium (ISM). The total mass loss rate of carbon stars in the Galaxy is obtained from

$$\dot{M}_{\mathrm{Gal}}=_0^{\mathrm{}}2\pi R\mathrm{\Sigma }<\dot{M}>𝑑R$$
(15)
where $`\mathrm{\Sigma }`$ is the surface density of carbon stars, $`<`$$`\dot{M}`$$`>`$ is the mean mass loss rate of the carbon stars. According to Guglielmo et al. Guglielmo et al. (1998) the infrared carbon stars have a roughly constant surface density out to the galactocentric distance of the Sun, after which it follows an exponential decline with a scale length of $``$2.5 kpc. We approximate this by assuming a constant value, $`\mathrm{\Sigma }_0`$ (the surface density in the solar neighborhood), to 10 kpc and zero beyond this galactocentric distance.

Based on the complete sample, i.e., stars within 500 pc from the Sun, we estimate $`\mathrm{\Sigma }_0`$$`<`$$`\dot{M}`$$`>`$ to be $``$1.7$`\times `$10-10 M yr-1 pc-2, where we have included also helium. This estimate is very sensitive to the number of high mass loss rate objects found within 500 pc, e.g., the high mass loss rate object CW~Leo contributes almost half of our estimate of the total mass returned to the ISM. Our estimate of the rate at which matter is returned to the ISM by carbon stars is consistent with previous estimates (considering the large uncertainties), e.g., Knapp & Morris Knapp and Morris (1985) and Jura & Kleinmann Jura and Kleinmann (1989) derive values of $``$2$`\times `$10-10 M yr-1 pc-2 and $``$1.5$`\times `$10-10 M yr-1 pc-2, respectively. Using Eq. 15 the annual return of matter to the ISM in the Galaxy is estimated to be $``$0.05 M for the carbon stars considered here. It is quite possible that a larger mass return from carbon stars is obtained during their final evolution on the AGB, a so called superwind phase, as is indicated by the estimated values for infrared carbon stars, $``$0.5 M yr-1 (Epchtein et al. 1990; Guglielmo et al. 1993 ), and PNe, $``$0.3 M yr-1 (Gustafsson et al. 1999). This confirms the importance of carbon stars for the cosmic gas cycle in galaxies. By comparison, high mass stars are estimated to contribute about 0.04 M yr-1 (Gustafsson et al. 1999).

We estimate the contribution from carbon stars to the carbon enrichment of the ISM to be $``$0.5$`\times `$10-4 M yr-1 during the AGB stage (if the subsequent superwind phase is included this value may increase to $``$5$`\times `$10-4 M yr-1). This corroborates the conclusions by Gustafsson et al. Gustafsson et al. (1999) that ‘normal’ carbon stars are not important for the total carbon budget of the Galaxy. According to these authors the main contributors should instead be high mass stars in the Wolf-Rayet stage, annually supplying the Galaxy with $``$0.01 M of carbon. This is roughly what is required to produce the $``$10-3$`\times `$1011 M of carbon present in the Milky Way over a period of 1010 years.

## 7 Conclusions

We have developed a detailed radiative transfer code to model circumstellar molecular line emission. The code also solves for the energy balance equation of the gas. It is found that the mass loss rate determination for low mass loss rate objects depends crucially on a number of assumptions in the CSE model, except on the dust properties, since the CO molecules are radiatively excited. On the other hand, the mass loss rate estimate for a high mass loss rate object depends essentially only on the temperature structure, and hence the uncertain dust parameters. In addition, for such an object the CO emission saturates, and becomes less useful as a mass loss rate measure. We also find that different lines respond differently to changes in the various parameters. Therefore, a reliable mass loss rate determination requires, in addition to a detailed radiative transfer analysis, good observational constraints in the form of multi-transition observations and radial brightness distributions.

This model has been applied to CO radio line observations of a large sample of optically bright N-type and J-type carbon stars on the AGB (69 objects). The sample is reasonably complete out to about 1 kpc, and all stars (41) within $``$500 pc of the Sun have been detected in circumstellar CO line emission. The derived mass loss rates span almost four orders of magnitude, $``$5$`\times `$10-9 Myr-1 to $``$2$`\times `$10-5 Myr-1, over which the physical conditions of the CSEs vary considerably. The fact that the model can be succesfully applied over such a wide range of environments gives us confidence in the results, even though we are aware of the fact that some assumptions are poorly constrained. The large majority of the stars have mass loss rates in a narrow range centered at $``$3$`\times `$10-7 M yr-1, and it appears that very few AGB carbon stars ($``$5%) lose matter at a rate less than $``$5$`\times `$10-8 M yr-1.

We find that the mass loss rate and the gas expansion velocity are relatively well correlated, but the scatter is large enough that we may conclude that the mass loss mechanism is able to produce a wide range of mass loss rates for a given expansion velocity. The mass loss rate is also well correlated with the pulsational period of the star, correlated with the stellar luminosity, and there is a weak trend with the stellar effective temperature, in the sense that the cooler stars tend to have higher mass loss rates. Also the gas expansion velocity is positively correlated with the period and the luminosity. We conclude that the mass loss rate increases with increased regular pulsation and/or luminosity, and that the expansion velocity increases with mass loss rate (for low mass loss rates) and luminosity. The observed trends are all supporting the common concensus that these winds are driven by radiation pressure on dust grains, and that pulsation may play an important role. Somewhat surprising there appears to be no dependence on the stellar C/O-ratio.

Our standard CSE model, assuming a single smooth expanding wind produced by a continuous mass loss, fails to reproduce the observational data for about 10% of the sample stars. Most notable among these are five stars with detached CSEs, presumably formed during a period of highly increased mass loss, and low present mass loss rates. We have found indications that the present mass loss rate is significantly higher for the star with the youngest dCSE, and that the shell mass may increase with shell age. Since our sample is reasonably complete, we can estimate that the time scale between mass ejections is about 105 years, if it is a repeatable phenomenon. An association with He-shell flashes is favoured. These objects have been the targets of a number of extensive, observational and theoretical studies.

For some of our sample stars there exist enough observational constraints to determine a combined dust parameter from the kinetic temperature structure. This result can be interpreted as a dust-to-gas mass ratio which increases with mass loss rate, but changes in the dust properties may also play a role. We also find that this means that the gas kinetic temperature in a carbon star CSE depends only weakly on the mass loss rate.

The size of the CO envelope is an important parameter in the mass loss rate determination, at least for the low mass loss rate objects. We have used published radial CO($`J`$$`=`$$`1`$$``$$`0`$) brightness distributions, in combination with the radiative transfer code, to show that the CO photodissociation calculation by Mamon et al. Mamon et al. (1988) gives reasonably accurate results. In a few cases the observed radial brightness distributions are clearly different than the model results, suggesting deviations from a simple $`r^2`$ density law, and hence time-variable mass loss.

We estimate that carbon stars, of the type studied here, return on the order of 0.05 M yr-1 of gas to the ISM making them marginally important for the gas cycle in galaxies. More extreme carbon stars may contribute an order of magnitude more. However, they are probably not important as regards the origin of carbon.

###### Acknowledgements.

We are grateful to Dr. F. Kerschbaum for generously providing estimates of some of the input parameters to the CO modelling, to Dr. R. Liseau for useful comments, and to an anonymous referee for constructive criticism which lead to an improved paper. Financial support from the Swedish Natural Science Research Council (NFR) is gratefully acknowledged.