© ESO, 2014

1. Introduction

Cool evolved stars are the main producers of dust in galaxies and are important drivers of the chemical evolution of matter in the Universe. Despite the importance of dust in a broad range of astrophysical phenomena, the conditions that lead to its formation in the outflows of evolved stars and its subsequent role in driving mass loss remain largely unknown. Oxygen-rich red supergiants (RSGs) are sources of inorganic dust (silicate and alumina) that is formed and plays a role in launching mass loss in the inner circumstellar environment of these stars. However, the dust properties in this region such as density, composition, and morphology have remained largely unknown because the spatial resolution is not adequate.

A prime target for studying the properties of dust around evolved stars is the enigmatic oxygen-rich RSG VY Canis Majoris (VY CMa). VY CMa is a mid-M spectral type (M5e Ia; Humphreys 1974) RSG with an extremely high mass-loss rate (Ṁ ~ 3 × 10^-4 M_⊙ yr^-1; Danchi et al. 1994) and can undergo periods of even more intense mass loss (Humphreys et al. 2007). This mass loss is about two orders of magnitude greater than the well known early-M spectral type RSGs, Betelgeuse and Antares. Consequently, VY CMa has a highly dense and dusty circumstellar envelope that obscures it and produces a reflection nebula at optical wavelengths (Humphreys et al. 2007). Thanks to its relative proximity (kpc; Zhang et al. 2012) and high intrinsic luminosity (L = 3 × 10⁵ L_⊙ using the photometry of Smith et al. 2001), the dust emission from the circumstellar envelope around VY CMa has been studied from optical to centimeter wavelengths (e.g., Lipscy et al. 2005; Muller et al. 2007). The dusty envelope consists of a diffuse and extended region with loops and arcs expanding over several arc seconds through a more uniform medium. At the sub-arc-second level, the envelope consists of a dense and dusty central core, but little is known about it. It is within this region that the dust condenses and radiation pressure on dust helps drive mass-loss. In this Letter we report the results of a sub-100 mas study of this region using ALMA Science Verification data.

2. Observations and results

VY CMa was observed at 321 GHz, 325 GHz, and 658 GHz as part of the ALMA Science Verification process on 2013 16–19 August, using 16 to 20 antennas of the main array with projected baselines ranging from 14 m to 2.7 km. Details of the observations and data processing are provided in Richards et al. (2014), hereafter R+14. For the purpose of the continuum analysis, we used a total of 1.74 GHz line free continuum channels around 321 GHz and 0.4 GHz centered on 658 GHz. As the data around 325 GHz provided similar results to the 321 GHz data, but were affected by the 325.15 GHz atmospheric water line resulting in increased noise, these data were not included in our subsequent analysis. The synthesized beam size, rms noise level, maximum recoverable scale (MRS), and total flux density of the two observing bands are given in Table 1. The rms noise is not uniform across the field, most likely because of low surface brightness and resolved-out emission (see below). In Table 1 we report the rms in the region of maximum noise level, while in R+14, the rms corresponds to the lower rms areas.

Fig. 1 HST/WFPC2 exposure of VY CMa using the F1042M filter (Smith et al. 2001) corrected for proper motion. The data has been truncated at 5000 counts per second and scaled with a power cycle of −1.5 to highlight the extended emission in the SW direction. The contours represent the ALMA data at 321 GHz (white) and 658 GHz (black) and the corresponding synthesized beams are located in the bottom left of the image. The 321 GHz contour levels are set at [5, 10, 30, 50, 100,....,300] × σrms, while the 658 GHz contour levels are set at [3,6,9,20,30,40,50] × σrms.

