© ESO, 2016

1. Introduction

The post-main-sequence evolution of stars with an initial mass in the range of 20−40 M_⊙ may lead through a cool red supergiant (RSG) phase with strong mass loss. These stars can potentially evolve through a Wolf-Rayet phase as well if mass loss is sufficient to remove the hydrogen envelope and, if not, they may explode as a supernova (Smith et al. 2009). The amount and nature of the mass loss dictates the future evolution of the star and determines their supernova (SN) appearance, be it as an SN of type II-P, II-L, IIb, IIn, or Ib/c. If a dense circumstellar envelope is still present at the moment of explosion, it will provide an obstacle for the SN blast wave to overtake (Pun et al. 2002).

Currently the driving mechanism of the mass loss of RSGs is not yet understood. Mechanisms based on radiation pressure on freshly synthesized dust grains, magneto-acoustic waves, and turbulent pressure owing to convection in combination with radiative pressure on molecular lines (Josselin & Plez 2007) are invoked, or a combination of the above in a so-called hybrid model featuring magneto-rotational effects with a dust driven-wind scenario (Thirumalai & Heyl 2012). Not being able to constrain the driving mechanism(s) represents a serious limitation for current stellar evolution models, which therefore adopt a simplified parametrized, time-averaged, mass-loss rate prescription (e.g. Reimers 1975; de Jager et al. 1988; Nieuwenhuijzen & de Jager 1990). Empirical studies, however, show that the mass loss rate, Ṁ(t), may vary significantly during the RSG evolution (e.g. Decin et al. 2006) and that these stars can display sudden episodic bursts of mass loss that is potentially linked to the magnetic and/or convective nature of the star (Smith et al. 2001; Adande et al. 2013).

The circumstellar envelopes (CSEs) created by the intense mass loss are rich chemical laboratories in the Universe. Because of the low temperatures in the envelopes and the long duration of the mass loss, molecules and solid-state species are formed in these CSEs. Through stellar winds, these elements are injected into the interstellar medium (ISM), thereby gradually enriching their surroundings.

Efforts to deepen our understanding of the origin of the mass-loss and of the thermodynamical and chemical behaviour of the CSE proceed on two separate tracks: more detailed studies of specific RSGs which are often selected by virtue of their proximity, and global studies of a larger sample of RSGs to derive general characteristics. The first of these approaches is pursued in this paper, in which we present ALMA high spatial resolution () observations of the close-by RSG VY CMa in thermal sodium chloride (NaCl) emission lines, tracing the inner wind region within ~1.5′′ (or ~280 R_⋆) from the central star.

1.1. The red supergiant VY CMa

VY CMa is one of the most intrinsically luminous red stars in the galaxy (L ~ 3 × 10⁵ L_⊙; Smith et al. 2001; Choi et al. 2008). Its current mass is estimated around 17 M_⊙, but initially it was some 25 M_⊙ (Wittkowski et al. 2012). VY CMa is also known to be one of the largest RSGs with a stellar radius, R_⋆ , of 1420 ± 120 R_⊙ (Wittkowski et al. 2012) or 5.53 mas at a distance of 1.2 ± 0.1 kpc (Choi et al. 2008; Zhang et al. 2012). The average recent mass-loss rate is high, around 1−4 × 10^-4 M_⊙ yr^-1 (Danchi et al. 1994; Monnier et al. 1999; Decin et al. 2006), producing a large circumstellar nebula some 10′′ (or 0.06 pc) across (Humphreys et al. 2007; Muller et al. 2007). The stellar velocity estimated from OH 1612 MHz and SiO maser observations is 22 ± 1 km s^-1 (Richards et al. 1998) (all velocities are with respect to the local standard of rest, LSR).

Optical and near-infrared high-resolution imaging and interferometry have revealed a complex distribution of knots and filamentary arcs in the asymmetric reflection nebula around the obscured central star (Kastner & Weintraub 1998; Monnier et al. 1999; Smith et al. 2001). Images, polarimetric data and radial velocity measurements (Smith 2004; Humphreys et al. 2007; Jones et al. 2007) show evidence that the arcs and knots in the ejecta are not only spatially distinct features, but are also kinematically separate from the surrounding diffuse material representing the steady wind (Humphreys et al. 2007; Smith et al. 2009). The current understanding is that VY CMa went through a series of episodic mass ejections with dense cloudlets moving at ~35 km s^-1 that produced a very high mass-loss rate of ~ 2 × 10^-3 M_⊙ yr^-1 between a few hundred and 1000 yr ago. The diffuse halo might expand at a slightly lower speed (~25 km s^-1) and corresponds to the current average mass-loss rate of ~ 3 × 10^-4 M_⊙ yr^-1 (Smith et al. 2009). Optical and infrared data also provide some evidence for a disk-like geometry near the central object (Herbig 1970, 1972; McCarthy 1979; Efstathiou & Rowan-Robinson 1990). The derived geometry is not always consistent between these studies. K-band interferometric measurements indicate a dusty disk oriented east-west (Monnier et al. 2004), but infrared images have been interpreted as evidence for a northeast-southwest axis of symmetry that represents the polar/rotation axis of the star and disk, which is tilted by 15° to 30° to our line of sight (for an overview, see Smith et al. 2001; Humphreys et al. 2007). Maser data have been interpreted as an expanding outflow, which at sub-arcsecond scale is more extended in the east-northeast and west-southwest direction (Morris & Bowers 1980; Bowers et al. 1983; Richards et al. 1998). The idea of a preferred northeast-soutwest axis is roughly supported by the Hubble Space Telescope (HST) data of Smith et al. (2001) of which the asymmetric emission in the South-West could be interpreted as being reflection along the lower density bipolar axis oriented northeast-southwest with the south-western lobe closer to us (see also the cartoon of the proposed geometry in Fig. 13 in Smith et al. 2009). The apparent random orientation of the arcs, filaments and knots in the wind of VY CMa suggest that they were produced by localized mass ejections, not necessarily aligned with either the star’s presumed northeast-southwest axis or its equator, and hence independent of any proposed bipolar geometry.

While data taken at optical and infrared continuum wavelength bands are prone to high extinction effects complicating an unambiguous interpretation of the inner wind structure of VY CMa, the scattering and extinction in the far-infrared and submillimeter domain are almost negligible. The retrieval of the morphological, kinematical and chemical structure of the inner wind region from prior observations was often hampered by a lack of spatial resolution to resolve the different kinematical components. Another difficulty was the association of molecular line emission to specific morphological features which was only based on a tentative interpretation of spatially unresolved velocity features in the broad asymmetric line profiles (Ziurys et al. 2007; Milam et al. 2007; Adande et al. 2013). Modern interferometers offer an angular resolution comparable to the angular scale of some arcs and knots in the CSE of VY CMa and interferometric far-infrared and submillimeter observations of molecular lines hence provide us with crucial diagnostics to unravel the physical and chemical processes shaping the inner envelope of VY CMa. With the aim to investigate the spatio-kinematical structure of the nebula around VY CMa, Kamiński et al. (2013b) have used the Submillimeter Array (SMA) to undertake an interferometric spectral imaging survey between 279 and 355 GHz at a spatial resolution of 0.9′′. Twenty-one molecules were detected, of which 11 show a “point-like” morphology, 4 were indicated as “double” with a second source located to the southwest of the star (at 0.̋8 to 1.̋2), one molecule (CN) has a “ring-like” distribution, and 5 showed an “extended” distribution with multiple components. One of the detected molecules is NaCl, showing a “point-like” angular distribution in the SMA data with the exception of the most intense NaCl line which shows weak emission in the southwest clump.

VY CMa was observed with ALMA during science verification in band 7 and band 9 (Richards et al. 2014). Four NaCl thermal emission lines are detected in the band 7 data, and are spatially resolved at an angular resolution of ~ (see Sect. 2). The benefits of studying the spatio-kinematical and chemical structure of the nebula around VY CMa via NaCl emission are manifold: (1) NaCl is a diatomic molecule with a simple Grotrian-diagram reducing the complexity of radiative transfer analysis and thus helping to retrieve the morphology of the nebula; (2) NaCl has a high dipole moment (μ = 9 D, Müller et al. 2005) enabling efficient radiative excitation and hence providing a link between the line emission and the radiation field throughout the nebula; (3) while the abundance of NaCl in the wind of VY CMa is estimated to be some 5 orders of magnitude lower than that of CO (Milam et al. 2007), its high dipole moment facilitates its detection; and (4) as most metal-bearing molecules NaCl is very refractory and is expected to quite easily form solid-state condensates so that one can investigate the gas-grain chemistry in the envelope. Understanding the gas-grain chemistry in this kind of oxygen-rich winds has a very broad application, since the abundances of the most chemically active species are quite close to the Galactic cosmic abundances and hence provides insight in the chemical processes at work on broad galactic scales.

