Evidence for localized onset of episodic mass loss in Mira

The interferometric data presented herein were obtained using the IOTA interferometer (Traub et al. 2003), which was a long baseline interferometer operating at near-infrared wavelengths until 2006. It consisted of three 0.45 m telescopes movable among 17 stations along two orthogonal linear arms. IOTA could synthesize a total aperture size of 35$\times$15 m, corresponding to an angular resolution of 10$\times$23 milliarcseconds at $1.65\,\mu m$. Visibility and closure phase measurements were obtained using the integrated optics combiner IONIC (Berger et al. 2003); light from the three telescopes is focused into single-mode fibers and injected into the planar integrated optics (IO) device. Six IO couplers allow recombinations between each pair of telescopes. Fringe detection was done using a Rockwell PICNIC detector (Pedretti et al. 2004). The interference fringes were temporally-modulated on the detector by scanning piezo mirrors placed in two of the three beams of the interferometer.

Observations were carried out in the H band ($1.5\,\mu$m $\leq\lambda\leq 1.8\,\mu$m) divided into seven spectral channels. The science target observations were interleaved with identical observations of unresolved or partially resolved stars, used to calibrate the interferometer’s instrumental response and effects of atmospheric seeing on the visibility amplitudes. The calibrator sources were chosen from the Borde et al. (2002) catalog. Mira was observed in October 2005 over nine nights and using five different configurations of the interferometer. Full observation information can be found in Table 1, including dates of observation, interferometer configurations, and visual phase computed from the Association Française des Observateurs d’Étoiles Variables (http://cdsarc.u-strasbg.fr/afoev/). The calibrators are listed in Table 2. Figure 1 shows the (u,v) coverage achieved during this observation run. The geometry of the IOTA interferometer and the position of the star on the sky constrained the frequency coverage, which is elongated and primarily located in the north-south direction.

Reduction of the IONIC visibility data was carried out using custom software by Monnier et al. (2006), which is similar in its main principles to those described by Coudé du Foresto et al. (1997). The premise is to measure the power spectrum of each interferogram (proportional to the target squared visibility, $V^{2}$) after correcting for intensity fluctuations and subtracting bias terms due to read noise, residual intensity fluctuations, and photon noise (Perrin, 2003). Next a correction is applied by the data pipeline for the variable flux ratios for each baseline by using a flux transfer matrix (Monnier, 2001). Finally, raw squared visibilities are calibrated using the visibilities obtained by the same means on the calibrator stars. As in Lacour et al. (2008) less than a 1% uncertainty in visibility is associated with the calibrators. This corresponds to a 2% calibration error for $V^{2}$. Therefore, we systematically added a 2% calibration error to all the calibrated squared visibilities in this paper. In order to measure the closure phase (CP), a real-time fringe tracking algorithm (Pedretti et al., 2005) ensured that interference occurred simultaneously for all baselines. We required that interferograms are detected for at least two of the three baselines in order to assure a good closure phase measurement. This technique, called baseline bootstrapping, allowed for precise visibility and closure phase measurements for a third baseline with very small coherence fringe amplitude. We followed the method of Baldwin et al. (1996) to calculate the complex triple amplitudes and derive the closure phases. Pair-wise combiners (such as IONIC) can have a large instrumental offset for the closure phase which can be calibrated by the closure phase of the calibrator stars.

The squared visibilities and the closure phases are presented in Fig. 2 and Fig. 3, respectively. The colors show different baseline orientations. Obviously the source is not circularly symmetric as the visibilities depend on azimuth. A zero crossing is visible in the visibility function near 50 cycles/arcsec, the signature of sharp edges for the main contributor to the flux. In addition, the closure phases depart from 0 modulo $\pi$ beyond the first zero of the visibility function, showing a departure from centro-symmetry at scales that are smaller or comparable to that of the central object.

An image was reconstructed using the monochromatic version of the MiRA software of Thiébaut (2008) (https://github.com/emmt/MiRA). As the (u,v) coverage was sparse, we used the same method as in Lacour et al. (2009) to regularize the image reconstruction process. The criterion minimized by MiRA in this section is a weighted sum of the squared distance to the data and of the squared distance to an initial image (SD2I in Appendix A). The initial image was built by fitting the data with a circularly symmetric star + shell model of Perrin et al. (2004), the parameters of which are given in Table 3. The total number of available data is 940 (705 squared visibilities and 235 closure phases). Images of 100x100 pixel were reconstructed. A hyper parameter $\mu$ controls the weight of the regularization in the reconstructed image with a total $\chi^{2}$ reading: $\chi^{2}_{\mathrm{tot}}=\chi^{2}_{\mathrm{data}}+\mu\chi^{2}_{\mathrm{rgl}}$. $\chi^{2}_{\mathrm{data}}$ is defined and discussed in Appendix A. The higher $\mu$ the higher the weight of the regularization and its impact on the reconstructed image. We tried several values of the hyper parameter and tested the impact on the image reconstructed with MiRA.
We tried other reconstructions with MiRA or MACIM using various types of initial images and priors imposed by the regularization. They are discussed in Appendix A and the images are presented in Fig. 8. For the level of information available at the angular resolution obtained with IOTA, all images share the same basic features: 1) a bright central object of size $\simeq 20\,mas$ elongated in the north-south direction; 2) surrounded by a fainter environment of $\simeq 40\,mas$, roughly circular; 3) a dark patch in the environment that is west of the bright central object. From among the available images, we selected the one with the lowest hyper parameter, $\mu=10^{8}$, as it is more constrained by the data relative to the others. The image is shown in Fig. 4.