Uncertainties in the phase transfer of the data resulted in ~110 mas offsets between the main continuum peaks in the 321 and 658 GHz images. To find the optimal (RA, Dec)-offset correction, we shifted the 658 GHz image relative the 321 GHz image and computed a 2D cross-correlation function for a range of shifts. The maximum correlation between the two images gave a (−80 mas, −30 mas)-offset correction for the 658 GHz image, similar to the correction found by R+14 when aligning the peak position of the second brightest continuum component. The 321 GHz image and the position corrected 658 GHz image are plotted in white and black contours in Fig. 1, respectively. Continuous emission was detected above the 3σ_rms level on scales as large as 1.4″ at 321 GHz and 0.8″ at 658 GHz.

There are two prominent components and a weaker extended N-NE feature that are reproduced at both frequencies, as seen in Fig. 1. The position around the peak of the secondary component has recently been shown to coincide with the center of expansion of many maser emission lines and has been deduced to be the location of the star (R+14). This component (VY) connects to the weaker emission extending for ~800 mas at a position angle (PA, measured east of north) of ~20°, which may contain some further substructure. Interestingly, the brightest component in both images (C) is not at the location of the star itself but is 334 mas SE of the star in the plane of the sky. Using a distance of 1.2 kpc to VY CMa, this angular distance corresponds to 400 AU, or 61 R_⋆ if a linear radius of 1420 R_⊙ (at 2 μm) is adopted (Wittkowski 2012). R+14 have found that both blue- and redshifted 325 GHz and 658 GHz masers, which are emitted out to ~240 AU, straddle this main continuum component, so the projected distance between this component and the star is likely to be close to the actual distance.

We fit 2D elliptical Gaussian components to the two main continuum features and found that both the C and the VY components were resolved at both frequencies. The sizes and flux densities of these components were estimated from the deconvolved fits and are listed in Table A.1. The main C component is highly elongated at both frequencies in the SE direction, with a major-to-minor axis ratio, θ_maj/θ_min = 1.6 ± 0.2 at 321 GHz and θ_maj/θ_min = 2.1 ± 0.1 at 658 GHz, with the axis of elongation pointing in the direction toward VY. We find that the size of component C is consistent at both frequencies within their errors. Comparing Tables 1 and A.1, we can see that 50% and 60% of the total continuum flux density emanates from the main C component at 321 GHz and 658 GHz, respectively. The other main continuum component VY, is also resolved at both frequencies. It is elongated toward the North with θ_maj/θ_min = 1.5 ± 0.4 at 321 GHz and θ_maj/θ_min = 2.1 ± 0.5 at 658 GHz, with emission extending to 27 R_⋆ (180 AU). It is centered on the location of the stellar photosphere, so it includes emission from the star itself at these sub-millimeter frequencies. The VY component contributes about 26% and 17% of the total ALMA continuum flux density at 321 GHz and 658 GHz, respectively. The significant difference in size at the two frequencies cannot be attributed to the MRS on these small angular scales, but is probably due to a combination of different sensitivities to low surface brightness emission and the fact that the higher frequencies probe emission from hotter regions of VY.

A spectral index map of the continuum emission was created after convolving the 658 GHz data with the synthesized beam from the 321 GHz image, and regridding the convolved image to match the pixel scale of the 321 GHz image. The resulting spectral index map is represented by the color image in Fig. 2, while the contour levels represent the convolved 658 GHz image. All data below the 6σ_rms noise level in the convolved 658 GHz image have been omitted from Fig. 2. The middle region of the main continuum component C has a spectral index of α ~ 1.9 (where S_ν ∝ ν^α). The spectral index value around the second continuum, or VY, and most of the extended emission in the N-NE direction has a spectral index of α ~ 2.7. The statistical error on the spectral index due to random noise errors is less than 0.2 in the 6σ_rms region and decreases to 0.05 in the peak region. The signal-to-noise ratio at 658 GHz is less than that at 321 GHz, so this error is dominated by the 658 GHz data. There is also a systematic absolute error of ~0.22 on the entire spectral index map based on the ±15% uncertainty of the absolute flux calibration (R+14). We note that for VY, the difference in spectral index between the pixel flux densities in the map and what can be derived from the integrated flux densities is due to the smaller size of VY at 658 GHz.