In Sect. 2 we describe the ALMA observations and some details of the data reduction. Section 3 provides a detailed description of the spectral line profiles, the spatial distribution of the NaCl emission seen in the channel maps and a non-exhaustive comparison to other tracers of the morphological structure not yet mentioned in the introduction. In Sect. 4 we discuss the NaCl excitation conditions, construct a simplified chemical network aiming at unravelling the cycling between gaseous and solid NaCl (Sect. 4.1), and finally embark on an in-depth discussion on the morphology of the inner envelope based on the ALMA NaCl channel maps (Sect. 4.2). We summarize our conclusions in Sect. 5.

2. Observations

VY CMa was observed for ALMA science verification on 2013 16–19 August using 16–20 12-m antennas on baselines from 0.014–2.7 km. The main objective was to map the H₂O maser lines at 321 and 325 GHz (Band 7) and 658 GHz (Band 9), but several thermal lines identified with various rotational transitions of NaCl, TiO₂, SO₂, and SiO were also present in the spectral setting in addition to the continuum data. Richards et al. (2014) summarises the observations, data reduction and H₂O maser results, while O’Gorman et al. (2015) describes the dust-dominated continuum in more detail. The TiO₂ thermal emission lines are presented by De Beck et al. (2015). Two prominent components are visible in the continuum dust emission around 321 GHz. The brightest dust component “C” has a peak position at α = 07h22m58.34293s and δ = −25°46′03.̋2227, while the position of the peak of the second continuum component has been shown to coincide with the location of the stellar position “VY” at α = 07h22m58.326s and δ = −25°46′03.̋042 (Richards et al. 2014; O’Gorman et al. 2015) (see Fig. 2, Fig. 7).

Four NaCl lines in the ground or first vibrational state are covered in the ALMA band 7 data¹. All of them are detected, albeit only two are unblended (see Table 1). The NaCl v = 1J = 25−24 is slightly blended with the TiO₂(33_8,26−33_7,27) line at 322.613 GHz (De Beck et al. 2015), while the NaCl v = 0J = 25−24 lines is blended with the 325.153 GHz H₂O maser line (Richards et al. 2014) and the SiS v = 1J = 18−17 line at 325.059 GHz (see Sect. 3.1). No rotational lines in the ground or first vibrational state of the minor isotopologue Na³⁷Cl were covered in the current ALMA spectral setup in band 7. The NaCl v = 3J = 51−50 line at 657.134 GHz was also detected, but is severely blended with the much stronger H₂O maser line at 658 GHz, making a quantitative analysis impossible and hence it is not further discussed.

The synthesized beamsize for the band 7 data is (see Table 1). The stellar continuum was detectable in the individual channels at all frequencies, so that the relative alignment of all lines is good to ~2 mas (Richards et al. 2014). The continuum was subtracted before making the final line cubes. All data were adjusted to constant velocity with respect to the Local Standard of Rest.

Image cubes were made at 1 km s^-1 spectral resolution around the NaCl frequencies listed in Table 1. The higher σ_rms noise around 325 GHz is due to the poor atmospheric transmission in this region. The noise in the brightest channels is slightly higher in all cases due to dynamic range limitations. The shortest baseline was 14 m and inspection of the continuum visibility amplitudes against baseline length shows that the flux density remains quite steady out to 170 m at 321−325 GHz, suggesting that we recover all the flux on scales < 13′′, with reliable imaging on scales up to at least 8′′. The absolute astrometric accuracy of the lines described here is 35 mas and the flux density scale is accurate to at least 10%.

3. Results

3.1. Spectral line profiles

Fig. 1 Flux densities for the four NaCl lines detected in the ALMA band 7 data. Flux densities are for an aperture with circular diameter of 1′′ centred at “VY”. The LSR velocity of 22 km s-1 is indicated by the vertical dotted line. Line blending causes the flux density for the 325 GHz NaCl v = 0J = 25−24 line to be skewed and to be too high with an expected contribution of the SiS v = 1J = 18−17 line of ~1.1 Jy (see text for more details).

Figure 1 shows the flux densities extracted for a 1′′ diameter region around the stellar position “VY” for the four NaCl lines; no NaCl v = 1 emission is detected beyond this aperture. The four NaCl lines are centred around the v_LSR of ~22 km s^-1 (Richards et al. 1998, 2014; Tenenbaum et al. 2010; Kamiński et al. 2013b) and show similar line shapes ranging from ~6 km s^-1 to ~40 km s^-1. Given the similar energy levels, quantum numbers, and Einstein-A coefficients, one would expect similar line intensities and spatial distributions for both lines in the same vibrational level. While this holds for the NaCl lines in the first vibrational level, the 325 GHz NaCl J = 25−24 line in the vibrational ground state is ~1.4 Jy higher than the 312 GHz NaCl v = 0J = 24−23 peak. This appears to be due to blending. The 325.153 GHz water maser (with a ~500 Jy peak) has a very broad profile, with a red-shifted wing >70 km s^-1 and contributes roughly 0.2 Jy to the low-velocity wing of the 325.069 GHz NaCl line. The 325.059 GHz SiS v = 1J = 18−17 line is even closer to the NaCl line (peaks separated by 9.2 km s^-1). To estimate the line strength of the SiS v = 1J = 18−17, we use the results of Kamiński et al. (2013b) who made a spectral line survey of VY CMa between 279 and 355 GHz using the SMA. They were unable to detect the 325.059 GHz SiS line due to the poor atmospheric transmission but, using a 1′′ radius aperture, they measured maximum flux densities of 0.96 Jy for the SiS v = 1J = 16−15 line and 0.97 Jy for the J = 17−16 line. The estimated maximum flux for the SiS v = 1J = 19−18 is indicated to be uncertain, with a value around 1.38 Jy. From these results, we estimate that the SiS v = 1J = 18−17 line can contribute ~1.1 Jy to the flux density profile of the 325 GHz NaCl line, mostly on the higher-velocity side. This explains why this line is skewed to the red-shifted side of the stellar velocity. If 1.1−1.3 Jy is subtracted from the apparent 325 GHz NaCl line profile, the remaining peak would be close to the 0.6 Jy peak of the 312 GHz NaCl line. After correction for line blending, both v = 1 lines have flux densities within a factor of ~2 of the v = 0 lines despite an energy separation of ~360 cm^-1, which is caused by the dominantly radiative excitation mechanism.

3.2. Spatial distribution

Fig. 2 Channel maps (colour scale) of the NaCl v = 0J = 24−23 emission at 312 GHz at a 1 km s-1 velocity resolution for the central 0.̋8 × 0.̋8 field around “VY” and “C”. Grey contours are the integrated NaCl line strength at (1,2,3,4) × 1 Jy/beam km s-1, red contours are the dust continuum measured with ALMA at 321 GHz at (1,2,3,4) × 18 mJy/beam (Richards et al. 2014; O’Gorman et al. 2015). The stellar position “VY” is indicated with a white cross, and the peak of the continuum emission “C” is indicated with a white plus-sign. The ordinate and co-ordinate axis give the offset of the right ascension and declination, respectively, in units of milli-arcseconds.

Fig. 3 Same as Fig. 2 but for the NaCl v = 1J = 24−23 line at 309 GHz.

Fig. 4 Same as Fig. 2, but showing the central 2.̋5 × 2.̋5 region. The white contours show the HST/WFPC2 data in the F1042M filter around 1.018 μm at (1,2,4,8,16) × 110 counts (Smith et al. 2001).

We created moment-zero maps (total intensity) for each of the NaCl lines. The NaCl lines in a given vibrational level show similar spatial distributions. Since the NaCl J = 25−24 lines are blended in both vibrational levels, we discuss the spatial distribution of NaCl on the basis of the J = 24−23 rotational lines; the channel maps for the NaCl J = 25−24 lines are shown in Appendix A. A zoom for the inner 0.8′′ around the stellar position “VY” is shown in Figs. 2 and 3. In this inner 0.8′′ region, the spatial distribution of both lines show a similar bimodal distribution with two clumps at small offsets from the stellar position. For most of the channel maps, the northern clump (“NB” from here onward) is brighter than the southern clump (“SB” from here onward); see also Fig. 5 for an indication of the nomenclature used in this paper. In addition, both rotational NaCl lines in the ground-vibrational state display a third emission component located at ~1.2′′ to the South-West (“SW” clump); see Fig. 4.

