Spectro-Imaging of the Asymmetric Inner Molecular and Dust Shell Region of the Mira Variable W Hya with MIDI/VLTI¹

The spectacular end of an intermediate-mass star (1 M_⊙ ≤ M_star ≤ 8 M_⊙) is its transformation from an AGB star into a planetary nebula developing a rich variety of shapes (e.g., Balick & Frank 2002; De Marco 2009). Recently, it has been shown that deviations from spherical symmetry in the circumstellar environment are already present in a large fraction of AGB stars (e.g., Castro-Carrizo et al. 2010; Cox et al. 2012; Johnson & Jones 1991; Weinberger & Aryal 2003, and references therein).

However, in most cases, it is not obvious if these asymmetries are intrinsic to the star, due to interactions with the surrounding environment or a companion. Various scenarios have been suggested to explain the origin of the observed asymmetries, but it is commonly agreed that understanding the mass-loss geometry is crucial for explaining the dust formation and wind acceleration mechanism itself (e.g., Dorfi & Hoefner 1996; Woitke & Niccolini 2005). In particular, in M-type (oxygen-rich) AGB stars, the pulsation-enhanced dust-driven wind scenario is still under intensive debate (e.g., Woitke 2006; Höfner 2008; Norris et al. 2012; Bladh et al. 2013) since this wind-driving mechanism has not been satisfactorily verified so far. The low-velocity winds can be easily affected by small perturbations influencing the dust production efficiency.

Stellar rotation and convection, cool stellar spots, or a weak magnetic dipole field likely trigger the dust formation in certain zones by increasing the density or reducing the temperature (e.g., Soker 2000; Matt et al. 2000; Szymczak et al. 1998). Dust also efficiently forms in a shocked environment and is believed to appear in clumps (e.g., Woitke & Niccolini 2005). A nonsphericity might also be caused by a companion (e.g., Mauron & Huggins 2006; Mayer et al. 2011; Kim & Taam 2012) or by interactions with the ambient interstellar medium (ISM; e.g., Ueta et al. 2010; Jorissen et al. 2011) on larger scales.

One way to study the asymmetric dust distribution and mass-loss is by mapping the inner region of the envelope where the mass loss originates. A promising method in optical/infrared interferometry is to use image reconstruction techniques instead of simple geometric modeling (Berger et al. 2012) to model more complicated source morphologies as, e.g., shown in the near-IR by Le Bouquin et al. (2009) and Ohnaka et al. (2013) for T Leporis and Antares, respectively. Both authors used AMBER data obtained within a few days to get a snapshot of these evolved stars in order to avoid the change of the appearance on time scales on the order of a few months. Compared to Le Bouquin et al. (2009) and Ohnaka et al. (2013), the MIDI data in this work are spread over several pulsation cycles. It is therefore expected that the reconstructed intensity maps represent the averaged appearance in the mid-IR spectral domain. This was not avoidable since a certain UV-coverage was needed. Images of W Hya were reconstructed with the MiRA (Multiaperture Image Reconstruction Algorithm) package (Thiébaut 2008; Cotton et al. 2008; Renard et al. 2011; Thiébaut 2013) in several spectral bins in the dust sensitive mid-IR spectral domain.