The image is generally consistent with the core-halo description of LPVs developed in Perrin et al. (2004). The two-component, circularly symmetric model correlates well with infrared spectroscopy, recorded in the same spectral region (Hinkle et al., 1984). The spectroscopy records molecular absorption against continuum emission from the effective photosphere. It is natural to associate the molecular absorption zone with the halo in Fig. 4 and to understand the core in the image as a representation of the deeper continuum-forming layer. This association is supported by the dynamical content of the spectra. The motions of the pulsating LPV atmosphere are followed in the Hinkle et al. (1984) study, showing in motions of gaseous CO the expansion, slowing, and falling back of the partially transparent halo. The radial velocity of this layer with respect to the Mira CM velocity (as inferred from the centroid of low excitation envelope emission lines) can be used to follow the halo expansion through a pulsational cycle. Using an average envelope speed of 7.5 km.s^-1 Hinkle et al. (1984) for half of the 321-day period, the total extent of the pulsational motion is found to be $1.1\times 10^{8}$ km. This corresponds well with the $1.4\times 10^{8}$ km radius difference between the top of the 10 mas model halo and the photosphere (based on the Hipparcos distance of $2.8\times 10^{15}$ km). Thus, to good approximation, the inner boundary of the imaged and modeled halo shows the H-band photosphere and the outer boundary of the halo maps, approximately, the extent of the periodically pulsating gaseous layers. The minimum observed CO excitation temperatures, $\approx$2700K, are in between the inferred model layer temperatures of 2000K and 3400K.

What remains for us to consider is the relation with the extended molecular layers described by Tsuji (2000). With a molecular excitation temperature of $\approx$ 1000K, it is, of course, outside the observed and modeled halo. The coupling of the extended layers to the star is not well observed. However, Hinkle & Lebzelter (2015) have shown in the LPV R Cas that gas at an intermediate temperature of 1300K shows a behavior that is intermediate with regard to the pulsating halo and the more stationary extended layers. The intermediate zone shows much-reduced mass motions, that are not coherent with the halo pulsation, and are consistent with ”puffs” of mass loss, maintaining line-of-sight coherence over two pulsation cycles, without yet attaining escape velocity and propelling mass loss. While longer spectroscopic time series of additional targets are needed to confirm this finding, it is appears to be a direct observation of the process of populating the extended layers with material raised above the pulsating halo. Expecting a comparable phenomenon in the similar LPV Mira, we suggest that the evidence favors dust formation, or at least critical grain growth, which is likely to occur in either the halo or the extended layers, where a combination of low ambient temperature and relatively enhanced density favors dust formation.

A core-halo structure similar to Fig. 4 has been reported in interferometric imaging of the LPV T Lep (Le Bouquin et al., 2009). An important difference is that the T Lep star is approximately circularly symmetric. The Mira image is not symmetric. In the following, we consider what the asymmetry may have to tell us about the envelope structure and, possibly, about the mass-loss process.

The bright core identified with the photosphere shows an elongated structure. A number of studies have reported strong azimuth dependence in the apparent size of LPV stars. Various models have been proposed, ranging from underlying non-radial oscillation to less dramatic but equally puzzling phenomena such as an azimuthally dependent envelope shape. While non-radial oscillation has not been ruled out for Mira by the new image, the strong asymmetry that is a primary feature of the observation indicates a more complex appearance than a simple elliptical shape of the photosphere or envelope. This is confirmed in the next section with the technical argument that the data are hard to fit solely with an elongated structure and, in fact, an asymmetrical brightness is required.