For the NaCl lines in the first vibrational state, all the emission is concentrated in “NB” and “SB”; the other two lines in the ground vibrational state show two (one) “SW” component(s) at 312 (325) GHz. The two “SW” components for the NaCl v = 0J = 24−23 line at 312 GHz are called “SW1” and “SW2” (see Fig. 5), with “SW1” being brighter than “SW2”. We measured the sizes and positions of each spatial region by fitting two dimensional elliptical Gaussian components (see Table 2). The position uncertainties given are estimated conservatively from the noise as (mean beam size)/(signal to noise ratio) but do not include the 2 mas line-to-line alignment uncertainty (the 35 mas absolute astrometric uncertainty does not affect alignment of these NaCl lines). The component sizes are the FWHM (full width at half maximum) deconvolved from the beam where possible; only “NB” is consistently fully resolved and the errors given take into account that the sizes of the other components are in some cases upper limits.

The brightest “NB” component is aligned to within 8 mas in both coordinates at 309, 312 and 322 GHz, which is within the 3 sigma position uncertainties. The fainter “SB” component is aligned to within 40 mas; although this could be significant it is more probably due to fitting uncertainties in the case of a non-Gaussian flux distribution. The average separations of “NB” and “SB” from the stellar position, “VY”, are 100 mas (or ~18 R_⋆) and 190 mas (or ~34 R_⋆), respectively. The separation between “NB” and “SB” is 265 mas (~48 R_⋆). The relative uncertainty in each case is 10 mas (~2 R_⋆). The brighter “SW” component, “SW1”, seen at 312 GHz is 1.22′′ (~220 R_⋆) from VY; the fainter “SW2” component is 0.3′′ further away. Formally, the angular separations are a lower limit to the physical separation, but since the NaCl clumps’ velocities span the stellar velocity, these are probably close to the actual separations.

“NB” and “SB” and a single “SW” clump are seen in the NaCl v = 0J = 25−24 line at 325 GHz. The uncertainties listed in Table 2 are a few times greater, due both to the noise level which is 6−10 times higher and blending with the other lines, in particular with the water maser which is widely distributed around the“NB” and “SB” clumps and beyond but not towards the “SW” clump.

The northern clump is visible in the channel maps for all observed NaCl transitions with velocities between 6 and 37 km s^-1. The flux density in “NB” shows two local maxima, one around 14 km s^-1 (“NB_blue”) and the other around 24 km s^-1 (“NB_red”); see Fig. 6. Emission from “SB” arises in channel maps with velocities between 12 and 26 km s^-1. The flux density of “SB” reaches a maximum around 19 km s^-1, i.e. coinciding with a local minimum in the flux density of “NB”. The ’SW’ extension in the ground vibrational lines is visible for velocity channels between 12−33 km s^-1 (see Fig. 4) with a maximum flux density around 25 km s^-1 (see Fig. 6).

3.3. Morphological comparison to other tracers of the envelope

In this section we compare the NaCl emission seen in the ALMA channel maps to some other results presented in the literature. We only focus on observational results already retrieved from ALMA data and on studies which discuss the NaCl emission and/or the extension seen at 1.2′′ to the south-west.

The maximum known angular extent of the molecular envelope around VY CMa as traced by CO (Muller et al. 2007; Kamiński et al. 2013b) and OH maser emission (Richards 1997) is ~10′′ and the typical velocity span is from −15 to + 60 km s^-1, with several species observed closer to the central star showing a more extended red wing up to + 105 km s^-1 (see, e.g. Kamiński et al. 2013b; De Beck et al. 2015). However, where imaging is available, these extreme velocities appear at only a few position angles, and most of the molecular shell appears to show radial acceleration out to a terminal velocity of 35–45 km s^-1 (Decin et al. 2006; Muller et al. 2007; Ziurys et al. 2007)². The NaCl emission detected by ALMA band 7 has hence a much more compact distribution and while its emission extends significantly beyond the estimated dust formation radius of 35−55 mas (Danchi et al. 1994; Monnier et al. 1999; Smith et al. 2001; Decin et al. 2006) it does not reach the terminal velocity.

Figures 2, 3, A.1, and A.2 compare the NaCl emission to the dust emission detected in the ALMA band 7 band. The emission of the brightest NaCl components (“NB” and “SB”) and the dust shows a similar central “waist” along a NE-SW axis, albeit the peak position of the NaCl clumps is to the north of “VY” and NW of “C” and are clearly offset from the dust continuum peaks.

As shown by Richards et al. (2014) the H₂O 658, 321, and 325 GHz maser emission cluster around the stellar position “VY”, with another accumulation in the “valley” between “C” and “VY” (see Fig. 7). Richards et al. (1998) found the 22 GHz H₂O masers to show a generally expanding envelope, with two “ears” outside the main 22 GHz shell, which tentatively was interpreted as being suggestive for a biconical outflow. The south clump (“SB”) of the NaCl emission lies in the “valley” between the two continuum peaks, to the south of some of the most intense water maser emission. The 22 GHz and 321 GHz H₂O maser features at the greatest angular separations lie approximately east and west or to the southwest.

Fig. 5 ALMA NaCl 312 GHz channel map at 25 km s-1 (colour scale) overlaid with the HST/WFPC2 F1042M data of Smith et al. (2001) (in contours at (1,2,4,8,16) × 110 counts). The cross sign indicates the position of the central star VY CMa, noted “VY” in the text, as derived from the ALMA H2O maser and continuum data (Richards et al. 2014; O’Gorman et al. 2015), the plus-sign indicates the maximum of the ALMA dust continuum data, noted as “C” in the text, (O’Gorman et al. 2015). The nomenclature of the different spatial regions refers to the northern clump ’NB’, the southern clump ’SB’, the brightest part of the south-western extension ’SW1’, and the faintest part of the north-western extension“SW2”. The ordinate and co-ordinate axis give the offset of the right ascension and declination, respectively, in units of arcseconds.

Fig. 6 Flux densities for the NaCl J = 24−23 lines in the ground-vibrational state v = 0 at 312 GHz (upper panel) and in the first vibrational state v = 1 at 309 GHz (lower panel). Flux densities are extracted in the north clump (NB, full line) at (07h22m58.326s, −25°46′02.̋9) for an aperture of (0.̋394 × 0.̋373), in the south clump (SB, dotted line) at (07h22m58.332s, −25°46′03.̋118) for an aperture of (0.̋36 × 0.̋36), and in case of the 312 GHz line in the south-west clump (SW, dashed line) at (07h22m58.273s, −25°46′04.̋028) for an aperture of (0.̋57 × 0.̋73). The LSR velocity of 22 km s-1 is indicated by the vertical dotted line.

Fig. 7 The H2O maser emission (Richards et al. 1998, 2014) is shown on top of the NaCl 312 GHz integrated line emission (gray contours, zeroth moment at 0.25, 0.5, 1, 2, 3, 4 Jy km s-1) and the ALMA dust continuum emission (colour scale; O’Gorman et al. 2015).

Milam et al. (2007) observed several rotational lines in the ground vibrational state of NaCl using single-dish telescopes with beam sizes between 29′′−44′′. The NaCl emission was spatially unresolved and modelling the observed spectra resulted in an abundance relative to H₂ of 5 × 10^-9 with a source size of 0.5′′. Kamiński et al. (2013b) observed several rotational lines of NaCl in vibrational states with v ≤ 3 at an angular resolution of ~0.9′′ using the SMA. The most intense lines were partially resolved, with the emission being slightly elongated with mean deconvolved FWHMs of and a position angle of −33 ± 2°. Weak emission in the SW clump at ~1′′ is seen only in the most intense line of NaCl and could not be mapped with the SMA. Thanks to the high spatial resolution of ALMA, the NaCl emission is for the first time well resolved and a clear distinction can be made between emission seen in the bright northern clump, the slightly weaker southern clump, and the weak south-western extension (see Table 2 in Sect. 3.2), the latter only detected for the NaCl lines in the ground-vibrational state. Based on the NaCl v = 0J = 24−23 transition – the only unblended line for which the emission traces all the different spatial components – the northern clump has a deconvolved size of , the upper limit for the size of the southern clump is , while the two components in the south-western extension have sizes with upper limits being and (see Table 2).