W Hya (13h 49 min 01.998 s, -28°22^'03.49'', M7 III e, μ_α ∗= -49.31 mas yr^-1, μ_δ = -59.71 mas yr^-1, v_rad = 42.3 km s^-1; van Leeuwen 2007; Wilson 1953) is one of the closest and best studied M-type AGB stars showing a departure from spherical symmetry (see below). The distance and mass-loss are estimated to be (Vlemmings et al. 2003) and (3.5 - 8) × 10^-8 M_⊙ yr^-1 (Justtanont et al. 2004), respectively. W Hya is listed in the Washington Visual Double Star Catalog (WDS; Mason et al. 2001). However, it is most likely only an optical pair rather then physically related as indicated by a different proper motion. The separation at epoch 1999 is given to be 68.6'' at a position angle (PA, measured from north to east) of 108°. No evidence for the presence of any additional close companions could be found in the Hipparcos astrometry data by Pourbaix et al. (2003) within the Hipparcos uncertainties.

Attempts to model asymmetries in the very close vicinity of W Hya were made with part of the data presented here before (Zhao-Geisler et al. 2011). With a simple geometric fully limb darkened disk, an elongation along a PA of (11.2 ± 6.2)° and an axis ratio of 0.87 ± 0.07 were found. The outer dust envelope, observed by Marengo et al. (2000) at 18 μm, shows a very similar N-S elongation on a 100 times larger scale of 8''. Tuthill et al. (1999) found a discrepant PA of (94 ± 33)° at visual wavelengths.

The SiO, H₂O, OH maser observations (Reid & Menten 1990, 2007; Imai et al. 2010; Szymczak et al. 1998), the radio photosphere at 43 GHz (Reid & Menten 2007), and the HCN thermal emission (Muller et al. 2008), tracing the gas at larger radii, are approximately oriented perpendicular to a N–S direction except of the 215.2 GHz SO line (Vlemmings et al. 2011) with a PA of the emission structure of (3 ± 10)°. In these images, large deviations from circularity are found. A velocity gradient is observed along the N–S direction (Szymczak et al. 1998; Vlemmings et al. 2011; Muller et al. 2008).

IRAS and Herschel-PACS images of the extended envelope of W Hya also show a complex spatial structure (Hawkins [1990] and Cox et al. [2012], respectively). Currently, it is not possible to establish a strict connection between the asymmetries seen in the small and large-scale structures, and between the distribution of dust and molecular gas.

2.1. Interferometric Observations with MIDI/VLTI

The data presented here were obtained with the mid-IR (8–13 μm) interferometric instrument MIDI (Leinert et al. 2003) at the Very Large Telescope Interferometer (VLTI) in service mode using the Auxiliary Telescopes (ATs). Some of the data were already presented by Zhao-Geisler et al. (2011). For this study, 26 additional measurements were included in order to conduct an image reconstruction. They were obtained between 2010 January and 2011 March under the program IDs 083.D-0294, 084.D-0334, and 085.D-0201 in GTO time. An overview of the new observations is given in Table 1.

Before or after each target observation a calibrator star was observed with the same set-up in order to calibrate the visibilities and fluxes. The properties of the calibrator stars are listed by Zhao-Geisler et al. (2012) in Table 3. 2 Cen (HD 120323) was predominantly used to calibrate the W Hya data. Observations were executed in SCI-PHOT mode, where the photometric and the interferometric spectra are recorded simultaneously. This has the advantage that the photometric spectrum and the fringe signal are observed under the same atmospheric conditions. The prism, with a spectral resolution of R = λ/Δλ = 30, was used to obtain spectrally dispersed fringes.

The UV-coverage of the newly obtained data is shown in Figure 1 together with the previous observations. It can be clearly seen that the UV-plane is not completely sampled. However, compared to most interferometric observations in the optical to mid-IR (e.g., Paladini et al. 2012; Klotz et al. 2012; Sacuto & Chesneau 2009; Ohnaka et al. 2008), the coverage is excellent for an image reconstruction effort. Since particularly the NW to SE direction is not well filled, the reconstructed images will be less well constrained along this axis.

Zoom In Zoom Out Reset image size

2.2. MIDI SCI-PHOT Data Reduction

The image reconstruction was performed with the freely available MiRA package⁹ (Thiébaut 2008; Cotton et al. 2008; Renard et al. 2011; Thiébaut 2013) (version 0.9.10). This algorithm was specifically developed to overcome the problems encountered with optical interferometric data like the sparseness of the UV-plane and missing Fourier phase information. MiRA is purely based on an inverse problem approach and can build an image without any Fourier phase information.

The basic way to account for the missing information is by applying a priori constraints in addition to the data fitting. The total criterion which needs to be minimized is then the sum of two terms:

with f_data(x) the likelihood term which measures the compatibility of the parameters with the available data, and f_prior(x) the regularization term which measures the compatibility with the prior information. One regularization term later used is the quadratic smoothness regularization which is given by:

where S is the smoothing operator implemented in MiRA via finite differences. The so-called hyperparameter μ is used to adjust the relative weight of the constraints set by the measurements and the ones set by the priors (Thiébaut 2008; Renard et al. 2011). Several different regularization terms are implemented in MiRA. Initial values for μ were obtained from the regularization benchmarking test by Renard et al. (2011).

The MIDI data were written into OI-fits files which are used as input for MiRA. MIDI does not provide absolute phase information. To create a complex visibility, we used the modulus of the visibility computed by MIA/EWS and a phase of zero. Setting the phase to 0° forces a point symmetric appearance of the object. For the modulus of the visibilities, the errors obtained by the method described in § 2.2 were applied. The error for the phase was arbitrarily set to 10°.

The image of each spectral bin was reconstructed independently. The "best" images (cf. discussion in the next section) could be obtained with the quadratic smoothness regularization (called "smoothness" in MiRA) which most closely describes a smooth image and tries to avoid unmeasured high spatial frequencies. A detailed example how this is done in MiRA is given in the appendix.

The quadratic smoothness regularization was compared against other regularization terms. The total variation regularization "totvar," which minimizes the total gradient in the image, and the compactness regularization "quadratic" with weights of power 2 and 3, which describes compactness in the image plane and hence smoothness in the Fourier plane, gave apparently less plausible results. In addition, images were reconstructed by using only parts of the data to study the influence of the UV-coverage.

The images which were reconstructed with the quadratic smoothness regularization term for all 25 wavelength bins are shown in Figure 2. The intensity is normalized to the highest value in each image. Since the image reconstruction does in general not conserve the flux, a flux comparison is not directly possible. The results for the additional three image reconstruction attempts (total variation and compactness regularization with weights of power 2 and 3) are shown in Figure 3 for one representative wavelength. Compared to the quadratic smoothness regularization, the images are less smooth, sometimes saturated, and show spikes and structures like a double peak. This is in particular apparent in the cuts along different position angles (PA) and believed not to be a source property.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

From the images in Figure 2, it can be clearly seen that W Hya is elongated and appears larger at longer wavelengths. Cuts along the major and minor axes as shown in Figure 4 confirm this. The PA of the major axis is (15 ± 10)°, taking for the uncertainty also the UV-coverage into consideration. Along that PA, the visibility sampling is sparse and only certain components can be traced. The traced structures along the major axis appear larger at longer wavelengths with a fitted Gaussian FWHM of 66 ± 5 mas (corresponding to 6.7 AU) between 8 and 10 μm, and 106 ± 3 mas (10.8 AU) beyond 11 μm. The errors are derived by taking the standard deviation of the diameters obtained with the mentioned four image reconstruction attempts. At the minor axis, with a PA of (105 ± 10)°, the overall fitted Gaussian FWHM diameter does not change as function of wavelength and has a value of about 42 ± 2 mas (4.1 AU). With the given UV-coverage along this axis, this could be well modeled.

Zoom In Zoom Out Reset image size

The derived Gaussian FWHM diameters are plotted in Figure 5 for the minor and major axis as function of wavelength as well as for several PAs around the minor axis. Altogether, this results in an axis ratio of 0.62 ± 0.06 between 8 and 10 μm, and 0.40 ± 0.02 beyond 11 μm. The axis ratio is therefore significantly smaller than previously estimated (Note, however, that the size of the major axis is not constrained well, as already mentioned.) At shorter wavelengths, from wavelength bin 8.12 μm to bin 9.26 μm (cf. Fig. 4), it seems that there is some extended emission along the minor axis. However, this could be traced back to remaining sidelobes from the image reconstruction process.

Zoom In Zoom Out Reset image size

The effect of using only certain parts of the data with only a limited UV-coverage as indicated in Figure 1 can be seen in Figure 6. By using only observations with a PA between 20° and 110° (modulo 180°) (case A), the appearance does not change considerably compared to the full data set image. If only observations at short baselines with spatial frequencies (SPF) less than 20 arcsec^-1 and long baselines with a PA between 110° and 200° (modulo 180°) are used (case B), the structures along the minor axis disappear. In the case that all long baselines (SPF≥20 arcsec^-1) are excluded (case C), all the fine structure vanishes as expected. In all three cases the PA of the major axis could be determined to be about (15 ± 10)° as well. These three tests show in particular that the long baseline observations with a PA between 110° and 200° (modulo 180°) do not influence the orientation of the elongation significantly.

Zoom In Zoom Out Reset image size

In all three cases, the major axis becomes shorter compared to the full data set image as can be derived from the intensity profiles, which confirms that the major axis is not well constrained. The image reconstruction largely interpolates the values here. The long baselines trace small scale structures as expected. However, their influence can be considered as not significant since the small scale structures are smoothed out by mixing data over several pulsation phases. A fourth case was tested using only long baselines but, as anticipated, no image could be reconstructed due to insufficient Fourier plane information. For the same reason, no images could be obtained if data only for discrete pulsation phase ranges were used.

5.1. The Image Reconstruction Process

5.2. The Diameter as Function Of Wavelength

5.3. The Asymmetries at Different Scales

The PA of the major axis could be relatively well determined to (15 ± 10)°. This is very similar to the value obtained with the previous purely geometrical visibility modeling by Zhao-Geisler et al. (2011). In that paper, it was already discussed that the elongation itself and its orientation are not an effect of the specific UV-coverage by investigating the behavior of the visibility as function of PA. This was reinvestigated and confirmed in this work by excluding parts of the data with certain PAs and spatial frequencies. An effect due to mixing data taken at different pulsation phases could been previously excluded as well. Due to this mixing, smaller scales and smaller asymmetries can not be probed and are smeared out.

There are basically five major shaping mechanisms which lead to asymmetries. These are wind-wind interactions, ISM interactions, interactions with a galactic magnetic field, binary interactions, and asymmetric mass loss. They are discussed in the following paragraph.

Figure 7 shows a comparison of the reconstructed MIDI image in the mid-IR with the Herschel-PACS image of W Hya at 70 μm (Cox et al. 2012). The Herschel image reveals a detached dust shell and a jet-like structure in the north–east direction at an approximate PA of 40°. This is the most disrupted detached shell seen in the Herschel sample. It is also the only one showing a bipolar-like outflow and a complex structure within the shell. It was discussed that the detached shell around W Hya is a typical phenomenon for a wind–wind and a wind–ISM interaction.

Zoom In Zoom Out Reset image size

A comparison with the tail of Mira AB (Wareing et al. 2007; Martin et al. 2007) seems at first glance to provide an interesting explanation for the NE jet-like structure. The PA of the space velocity is 211°, thus opposite to the direction of the prominent NE outflow as indicated in Figure 7. However, there are several problems with this analogy. GALEX ultraviolet observations of W Hya in the far and near UV (Martin et al. 2005), tracing possible shock regions with the ISM, show no sign of a tail-like structure and only some irregular emission at larger scales as shown in Figure 8. It would also not be plausible to have a symmetric outer shell and a tail due to interactions with the ISM within.

Zoom In Zoom Out Reset image size

Another way of forming an asymmetric outflow is by a galactic magnetic field which would shape an outflow in an eye-like structure (van Marle et al. 2014). This was found several times in the Herschel sample (Cox et al. 2012). Such an appearance is caused by the fact that the wind velocities perpendicular to the field lines are lower. This is however not seen here.

More interesting is that the MIDI images reveal an axisymmetric structure at a similar PA as the Herschel-PACS image on a scale which is three orders of magnitudes smaller. Even though the PAs are somewhat different (but not by a large amount), this might be a hint for a connection between small and large scales. This therefore might show that the large-scale asymmetric mass-loss originates in the very close vicinity of the star. Both structures, however, probe different epochs and time scales: a few years of wind expansion on the smallest scales (50 mas; 4.9 AU) versus up to 4600 years at the largest scales (80''; 7800 AU), assuming a terminal velocity of 8.0 km s^-1 (Szymczak et al. 1998).

The small PA discrepancies between the small scale asymmetries and those on large-scale might be due to a precession effect as, e.g., seen in planetary nebulae (e.g., Sahai et al. 2005; Volk et al. 2007). For this kind of precession of a collimated or nonspherical outflow an accretion disk around a companion is needed. According to Soker & Harpaz (1992), Soker (2006), and Nordhaus & Blackman (2006), a single star cannot supply the energy and angular momentum to shape the circumstellar wind since magnetic fields are not efficient enough. However, as mentioned in the introduction, there is no sign for a close companion around W Hya (Pourbaix et al. 2003). On the other hand, it would also be very difficult to observe a secondary as faint as a brown dwarf or even a Jupiter-like object which would be suitable to trigger a collimated outflow (cf., e.g., Nordhaus & Blackman 2006; Hrivnak et al. 2011).

In the end, no clear shaping mechanism can be identified which would also explain the prominent NE dust clump seen in the Herschel-PACS image. It can only be said that there might be a possible connection between small and large scales. In order to more clearly disentangle effects which are due to the star itself and due to its close environment, the inconclusive asymmetries seen in the molecular line/maser emissions of e.g. SiO and H₂O very close to the star (as discussed by Zhao-Geisler et al. [2011]) as well as the onset of a possible asymmetric mass-loss, a kinematic study of the innermost part of the circumstellar environment is required. In addition, a search for a companion and detailed modeling of the outflow structure will help to get a better insight as well.

Asymmetries in the very close vicinity of W Hya were investigated in the dust sensitive mid-IR domain between 8 and 12 μm. To our knowledge, this is the first time that images could be obtained with the image reconstruction package MiRA using data from the MIDI instrument. The apparently best images were reconstructed with the smoothness regularization term.

W Hya appears elongated with a position angle of the major axis of about (15 ± 10)° and a less well-defined axis ratio between 0.4 and 0.6. In contrast to the major axis, the minor axis is constant over the wavelength range with a Gaussian FWHM diameter of about 42 ± 2 mas (4.1 AU) which is approximately equivalent to the photospheric diameter of the star. The spectral dependence of the major axis diameter points again to the existence of presumably amorphous aluminum oxide (Al₂O₃) as one of the first dust condensates at around 2 photospheric radii, which is concentrated along the major axis. Smaller scales could not been probed due to the combination of observations made at different pulsation phases, which was necessary to preserve a good UV-coverage.

A comparison of the elongated structure seen with MIDI with the Herschel-PACS 70 μm image at much larger scales indicates that asymmetric mass-loss likely originates in the very close vicinity of the star. Further investigations are necessary to better understand the effects which are due to the star itself or in its close vicinity, including the possible existence of a companion.

The final images were produced in two steps. In the first run, with 500 iterations, the hyperparameter μ (cf. § 3) was set to 10⁹ while for the second run, also with 500 iterations, μ was set to 10⁶ for the quadratic smoothness regularization "smoothness." A large μ in the first round gives more weight to the regularization itself in order to get a good initial image. A lower μ in the second run adjusts the image to represent the observational data more accurately.

The basic steps done in MiRA are as following:

> mira_config, db, dim=256, xform="exact", pixelsize=5.0*MIRA_MILLIARCSECOND

> rgl = rgl_new("smoothness")

> img1 = mira_solve(db, maxeval=500, verb=1, xmin=0.0, normalization=1, regul=rgl, mu=1e9)

> img2 = mira_solve(db, img1, maxeval=500, verb=1, xmin=0.0, normalization=1, regul=rgl, mu=1e6)

with db being the input data set.

For the total variation regularization "totvar", μ was set to 10³ and 10² in the first and second run, respectively, and for the compactness regularization "quadratic" with weights of power 2 and 3, μ was set to 10⁷ and 10⁵ in the first and second run, respectively.¹⁰