A natural question is to ask whether this differential effect represents a dimming or a brightening. This can be answered with reasonable confidence by referring to the AAVSO light curve for the last 15 years reproduced in Fig 5. This shows that the 2005 maximum of Mira was one of the faintest in recent years and more than a magnitude fainter than the second following maximum. At the same time, the image appears to show that most of the stellar disk is rather uniformly bright. We tentatively identify the asymmetry with a differential darkening phenomenon.

A related question considers whether the dimming represents a zone of reduced brightness on the photosphere or a region of increased opacity in the envelope. If the dimming is due to a dark region on the star, a possible explanation would be a long-lived convective structure. Convection is expected to be important in LPVs and present at the same time as pulsation. We note that pulsation does not transport energy effectively and convection is still required for energy transport (Wood, 2007). The relative apparent brightness differential across the photosphere (approximately $6\times$), which implies a similar differential in total surface emittance, also seems too large in the H band for a temperature differential. The dark region covers a significant fraction of the hemisphere and even though large convective cells are predicted, this scale seems unprecedented for a convective effect. The 3D hydrodynamical models of AGB stars combine the effects of pulsation and convection, providing evidence that convection plays a role to produce surface brightness asymmetries (Freytag et al., 2017). These simulations show that dark patches cannot be present alone without bright patches. The spatial convective signature is better characterized as localized, high-contrast bright spots associated with hot, rising material, rather than localized dark spots of falling, cool material. Since the present data set do not demonstrate or require a bright patch, we take this as an additional confirmation for a dark region.

The darkening could be due to an increased opacity in the halo or above. This appears more consistent with the apparent large spatial scale - the dark spot appears to significantly obscure a region of the halo. It is also consistent with previous evidence for irregularities in the near circumstellar region of LPVs, including Mira (Lopez et al., 1997). Therefore, we suggest that the Mira image is showing the effects of a localized enhancement in opacity above the halo. This could be consistent with localization within the intermediate layer feeding the extended molecular layers or with that material itself, with coherent structure lasting over several pulsational periods, but not indefinitely consistent with the observed shell expansion velocity $\approx$few km/sec - Hinkle & Lebzelter (2015). It is reasonable to find that the cause is absorption by dust located in or near the extended molecular layers as the dust condensation radius is expected in this region, as discussed in, for example, Reid & Menten 1997, Perrin et al. 2004 & 2015, and Wittkowski et al. 2007. In the next section, an estimate of the mass of the clump with respect to the net mass loss rate leads to another plausible explanation for a location near or somewhat interior to this region, which Tsuji (2000) describe as a MOLsphere.

Previous images of LPVs with large telescope aperture masking (Woodruff et al., 2008) show a rich variety of information about the upper layers of LPVs, but they do not probe asymmetries in the apparent shape of the photosphere well owing to the lack of (u,v) sampling beyond the first visibility null. Multi-telescope interferometry with COAST has reported some significant non-zero closure phases, but did not distinguish the spatial location of the responsible asymmetries (Burns et al., 1998). Hence the inference that asymmetry and inhomogeneity in the circumstellar dust opacity is consistent with previous imaging work, but not well constrained by it. The objective of the next section is to study what constraints on the opacity distribution are provided by the Mira interferometric image.

The image in Fig. 4 clearly shows asymmetries that require a more complex model than the one used in Sect. 3. We have tried several approaches to explain the data with a central source elongated in the north-south direction and various offsets with respect to the more circularly symmetric shell. These did not lead to satisfactory fits as, in particular, they were inconsistent with the closure phases.

Looking at the reconstructed image, however, it appears generally compatible with the simple circularly symmetric model of Sect. 3 in the sense that it is globally symmetric at low spatial frequency scale and with characteristic sizes of 20 and 40 mas for the star and shell around it. With the model temperatures fixed as described above, the inner and outer diameters and the optical depth of the shell are essentially given by the first visibility lobe and the zero crossing around 50 cycles/arcsec. The parameters of the star+shell model are therefore those of Table 3. The star+shell model has five free parameters. It defines a featureless surface brightness distribution $B_{0}(x,y)$. After determining the minimal impact of adding high spatial resolution features to the star+shell model, we decided to separately fit the parameters of Table 3 and the parameters of the more complex model described below. This method with two separate steps led to satisfactory results.