A similar “SW” extension/clump at ~1′′ for a position angle (PA) of ~220° w.r.t. “VY” is seen in some ALMA TiO₂ channel maps (De Beck et al. 2015) and is also reported by Kamiński et al. (2013b) for some sulfur bearing species (CS, NS, SiS, and H₂S) and the most intense NaCl (J = 26−25) line. The line profile of the SMA H₂S and CS emission in the “SW” clump seen is rectangular with FWHM of 15 km s^-1. The emission of this “SW” component very well matches the south-west extension seen in the Hubble Space Telescope (HST) data at wavelengths between 1–2.14 μm that image the reflection nebula around VY CMa (see Fig. 5; Smith et al. 2001) and with other near-infrared and polarimetric data as for instance presented by Monnier et al. (1999) and Jones et al. (2007). The direction of this southwest extension/knot is roughly perpendicular to the possible density-enhanced equatorial plane as deduced from the one-sided reflection nebula seen in the visible and the maser elongation (see Sect. 1). More recently, Shenoy et al. (2013) re-confirmed the existence of a peculiar south-west clump at approximately 1.4′′ to the southwest of the star using LBT/LMIRCam data. Humphreys et al. (2007) studied the total space motions and directions of several features in the envelope of VY CMa. They deduced that the southwest clump is slightly redshifted w.r.t. the systemic velocity and is moving away from us at an angle of + 8° with a total velocity of only ~18 km s^-1, of which the transverse component is ~17.7 km s^-1 and with a negligible radial component of ~2.5 km s^-1. This derived redshifted behaviour is in accord with the ALMA NaCl flux density shown in Fig. 6 which shows a peak around 25 km s^-1, or 3 km s^-1 greater than the stellar velocity of 22 km s^-1. However, this is not in agreement with the ALMA TiO₂ emission showing up blue shifted (De Beck et al. 2015) (see Sect. 4.2.2).

4. Discussion

Metal-bearing molecules are rarely observed in the gas-phase in the ISM due to their large refractory character that makes them to easily form solid condensates. Metal halides are among the first metal-containing molecules detected in space. The alkali halide salt molecule, NaCl, was already detected in VY CMa and in the oxygen-rich AGB star IK Tau (Milam et al. 2007; Kamiński et al. 2013b). The high dipole moment of NaCl is favourable for radiative excitation to play a significant role in populating the (vibra-)rotational levels. For Einstein A-values in the range of 0.013−0.019 s^-1 (see Table 1) and assuming collision rates within v = 0 around 10^-10 cm³ s^-1 (based on the collisional rates of SiS, Agúndez et al. 2012), one gets a critical density n_crit around 2 × 10⁸ cm^-3. Using the output from the thermodynamical model for the wind of VY CMa as presented by Decin et al. (2006), this implies that there is only a small region for r ≲ 9 R_⋆ (~0.05′′) where the gas number density is higher than the critical density and hence where LTE would prevail.

Milam et al. (2007) determined the formation locus of gaseous sodium chloride to be close to the central star, between ~5−40 R_⋆, with a fractional abundance relative to H₂ of 5 × 10^-9. This is in accord with thermochemical equilibrium predictions, not taking phase changes into account, predicting that the gaseous NaCl abundance peaks around 500 to 1200 K at densities n ~ 10¹¹ cm^-3 (see Fig. 3 in Milam et al. 2007), with NaCl being the main sodium carrier in oxygen-rich environments³. For temperatures above 1500 K, the fractional abundance of gaseous NaCl is below 10^-12. Kamiński et al. (2013b) derived the excitation temperatures and column densities for the molecules detected in their SMA scan of VY CMa from a population diagram analysis, i.e. assuming an isothermal gas in local thermodynamic equilibrium. They determined a column density for NaCl of 1.3 ± 0.2 × 10¹⁵ cm^-2 for a rotational excitation temperature of 210 ± 25 K; the temperature derived from transitions including vibrationally excited levels is a factor of a few larger and the column density a factor of 2–3. Our source size implies that column densities derived by Kamiński et al. (2013b) should be increased by a factor of ~3, i.e. N_rot ≈ 5 × 10¹⁵ cm^-2. A mass loss rate of 10^-4 M_⊙/yr and a velocity of ~15 km s^-1 (as derived in the NB) yields a fractional NaCl abundance relative to H₂ of ~ 4 × 10^-8 in the northern blob. We note, however, the very provisional nature of this number. The lack of information on cross sections, physical conditions, radiative pumping efficiency, non-LTE effects, etc. means that extreme caution should be taken when using this value of the column density and NaCl abundance which is likely uncertain by about a factor of 10.

The confined distribution of gaseous sodium chloride, the fact that the observed velocities do not reach the terminal velocity, and realising the refractory nature of sodium chloride hint toward the fact that NaCl condenses easily onto grains. Thermochemical calculations including gas chemistry and the formation of possible liquid and solid condensates predict that sodalite, Na₄[AlSiO₄]₃Cl, is the major silicated chlorine containing condensate in an oxygen-rich environment, while halite, crystalline NaCl salt, might be formed in a carbon-rich wind (Lodders & Fegley 1999). However, the fact that the ALMA data prove the existence of gaseous NaCl at distances far beyond the major dust condensation region, specifically around 1.22′′ (~220 R_⋆) in the “SW” extension, is an incentive to reconsider the different gas-grain chemical processes at work that might prevent all gaseous NaCl species to be condensed.

A simplified NaCl chemical network is developed in Sect. 4.1 with the aim to constrain the main chemical processes that control the gaseous and solid NaCl content. Potential geometrical structures as derived from the multi-component NaCl emission are discussed in Sect. 4.2.

4.1. NaCl gas-grain chemistry

Our aim is to understand the cycling between gaseous and solid NaCl through the envelope of VY CMa. Since many parameters such as the sticking coefficient, surface binding energy, and the surface density of adsorbing sites are not known, we have constructed a simplified chemical network that lets us control several parameters that directly influence the NaCl abundance profile. Analogous simplified networks are constructed to study, for instance, the water content in low-mass forming regions (Schmalzl et al. 2014). The main chemical processes to be considered are accretion of gaseous NaCl onto grains, thermal desorption and photodesorption.

4.1.1. Simplified NaCl network

The rate of accretion of gaseous NaCl onto existing dust grains is given by ${\mathit{R}}_{\mathrm{ACC}}\mathrm{=}\mathit{S}{\mathit{\sigma }}_{\mathrm{d}}\mathrm{⟨}{\mathit{v}}_{\mathit{X}}\mathrm{⟩}{\mathit{n}}_{\mathrm{g}}\mathrm{\left(}\mathrm{NaCl}\mathrm{\right)}$(1)in units of [cm^-3 s^-1], where S is the sticking probability, σ_d the total grain surface cross-section per unit volume in units of [cm² cm^-3], ⟨v_X⟩ is the average thermal velocity $\sqrt{\mathrm{(}\mathrm{8}\mathit{k}{\mathit{T}}_{\mathrm{g}}\mathit{/}\mathit{m}\mathrm{)}}$ in units of [cm s^-1] with T_g the gas kinetic temperature, k the Boltzmann constant, m the mass of gaseous NaCl (58.5 atomic mass units), and n_g(NaCl) the number density of gaseous NaCl in units of cm^-3.

