K. M. Menten¹ - M. J. Reid² - E. Krügel¹ - M. J. Claussen³ - R. Sahai⁴

1 - Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany
2 - Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA
3 - National Radio Astronomy Obsrvatory, Array Operations Center, PO Box O, Socorro, NM 87801, USA
4 - Jet Propulsion Laboratory, MS 183-900, 4800 Oak Grove Drive, Pasadena, CA 91109, USA

Abstract
We describe Very Large Array observations of the extreme carbon star IRC+10216 at 8.4, 14.9, and 22.5 GHz made over a two year period. We find possible variability correlated with the infrared phase and a cm- to sub-millimeter wavelength spectral index very close to 2. The variability, observed flux densities, and upper limit on the size are consistent with the emission arising from the stellar photosphere or a slightly larger radio photosphere.

Key words: stars: carbon - stars: AGB and post-AGB - radio continuum: stars - techniques: interferometic

1 Introduction

The "extreme'' carbon star IRC+10216 (CW Leonis) was discovered in the early days of modern infrared astronomy (Becklin et al. 1969). The star is one of the brightest near- and mid-infrared sources (LeBertre 1987 and references therein). IRC+10216 has an extremely rich molecular spectrum arising from a dense envelope created by a powerful mass outflow (Cernicharo et al. 2000). These attributes, together with its proximity (distance of $\approx$ 130 pc) and its prodigious mass-loss rate ( $1.5{-}3.3\times10^{-5}$ $M_{\odot}$ yr^-1), make IRC+10216 a "keystone'' object for which many aspects of (carbon-star) asymptotic giant branch (AGB) evolution can be studied in detail. IRC+10216's properties are summarized in Table 1.

IRC+10216 has been detected at wavelengths of 1.5 and 2 cm (Sahai et al. 1989; Drake et al. 1991). From Very Large Array (VLA) observations it should be possible to measure its size, position (which is extremely difficult at other wavelengths owing to absorption and scattering), proper motion, and the spectral index of its radio emission. Very Large Array data can strongly constrain our information on the star's nature.

This paper has the following structure: in Sect. 2, we describe six-epoch, three-wavelength VLA observations. These, and an additional dataset (published by Drake et al. 1991), were retrieved from the VLA archive. The properties of the continuum emission are discussed in Sect. 3. In particular, we address the variability, spectral index, absolute position, and source morphology. In Sect. 4 we synthesize these results to reach conclusions on the nature of the radio emission.

2 VLA observations and data reduction

2.1 1991-1993 observations

The data reported here were taken by two of us (RS & MC) with the NRAO Very Large Array (VLA). We analyzed observations made between 1991 June 1 and 1993 June 9, spaced typically by 3 to 4 months. Observations were made at 8.4, 14.9, and 22.5 GHz (X-, U-, and K-band, respectively). Details are given in Table 2. A typical observing run consisted two $\sim$ 20-28 min long X- and U-band observations, preceded and followed by observations of the calibrator 0953+254. Because of the generally short coherence times at the highest frequency, the K-band scans were shorter ( $\approx$ 15 min duration) and more scans were made in this band. At the end of each run, we observed the absolute flux density calibrator 3C 286.

The flux density of 0953+254 was determined by comparing its visibility amplitude with that of 3C 286. The elevations of 3C 286 were comparable to that of 0953+254, minimizing systematic amplitude calibration errors due to elevation-dependent antenna gain. By comparing the variations of the flux densities for several scans of 0953+254, we estimate our flux density scale is accurate to within 10% at X-band, 15% at U-band, and 20% at K-band. Data from one epoch (1992 October 2) suffered serious instrumental difficulties and were discarded. Table 2 contains the measured flux densities of 0953+254. Complex gain correction factors were interpolated from the amplitude and phase of 0953+254 to the times of the IRC+10216 data.

We made maps using the AIPS task IMAGR and fit elliptical Gaussian models to the images. The integrated flux densities are listed in Table 2. Because IRC+10216 is weak, self calibration was not possible, and the "fast-switching'' technique, routinely employed now, had not been developed at the time of our observations. Nevertheless, during one epoch (1993 January 26) when the array was in the A-configuration, superb phase coherence persisted throughout the observations. We use those data to establish the source position (in Sect. 3.3).