Fig. 2 Spectral index map derived from the 321 GHz and 658 GHz ALMA continuum maps. The strongest continuum component has a spectral index α = 1.9 ± 0.05 at the peak, and the weaker elongated emission has a typical spectral index α = 2.7 ± 0.1. The error values here do not include the systematic error based on the ± 15% uncertainty of the absolute flux calibration, which is in the same direction for both regions of the image. The contour map represents the 658 GHz map, which was convolved with the 321 GHz restoring beam to create the spectral index map. Contour levels are set at [6,9,...,39] × σrms.

3. Discussion and conclusions

3.1. Comparison with HST

The color image in Fig. 1 is an HST/WFPC2 exposure of VY CMa taken with the F1042M filter on 1999 March 22 (Smith et al. 2001). The data is adjusted for proper motion using the values of Zhang et al. (2012), and it can be seen that the peak emission fits the second main continuum component VY very well. The typical positional uncertainty rms in this HST image is 48 mas in RA and 146 mas in Dec, while the ALMA positional error is 35 mas, ruling out a correspondence between the C component in the ALMA images and the peak of the HST image. We can conclude that most of the escaping scattered light at optical wavelengths is emitted along the line of sight to the star itself. There is no evidence of the main C component in the HST image, probably because the HST emission does not trace dense dusty emission, but rather more tenuous regions in the circumstellar environment. Scattered light polarimetry with the HST does show a pronounced lower polarization toward the SE of the star (Jones et al. 2007), potentially indicating enhanced obscuration along the direction between the star and component C.

The extended SW emission feature in the HST image, which may be a density cavity in the circumstellar environment formed by a previous episodic mass-loss event, has no detectable signature in our ALMA continuum images. This extended feature ends in what is named the SW clump in Smith et al. (2001) and Humphreys et al. (2007). This clump is thought to be related to an ejection event from the star and has been estimated to have a dust-mass lower limit of 5 × 10^-5 M_⊙ and a temperature between 80 K and 210 K (Shenoy et al. 2013). However, we see no sign of this clump in the ALMA data. With our 3σ_rms limits, we rule out a mass >7 × 10^-6 M_⊙ (for T = 80 K) or > 3 × 10^-6 M_⊙ (for T = 210 K). For the ALMA observations to be consistent with the lower limit estimates of Shenoy et al. (2013), the dust temperature would have to be <20 K. Alternatively, a different dust composition with different dust scattering and emissivity indices needs to be invoked. The exact nature of the SW clump is thus not clear. We also find no other compact submillimeter emission in the larger ~8″ region of dust seen with the HST.

3.2. Properties of the dust

Dust is the main source of the observed thermal submillimeter emission around VY CMa. If the dust is optically thin at these frequencies, and we assume the Rayleigh-Jeans approximation, then the observed flux density minus the stellar flux density contribution S_⋆ is proportional to the mass of the emitting dust M_d, such that (1)where Q_ν is the grain emissivity, a_g and ρ_g are the radius and mass density of the dust grains, respectively, d is the distance to the star, and T_d is the dust temperature (Knapp et al. 1993). Assuming that the grain emissivity has a power law dependence on frequency, Q(ν) ∝ ν^β, then the flux density is S_ν ∝ ν^α, where α = 2 + β. Here, β is the dust emissivity spectral index and has typical values of β = 1.8 ± 0.2 for both diffuse and dense clouds in the interstellar medium (Draine 2006). For evolved stars, the dust emissivity index is generally lower (Knapp et al. 1993), and for VY CMa previous observations have found values between β ~ 0 − 1.0 (Shinnaga et al. 2004, Kamiński et al. 2013), while its circumstellar environment has been modeled using values of β = 0.9 (Knapp et al. 1993).