To fit the closure phases and account for the azimuth dependency of the squared visibilities, we considered a dark patch to the west combined with the star+shell model. We tried a single patch as increasing the number of the patches made the fit more complex but did not bring improvement.

Experimenting with several patch options, a smooth edge better matches the data. We used two types of smoothing functions: exponential and Gaussian. Technically, the star+shell model is multiplied by an attenuation function:

The dark patch model then has four parameters: 1) the patch offset $\alpha$ and $\delta$ with respect to the star center; 2) the attenuation strength $a$ (between 0 and 1); 3) the patch FWHM. The total number of parameters of the star+shell+patch model is 9 and the corresponding surface brightness distribution is $B(x,y)$ with $B(x,y)=B_{0}(x,y)\times f(x,y)$.

The closure phases outside the $\pm 150\degree$ range were found to have a diverging impact on the fit as they are very sensitive to noise. We therefore decided to crop them to search for the most robust solution. The best parameters we found are listed in Table 4. The error bars were derived by varying the $\chi^{2}$ by 1. It is clear, however, that these are lower limits since the model does not perfectly account for the data. We did not consider it useful to refine error bars as the goal is to obtain orders of magnitude for physical parameters of the dust clump later in Sect. 5.3 and as we used two different models for the patch which give an idea of the scatter for these parameters as reported in Table 5.

The surface brightness distributions, residues, and closure phase fits are presented in Fig. 6. The fit to the closure phases is quite good, even at high spatial frequencies, with residuals smaller than $25\degree$. The shape of the closure phases is consistent with a smooth spot as a uniform spot did not provide such nice fit. The $V^{2}$ residuals at high visibility are quite large. We investigated this issue by fitting the low spatial frequency data (below 40 cycles/arcsec) with a uniform disk model with an azimuth-dependent diameter ($10\degree$ bins). The amplitude of the scatter remained the same albeit more symmetric around 0. We concluded that the scatter is real and does not depend on the model we chose to interpret the data.

The two models yield qualitatively similar results and although the $\chi^{2}$ for the Gaussian model is a bit lower, both were retained for further consideration. The total $\chi^{2}$ may look large but the two branches of models capture most of the visibility and closure phase features that are in very good agreement with the reconstructed image so that we consider that it is a satisfactory description of the object.

We concluded that the dark patch is probably located above the pulsation/convective zone and is accounted for by absorption from a patch of material. It therefore equally absorbs the radiation from the star and from the inner envelope. We neglect the H-band emission by the patch as it is at a higher, cooler elevation. The attenuation function can be locally equated to an optical depth effect:

The surface brightness distribution then is written as:

We can therefore derive the local optical depth from the ratio of the obscured to unobscured surface brightnesses. The local optical depth can then be converted into a mass of absorbing material. Let $\rho$ [g.cm^-3] be the density of the absorbing material, $\kappa$ [g^-1.cm²] the opacity and $z(x,y)$ the local depth of absorbing material. In the following relation:

where $\rho z(x,y)$ is the mass column density in g.cm^-2. The integration of this quantity over the stellar disk, therefore, yields the total mass $m$ of absorbing material, $m$:

The integral can be evaluating using the surface brightnesses of the star+shell and star+shell+patch models:

We can also evaluate the average attenuation by the patch with the ratio of the integrated surface brightnesses:

The results are given in Table 5. Uncertainties were derived with the Monte-Carlo technique using the values of Table 4. The two patch models yield similar results. The patch causes an average attenuation of the H-band flux of $0.925\pm 0.007$ (minimum relative surface brightnesses of $2\times 10^{-3}$ and $2\times 10^{-1}$ for the exponential and Gaussian respectively, compatible with the dynamic range of the reconstructed image). We use the average value of $(6.31\pm 0.72)\times 10^{25}\,$cm² for the mass-opacity product (optical depth integrated over the star surface), denoted as $m.\kappa_{\mathrm{H}}$ from this point forward.