2.2 Other VLA data

VLA A-configuration observations of IRC+10216 have been reported by Drake et al. (1991). We re-analyzed their 1987 June 2 and 3 U-band data to compare with our U-band data. Fitting an elliptical Gaussian using the AIPS task JMFIT, we find the source to be nearly unresolved, with upper limits on the deconvolved size of 95 milliarcesconds (mas) and a peak brightness and integrated flux density of $1.29\pm0.15$ mJy beam^-1 and $1.40\pm0.27$ mJy, respectively. The latter value can be compared with our observed range of integrated U-band flux densities of 1.77 to 2.67 mJy (see Table 2). Since we used identical calibration procedures for both datasets and the same absolute flux calibrator, this discrepancy is probably real. Since the 1987 data employed a secondary gain calibrator nearer to IRC+10216 than used for our data, one would expect that the 1987 data should yield a greater flux density than the 1993 data, if random atmospheric phase fluctuations degrade the imaging.

In order to estimate of the absolute position accuracy, we made images from the the June 2 and 3 data separately. The positions of the fitted Gaussians differed by 18 mas in right ascension and 10 mas in declination, testifying to the excellent phase coherence during these observations. Assuming that large-scale atmospheric effects, which usually limit astrometric accuracy, were independent for these two days, this indicates an uncertainty of $\approx$ 20 mas in each coordinate.

3 Properties of the radio emission

3.1 Variability

$\begin{figure} \par\includegraphics[width=7cm,angle=0]{4463fig1.eps}\end{figure}$