Our spectral index map in Fig. 2 shows that the two spatially resolved continuum components in the inner circumstellar envelope have different spectral indices. The second brightest continuum component, VY, has a typical dust emissivity index of β = 0.7 ± 0.1 and is consistent with the previous single-dish observations outlined in Appendix B. However, much of the main continuum component, C, has a dust emissivity index of β = −0.1 around the location of the peak emission, meaning that the submillimeter emission from this component is, or is close to, becoming optically thick. Following Herman et al. (1986), optically thin dust at a radius r_d is heated by the incident stellar radiation field to a dust temperature (2)where L_⋆ and T_⋆ are the stellar luminosity and temperature, respectively. Assuming T_⋆ = 3490 K (Wittkowski et al. 2012), the main continuum component, C, would have an isothermal temperature of 450 K at 400 AU from the star, while the VY component would have a value of 970 K at a mean dust radius of 67 AU. We stress that in deriving the temperature for the C component, we have assumed that the dust is optically thin even though the spectral index map indicates otherwise. We thus use 450 K as an upper limit to the temperature of C. These rough estimates for the dust temperature allow us to empirically calculate the optical depth at both observing frequencies, the values of which are given in Table A.1. The VY component is optically thin at both frequencies with τ_{321 − 658 GHz} = 0.02 − 0.06, while the main continuum component C, has an optical depth of τ_{321 − 658 GHz} ~ 0.2, and would become optically thick at 110 − 120 K at both frequencies. Therefore, a cooler dust component in the center of C would explain the spectral index having a lower value than expected for optically thin dust. Alternatively, larger dust grains could also result in low β values although our observing frequencies would then suggest that the dust grains would need to be much larger than 100 μm. Interestingly, we find that the thermal molecular emission lines from the ALMA data avoids the C component, as is consistent with the molecular emission offsets found in SMA observations (Muller et al. 2007; Kamiński et al. 2013). This could potentially support the low temperature and high optical depth suggestion because of low excitation and/or depletion onto dust grains.

Equation (1) can be re-arranged to find the total mass of dust in each of the two main continuum components. The main component C may not be fully optically thin at both frequencies and so its dust mass estimates will be lower limits. Following Knapp et al. (1993) we assume typical oxygen-rich star values for the radius and mass density of the dust grains, a_g = 0.2 μm and ρ_g = 3.5 g cm^-3. We slightly alter their grain emissivity function to include our derived dust emissivity spectral index so that Q_ν = 5.65 × 10^-4 (ν/ 274.6 GHz)^0.7. We also assume that the main C component has the same dust composition as VY, so the dust emissivity spectral index of C would also be β = 0.7 if fully optically thin. We then derive dust masses of 2.5 × 10^-4 M_⊙ at 321 GHz and 1.6 × 10^-4 M_⊙ at 658 GHz for the main dust component C, having assumed no photospheric contribution. The lower mass value at 658 GHz is probably due to a combination of optical depth effects and an inaccurate grain emissivity law. We derive dust masses of 4.0 × 10^-5 M_⊙ at 321 GHz and 1.8 × 10^-5 M_⊙ at 658 GHz for the VY component having subtracting the stellar contribution, which is calculated in Appendix C. The discrepancy in the dust mass values here could be due to an inaccurate grain emissivity law or differences in sensitivities to low surface brightness emission.

3.3. Dust mass-loss history and mechanisms