As discussed in Sect. 4, our favored hypothesis to explain the lower surface brightness at the spot is dimming by a dust patch. Al₂O₃ is a good candidate for dust formation at high temperature in the molecular shell of Mira stars (see e.g. Perrin et al. (2015)). But the condensation of Al₂O₃ at distances and concentrations implied by observations requires high transparency of the grains in the visual and near-IR region to avoid destruction by radiative heating (Höfner et al., 2016). We therefore consider silicate grains to explain the obscuration. Amorphous Mg₂SiO₄ grains are the main candidate for driving the winds of O-rich AGBs according to Bladh & Höfner (2012). We interpret the dimming in H measured with IOTA/IONIC and in V from the AAVSO as being due to this single dust component. Within the scope of the single-grain type hypothesis, the relative V/H extinction constrains the dust grain characteristics. We have decided to consider two types of magnesium-silicate dust grains: MgSiO₃ and Mg₂SiO₄. As a remark, we note that Khouri et al. (2018) did their modeling with aluminum oxides, adding that magnesium-silicate dust grains have scattering properties that are similar to those of Al₂O₃.

Opacities can be obtained from the literature for specific grain sizes and types. We performed calculations of the opacity of dust in the Mie approximation using the Fortran code, bhmie, from Bohren & Huffman (1983), which was modified by Bruce Draine (http://scatterlib.wikidot.com/mie). Optical constants were taken from Jäger et al. (2003). The numerical values are listed at the Jena group web page: http://www.astro.uni-jena.de/Laboratory/OCDB/amsilicates.html (OCDB database).

The extinction efficiency $Q_{\mathrm{ext}}$ was computed at 0.55 and 1.65 $\mu$m wavelengths. In the Mie approximation, the extinction efficiency can be translated to an opacity via the equation:

where $a$ is the radius of the dust sphere and $\rho$ is the density, all in CGS units. The typical density of 3 g.cm^-3 from the OCDB database was used. The opacities were computed for radii ranging from 0.5 nm to 100 $\mu$m to explore the different regimes of the Mie scattering. The plots of Fig 7 show the opacities for the two species of dust as a function of radius at 0.55 and 1.65 $\mu$m as well as the $\kappa_{\mathrm{V}}/\kappa_{\mathrm{H}}$ ratio. The opacities are quite similar for the two species of dust except for radii below 0.01 $\mu$m.

The AAVSO data of Fig 5 show that the two V-band maxima before and after the median date of the observation (JD2453654) are brighter, respectively, by a 0.4 and 0.0 magnitude. Assuming the difference in maximum magnitude is entirely due to the obscuration of a persisting dust cloud, the V-band data give an upper limit on the attenuation due to dust of 0.4 magnitude. We made the following assumptions to compute the flux attenuation by the dust patch in V: 1) the stellar brightness in the visible is dominated by the brightness of the molecular shell, the photosphere is not visible: this is supported by the large optical depth of the shell at visible wavelengths derived both from observations (diameters in the visible are much larger, see for example the first observational evidence by Labeyrie et al. (1977)) and by the 3D hydrodynamical simulations; 2) the diameter of the molecular shell is the same as in H; 3) the dark patch characteristics are the same as at H except for the optical depth.

With these assumptions, the maximum global flux attenuation can be computed as a function of the strength parameter $a$ in Eq. 1. The maximum attenuation is 5% with the Gaussian model, whereas it reaches 10% with the exponential model. In consequence, we used the exponential model to explore the variations of the attenuation and of the $m.\kappa_{\mathrm{V}}$ product versus a. The maximum for the geometrical characteristics of the patch derived in the H band is an attenuation of 10% of the flux corresponding to a difference in magnitude of $\Delta V=0.11$. This is compatible with the maximum brightness decrease with respect to the two neighboring V-band maxima. This cannot, however, be the only source of the dimming in V as the AAVSO data show cycle-to-cycle variations up to a magnitude or more. Unless several such patches are episodically produced. Here, we acknowledge the simplicity of our approach as the reasons for the cycle-to-cycle V magnitude variations of Mira are likely to be multiple and not strictly due to dust patches occurring at the surface of the star. Our simple hypothesis, therefore, provides an upper value for the opacity of the dust patch in V. Following this approach, it yields a maximum mass-opacity product at visual wavelengths of $m.\kappa_{\mathrm{V}}=1.49\,10^{26}\,$cm² and, thus, a maximum ratio for the dust opacity between V and H of:

The relative fraction of V/H extinction assumed does not strongly impact the deductions concerning dust mass. In any event, the inferred V-band obscuration cannot explain the large variations of the maximum visual magnitude shown by the AAVSO measurements.

According to Fig. 7, this is compatible with both MgSiO₃ and Mg₂SiO₄ dust. The particle radii ranges are shown in Table 6, along with the implied total mass of the dust in the patch, and also the total mass of the associated gas - assuming a gas/dust ratio of 89.5 (Knapp, 1985).