Figure 1: Radio flux densities of IRC+10216 from 1991 June 1 to 1993 June 8 (second panel from top) compared to that of the calibrator 0953+254 (top panel). For both sources the X-, U-, and K-band integrated flux densities from Gaussian fits are plotted bottom to top. Note the different flux density scales. Error bars represent the quadratic sum of formal fitting errors and estimates of the errors due to absolute flux density calibration (10, 15, and 20% for the X-, U-, and K-band data points). In the IRC+10216 panel, the solid lines represent the results of the model calculations described in the text. The third to fifth panels from top show K-, U- and X-band residuals in percent-for each epoch the flux density predicted by the model was subtracted from the measured value and divided by the model value. X-, U-, and K-band data points are represented as crosses, squares, and circles, respectively. In all panels the same time range is shown. The abscissa gives Julian day number in the top panel and infrared (L'-band) stellar phase, $\Phi$ , in the others. The date of maximum infrared light ( $\Phi =0$ ) was calculated from the value determined by LeBertre 1992 (JD 2447483), using his period of 649 d. The lowermost panel represents the L'-band ( $3.76~\mu$ m) light curve converted to kJy units (see Fig. 2k of Le Bertre 1992).

Open with DEXTER

A weighted least-squares fit to the data yielded the following results: $S(\nu_0) = 1.02\pm0.14$ mJy, $\Delta S(\nu_0) = 0.25\pm0.09$ mJy, $P = 535\pm51$ d, $t_{\rm max} = ({\rm JD}~2440000+) 8841\pm24$ d, and $\alpha = 1.88\pm0.14$ . This fit has a reduced $\chi^2$ per degree of freedom of 1.6, and the parameter uncertainties are the formal $1\sigma$ values scaled by $\sqrt{1.6}$ . Thus, IRC+10216 appears to display flux density variations of $\pm$ 25% at a significance level of nearly $3\sigma$ . The variability will be further discussed in the Appendix.

IRC+10216 shows strong periodic variability at infrared wavelengths (Le Bertre 1992; Dyck et al. 1991). Dyck et al. determine a period of 638 d at (infrared) K-band, while Le Bertre's measurements in the J-, H-, K-, L'-, and M-bands yield values between 636 and 670 d. If we assume the infrared period has a 5% uncertainty, then the cm-wave period of $535\pm51$ d differs from the infrared period by about $2\sigma$ . The cm-wave maximum occurs at $8841\pm24$ , whereas the infrared maximum (extrapolated by about 1200 d) would be predicted to occur near 8770 with an estimated uncertainty of about $\pm$ 50. Thus, the cm-wave minus infrared maxima differ by $71\pm71$ d and are consistent within their joint uncertainty.

3.2 Spectral index

As discussed above, the radio emission from IRC+10216 is probably variable. Therefore, we calculate the X-U-K-band spectral index (SI) for each epoch and find values between 1.6 and 2.2, with uncertainties of up to 20%. These daily SIs are consistent with a constant value, and a variance-weighted average of these indices is $1.92\pm 0.11$ . We also listed in Table 2 variance-weighted, time-averaged flux densities for each band. Using these to determine a SI yields $1.91\pm 0.11$ . The close correspondence between the daily and time averaged SIs suggests that the spectral index is not strongly variable.

The cm-wave spectral index obtained by least-squares fitting of all of the data simultaneously with the model described above is $1.88\pm0.14$ . This value is consistent with those observed for Mira variables and with the value of 1.86 for a model of Miras Reid & Menten (1997a). However, the cm-wave spectral index is also consistent with a value of 2.0 for optically-thick black-body emission.

At millimeter- and sub-millimeter wavelengths, flux densities measured with single-dish telescopes (with >10'' FWHM beams) contain contributions from the star as well as from the extended dust envelope. If the measurements are made with wideband bolometer detectors, circumstellar molecular emission may also contribute $\approx$ 30% to the total flux density (Groenewegen 1997; Walmsley et al. 1991). Whereas emission by dust dominates in the submillimeter region on large scales (see Fig. 4), it is negligible at cm-wavelengths. Based on the model discussed in the Appendix, we expect that the dust envelope contributes only $\sim$ 0.1 mJy to the total flux density of $\sim$ 2 mJy in our VLA K-band measurements.

A few interferometric observations exist that should deliver accurate measurements of compact emission free of dust contamination. Also by selecting portions of the band largely free of molecular lines, nearly pure continuum emission can be measured. Lucas & Guélin (1999), using the IRAM Plateau de Bure Interferometer (PdBI), find a "point source'' at 89 and 242 GHz with flux densities of $65\pm7$ and $487\pm70$ mJy, respectively. (Note these uncertainties are much larger than the formal errors given in the cited paper; R. Lucas personal communication.) These authors also detect "extended'' emission with flux densities of 12 and 57% of the point source values at 89 and 242 GHz, respectively, which they ascribe to dust emission from the inner envelope. Young et al. (2004) used the Smithsonian Sub-Millimeter Array (Ho et al. 2004) to image the 680 GHz continuum emission of IRC+10216. They also find a "compact, unresolved component'' in their $\approx$ 2'' FWHM beam with a flux density of $3.9\pm1.2$ Jy. The flux densities discussed above are plotted in Fig. 2. A least-squares fit to all the data (cm through sub-millimeter) delivers an SI of $1.96\pm 0.04$ .

$\begin{figure} \par\includegraphics[width=7.5cm,clip]{4463fig2.eps} \end{figure}$

Figure 2: Interferometric flux densities of IRC+10216 from cm to sub-mm wavelengths. The X-, U-, and K-band cm-wave data points are the weighted averages of data from six epochs (see Table 1). These error bars are quadratic sums of formal uncertainties and absolute calibration errors, which are assumed to be 10, 15, and 20% for the X-, U-, and K-band data, respectively. The error bars of the mm and submm wavelength data are discussed in the text. The straight line has a slope of 2.0, the blackbody value. A least squares fit to all the shown data delivers a slope of $1.96\pm 0.04$ .

Open with DEXTER

3.3 Absolute position of the radio emission distribution

During our 1993 January 26 observations, the VLA was in its most extended configuration and the phase coherence was excellent. Data from that epoch can be used to obtain meaningful images of the emission distribution and to make an astrometric position determination. To check the validity of the size information and the position accuracy, we determined, for the K-band data, phase and amplitude calibration of only the first scan of 0953+254 and then applied extrapolations of these solutions to subsequent scans of that source, which were spaced by $\approx$ 16 min. As expected, the phase coherence degraded with elapsed time between the first and subsequent scans. However, in maps made from the uv-data of the second scan, taken 16 min after the first, the source showed <3% amplitude degradation, a (formal) source size <31 mas and an offset of <5 mas in both coordinates from the nominal position. Since IRC+10216 scans were placed between 0953+254 scans, we conclude that some, but probably little, blurring occurs due to deficient phase calibration. The test described above does not account for phase calibration errors caused by the fact that 0953+254 is > $12\deg$ away from IRC+10216. To take that into account, we increase our value for the systematic position uncertainty to 20 mas.

The position uncertainty may be dominated by statistical errors due to the modest signal to noise ratio of our continuum images. JMFIT delivered positions for the X-, U-, and K-band images differing by maximally 20 mas from each other, with formal errors of between 5 and 16 mas. The total uncertainty, obtained by quadratically adding the systematic and statistical uncertainty is 15 mas in each coordinate. The variance-weighted mean position and its uncertainty is given in Table 1. The position derived from the data of Drake et al. (1991) is $(\alpha,\delta)_{{\rm J}2000}= 09^{\rm h}47^{\rm m}57\hbox{$.\!\!^{\rm s}$ }382, +13^{\circ}16'43\hbox{$.\!\!^{\prime\prime}$ }61,$ for which we estimate an uncertainty of 20 mas in each coordinate (see Sect. 2.2). The latter position (measured on 1987 June 2.5) is offset from the 1993 January 26 position by (-146, -20) mas in the (East, North) directions. Taking these offsets at face value, we derive proper motion components of IRC+10216 of ( $+26\pm6$ , $+4\pm6$ ) mas yr^-1, respectively. A motion of that magnitude is plausible. For example, Mira (o Ceti), which is at a similar distance as IRC+10216, has a larger proper motion of (+33,-239) mas yr^-1 measured by Hipparcos (Perryman et al. 1997).

The proper motion derived above is consistent with the upper limit of 30 mas yr^-1 that Becklin et al. (1969) derive from a comparison of a Palomar Sky Survey plate taken in 1954 with a plate they took in 1969.

3.4 Morphology

While an elongated source would be tantalizing, we have greater confidence in the 1987 than the 1993 data. The separation between IRC+10216 and 0952+179, the phase calibrator used by Drake et al., is only $1\hbox{$.\!\!^\circ$ }7$ in right ascension and $4\hbox{$.\!\!^\circ$ }4$ in declination. In contrast, the phase calibrator used in 1993 is displaced by $2\hbox{$.\!\!^\circ$ }1$ in right ascension and a much larger distance, $-12\hbox{$.\!\!^\circ$ }2$ , in declination from IRC+10216. While it is difficult to tell conclusively, we believe that the source size observed in 1993 January may be affected by uncorrected residual phase errors stemming from the large angular separation of IRC+10216 and the calibrator.

Groenewegen (1996) takes the elongation found by Drake et al. (1991) at face value and argues for asphericity. We note that on scales of a few stellar radii, IRC+10216's morphology is highly asymmetric and changes with time. Multi-epoch 2.2 $\mu$ m speckle interferometry between 1995 and 1998 shows an asymmetric clumpy shell with two dominant components, A and B, whose separation grows from 191 to 265 milliarcseconds (mas) between 1995 and 1998 (Osterbart et al. 2000). 2.2 $\mu$ m Keck I aperture masking interferometry between 1997 and 1999 confirm the observed morphology (Tuthill et al. 2000). While both studies find motions in the material, these are inconsistent with simple spherical expansion and a bipolar flow. Osterbart et al. (2000) suggest that the star itself is identical to, or located near, component B. On the other hand, (Richichi et al. 2003), modeling their infrared lunar occultation data, favor an identification of the star with component A.

Future VLA observations outlined in Sect. 5, together with absolute infrared astrometry, may unambiguously establish the star's position relative to the clumpy components of its envelope. Also, future VLA array observations at 1.3 and 0.7 cm can definitively settle the question of IRC+10216's radio source size and morphology.

4 The nature of the radio emission

Groenewegen (1997) performed model calculations of IRC+10216's infrared through radio continuum emission. His models and assumptions under-predict our measured radio flux densities by a factor of $\approx$ 1.8. A nearer distance of $\approx$ 100 pc (instead of the 135 pc assumed) would provide a better fit to the data. For X- and U-band, his models predict the emission to arise solely from the stellar photosphere, and at K-band the dust contribution from the whole envelope would add about 25% to the photospheric emission. Assuming an effective temperature of 2000 K, a photospheric diameter of 70.2 mas is derived, which is consistent with our upper limit. As shown in the Appendix, assuming these values for temperature and diameter, IRC+10216's spectral energy distribution can be modeled successfully.

A basic question regarding the cm-wave emission is whether it comes from the stellar photosphere or a somewhat larger radio photosphere. The radio emission from Mira variables comes from a radio photosphere, which is controlled by H^- free-free opacity and is approximately twice the size of the stellar photosphere (defined in line-free regions of the IR spectrum) (Reid & Menten 1997a). If a size (or stricter upper limit) could be measured with the VLA, that would directly yield a brightness temperature. A temperature below 2000 K would be expected for a radio photosphere, whereas a temperature in excess of 2000 would be expected for the stellar photosphere.

The amplitude of the radio variations of IRC+10216 of $\pm$ 25% is larger than for Mira variables, where variations of < $\pm$ 15% are observed by Reid & Menten (1997a). This is in contrast to variations of $\pm$ 40% (see appendix) that would be expected from the stellar photosphere of IRC+10216, based on IR variations of nearly 2 mag. For Miras, simple models of radio photospheres with shocks propagating outward with speeds less than $\approx$ 7 km s^-1 can explain the absence of strong radio variability. However, shocks propagating outward with speeds of $\approx$ 10-15 km s^-1 can produce variation of the magnitude we find in IRC+10216 (Reid & Menten 1997b).

5 Conclusions and outlook

Using the VLA, observations of IRC+10216 at 8.4, 14.9, and 22.5 GHz were made over a three year period. We find probable variability and a spectral index of about 1.9, i.e. near that of optically thick blackbody emission, but consistent with a radio photosphere as observed in Mira variables.

Future VLA A-configurations observations at Q-band (7 mm wavelength) with a resolution of 35 mas should definitely resolve the emission and yield a measurement of the brightness temperature. This could discriminate between a stellar and radio photosphere. Such observations may even lead to the detection of surface structure, as observed for $\alpha$ Orionis by Lim et al. (1998). These observations will benefit greatly by employing the "fast-switching'' technique.

Appendix: Modeling IRC+10216's spectral energy distribution

We add to the existing radiative transfer calculations another one, on a much more modest level, which has the goal of clarifying the relative contributions of the star and its dusty envelope to the long wavelength emission. In light of the huge amount of observational data, simplicity has, despite its obvious shortcomings, the advantage of transparency and of not being entangled by occasionally secondary details. The numerical code is based on a standard ray tracing method assuming radial symmetry. It can handle an arbitrary number of grain types, includes scattering and is described in Sect. 13.2 of Krügel (2003).

Our model parameters are as follows: The star has at maximum a luminosity $L = 1.1\times 10^4~L_\odot$ and emits as a blackbody of temperature T_*=2000 K. The density distribution of the dust is radially symmetric with $\rho(r) \propto r^{-2}$ . The outer radius of the circumstellar envelope equals $3\times 10^{17}$ cm, corresponding to a source diameter of 6' at a distance of 110 pc, the inner radius is $3\times 10^{14}$ cm, or about five stellar radii. There the dust temperature is around 1000 K, which is close to the expected condensation temperature. The optical depth in the V-band toward the star amounts to $A_{\rm V}= 48$ mag. Neglecting minor mineralogical components, the dust consists only of spherical amorphous carbon grains with a $n(a) \propto a^{-3.5}$ size distribution and lower and upper limits a_-=80 Å, a₊=640 Å, respectively, and optical constants as compiled by Zubko et al. (1996) for their type BE. At 1 mm wavelength, the mass absorption coefficient, $K_\nu$ , is about 0.5 cm² per g of dust with a $\nu^{1.6}$ frequency dependence. As the grains are much smaller than the wavelength, their size distribution does not affect the value of $K_\nu$ at FIR or longer wavelengths. Being a radiative transfer model for a dusty medium, the dust-to-gas ratio does not enter our calculations.

$\begin{figure} \par\includegraphics[width=7.5cm,clip]{4463fig3.eps} \end{figure}$

Figure A.1: Spectral energy distribution for IRC+10216. Data (squares): J- to M-band from Le Bertre (1987); 12 to 100 $\mu$ m from IRAS; 450 $\mu$ m from Young et al. (2004) for 2'' beam; 1.24 mm fluxes from Lucas & Guelin (1999), lower point for $1\hbox{$.\!\!^{\prime\prime}$ }6$ beam, upper point for total envelope; 88% of the total 3 mm flux (same authors) come from the inner 2 $\hbox{$.\!\!^{\prime\prime}$ }$ 4; 1.3 cm point, this paper. Models: the solid curve is the model described in the Appendix for maximum luminosity, the dotted line is for minimum light, when the stellar luminosity is 2.5 times weaker and the surface temperature is 1700 K. The dash-dots, which overshoot the near IR data points, are flux densities from the star alone. The dashed line gives only the stellar flux, after extinction by the envelope of $A_{\rm V}= 48$ mag.

Figure A.1 displays model results together with observational data appropriate near maximum luminosity. The solid curve is our fit near maximum luminosity within a 6' beam. It was obtained by adjusting the total optical depth, the inner radius, and the grain sizes, the other above mentioned parameters were kept fixed. Overall, the fit is satisfactory. The squares which lie significantly below the model refer to smaller angular sizes (see the caption of Fig. A.1).

Figure A.2 shows how the flux increases with beam size. As the curves are flat in the very left part of the figure, the emission from diameters less than 0 $\hbox{$.\!\!^{\prime\prime}$ }$ 2 is due to the star, not to dust, despite the central dust density peak. One can read from the figure that the fraction of the total flux contributed by the stellar photosphere increases from $\sim$ 4% at 100 $\mu$ m to $\sim$ 96% at 1.3 cm.

$\begin{figure} \par\includegraphics[width=7cm,clip]{4463fig4.eps} \end{figure}$

Figure A.2: Normalized model flux density as a function of aperture at five wavelengths. The star has an assumed angular diameter of 0 $\hbox{$.\!\!^{\prime\prime}$ }$ 075; at a comparable resolution essentially all the emission comes from the star. Substantial dust emission is indicated for beam sizes >1'' at IR wavelengths.

One can qualitatively understand why the brightness variations over a pulsational period are considerably smaller in the radio regime (Fig. 1) than at near infrared wavelengths. The monochromatic flux density which we would receive without foreground extinction, $F_\nu$ , from a star with a blackbody photosphere of temperature T and luminosity L is

How do L and T vary during a cycle? Reid & Menten (1997a) in their Eq. (1), which is based on observations by Pettit & Nicholson (1933), approximate the stellar radius and temperature of a Mira variable as a function of phase $\phi$ by $T(\phi) = 2300 + 300 \cos\phi$ K and $R(\phi)= 2 + 0.4 \sin\phi$ AU. Oscillations in luminosity, $L(\phi)$ , follow because $L=4\pi R^2\sigma T^4$ . According to these formulae, temperature and luminosity are roughly in phase, and therefore the observed radio flux $F_\nu \propto L / T^3$ changes less than L itself. The radio $F_\nu$ varies from maximum to minimum by a factor of 2 (corresponding to excursions from the mean by $\pm$ 40%) and L varies by a factor 3.76. We note that at infrared wavelengths, assuming a constant optical depth over the cycle and no dust emission, the monochromatic flux mimics L closely. The dotted curve in Fig. A.1 is our model for IRC+10216 at minimum (2.5 times less luminous, T=1700 K). At minimum, the 1.3 cm radio brightness is 25% lower and the J- and H-band flux densities are 2.5 times lower compared to maximum light.

The radio variations around the mean suggested by Fig.1 are smaller than $\pm40$ %. This suggests that the cm-wave emission may not come directly from the stellar photosphere. However, other possible reasons for this discrepancy might be a smaller variability in the bolometric luminosity (by a factor 2.5, as Men'shchikov et al. 2001 assume, and not by 3.76) or inaccuracies in the simple expressions for $R(\phi)$ and $T(\phi)$ .

Radio continuum monitoring of the extreme carbon star IRC+10216