Considering that we find that part of the dust emission around VY CMa arises in compact clumps, we can obtain a new estimate of the average dust mass-loss rate by determining the dust mass that is resolved out in the ALMA observations. We assume the dust to be extended over the ~8″ size of the HST image and use an average dust velocity V_d that depends on grain size and that for VY CMa is ~60 km s^-1 for grains of ~0.25 μm (Decin et al. 2006). We assume an average temperature of 200 K throughout the shell. The observed VY component arises from the dust lost in the last ~20(60 km s^-1/V_d) yr and the resolved-out emission (~1.1 Jy at 321 GHz and ~7.8 Jy at 658 GHz as shown in Fig. B.1) corresponds to the dust lost over ~360(60 km s^-1/V_d) yr before. For β = 0.7 using Eq. (1), we find a dust mass of ~1.8 × 10^-3 M_⊙ which is consistent at both frequencies, and a dust mass-loss rate of 5 × 10^-6 (V_d/ 60 km s^-1)  M_⊙yr^-1, which is four times higher than previously estimated by Sopka et al. (1985) and slightly higher than Knapp et al. (1993), taking the different adopted velocities and distances into account. Then, taking all compact emission from the ALMA observations, with the exception of what is around the star, to be in optically thin clumps at ~450 K, those clumps contain ~3 × 10^-4 M_⊙ of dust. This is a lower limit, because a significant component of cold dust can potentially be present in C and in the SW and other HST features not detected by ALMA. We can thus conclude that >17% of the dust mass around VY CMa is located in clumps ejected within a more quiescent, roughly spherical stellar wind. This picture of a quiescent mass loss disrupted with periods of more localized intense ejections of mass is consistent with the findings of Humphreys et al. (2007) for VY CMa and Ohnaka (2014) for Antares (M1 Iab).

It has been hypothesized that large convection cells initiate mass loss in RSGs by lifting material out to distances beyond the photosphere where dust can form and drive mass loss (e.g., Lim et al. 1998). These convection cells are predicted to have time scales of ~150 days (Schwarzschild 1975) and are expected to eject material at random PAs from the star, generating supersonic motions and shocks in the process (Lion et al. 2013). Our observations have revealed a massive dust clump centered at 60 R_⋆ and a weaker continuous dust feature extending 90 R_⋆ to the North. These features imply a continuous directed dust mass-loss from the star over the past 30–50 yr. This contrasts with what would be expected from mass loss that is only driven by convection, which on these spatial scales would show many dust clumps at random PAs. A mass-loss mechanism that is localized but much more stable over time (i.e., over a few decades) is required to explain our observations. This would suggest a magnetic origin. Indeed a surface magnetic field strength of 10³ G has been extrapolated from circumstellar maser measurements for VY CMa (Vlemmings et al. 2002). Cool spots on the photosphere due to localized long-lived MHD disturbances could then enhance local dust formation, and hence drive mass loss in the localized directions that are evident in our ALMA data. Such a mass-loss mechanism would still appear as random mass ejections on larger spatial scales and could be the source of the arcs and clumps seen in VY CMa’s outer circumstellar environment by Smith et al. (2001) and Humphreys et al. (2007).

Acknowledgments

This paper makes use of the following ALMA data: ADS/JAO.ALMA#2011.0.00011.SV. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada) and NSC and ASIAA (Taiwan), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. EOG and WV acknowledge support from Marie Curie Career Integration Grant 321691 and ERC consolidator grant 614264.

References

Online material

Appendix A: Properties of the two main continuum components

In Table A.1 we list the properties of the two main continuum components, C and VY.

Appendix B: Comparison with other (sub-)mm observations

There have been many previous (sub-)millimeter continuum flux density measurements of VY CMa. In Fig. B.1, we present both the single-dish observations, using the compilations of Knapp et al. (1993) and Ladjal et al. (2010), and the interferometric measurements from this paper and from SMA observations (Shinnaga et al. 2004; Muller et al. 2007; Kamiński et al. 2013). It is immediately apparent that the interferometric observations systematically underestimate the total flux density and can be attributed to flux density that is resolved out by missing short baselines, as well as by having low surface brightness. The MRS at higher frequencies is smaller so more flux density is resolved out, resulting in a larger discrepancy between the single-dish and interferometric measurements. The spectral index determined from the single-dish observations is α = 2.5 ± 0.2 and is consistent with our findings in Sect. 3.2 for the optically thin dust. We note that the single-dish data points in Fig. B.1 will also contain flux density from molecular emission lines. Kamiński et al. (2013) find that this accounts for about 25% of their measured flux density between 279 and 355 GHz. This level will vary at the different frequencies and will mainly result in an added upward scatter but should not affect the overall spectral index. Assuming the true spectral index of the interferometric observations to be the same as that of the single-dish observations and assuming there are no significant changes in dust properties on the different scales, we can conclude that the ALMA observations lose ~65% of the emission at 321 GHz and ~72% of the emission at 658 GHz. Considering the MRS and sensitivity of the ALMA observations, most of the submillimeter dust continuum flux density is thus located in a smooth low surface brightness distribution stretching beyond ~4″, which is consistent with the dust distribution seen over ~8″ with the HST (e.g., Smith et al. 2001).

