ALMA Reveals Hidden Morphologies in the Molecular Envelope of VY Canis Majoris

Evolved stars with initial masses between 9 and 40 M_⊙ are considered to end their lives as Type II-P or Type II-L supernovae (SNe) on the red supergiant (RSG) track. However, 20 yr of observations have shown a lack of RSGs above 18 solar masses among SNe progenitors (Smartt et al. 2009; Smartt 2015), which suggest a different outcome for these stars. With their enhanced mass loss, such RSGs could collapse directly to black holes or evolve back to warmer temperatures. With its uniquely explored history of high mass-loss events, the red hypergiant VY Canis Majoris (VY CMa) is the ideal candidate for understanding this puzzle.

VY CMa is one of the most massive (∼40 M_⊙) and most luminous (5 × 10⁵ L _⊙) RSGs known with a variable mass-loss rate of ∼ 2–4 × 10⁻⁴ M_⊙ yr ⁻¹ (Shenoy et al. 2016; Dinh-V-Trung et al. 2022). Mass-loss events occurring over the course of several hundred years (e.g., Humphreys et al. 2007, 2021) have created an irregular circumstellar envelope consisting of numerous clumps, knots, arcs, and filaments (e.g., Monnier et al. 1999; Smith et al. 2001; Humphreys et al. 2005,2007, 2021). Recent Hubble Space Telescope (HST) observations probed scattered light from dust and KI emission (e.g., Humphreys et al. 2007, 2019, 2021). Doppler shifts of KI lines and proper motions of the dust knots were critical in assigning velocities to the structures and demonstrated that the various winds in VY CMa are spatially and kinematically distinct, tracing a complex mass-loss history.