We can write the total grain surface cross section per unit volume as ${\mathit{\sigma }}_{\mathrm{d}}\mathrm{=}\mathrm{\int }\mathrm{\left(}\mathit{\pi }{\mathit{a}}^{\mathrm{2}}\mathrm{\right)}\hspace{0.17em}{\mathit{n}}_{\mathrm{d}}\mathrm{\left(}\mathit{a,r}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{a,}$(2)with n_d(a,r) being the grain volume density for a grain with radius a. We adopt a grain size distribution n_d(a,r)da = Ka^-3.5n_H(r) da, with a minimum grain size, a_min, of 0.005 μm and a maximum grain size, a_max, of 0.25 μm. The slope of −3.5 is typical for interstellar grains (Mathis et al. 1977). n_H is the total hydrogen number density in units of [cm^-3] and K represents an abundance scaling factor giving the number of dust particles in units of particles per H atom. For the interstellar medium, K is estimated to be around 7.9 × 10^-26 cm^2.5/H for silicate grains (Draine & Lee 1984). For a dust-to-gas mass ratio of 0.002 (Decin et al. 2006), K is ~2.8 × 10^-26 cm^2.5 per H atom (see Appendix B). For the assumed grain size distribution, the total grain surface cross section σ_d can hence be written as σ_d ~ 2500 πKn_H(r) or σ_d ~ 2.2 × 10^-22n_H(r).

The number density of gaseous NaCl, n_g(NaCl), is written as $\begin{array}{ccc}{\mathit{n}}_{\mathrm{g}}\mathrm{\left(}\mathrm{NaCl}\mathrm{\right)}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\mathit{X}\mathrm{\left(}\mathrm{NaCl}\mathrm{\right)}\hspace{0.17em}{\mathit{n}}_{\mathrm{H}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathit{.}& & \end{array}$Assuming solar elemental abundance of chlorine (~3.2 × 10^-7) and sodium (~2.1 × 10^-6) the abundance of NaCl, X(NaCl), is controlled by the chlorine content.

The rate of thermal desorption is given by (Tielens & Hagen 1982) ${\mathit{R}}_{\mathrm{TD}}\mathrm{=}{\mathit{n}}_{\mathrm{s}}\mathrm{\left(}\mathrm{NaCl}\mathrm{\right)}\hspace{0.17em}{\mathit{\nu }}_{\mathrm{0}}\hspace{0.17em}\exp \mathrm{(}\mathrm{-}{\mathrm{E}}_{\mathrm{b}}\mathit{/}\mathrm{k}\hspace{0.17em}{\mathrm{T}}_{\mathrm{d}}\mathrm{)}$(3)in units of [cm^-3 s^-1], where n_s(NaCl) is the number density of NaCl molecules on the grains per cm³, ν₀ a characteristic vibration frequency in units of [s^-1], E_b the surface binding energy, and T_d the dust grain temperature.

To get the characteristic vibration frequency, ν₀, we adopt the harmonic oscillator approach (Hasegawa et al. 1992) ${\mathit{\nu }}_{\mathrm{0}}\mathrm{=}\sqrt{\mathrm{2}\hspace{0.17em}{\mathit{N}}_{\mathrm{s}}\hspace{0.17em}{\mathit{E}}_{\mathrm{b}}\mathit{/}{\mathit{\pi }}^{\mathrm{2}}\hspace{0.17em}\mathit{m}}\mathit{,}$(4)where N_s is the surface density of adsorbing sites on the grains and m is the mass of the adsorbed NaCl particle. The number of adsorbing sites per unit area, N_s, is assumed to be ~ 1.5 × 10^-15 cm^-2 (Hasegawa et al. 1992). Hence ν₀ is given by ${\mathit{\nu }}_{\mathrm{0}}\mathrm{=}\mathrm{2.07}\mathrm{\times }{\mathrm{10}}^{\mathrm{12}}\sqrt{\frac{{\mathit{T}}_{\mathrm{b}}}{\mathrm{10}\hspace{0.17em}\mathrm{000}}}\mathit{,}$(5)with T_b the surface binding temperature in Kelvin. This is in good accord with the observed frequency of 336 cm^-1 (or 1 × 10¹³ s^-1) listed by Martin (1983) or the value of 380 cm^-1 (or 1.13 × 10¹³ s^-1) calculated using the Born-Mayer lattice theory by Rittner (1951).

One has to realize that even if all Cl is in NaCl and if all NaCl molecules start off in the gas and accretes onto already formed dust grains, NaCl will only occupy at most 25% of the available surface sites, i.e. NaCl is in the sub-monolayer regime. The reason for that is that the number of sites per unit volume is given by 4σ_dN_s or 1.3 × 10^-6n_H(r) with the cosmic abundance of Cl being 3.2 × 10^-7n_H(r). As shown by Tielens & Hagen (1982) quantum mechanical tunneling and thermal hopping of the adsorbed species on the grain surface are not relevant for NaCl.

For sub-monolayer coverage NaCl can photodesorb from the grain surface through far-ultraviolet (FUV) photons at a rate ${\mathit{R}}_{\mathrm{PD}}\mathrm{=}\mathrm{(}{\mathit{F}}_{\mathrm{FUV}}{\mathit{Y}}_{\mathrm{PD}}\mathit{/}\mathrm{4}{\mathit{N}}_{\mathrm{s}}\mathrm{)}{\mathit{n}}_{\mathrm{s}}\mathrm{\left(}\mathrm{NaCl}\mathrm{\right)}\mathit{,}$(6)in units of [cm^-3 s^-1]. F_FUV is the flux of FUV photons reaching the grain and Y_PD the photodesorption yield.

The photodesorption yield Y_PD is difficult to assess. For grains for which water ice is the dominant mantle component, water molecules are photodesorbed either through the excess energy caused by the recombination of H and OH radicals following direct photodissociation or via the “kick-out” mechanism where an excited H atom ejects an H₂O molecule upon collision (Arasa et al. 2015). In VY CMa Shenoy et al. (2013) derive a grain temperature of ~130−170 K at the “SW” clump, making it marginal whether ice mantles can be present – thermal desorption of water occurs at 160 K (Collings et al. 2004). If no water ice is present, the yield is likely to be much less than 3 × 10^-3 which is a typical value for molecules in water ice (Öberg et al. 2009a,b; Westley et al. 1995)

Although it is commonly thought that the internal FUV radiation field of RSG stars is negligible, this is not necessarily the case for material sitting close to the stellar photosphere. To calculate the flux at the “SW” clump, we use an effective stellar temperature, T_⋆, of 3480 K, the photodissociation cross-sections of NaCl measured by Silver et al. (1986), and a distance of 250 R_⋆ to the “SW” clump. For wavelengths shorter than 2910 Å, the threshold wavelength for NaCl photodissociation, we find that the total FUV flux is roughly 10⁴ times larger than the general ISM flux (Draine 1978) and calculate the unshielded photodissociation rate to be 3.2 × 10^-5 s^-1, which is a lower limit to the rate given that the cross-section has not been measured below 1890 Å and some 10⁴ times larger than the unshielded ISM rate (van Dishoeck 1988).

The change of the gaseous and solid NaCl number densities is then given by $\frac{\mathrm{d}{\mathit{n}}_{\mathrm{g}}\mathrm{\left(}\mathrm{NaCl}\mathrm{\right)}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathrm{-}\frac{\mathrm{d}{\mathit{n}}_{\mathrm{s}}\mathrm{\left(}\mathrm{NaCl}\mathrm{\right)}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathrm{-}{\mathit{R}}_{\mathrm{ACC}}\mathrm{+}{\mathit{R}}_{\mathrm{TD}}\mathrm{+}{\mathit{R}}_{\mathrm{PD}}\mathit{.}$(7)and we make the assumption that $\mathit{X}\mathrm{\left(}\mathrm{Cl}\mathrm{\right)}{\mathit{n}}_{\mathrm{H}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathit{n}}_{\mathrm{g}}\mathrm{\left(}\mathrm{NaCl}\mathrm{\right)}\mathrm{+}{\mathit{n}}_{\mathrm{s}}\mathrm{\left(}\mathrm{NaCl}\mathrm{\right)}$(8)at any time t with X(Cl) the cosmic abundance of chlorine relative to hydrogen.

4.1.2. Chemical network model results

Fig. 8 Steady state solutions for the gaseous and solid NaCl abundances (relative to the hydrogen number density) as a function of the dust temperature, Td, for the case in which accretion, thermal desorption and photodesorption are included in the chemical network. The different panels are for different values of the binding temperature (5000, 10 000, 15 000 K) at a hydrogen number density of 1 × 108 cm-3. Each panel also shows the inverse timescale for accretion, thermal desorption and photodesorption in units of s-1, calculated with FFUVYPD = 3 × 105.

Fig. 9 Fractional abundances of gaseous and solid NaCl in function of time for a binding temperature of 15 000 K at three different values of the dust temperature Td (290, 300, 310 K).

The value for the surface binding energy, T_b, is not known. It depends on the polarizability of NaCl and its ability to form specific bonds (Tielens & Hagen 1982) to the surface. In the case of H₂O with its strong dipole moment and its ability to from hydrogen bonds, the surface binding energy on an ice surface has been calculated to be about 4000 K (Hale et al. 1981). On the other hand, the surface binding energy of Na to a NaCl polycrystalline target is about 1.3 eV (or ~15 000 K). Considering this, we have decided to study the evolution of the NaCl abundance for three different values of the surface binding energies: 5000 K, 10 000 K, and 15 000 K.

It is instructive to follow the abundance of gaseous and solid NaCl by solving Eq. (7) for n_s(NaCl). Using conservation of NaCl at any time t (Eq. (8)) and writing σ_d(r) = 2.2 × 10^-22n_H(r) = σn_H(r), the rate equation can be rewritten as $\frac{\mathrm{d}{\mathit{n}}_{\mathrm{s}}\mathrm{\left(}\mathrm{NaCl}\mathrm{\right)}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathit{a}\mathrm{-}\mathit{b}\hspace{0.17em}{\mathit{n}}_{\mathrm{s}}\mathrm{\left(}\mathrm{NaCl}\mathrm{\right)}\mathit{,}$(9)where $\mathit{a}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathit{S\sigma }{\mathit{n}}_{\mathrm{H}}^{\mathrm{2}}{\mathit{v}}_{\mathit{X}}\mathit{X}\mathrm{\left(}\mathrm{Cl}\mathrm{\right)}$ and b = Sσn_Hv_X + ν₀exp(−T_b/T_d) + (F_FUVY_PD/ 4N_s). This has a simple time-dependent solution ${\mathit{n}}_{\mathrm{s}}\mathrm{\left(}\mathrm{NaCl}\mathit{,t}\mathrm{\right)}\mathrm{=}\frac{\mathit{a}}{\mathit{b}}\left(\mathrm{1}\mathrm{-}\exp \mathrm{(}\mathrm{-}\mathit{bt}\mathrm{)}\right)\mathit{.}$(10)n_s(NaCl,t = 0) = 0 and the time taken to reach steady state is a few times 1 /b. The steady-state solution is n_s(NaCl, ss) = a/b and is relatively insensitive to the value of S for S ~ 1.

It is difficult to determine the product F_FUVY_PD with any degree of certainty. The flux of FUV photons at the “SW” clump is very large if there is no extinction due to dust along the line-of-sight to the star. If this really is the case, then NaCl is photodesorbed completely on a very fast time-scale and the NaCl released from the dust will be photodissociated rapidly, on a time-scale of a few days. In this case, NaCl is very unlikely to reform in the gas phase given that both Na and Cl are very minor elements. This scenario is unlikely, however, as the dust condensation radius lies much closer to the star, on the order of 35–55 mas, see Sect. 3.3. Using our derived dust-to-gas mass ratio, we calculate that the photosphere is hidden by some 10⁴ mag of visual extinction. We note, however, that even 10 mag of visual extinction is enough to reduce the unshielded FUV flux to interstellar values. Thus stellar FUV photons are unlikely to play a role in NaCl chemistry.

It is instructive, though, to consider a calculation for which the value of F_FUVY_PD = 3 × 10⁵, the typical interstellar value. Assuming a sticking coefficient, S, of 1, the steady-state solution for the gaseous and solid NaCl fractional abundances (with respect to the hydrogen number density) are shown in Fig. 8 for different values of the binding energy. The hydrogen number density in Fig. 8 is taken to be 1 × 10⁸ cm^-3. This value is derived from the equation of mass conservation Ṁ = 4 πr²ρ_g(r) v_g(r), with ρ_g the gas number density and v_g(r) the gas velocity. Using a mass-loss rate Ṁ of 1 × 10^-4 M_⊙ yr^-1 and the gas velocity structure as derived from solving the momentum equation (see Decin et al. 2006), a hydrogen number density of 1 × 10⁸ cm^-3 is reached at ~10 R_⋆, i.e. the dust condensation locus of interest to this study. Figure 8 also displays the inverse timescale for accretion and thermal desorption. Figure 8 shows that photodesorption occurs at a slower rate than accretion and that only about 10% of the NaCl abundance resides in the gas phase at low dust temperatures (relative to the surface binding temperature). For example, for a binding energy of 10 000 K, once the grain becomes warm enough, around 200 K, all solid-state NaCl is thermally desorbed into the gas phase. It is clear from these simulations that the crucial factor determining the amount of solid versus gaseous NaCl is the ratio of the binding energy to the dust temperature, which appears in Eq. (3) for the thermal desorption. For a binding temperature of 5000 K all NaCl will be gaseous unless the dust temperature is below ~120 K. A higher binding energy implies that the dust grains may attain higher temperatures before thermal desorption starts to dominate over the accretion rate. The turn-over point is around T_b/T_d ~ 50 (see also Fig. 9). Close to this value, steady state solutions occur on timescales of 50−100 yr as shown in Fig. 9.

At the dust condensation locus of ~10 R_⋆, the grain temperatures might be in the range of 500−1000 K, so that all NaCl will be in gaseous form independent of the assumed binding energy. This is in line with the ALMA observations presented in this paper. A steady state solution of ~100 yr for a velocity of ~15 km s^-1 (as derived from the ALMA NaCl data) yields a travel distance of 30 R_⋆ (or 160 mas). When the grain temperatures becomes lower than 300 K, solid NaCl starts to form from accretion onto existing grains. Whenever T_d< 100 K, another chemical process than thermal desorption needs to be invoked for NaCl to remain (partly) in gaseous form. The ALMA data clearly show gaseous NaCl emission around 1.22′′ in the south-western extension, at a place where also scattering by dust grains peaks (see Fig. 5). Assuming only direct stellar radiation, Shenoy et al. (2013) derived a radiative equilibrium temperature for the grains at that distance of ~130−170 K. If the grains indeed have that temperature at 1.22′′, this would imply that either the binding temperature is around 5000 K so that NaCl stays gaseous that far away from the star or that, in the case of higher values for the binding temperature, another process is important. One possibility is shock-induced sputtering driven by localised mass ejections (see Sect. 4.2.2. Such a process could lead, for example, to the different spatial distributions of refractory molecules, such as observed for NaCl and TiO₂, since their sputtering rates likely depend on their binding energies.

4.2. Morphology

The NaCl emission does not peak at the stellar position “VY”, which can be explained by the presence of a density cavity close to the star, where the temperature is too high for gaseous NaCl to exist in high enough amounts (see previous sections). The morphology in the channel maps is clearly inconsistent with a spherical, uniformly expanding, envelope, which would exhibit a circular or ring morphology with decreasing radius for increasing velocity difference w.r.t. systematic velocity. Previous studies hint toward a preferred axis of symmetry oriented northeast-southwest at a position angle of ~50° (see Sect. 1). The central waist in the integrated NaCl intensity map aligns with the central waist seen in the dust continuum map and, as seen in Figs. 4 and 5, the ALMA data also roughly support the idea of a preferred axis oriented northeast-southwest.

Deducing the morphology of the circumstellar environment of VY CMa is not a straightforward task. However, the ALMA NaCl data seem to favour two types of geometries: (1) an axisymmetric geometry, be it a disk, a bipolar outflow, or in general a kind of toroidal overdensity or (2) localized mass ejections. We consider the merits of both geometries regarding the morphology seen in the ALMA data and other optical, infrared and maser data published in the literature.

4.2.1. Axisymmetric geometry

Within 1′′ of the central star, the ALMA NaCl data suggest an axial symmetry with a preferred axis oriented NE/SW at a PA of ~50° (defined North-East), showing both a bright NaCl and dust continuum peak to either side of the symmetry axis. This kind of spatial-kinematic structure can be interpreted as an equatorially enhanced density structure with the polar axis, hence lower density region, along the NE/SW axis at PA of ~50°.

Different studies (see Sect. 1 and below) support the axisymmetry revealed by the ALMA data. Herbig (1970) invoked a flared disk or ring in latitudes 10° to 30° and seen almost edge-on to explain the infrared energy distributions. Using optical and near- and mid-infrared data Smith et al. (2001) gave support to that interpretation and qualitatively interpreted the morphological structure as being caused by an equatorial density enhancement with a bipolar axis projected roughly northeast-southwest, and with the southwest lobe closer to us. The lack of optical and near-infrared emission in the northeast region of their data could be explained by high extinction along the line-of-sight and inefficient backscattering of the light by the dust particles. The SMA data of ¹²CO and ¹³CO at a spatial resolution of ~2′′ were modelled by Muller et al. (2007) using a dense, compact and dusty central component, embedded in a more diffuse extended envelope and with a high-velocity bipolar outflow oriented in the east-west direction with a wide opening angle (~120°) viewed close to the line-of-sight (i = 15°). The inner radius of the bipolar lobes was determined at ~0.6′′ and the velocity in the bipolar lobes increased linearly from 15 km s^-1 to 45 km s^-1 at the outer radius. An enhanced density shell within the bipolar lobes was needed to fit the observed CO data.

Fig. 10 Position-velocity diagram of the NaCl 312 GHz for an axis length of 1′′ and an averaging width of 650 mas centred at α = 07h22m58.332s and δ = −25°46′03.̋053. In the left and middle panels, the slit has a PA of 50° (i.e. along the presumed polar/rotation axis; with north-east positions having negative offsets and south-west positions positive offsets), in the right panel the slit PA is 140° (i.e. along the presumed equatorial plane; with south-east positions having negative offsets and north-west positions positive offsets). In the left and right panels the different velocity components are indicated, while in the middle and right panels the main contributors to the different intensity peaks are indicated. The contours in the left panel are at [1,2,...10] × 0.01 Jy/beam.

Two position-velocity (PV) diagrams of the NaCl 312 GHz ALMA observations are shown in Fig. 10; the panel to the left (PV1) gives the PV diagram for a slit oriented along the presumed polar/rotation axis at a PA of 50° to the north-east and the panel to the right shows the PV diagram (PV2) for a slit perpendicular to the polar/rotation axis, i.e. along the presumed equatorial plane.

PV1 diagram is characterized by a S-type signature. The northern clump is the main contributor to the bright peak around 24 km s^-1 (see middle panel in Fig. 10), while both the northern (NB_blue) and southern clump contribute to the blue-shifted peak around 19 km s^-1. The blue-shifted extension around 14 km s^-1 arises from NB_blue, while NB_red explains the extension toward velocities of ~35 km s^-1. It is tempting to interpret the two horizontal ridges at Δv ~ ± 2.5 km s^-1 from the systemic velocity as arising from a disk (or in general a density enhanced equatorial region) inclined to the line-of-sight (hereafter referred to as the low-velocity component; see left panel in Fig. 10). The disk has a projected velocity around 2.5 km s^-1. The blue-shifted part of the disk is to the south-west, and the red-shifted part to the north-east. There is also a more or less vertical ridge extending to ±15 km s^-1 w.r.t. the systemic velocity, with 2 other horizontal ridges at Δv ~ ± 9 km s^-1 from the systemic velocity (i.e. at 14 and 31 km s^-1). This medium-velocity component⁴ is only seen in the northern clump, with the blue-shifted part slightly more to the north-east than the red-shifted part, i.e. an orientation opposite to the low-velocity disk. This tentatively hints towards an interpretation that the medium-velocity component is a biconical polar outflow, with wide opening angle, perpendicular to the equatorial disk plane. An analogous morphology was also deduced for the carbon star V Hya by (Hirano et al. 2004) and the S-type AGB star π¹ Gru by Chiu et al. (2006) from the analysis of ¹²CO J = 2−1 SMA data, showing similar PV diagrams as PV1. With the intensity of “NB” being brighter around 31 km s^-1 than at 14 km s^-1, this suggests that the far-side of the bipolar outflow is situated to the north. This orientation supports the analysis by Smith et al. (2001), who suggested that the far side of the optical/near-infrared nebula is situated to the north-east, while the southwest lobe is closer to us. A geometry with an equatorial density enhancement roughly northwest-southeast and with the polar lower density region oriented along a NE/SW axis can also explain the morphology of the dust continuum as seen in the ALMA data and shown in Fig. 2 (O’Gorman et al. 2015).

However, one then might wonder why the medium-velocity wind (aka bipolar outflow) is only seen in the northern clump (see also PV2), if the southwest lobe is closer to us. In the case of optically thin emission, it might be that the accumulation of optical thickness along the line sight is too low when reaching the southern clump and that the sensitivity of the present ALMA data is not high enough to trace these low intensities. Another suggestion is that for one reason or another (see Sect. 4.1) gaseous NaCl condenses rapidly in the southern lobe and that its abundance is too low to be traced.

However, this axisymmetric geometry can not explain the full morphological picture of the envelope around VY CMa. As clearly shown by the ALMA detection of NaCl emission in the “SW” clump and as already mentioned in this paper, there are convincing indications of localised density structures that might be caused by randomized mass ejections from the ill-defined convective surface of the star. However, even in that case, the proposed lower density polar regions in the northeast-southwest direction might explain why we see the emission from the NaCl clump at 1.2′′ to the southwest: this clump, that might be caused by a higher density mass-loss event, can receive more illumination from the central star through the low density region which will enhance its excitation. This effect was also invoked by De Beck et al. (2015) to explain the detection of TiO₂ in the “SW” clump.

4.2.2. Localized mass ejections

The present ALMA NaCl data can also be geometrically interpreted as showing evidence for localized mass ejections. Episodic mass loss events were already suggested by, e.g. Ziurys et al. (2007) and Smith et al. (2009); see Sect. 1. As described in Sect. 3.2, the flux density of the “NB” shows two local maxima, one around 14 km s^-1 (“NB_blue”) and one around 24 km s^-1 (“NB_red”). The minimum around 18 km s^-1 coincides with the maximum in flux density of the southern clump (“SB”). Careful inspection of the channel maps in Fig. 2 also shows that (albeit barely spatially resolved) the morphology of “NB_blue” has a more pronounced north-south extension than “NB_red”. While the TiO₂ emission coinciding with the “SW” clump is blue-shifted (De Beck et al. 2015), the NaCl channel maps clearly show that the main emission component in this case is slightly red-shifted. The spatially distinct clumps at different projected velocities might hence be interpreted as tracers of different localized mass ejection events, and it is even plausible that “NB_blue” marks an event different from “NB_red”. Some of these events are heading away from the observer, other ones are approaching us. This type of 3-dimensional geometry has already been suggested in the literature (see Sect. 1). Combining visual and near-IR images and 4.6 μm CO spectra Smith et al. (2009) proposed a geometry where a halo is responsible for an ubiquitous 25 km s^-1 outflow with potentially mild axisymmetry. Faster (~40 km s^-1) dense CO cloudlets move through this halo and create a very asymmetric envelope. The NaCl clumps clearly trace much lower velocities, albeit their emission is detected far beyond the dust formation radius estimated to be around 10 R_⋆ (or 55 mas Monnier et al. 1999). We can not exclude the coincidental fact that the different NaCl clumps have space vectors quite close to the plane of the sky. However, it might be that these clumps indeed trace lower velocity structures. We note that Kamiński et al. (2013a) also deduced a velocity around 15 km s^-1 for the “SW” clump, but other molecules in the SMA survey which emit closer to the central target trace velocities larger than the ALMA NaCl “NB” and “SB” clumps. In the case that the NaCl clumps and CO cloudlets are part of the same physical structure, the interaction interface between the cloudlet and the surrounding material could enhance the gas temperature thus facilitating the excitation of both CO and NaCl.

In the case of localized mass ejections, each event will have more or less the same chemical content at the moment of its ejection and the same laws of thermodynamics and chemistry will apply to each of them. One then might wonder why the morphology as seen in NaCl is so different compared to what is traced by the ALMA TiO₂ channel maps presented by De Beck et al. (2015), a molecule that also has a high dipole moment (6.33 Debye) and with an expected abundance slightly higher than NaCl in the inner wind region (Gobrecht et al. 2016). We simply note here that the Grotrian diagram of NaCl is less complex than for TiO₂ and that even different TiO₂ lines with similar excitation energies show quite different channel maps. This leads us to propose that the different mass ejections have very short dynamical times and different molecules/excitations might (partly) trace different events where the local thermodynamical, radiative and chemical properties are different. The nature of the star led Smith et al. (2001) to suggest that these ejection events could be of convective or magnetic origin, alike the Sun with its solar flares and eruptions.

There is, however, some similarity between the ALMA NaCl and TiO₂ emission, i.e. both molecules trace the “SW” clump out to ~1.2′′, a clump also detected in SiS, CS, NS, and H₂S in the SMA data of Kamiński et al. (2013a) and in dust scattered light out to 1.4′′ (Smith et al. 2001; Shenoy et al. 2013). Some of these molecules are expected to be quite abundant (as SiS) and others not, and they have very different dipole moments (from 1 to 9 Debye). This might indicate that they trace the same physical event along a northeast-southwest axis. The “SW” clump is potentially embedded within the lower density southwest bipolar lobe (see Sect. 4.2.1) or it might be that the dust clump is a high mass-loss event leaving a lower density wake behind it. In both cases, the lower density facilitates the radiative excitation of the molecules by direct stellar light. And while the nature of the bright “C” component in the ALMA dust continuum map is not yet clear (O’Gorman et al. 2015), in retrospective it might be that this dust clump also swept up material while moving away from the star leaving a lower density region behind it that facilitates the excitation of NaCl seen in the “SB”, this under the assumption that both are in the same direction w.r.t. the central star.

5. Conclusions

In this paper we have analysed gaseous NaCl in the inner wind of the luminous red supergiant VY CMa. The high sensitivity and high spatial resolution capabilities of ALMA allow a detection of this molecule within 250 R_⋆ from the central star. Four NaCl lines are detected, being the v = 0, 1 J = 24−23 and J = 25−24 rotational lines in ALMA band 7 at a spatial resolution of . The ALMA data permit the first spatial decomposition of the NaCl emission in up to four different spatial regions. We show that all NaCl lines show a bimodal spatial distribution, with a bright emission component to the north of the central star and a slightly weaker emission component to the south at an average separation of 100 mas (or ~18 R_⋆) and 190 mas (or ~34 R_⋆), respectively. The rotational NaCl lines in the ground vibrational state show in addition an extra emission component at ~1.22′′ to the southwest of the star. This southwest clump is already detected in other molecules (CS, SiS, NS, and H₂S) with the SMA (Kamiński et al. 2013a) and in TiO₂ with ALMA (De Beck et al. 2015) and is aligned with the southwest extension seen in dust scattered light (Smith et al. 2001) and other near-infrared and polarimetric data (e.g. Monnier et al. 1999; Jones et al. 2007).

The NaCl lines are detected up to a velocity of 17 km s^-1 w.r.t. the local standard of rest velocity, v_LSR, of 22 km s^-1. This is significantly lower than the terminal velocity, which is around 35 km s^-1 as estimated from low-J rotational CO lines. This low velocity was already noticed by Milam et al. (2007), who explained this as being due to the very efficient condensation of NaCl, being a highly refractory species. However, the fact that we detect gaseous NaCl at distance of 1.2′′, far beyond the main dust formation region estimated around 35−55 mas (~10 R_⋆), points toward a chemical process preventing all gaseous NaCl from nucleating. We construct a simplified chemical network aiming at understanding the cycling between gaseous and solid NaCl through the envelope of VY CMa. We show that at the dust condensation locus where the grain temperatures can be above 300 K, thermal desorption dominates the rate of accretion of gaseous NaCl onto existing dust grains. For grain temperatures lower than 100 K, photodesorption is unlikely to be important and removal of NaCl from the dust is likely due to shock sputtering driven by localized mass ejection events.

It is not straightforward to interpret the geometrical signatures in the inner wind. Based on the ALMA data, and previous optical, near-IR, and mid-IR data two scenarios are favoured. In the first scenario, the overall density structure in the inner wind is described by an equatorial density-enhanced region with the lower density bipolar regions oriented northeast-southwest at a position angle of ~50°. The density-enhanced region, be it a disk, has a projected velocity of 2.5 km s^-1. The far side of the bipolar lobes is situated to the north, in agreement with the non-detection of the optical reflection nebula to the northeast of the star by Smith et al. (2001). However, this picture alone can not explain all morphological signatures. The south-west clump at 1.2′′ and the presence of arcs and knots in optical images point toward the action of localised mass ejections creating higher density structures within the diffuse halo. The present data, obtained with only 16–20 ALMA antennas, can not prove the validity of both scenarios, albeit it seems likely that both play a role in shaping the wind of VY CMa.

Acknowledgments

This paper makes use of the following ALMA data: ADS/JAO.ALMA2011.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.

References

Appendix A: Channel maps of the NaCl v = 0, 1 J = 25–24 lines

Figures A.1 and A.2 show the channel maps for the NaCl v = 0,1J = 25−24 emission lines.

Fig. A.1 Same as Fig. 2 but for the NaCl v = 1J = 25−24 line at 322 GHz.

Fig. A.2 Same as Fig. 2 but for the NaCl v = 0J = 25−24 line at 325 GHz. Since this line is severely blended with the 325 GHz H2O maser line and the SiS v = 1J = 18−17 line the integrated line strength is not shown.

Appendix B: Additional equations

The dust-to-gas mass ratio, ψ, is defined as $\mathit{\psi }\mathrm{=}\frac{{\mathit{n}}_{\mathrm{d}}\hspace{0.17em}{\mathit{m}}_{\mathrm{d}}}{{\mathit{\rho }}_{\mathrm{g}}}\mathit{,}$(B.1)with n_d(a,r) being the grain volume density (see Sect. 4.1), m_d the mass of the grain being 4 / 3 πa³ρ_s, with ρ_s the specific

density of the dust grain taken to 3.3 g cm^-3 for silicate grains (Draine & Lee 1984), and ρ_g the gas mass density being n_H(r) m_H, with m_H the mass of a hydrogen atom. Applying the same grain size distribution as described in Sect. 4.1, it is easily shown that $\mathit{\psi }\mathrm{=}\mathrm{7}\mathrm{\times }{\mathrm{10}}^{\mathrm{22}}\hspace{0.17em}\mathit{K}\mathit{.}$(B.2)For a dust-to-gas mass ratio of 0.002, this implies that K is equal to 2.8 × 10^-26 cm^2.5 per H atom.

The grain density n_d(r) = ^∫n(a,r) da is then given by 2.2 × 10¹⁵Kn_H(r), which for the given value of K results in n_d(r) = 6 × 10^-11n_H(r).

All Tables

All Figures

Fig. 1 Flux densities for the four NaCl lines detected in the ALMA band 7 data. Flux densities are for an aperture with circular diameter of 1′′ centred at “VY”. The LSR velocity of 22 km s-1 is indicated by the vertical dotted line. Line blending causes the flux density for the 325 GHz NaCl v = 0J = 25−24 line to be skewed and to be too high with an expected contribution of the SiS v = 1J = 18−17 line of ~1.1 Jy (see text for more details).

In the text

Fig. 2 Channel maps (colour scale) of the NaCl v = 0J = 24−23 emission at 312 GHz at a 1 km s-1 velocity resolution for the central 0.̋8 × 0.̋8 field around “VY” and “C”. Grey contours are the integrated NaCl line strength at (1,2,3,4) × 1 Jy/beam km s-1, red contours are the dust continuum measured with ALMA at 321 GHz at (1,2,3,4) × 18 mJy/beam (Richards et al. 2014; O’Gorman et al. 2015). The stellar position “VY” is indicated with a white cross, and the peak of the continuum emission “C” is indicated with a white plus-sign. The ordinate and co-ordinate axis give the offset of the right ascension and declination, respectively, in units of milli-arcseconds.

In the text

	Fig. 3 Same as Fig. 2 but for the NaCl v = 1J = 24−23 line at 309 GHz.
In the text

	Fig. 4 Same as Fig. 2, but showing the central 2.̋5 × 2.̋5 region. The white contours show the HST/WFPC2 data in the F1042M filter around 1.018 μm at (1,2,4,8,16) × 110 counts (Smith et al. 2001).
In the text

Fig. 5 ALMA NaCl 312 GHz channel map at 25 km s-1 (colour scale) overlaid with the HST/WFPC2 F1042M data of Smith et al. (2001) (in contours at (1,2,4,8,16) × 110 counts). The cross sign indicates the position of the central star VY CMa, noted “VY” in the text, as derived from the ALMA H2O maser and continuum data (Richards et al. 2014; O’Gorman et al. 2015), the plus-sign indicates the maximum of the ALMA dust continuum data, noted as “C” in the text, (O’Gorman et al. 2015). The nomenclature of the different spatial regions refers to the northern clump ’NB’, the southern clump ’SB’, the brightest part of the south-western extension ’SW1’, and the faintest part of the north-western extension“SW2”. The ordinate and co-ordinate axis give the offset of the right ascension and declination, respectively, in units of arcseconds.

In the text

Fig. 6 Flux densities for the NaCl J = 24−23 lines in the ground-vibrational state v = 0 at 312 GHz (upper panel) and in the first vibrational state v = 1 at 309 GHz (lower panel). Flux densities are extracted in the north clump (NB, full line) at (07h22m58.326s, −25°46′02.̋9) for an aperture of (0.̋394 × 0.̋373), in the south clump (SB, dotted line) at (07h22m58.332s, −25°46′03.̋118) for an aperture of (0.̋36 × 0.̋36), and in case of the 312 GHz line in the south-west clump (SW, dashed line) at (07h22m58.273s, −25°46′04.̋028) for an aperture of (0.̋57 × 0.̋73). The LSR velocity of 22 km s-1 is indicated by the vertical dotted line.

In the text

	Fig. 7 The H2O maser emission (Richards et al. 1998, 2014) is shown on top of the NaCl 312 GHz integrated line emission (gray contours, zeroth moment at 0.25, 0.5, 1, 2, 3, 4 Jy km s-1) and the ALMA dust continuum emission (colour scale; O’Gorman et al. 2015).
In the text

Fig. 8 Steady state solutions for the gaseous and solid NaCl abundances (relative to the hydrogen number density) as a function of the dust temperature, Td, for the case in which accretion, thermal desorption and photodesorption are included in the chemical network. The different panels are for different values of the binding temperature (5000, 10 000, 15 000 K) at a hydrogen number density of 1 × 108 cm-3. Each panel also shows the inverse timescale for accretion, thermal desorption and photodesorption in units of s-1, calculated with FFUVYPD = 3 × 105.

In the text

	Fig. 9 Fractional abundances of gaseous and solid NaCl in function of time for a binding temperature of 15 000 K at three different values of the dust temperature Td (290, 300, 310 K).
In the text

Fig. 10 Position-velocity diagram of the NaCl 312 GHz for an axis length of 1′′ and an averaging width of 650 mas centred at α = 07h22m58.332s and δ = −25°46′03.̋053. In the left and middle panels, the slit has a PA of 50° (i.e. along the presumed polar/rotation axis; with north-east positions having negative offsets and south-west positions positive offsets), in the right panel the slit PA is 140° (i.e. along the presumed equatorial plane; with south-east positions having negative offsets and north-west positions positive offsets). In the left and right panels the different velocity components are indicated, while in the middle and right panels the main contributors to the different intensity peaks are indicated. The contours in the left panel are at [1,2,...10] × 0.01 Jy/beam.

In the text

	Fig. A.1 Same as Fig. 2 but for the NaCl v = 1J = 25−24 line at 322 GHz.
In the text

	Fig. A.2 Same as Fig. 2 but for the NaCl v = 0J = 25−24 line at 325 GHz. Since this line is severely blended with the 325 GHz H2O maser line and the SiS v = 1J = 18−17 line the integrated line strength is not shown.
In the text