Fig. B.1 Compilation of the single dish bolometer and interferometric mm/sub-mm continuum observations of VY CMa. The gray filled circles are single-dish flux density measurements, while the blue filled diamonds and the red filled squares are interferometric measurements. The solid line indicates the spectral index fit to the single-dish observations, with α = 2.5 ± 0.2, while the dotted line is a linear fit to the interferometric observations.

Appendix C: Stellar flux density contribution

The stellar flux density contribution S_⋆, needs to be considered when calculating the dust mass from the VY component at these ALMA frequencies. Previous studies have estimated this contribution by calculating the flux density of an optically thick blackbody of radius R_⋆ and temperature T_eff. Estimating the contribution from this method yields S_⋆ = 26.5 mJy at 321 GHz and S_⋆ = 111 mJy at 658 GHz, where we have assumed a stellar radius of R_⋆ = 1420 R_⊙ (6.84 AU or 5.7 mas) and an effective temperature T_eff = 3490 K (Wittkowski et al. 2012). However, RSGs have weakly ionized extended atmospheres, which will become opaque at (sub-)mm frequencies. To estimate this contribution at these ALMA frequencies, we scaled the Harper et al. (2001) semi-empirical model for the M2 Iab RSG Betelgeuse (without the silicate dust) to the angular diameter of VY CMa. The main source of opacity in this model is the H⁻ and H free-free opacity from electrons produced by photoionized metals. We then assumed the same ionization fraction and scaled the particle densities so that good agreement was obtained with the Very Large Array centimetre observations of VY CMa from Lipscy et al. (2005). The stellar contribution from this weakly ionized atmosphere is 36 mJy at 321 GHz and 124 mJy at 658 GHz.

All Tables

All Figures

Fig. 1 HST/WFPC2 exposure of VY CMa using the F1042M filter (Smith et al. 2001) corrected for proper motion. The data has been truncated at 5000 counts per second and scaled with a power cycle of −1.5 to highlight the extended emission in the SW direction. The contours represent the ALMA data at 321 GHz (white) and 658 GHz (black) and the corresponding synthesized beams are located in the bottom left of the image. The 321 GHz contour levels are set at [5, 10, 30, 50, 100,....,300] × σrms, while the 658 GHz contour levels are set at [3,6,9,20,30,40,50] × σrms.

In the text

Fig. 2 Spectral index map derived from the 321 GHz and 658 GHz ALMA continuum maps. The strongest continuum component has a spectral index α = 1.9 ± 0.05 at the peak, and the weaker elongated emission has a typical spectral index α = 2.7 ± 0.1. The error values here do not include the systematic error based on the ± 15% uncertainty of the absolute flux calibration, which is in the same direction for both regions of the image. The contour map represents the 658 GHz map, which was convolved with the 321 GHz restoring beam to create the spectral index map. Contour levels are set at [6,9,...,39] × σrms.

In the text

Fig. B.1 Compilation of the single dish bolometer and interferometric mm/sub-mm continuum observations of VY CMa. The gray filled circles are single-dish flux density measurements, while the blue filled diamonds and the red filled squares are interferometric measurements. The solid line indicates the spectral index fit to the single-dish observations, with α = 2.5 ± 0.2, while the dotted line is a linear fit to the interferometric observations.

In the text