The envelope of VY CMa has also been explored at millimeter/submillimeter wavelengths. A 1 mm spectral survey (214–285 GHz) of VY CMa was conducted using the Submillimeter Telescope (SMT) of Arizona Radio Observatory (ARO) by Tenenbaum et al. (2010). This and other studies revealed a rich chemical inventory in the envelope, with exotic molecules TiO, VO, and AlOH (e.g., Tenenbaum & Ziurys 2010; Kamiński et al. 2013b). Atacama Large Millimeter Array (ALMA) observations at an angular resolution of 01–02 revealed clumpy structures in TiO₂ and NaCl in the inner wind close to the dust formation zone and emission ∼ 220 R_* from the star, in the "SW clump" (De Beck et al. 2015; Decin et al. 2016). Images of continuum emission at Band 6 and Band 7 (O'Gorman et al. 2015; Kamiński 2019) show five dusty clumps forming an elongated structure near the star, including bright "Clump C" (Richards et al. 2014). More recently, Quintana-Lacaci et al. (2023) imaged lines of SiO, H₂S, SO, and NaCl with 0025 resolution, showing collimated ejecta in different directions disrupting a slower wind, as suggested by Ziurys et al. (2007).

Single-dish observations of CO, on the other hand, suggested that the outflows in VY CMa extend to much larger spatial scales ∼12'' in extent (Ziurys et al. 2007). Furthermore, molecular images of VY CMa, made with the Submillimeter Array (SMA) showed emission extending well beyond 2'' but resolved out about half the flux on larger scales (Kamiński et al. 2013a). The 37 μm image of Shenoy et al. (2016) also revealed a larger envelope size of ∼15''.

To further investigate the extent of the envelope, ALMA has been used to image VY CMa in 1 mm molecular emission with resolutions between 02 and 15. Combined with previous single-dish spectra, the images are sensitive to scales as large as ∼30''. For the first time, the complete envelope of VY CMa has been revealed. New winds have been discovered in CO that extend northeast and east of the star as far as ∼9'', increasing the extent imaged in detail by a factor of ∼2. In the Letter, we present these new, comprehensive images of VY CMa in CO and HCN.

VY CMa was observed using ALMA at Band 6 for the project 2021.1.01054.S between 2021 December and 2022 June. A total of 43 12 m antennas were used, covering frequencies 216–270 GHz in two array configurations, yielding high- and low-angular resolutions around 02 and 15, which were used separately and in combination to produce continuum-subtracted image cubes at ∼1.3 km s⁻¹ spectral resolution. The maximum recoverable scale was 11'' with a full width half maximum field of view (FOV) of 26''. Data taken using the ARO SMT (Tenenbaum et al. 2010) were added to provide sensitivity to angular scales up to 30''. The rms noise per channel was <1 mJy. Velocities are in the kinematic local standard of rest (LSR) convention relative to the rest frequency of the given line, assuming V_LSR= 19.0 km s⁻¹ for the star. The astrometric accuracy is 15 mas for high-resolution cubes and ≤30 mas for the combined data; the absolute flux scale accuracy is better than ∼10%. The position used for VY CMa (α = 07^h22^m58324; δ = −25°46$^{\prime}$0301 J2000 ICRS) was estimated from our highest-resolution 1 mm continuum image by fitting a 2D Gaussian component, with an uncertainty of 002 in each coordinate (see the ALMA Technical Handbook and ALMA Memo 620). More information on data reduction, resolution, sensitivity, and other parameters of the cubes used are given in Appendix A.

While many spectral lines were observed in the ALMA bandpass (see Appendix B, Table B1), this work focuses on CO (J = 2 → 1: 230, 538.0 MHz) and HCN (J = 3 → 2: 265, 886.4 MHz). Single-dish data from ARO SMT were included using feathering. ALMA images presented are moment 8 maps, which use the peak intensity of the spectrum, as the best choice to decipher the complex features in VY CMa (rather than the commonly used moment 0 channel maps that integrate over the line profile.)

Spectra arising from certain features in the CO maps were fitted with the non-LTE radiative transfer code ESCAPADE (Adande et al. 2013; also, Singh et al. 2021, 2022). The code, based on the Sobolev approximation, employs the escape probability formalism, considering collisional and infrared dust excitation. The code was used to model the orientation of collimated outflows with respect to the plane of the sky [(x, y) plane], and the outflow velocity perpendicular to this plane or V_xz . Model fits are shown in Figure 3 as overlays (in red or gray scale).

Images for the ¹²CO, J = 2 → 1 and HCN, J = 3 → 2 lines are shown in Figures 1 and 2, respectively (left panel: high-resolution data; right panel: composite image). The red contours in the left panels indicate the approximate position of the star. In the combined image, all the flux has been recovered, as verified by the close resemblance of the ALMA-combined CO and HCN line profiles with those from the SMT (see Figure 3 versus Tenenbaum et al. 2010).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 1 shows the distinctive new structure in the molecular envelope in VY CMa. There are clear features to the southwest (Arc 1 and Arc 2) and northwest (NW Arc), observed previously in HST images of VY CMa (Humphreys et al. 2007). In addition, there are newly revealed outflows to the east and northeast, one extending as far as ∼9'' from the star. The total extent of the CO emission is ∼14'' in diameter, as noted previously by Ziurys et al. (2007). These regions were not detected earlier by HST because of high obscuration in this direction (Smith et al. 2001), while past ALMA measurements were insensitive to large-scale structure (Decin et al. 2016; Kamiński 2019; Quintana-Lacaci et al. 2023). HCN emission (Figure 2) is primarily located near the star and appears only faintly in some extended features (Arc 2, NW Arc and NE Arc). It is absent from Arc 1 and other E structures.

Table 1 summarizes the features observed in the ALMA images, approximate average positions, and LSR velocities in CO, and when observed, in HCN, assuming V_LSR = 19 km s⁻¹ for the star. These velocities are determined from the strongest emission peaks in composite spectra generated for each region; see Figure 3. In certain cases, a strong secondary peak was also visible in the data and is listed in the table. Other, less intense emission appear in some of the spectra, which could also arise in part from contaminating regions. From the ESCAPADE modeling of line profiles, rough estimates of expansion velocities V_xz and their approximate angle with respect to the plane of the sky for Arc 1, NE Extension, and NW Arc were derived. These values are listed in Table 1. The table also contains preliminary ¹²C/¹³C ratios for specific outflows.

4.1. ALMA View of Previously Observed Features

The features observed by ALMA in CO almost have a one-to-one correspondence to those seen in the HST image of dust-scattered starlight, as shown in Figure 4. Here the high-resolution ALMA map and the HST image (Humphreys et al. 2007) are shown on the same spatial scale. Individual structures are indicated by letters. Arc 1, Arc 2, and the NW Arc (a, b, e) appear to have virtually identical structures in both CO and HST images. This result suggests that the dust and gas are well mixed in these outflows. CO is likely to be excited by collisions, such that the rims of these arcs may represent shocked material. The less conspicuous S Arc (c), SW Clump (d), S Knots (f), and SE Loop (g) also are present in both images. The NW Arc and Arc 2 are visible in the HCN image, while weaker features are less apparent, and Arc 1 is absent. Such variation likely occurs because HCN requires higher densities for collisional excitation than CO (10⁷ versus 10⁴ cm⁻³ for the observed lines). Individual features are discussed below.

Zoom In Zoom Out Reset image size

Arcs 1 and 2

Initially detected in the near-IR by Monnier et al. (1999), Arcs 1 and 2 lie ∼7'' and ∼ 5'' southwest from the star. Arc 1 was thought to trace blueshifted gas (Humphreys et al. 2007; Ziurys et al. 2007). The ALMA CO spectrum shows that Arc 1 has a prominent, blueshifted velocity component near V_LSR = −9 km s⁻¹ (Figure 3 and Table 1), strongest near the tip of the arc. A second blueshifted component is also present near 11 km s⁻¹, more obvious in the arc "wings." HCN was not observed in Arc 1 as the density is too low, as noted. Modeling the strongest peak in the profile suggests that the most blueshifted material is moving at an angle ∼ −50° relative to the plane of the sky, with an expansion velocity V_xz ∼ 40 km s⁻¹. These estimates are in good agreement with the HST results, which suggested an angle of ∼ −33° and V_xz ∼ 48 km s⁻¹, averaging over Arc 1 positions A, C, and D. Note that the total expansion velocity for Arc 1 from the HST analysis is V_tot = 68.2 ± 2.5 km s⁻¹ (Humphreys et al. 2007).

Arc 2 has a complicated line profile of six "sub-peaks" with the maximum near V_LSR ≈ 24–25 km s⁻¹, seen both in CO (Figure 3) and HCN. The redshifted gas is observed toward the eastern part of the arc, while the blueshifted material lies to the west. This complex velocity structure suggests that Arc 2 is composed of multiple outflows or is even a bubble. The KI profile for Arc 2 traces the more western region and exhibits two velocity components near V_LSR ∼ 5 and 50 km s⁻¹ (Humphreys et al. 2005), also seen in CO at this position.

NW Arc

One of the brightest features in the HST images, the NW Arc is observed in both CO and HCN. The CO and HCN peak velocities of 48 and 40 km s⁻¹ do not entirely agree and suggest a density gradient. The NW Arc is also known to have a velocity gradient (Humphreys et al. 2005), which is also apparent in the broad CO line profile. This spectrum extends from −10 to 60 km s⁻¹, east to west, and the intensity decreases to negative velocities (Figure 3). The primary peak in the CO line profile suggests that the NW Arc is behind the star, as previously surmised (Humphreys et al. 2005), but is tilted ∼ 80° with respect to the plane of the sky, with V_xz ∼ 40 km s⁻¹ (Table 1). Humphreys et al. (2007) placed the arc at an angle of 22° with respect to the plane of the sky, based on sampling the eastern half of the structure. (There are no proper motion data for the western edge). Polarimetry data suggest an angle that varies from 35° to 45° from east to west (Jones et al. 2007). The CO spectrum, which preferentially samples the western material, indicates that the NW Arc may curve even further away from the plane. Humphreys et al. (2007) suggested the tip may eventually curve back toward the star. Humphreys et al. (2007) found a total expansion velocity of V_tot = 46 ± 4 km s⁻¹ and V_xz = 39.5 km s⁻¹, in good agreement with the CO-derived velocity.

South Features

While the SW Clump is prominent in CO (d), the SE Loop, S Knots, and S Arc ((g,), (f), and (c)) appear as weak features. None of these regions are clearly defined by HCN emission. These structures will be discussed in a later paper; only a cursory summary will be given here.

The "SW Clump," which is actually multiple features, has been the subject of several studies (e.g., Decin et al. 2016; Kamiński 2019; Humphreys et al. 2021). This composite clump has a complex velocity structure in CO with three distinct peaks at 19, 50, and −11 km s⁻¹, the first being the most intense. This result is consistent with the multiple component "bubble" structure identified by Kamiński (2019). The SE Loop has notable velocity structure, with strongest CO emission near 5 km s⁻¹, most apparent toward the ends of the loop. Another weaker component is near 25 km s⁻¹, lying at the loop center. The S Knots, the brightest structures in the SE envelope, are blueshifted, with a peak velocity in CO near −7 km s⁻¹. A weaker, redshifted tail is also apparent in the spectrum. The S Arc is kinematically complex with three components: the strongest near 23 km s⁻¹; a major, blueshifted feature at −10 km s⁻¹; and a less distinct, redshifted peak at ∼50 km s⁻¹.

4.2. Newly Discovered Outflows

4.3. Preliminary ¹²C/¹³C Ratios in Select Features

This work is supported by NSF Grant AST-1907910. The following data were used: ADS/JAO.ALMA#2021.1.01054.S. ALMA is a partnership of ESO, NSF, and NINS, with NRC, MOST, ASIAA, and KASI, in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. NRAO is a facility of the NSF operated under cooperative agreement by Associated Universities, Inc.

Three separate tunings were observed (each covering four spectral windows, SPWs, of 1.875 GHz bandwidth) between 216 and 270 GHz, at 1.25–1.45 km s⁻¹ velocity resolution (depending on frequency). The calibration for the bandpass, flux, and pointing was done using J0538-4405, while the phase calibrator was J0725-2640. Two configurations individually provided ∼02–04 or ∼10–15 resolutions and ∼34 or ∼11'' maximum recoverable scales, hereafter high and low resolution. Observations followed standard ALMA procedures (ALMA Technical Handbook ⁶ ), and details of conditions and calibrators are given in the ALMA archive. All data reduction was performed in CASA 6 (2022 PASP, 134k4501C). The data were pipelined, followed by extraction of the calibrated VY CMa data and identification of line-free channels for each tuning and configuration to allow phase and amplitude self-calibration of the continuum. The solutions were applied to all data. We followed procedures similar to those outlined in Gottlieb et al. (2022) to ensure a consistent stellar position and flux scale for all ALMA data.

From the ALMA Technical Handbook, the expected astrometric accuracy at least ∼18 (60) mas for the high- (low-) resolution data for features with a signal-to-noise ratio > 20 (as for the continuum peaks). The recovered positions of the check source were accurate to <∼15 mas in the three high-resolution tunings, and the positions of the low-resolution continuum peaks agreed with each high-resolution peak within 30 mas so we conservatively take this as the astrometric accuracy of the combined images.

The flux scale of ALMA standard calibrators is accurate to ∼5% in Band 6, but flux density errors increased in transferring amplitude corrections to the target and combining data from different epochs. We checked by comparing the high- and low-resolution flux densities for continuum emission on the same scales at each frequency. This showed that the flux scale uncertainties were <10% except for the lowest-frequency low-resolution data which were discrepant by 11%; we scaled these to make the flux scale more consistent.

From O'Gorman et al. (2015) Table A1, the flux-weighted mean of spectral indices for the two main continuum sources (clump C and the star VY) is α = 2.15. We measured the flux densities for the continuum in each separate tuning of the combined ALMA data in an ellipse of full-width axes 21 × 24 (E × N), giving α = 2.1 (0.1) showing that the relative flux density across our band is accurate to ∼5%. The predicted combined flux density at our intermediate frequency is 293 mJy, compared with an observed 263 mJy, or 11% discrepancy, which includes the earlier observations uncertainty and possibly genuine variability.

After subtracting the continuum, image cubes were made for each SPW for separate configurations and for combined configurations. Spectra taken using the ARO SMT were added to the combined ALMA data using feathering to provide sensitivity to angular scales up to 30''. All image cubes used here have had the primary beam correction applied. Velocities are in the kinematic LSR convention relative to the rest frequency of the line shown.

The ¹²CO/¹³CO ratios were estimated after smoothing the data cubes to the lower resolution of the ¹³CO images and blanking both at the level where the weaker ¹³CO emission fell below 3× local rms noise.

A table of molecular lines, frequencies, and intensities, observed in the dataset, is given in Table B1.