# Evidence for localized onset of episodic mass loss in Mira

G. Perrin, S. T. Ridgway, S. Lacour, X. Haubois, É. Thiébaut, J. P. Berger, M. G. Lacasse, R. Millan-Gabet, J. D. Monnier, E. Pedretti, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

G. Perrin, S. T. Ridgway, S. Lacour, X. Haubois, É. Thiébaut, et al., Evidence for localized onset of episodic mass loss in Mira. Astronomy and Astrophysics - A&A, 2020, 642, pp.A82. 10.1051/0004-6361/202037443. hal-02962090

HAL Id: hal-02962090

https://hal.science/hal-02962090

Submitted on 8 Oct 2020

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Evidence for localized onset of episodic mass loss in Mira

G. Perrin<fnref n="1" />, S. T. Ridgway<fnref n="1" />,<fnref n="2" />, S. Lacour<fnref n="1" />, X. Haubois<fnref n="3" />, É. Thiébaut<fnref n="4" />, J. P. Berger<fnref n="5" />, M. G. Lacasse<fnref n="6" />, R. Millan-Gabet<fnref n="7" />, J. D. Monnier<fnref n="8" />, E. Pedretti<fnref n="9" />, S. Ragland<fnref n="10" />, and W. Traub<fnref n="11" />

Received 3 January 2020 / Accepted 14 August 2020

## ABSTRACT

**Context.** Mass loss from long-period variable stars (LPV) is an important contributor to the evolution of galactic abundances. Dust formation is understood to play an essential role in mass loss. It has, however, proven difficult to develop measurements that strongly constrain the location and timing of dust nucleation and acceleration.

**Aims.** Interferometric imaging has the potential to constrain the geometry and dynamics of mass loss. High angular resolution studies of various types have shown that LPVs have a distinct core-halo structure. These have also shown that LPV images commonly exhibit a non-circular shape. The nature of this shape and its implications are yet to be understood.

**Methods.** Multi-telescope interferometric measurements taken with the Interferometric Optical Telescope Array (IOTA) provide imagery of the LPV Mira in the *H*-band. This wavelength region is well suited to studying mass loss given the low continuum opacity, which allows for emission to be observed over a very long path in the stellar atmosphere and envelope.

**Results.** The observed visibilities are consistent with a simple core-halo model to represent the central object and the extended molecular layers but, in addition, they demonstrate a substantial asymmetry. An analysis with image reconstruction software shows that the asymmetry is consistent with a localized absorbing patch. The observed opacity is tentatively associated with small dust grains, which will grow substantially during a multi-year ejection process. Spatial information along with a deduced dust content of the cloud, known mass loss rates, and ejection velocities provide evidence for the pulsational pumping of the extended molecular layers. The cloud may be understood as a spatially local zone of enhanced dust formation, very near to the pulsating halo. The observed mass loss could be provided by several such active regions around the star.

**Conclusions.** This result provides an additional clue for better understanding the clumpiness of dust production in the atmosphere of AGB stars. It is compatible with scenarios where the combination of pulsation and convection play a key role in the process of mass loss.

**Key words.** stars: atmospheres – stars: AGB and post-AGB – stars: mass-loss – techniques: interferometric – infrared: stars – stars: individual: Mira

## 1. Introduction

All stars of low to intermediate mass evolve through the AGB phase, with a dredging up of synthesized elements, and undergo extensive mass loss as long-period variables (LPVs). This return of heavy elements to the interstellar medium constitutes a critical step in the evolution of stars and galaxies. The LPV is a very complex phenomenon and our understanding of it remains incomplete. Both modeling and observations continue to guide our progress in this area.

LPVs vary in brightness cyclically, with periods on the order of a few hundred days. Photometric monitoring shows that these periodic variations are often slightly or even severely irregular.

The improving capabilities at present for direct imagery of LPVs show that the LPV atmosphere is asymmetric, at least in some cases and potentially for all. Speckle and other imaging techniques have been used to obtain high angular resolution information about LPVs with ground-based telescopes and interferometers, as well as with HST. Scholz (2003) gives an extensive list of such observations prior to 2003. Creech-Eakman & Thompson (2009) summarize a seven-year study of 100 LPVs using the Palomar Testbed Interferometer and report extensive results for apparent sizes and their variation with position angle and time. Ragland et al. (2006), using the Interferometric Optical Telescope Array (IOTA), found that 1/3 of 56 nearby LPV stars surveyed showed asymmetry. The authors speculated that

<footnote n="1">¹ LESIA, Observatoire de Paris, Université PSL, CNRS, Sorbonne Université, Université de Paris, 5 place Jules Janssen, 92195 Meudon, France
e-mail: guy.perrin@obspm.fr</footnote>

<footnote n="2">² NSF's NOIRLab, 950 N. Cherry Ave., Tucson, AZ 85719, USA</footnote>

<footnote n="3">³ European Southern Observatory, Alonso de Córdova 3107, Vitacura, Casilla 19001, Santiago de Chile, Chile</footnote>

<footnote n="4">⁴ Université de Lyon, Université Lyon1, ENS de Lyon, CNRS, Centre de Recherche Astrophysique de Lyon UMR 5574, 69230 Saint-Genis-Laval, France</footnote>

<footnote n="5">⁵ Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France</footnote>

<footnote n="6">⁶ Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA</footnote>

<footnote n="7">⁷ Infrared Processing and Analysis Center, California Institute of Technology, Pasadena, CA 91125, USA</footnote>

<footnote n="8">⁸ University of Michigan, 941 Dennison Building, 500 Church Street, Ann Arbor, MI 48109-1090, USA</footnote>

<footnote n="9">⁹ RAL Space, STFC Rutherford Appleton Laboratory, Harwell, Didcot, OX11 0QX, UK</footnote>

<footnote n="10">¹⁰ W. M. Keck Observatory, 65-1120 Mamalahoa Hwy, Kamuela, HI 96743, USA</footnote>

<footnote n="11">¹¹ Jet Propulsion Laboratory, California Institute of Technology, M/S 301-451, 4800 Oak Grove Dr., Pasadena, CA 91109, USA</footnote>
all of them would appear asymmetric if studied with sufficient resolution. This trend was confirmed by the observations of Wittkowski et al. (2011) at VLTI with AMBER which measured deviations from centro-symmetry on four Mira variables. Later, Wittkowski et al. (2016) used AMBER to observe large-scale inhomogeneities that contribute a few percent of the total flux in six Mira variables. The detailed structure of the inhomogeneities show a variability on time scales of three months and above. This is compatible with recent 3D dynamical simulations of AGB stars which deal self-consistently with convection and pulsations (Freytag et al. 2017). Convection, together with pulsational shocks, contributes to the raising of material up to altitudes where dust can condense. This knowledge leads us to consider whether the spatial inhomogeneities revealed by observations and predicted by models could also lead to inhomogeneous zones of dust formation.

Ohnaka et al. (2016) used visible polarimetric imaging observations of W Hya taken with VLT/SPHERE-ZIMPOL to detect three clumpy dust clouds at two stellar radii. The authors deduced large dust grains (0.4–0.5 µm) of Al₂O₃, Mg₂SiO₄, or MgSiO₃. The grain size is consistent with the prediction of hydrodynamical models for mass loss driven by scattering due to Fe-free micron-sized grains (Höfner 2008). The detection of clumpy dust clouds close to the star lends support to dust formation induced by pulsation and large convective cells. Observations at similar spatial scales but at longer wavelengths in the millimeter range have also revealed asymmetries. Wong et al. (2016) detected bright SiO-emitting clumps with ALMA in Mira at distances around two stellar radii. Also with ALMA on Mira, Kamiński et al. (2017) revealed the anisotropy of the distribution of TiO₂ in the atmosphere of Mira up to 5.5 stellar radii. This may suggest inhomogeneous distributions for both dust and its precursors near the photosphere, which is compatible with clumpy nucleation sites.

The shapes suggest a structure that may eventually contribute to the form of the planetary nebula (PN) expected to follow the AGB phase. However, PNe almost certainly represent the shape of the outer envelope and beyond, arising perhaps at some distance from the asymmetrical physical phenomena occurring at yet deeper layers which give rise to them. In addition, part of the shape may be due to binarity. Measurements with higher angular resolution at wavelengths that can penetrate deeply into the envelope promise to stand as the best opportunity to tie together the most reliable part of the LPV theory – the unobservable interior processes – with the complexity of the observed exterior.

Here, we report the continuation of a series of LPV studies carried out at IOTA (Traub 1998), which offers excellent (u, v) coverage, which can be enhanced with multi-spectral measurements. In Perrin et al. (2004), we showed strong evidence in favor of the presence of a very deep, but surprisingly transparent, extended atmosphere that is analogous to, and perhaps similar to, the MOLsphere, which Tsuji ascribed to luminous, non-Mira stars. This interpretation was supported by other observations, for example Lacour et al. (2009) and Le Bouquin et al. (2009). Since dynamical models reproduce the characteristics of Tsuji's MOLsphere, without, however, an actual discrete layer (Bladh et al. 2015), the historical term is no longer understood as descriptively correct. We follow our referee's recommendation to describe this region as "extended molecular layers".

In our previous studies, we utilized pair-wise combination of telescopes and only reported the analysis based on fringe visibilities and spherical models. In this study, we use the IONIC beam combiner instrument of IOTA (Berger et al. 2003), which combines three telescopes simultaneously and measures closure

phases. Squared visibilities and closure phases can be combined to reconstruct images. We present the observations and data reduction in Sect. 2. The data are first fit with a circularly symmetric model in Sect. 3. These data have allowed us to synthesize an image of the LPV Mira (omicron Ceti), presented in Sect. 4, which is sufficiently indicative of the asymmetric structure to motivate inclusion of an asymmetric component in our models and to derive some related parameters in Sect. 5. We follow with a discussion of how this model relates to other measurements and how it may fit in the context of expected phenomena in LPV oscillation and mass loss in Sect. 6. Finally, we comment on the future of LPV imagery in Sect. 7.

## 2. Observations and data reduction

### 2.1. Description of IOTA observations

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 × 15 m, corresponding to an angular resolution of 10 × 23 milliarcseconds at 1.65 µ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 µm ≤ λ ≤ 1.8 µ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<fnref n="1" />. 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.

### 2.2. Data reduction

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

<footnote n="1">¹ http://cdsarc.u-strasbg.fr/afoev/</footnote>
<table><caption>**Table 1.** Log of the observations of Mira with IONIC on IOTA.</caption><thead><tr><th>Date (UT)</th><th>Visual phase</th><th>Configuration</th><th>Baseline lengths (m)</th></tr></thead><tbody><tr><td>2005 October 5, 6, 7</td><td>0.41</td><td>A5-B5-C0</td><td>5.00, 5.00, 7.07</td></tr><tr><td>2005 October 8</td><td>0.42</td><td>A5-B15-C0</td><td>5.00, 15.00, 15.81</td></tr><tr><td>2005 October 10, 11</td><td>0.43</td><td>A15-B15-C0</td><td>15.00, 15.00, 21.21</td></tr><tr><td>2005 October 12, 13</td><td>0.44</td><td>A25-B15-C0</td><td>15.00, 25.00, 29.15</td></tr><tr><td>2005 October 15</td><td>0.45</td><td>A30-B15-C15</td><td>15.00, 21.21, 33.54</td></tr></tbody></table>

**Notes.** The visual phase was computed from data from the AFOEV. The configurations give the positions of the three telescopes in meter with respect to the corner of the array.

<table><caption>**Table 2.** Reference stars used to calibrated the Mira data.</caption><thead><tr><th>Star</th><th>Spectral type</th><th>Uniform disk diameter (mas)</th></tr></thead><tbody><tr><td>HD 8512</td><td>K0 IIb</td><td>2.69 ± 0.030 mas</td></tr><tr><td>HD 16212</td><td>M0 III</td><td>3.02 ± 0.032 mas</td></tr></tbody></table>

**Notes.** Diameters are from Borde et al. (2002).

<figure><img src="image_3.png" /><figcaption>**Fig. 1.** $(u, v)$ coverage achieved on Mira during the observations with IOTA/IONIC. Each point marks an exposure.</figcaption></figure>

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 Figs. 2 and 3, respectively. The colors show different baseline

<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Squared visibilities of Mira measured with IOTA/IONIC vs. spatial frequency. The color code gives the orientation of the baselines arranged in four bins: [−88°,−64°] (blue), [−59°,−34°] (red), [40°,63°] (green), [65°,59°] (yellow). The color code directly shows that the source is not circularly symmetric.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Closure phases of Mira measured with IOTA/IONIC. The closure phases are plotted against the maximum spatial frequency of each triangle. The closure phases are compatible with $0(\pi)$ up to the first zero crossing of the squared visibilities near 50 cycles arcsec$^{-1}$. Beyond that spatial frequency, closure phases depart from $0(\pi)$ meaning the source is not centro-symmetric at scales smaller than the typical radius of the star.</figcaption></figure>

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$^{-1}$, 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.

### 3. Model fitting with a circularly symmetric model

The $V^2$ data were first compared to a model star photosphere + partially transparent shell (Perrin et al. 2004) to derive the main parameters of the model: photosphere diameter, shell diameter,
<table><caption>**Table 3.** Parameters of the star+shell model fit to the V<fnref n="2" /> data.</caption><thead><tr><th>D<sub>★</sub></th><th>T<sub>★</sub></th><th>D<sub>shell</sub></th><th>T<sub>shell</sub></th><th>τ</th></tr></thead><tbody><tr><td>21.1 mas</td><td>3400 K</td><td>40.00 mas</td><td>2000 K</td><td>0.70</td></tr></tbody></table>

shell optical depth, shell temperature, and photosphere temperature. The degeneracy of the temperature fit is solved by fitting the flux. For this purpose, we used the phases from the AFOEV data, as shown in Table 1. For the average photometric phase of 0.43 the magnitudes interpolated from Whitelock et al. (2000) are *H* (−1.7 to −1.8) and *K* (−2.2 to −2.3). We note that the infrared fluxes are much less variable than the visible fluxes. The *H* and *K* bands gave weak constraints on the photospheric temperature, probably because the photospheric diameter is not uniquely determined. The flux gave a minimum value of 3000 K, and it was set to 3400 K, which is a rough average of the temperatures measured by Perrin et al. (2004) at earlier phases. This yields the global parameters for the model in Table 3.

# 4. Imaging of Mira

## 4.1. Image reconstruction

An image was reconstructed using the monochromatic version of the MiRA software of Thiébaut (2008)<fnref n="2" />. 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 100 × 100 pixel were reconstructed. A hyper parameter μ controls the weight of the regularization in the reconstructed image with a total χ<fnref n="2" /> reading: $\chi^2_{\text{tot}} = \chi^2_{\text{data}} + \mu\chi^2_{\text{rgl}}$. $\chi^2_{\text{data}}$ is defined and discussed in Appendix A. The higher μ 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. A.1. 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 ≈ 20 mas elongated in the north-south direction; (2) surrounded by a fainter environment of ≈ 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, μ = 10⁸, as it is more constrained by the data relative to the others. The image is shown in Fig. 4.

## 4.2. The inner envelope

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

<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Image of Mira reconstructed with the MiRA software. A dark patch is visible west of the bright central peak. The scale is linear with 0 intensity for the background and 5 for the peak. This image is one of six images presented in Appendix A.</figcaption></figure>

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⁻¹ Hinkle et al. (1984) for half of the 321-day period, the total extent of the pulsational motion is found to be 1.1 × 10⁸ km. This corresponds well with the 1.4 × 10⁸ km radius difference between the top of the 10 mas model halo and the photosphere (based on the HIPPARCOS distance of 2.8 × 10¹⁵ 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, ≈2700 K, are in between the inferred model layer temperatures of 2000 K and 3400 K.

What remains for us to consider is the relation with the extended molecular layers described by Tsuji (2000). With a molecular excitation temperature of ≈1000 K, 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 1300 K 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

<footnote n="2">² https://github.com/emmt/MiRA</footnote>
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** AAVSO V-band data for the period around the dates of the IOTA/IONIC observations.</figcaption></figure>

approximately circularly symmetric. The Mira image is not sym-
metric. In the following, we consider what the asymmetry may
have to tell us about the envelope structure and, possibly, about
the mass-loss process.

4.3. Interpretation of the image

The bright core identified with the photosphere shows an elon-
gated 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 yr 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 fol-
lowing 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 (approx-
imately 6×), which implies a similar differential in total surface
emittance, also seems too large in the H band for a tempera-
ture differential. The dark region covers a significant fraction of
the hemisphere and even though large convective cells are pre-
dicted, 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 con-
vection 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 evi-
dence 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 localiza-
tion 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 ≈few km s⁻¹ – Hinkle
& Lebzelter (2015). It is reasonable to find that the cause is
absorption by dust located in or near the extended molecular lay-
ers 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 mask-
ing (Woodruff et al. 2008) show a rich variety of information
about the upper layers of LPVs, but they do not probe asym-
metries 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 inhomogene-
ity 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.

**5. Obscuring dust patch**

5.1. Dust patch model

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 gen-
erally compatible with the simple circularly symmetric model of
Sect. 3 in the sense that it is globally symmetric at low spatial fre-
quency 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 opti-
cal depth of the shell are essentially given by the first visibility
lobe and the zero crossing around 50 cycles arcsec⁻¹. The param-
eters 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₀(x, y). After determining the
<table><caption>**Table 4.** Best parameters for the dark patch model fit of the data with an exponential and a Gaussian attenuation function.</caption><thead><tr><th rowspan="2"></th><th colspan="2">Dark patch model</th></tr><tr><th>Gaussian</th><th>Exponential</th></tr></thead><tbody><tr><td>FWHM (mas)</td><td>8.57 ± 0.40</td><td>6.89 ± 0.65</td></tr><tr><td>α (mas)</td><td>6.94 ± 0.23</td><td>7.84 ± 0.30</td></tr><tr><td>δ (mas)</td><td>-1.69 ± 0.15</td><td>-2.04 ± 0.21</td></tr><tr><td>a</td><td>0.783 ± 0.005</td><td>1.00<sup>+0.00</sup><sub>-0.008</sub></td></tr><tr><td>χ<sup>2</sup><sub>V2</sub></td><td>134</td><td>165</td></tr><tr><td>χ<sup>2</sup><sub>CP</sub></td><td>83</td><td>69</td></tr><tr><td>χ<sup>2</sup><sub>tot</sub></td><td>124</td><td>146</td></tr></tbody></table>

<table><caption>**Table 5.** Comparison of the average attenuations $\bar{f}$ and integrated optical depths $m.\kappa$ for the two models of patches.</caption><thead><tr><th>Patch model</th><th>Exponential</th><th>Gaussian</th></tr></thead><tbody><tr><td>$\bar{f}$</td><td>0.927 ± 0.007</td><td>0.924 ± 0.005</td></tr><tr><td>$m.\kappa$ (cm<sup>2</sup>)</td><td>(6.05 ± 0.72) × 10<sup>25</sup></td><td>(6.57 ± 0.67) × 10<sup>25</sup></td></tr></tbody></table>

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:

$$f(x + \alpha, y + \delta) = \begin{cases} 1 - a.e^{-\frac{2\ln2 \sqrt{2+y^2}}{FWHM}}, & \text{if exponential patch} \\ 1 - a.e^{-\frac{4\ln2(x^2+y^2)}{FWHM^2}}, & \text{if Gaussian patch} \end{cases}. \quad (1)$$

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 ±150° 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°. 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$^{-1}$) with a uniform disk model with an azimuth-dependent diameter (10° 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.

## 5.2. Interpretation of the attenuation function

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:

$$f(x, y) = e^{-\tau(x, y)}. \quad (2)$$

The surface brightness distribution then is written as:

$$B(x, y) = B_0(x, y) \times e^{-\tau(x, y)}. \quad (3)$$

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$^2$] the opacity and $z(x, y)$ the local depth of absorbing material. In the following relation:

$$\frac{\tau(x, y)}{\kappa} = \rho z(x, y), \quad (4)$$

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$:

$$m = \iint \rho z(x, y) \, dx dy. \quad (5)$$

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

$$m = \frac{1}{\kappa} \iint -\ln \left[ \frac{B(x, y)}{B_0(x, y)} \right] dx dy. \quad (6)$$

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

$$\bar{f} = \frac{\iint B(x, y) \, dx dy}{\iint B_0(x, y) \, dx dy}. \quad (7)$$

The results are given in Table 5. Uncertainties were derived with the Monte-Carlo technique using the values of Table 4.
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Result of the fits with the patchy model. Left: exponential attenuation model. Right: Gaussian attenuation model. From top to bottom: surface brightness distribution, V² residues, closure phases (red), model closure phases (black), and closure phases residues. The closure phase and closure phase residues are displayed as a function of the maximum spatial frequency of the triangle of baselines.</figcaption></figure>

The two patch models yield similar results. The patch causes an average attenuation of the H-band flux of 0.925 ± 0.007 (minimum relative surface brightnesses of 2 × 10⁻<fnref n="3" /> and 2 × 10⁻¹ for the exponential and Gaussian respectively, compatible with the dynamic range of the reconstructed image). We use the average value of (6.31 ± 0.72) × 10²⁵ cm² for the mass-opacity product (optical depth integrated over the star surface), denoted as m.κH from this point forward.

### 5.3. Interpretation of the attenuation as a dust opacity effect

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<fnref n="3" />. Optical constants were taken from Jäger et al. (2003). The numerical values are listed at the Jena group web page: (OCDB database<fnref n="4" />).

The extinction efficiency Qext was computed at 0.55 and 1.65 µm wavelengths. In the Mie approximation, the extinction efficiency can be translated to an opacity via the equation:

$$ \kappa = \frac{3Q_{\text{ext}}}{4\alpha\rho}, \qquad (8) $$

where a is the radius of the dust sphere and ρ is the density, all in CGS units. The typical density of 3 g cm⁻<fnref n="3" /> from the OCDB database was used. The opacities were computed for radii ranging from 0.5 nm to 100 µ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 µm

<footnote n="3">³ http://scatterlib.wikidot.com/mie</footnote>

<footnote n="4">⁴ http://www.astro.uni-jena.de/Laboratory/OCDB/amsilicates.html</footnote>
<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Result of computations done with the *bhmie* code for MgSiO₃ and Mg₂SiO₄ dust with grain sizes ranging from 0.5 nm and 100 µm. *Left:* opacities at V and H. *Right:* $κ_V/κ_H$ ratio.</figcaption></figure>

<table><caption>**Table 6.** Dust parameters and mass loss rates for silicate species for small or large grains.</caption><thead><tr><th>Particle size range</th><th>Opacity range<br>κ<sub>H</sub> (cm<sup>2</sup> g<sup>-1</sup>)</th><th>Dust mass<br>(M<sub>☉</sub>)</th><th>Gas mass<br>(M<sub>☉</sub>)</th></tr></thead><tbody><tr><td>1–5.5 nm</td><td>0.6–1.45</td><td>2.2–5.3 × 10<sup>-8</sup></td><td>2.0–4.7 × 10<sup>-6</sup></td></tr><tr><td>0.45–0.47 µm</td><td>1.07–1.27 × 10<sup>6</sup></td><td>2.5–3.0 × 10<sup>-14</sup></td><td>2.2–2.7 × 10<sup>-12</sup></td></tr></tbody></table>

as well as the $κ_V/κ_H$ ratio. The opacities are quite similar for the two species of dust except for radii below 0.01 µ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·κ_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 $Δ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·κ_V = 1.49 × 10^{26}$ cm² and, thus, a maximum ratio for the dust opacity between V and H of:

$$ \frac{\kappa_V}{\kappa_H} \le \frac{1.49 \times 10^{26}}{(6.31 \pm 0.72) \times 10^{25}} \approx 2.36 \pm 0.27. \quad (9) $$

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).

## 6. Discussion

In the previous sections, we modeled the IOTA/IONIC data on Mira in the H band using a patch of absorbing material obscuring the core-halo, that is, the modeled H-band photosphere and pulsating envelope. We show, based on a combination of these data with visual photometric data, that this material is compatible with silicate dust. As summarized in Table 6, the inferred dust mass associated with the absorbing material is strongly dependent on the dust physical characteristics. Our data set is limited and the characteristics of the dust remain rather speculative in this analysis and further interpretation must rely on indirect evidence.

Since micron size grains are commonly inferred from models of the dust at several AU from an LPV (Ohnaka et al. 2016), such a grain size might be favored for the observed patch. However, from Table 6, the associated dust/gas mass would be only $\approx 10^{-12} M_☉$. This is much less than the known yearly mass loss rate for Mira of $3 \times 10^{-7} M_☉ yr^{-1}$ (Ryde et al. 2000). In order for such a patch to contribute significantly to the observed mass loss
rate, it would be necessary to suggest an implausibly large num-
ber of patches or high frequency and speed of ejection. So we
tentatively reject the large grain option for the observed patch.

There is a precedent for considering small grains. Recent
studies of the mira W Hya by Ohnaka et al. (2017) show a shift
toward smaller grain sizes near photometric minimum (0.54),
compared to their earlier study at 0.77 (Ohnaka et al. 2016).
This is not surprising since the minimum brightness corresponds
approximately to maximum expansion and greatest atmospheric
cooling. Dynamical models reported by Höfner et al. (2016),
Aronson et al. (2017) and Liljegren et al. (2017) indeed show
dust formation near the brightness minimum. The observations
reported here are also near minimum.

From Table 6, a small grain size would be consistent with
a reservoir of gas/dust mass in the patch ≈10 times larger than
the annual stellar mass loss if the gas/dust ratio can be applied at
short distance from the photosphere and comparable to or a bit
below the yearly mass loss rate otherwise. Therefore, it is reason-
able to consider the possibility that the patch can be linked with
the mass ejection process and what this implies for the ejection
efficiency.

Only a single patch was observed. However, from the model
suggesting that the patch is near but outside the halo and from
the knowledge that approximately 1/3 of the progeny of LPVs
show an asymmetry (Bobrowsky 2007, for Mira, here, attributed
to dust patches), it can be estimated from geometry and assumed
random distribution that at any given time there may be on the
order of five such patches around Mira, of which zero or one
would typically be in the line of sight to the halo at any given
time. Thus, the elevation of a few percent of each of the mass
of each of these “puffs” of material over each pulsation cycle,
on average, could support the observed Mira mass loss rate.
Chandler et al. (2007) reported their detection of episodic mass
loss around Mira at 11.15 µm with the ISI interferometer. Our
detection in H band is compatible with their finding.

With low envelope expansion velocities, compared to the
escape velocity, and being on the same order as those observed
in the extended molecular layers, puffs of rising dust-enhanced
material would have ample time for grain growth during the
elevation to the vicinity of 2+ stellar radii, where much larger
grain sizes (with much reduced opacity) are found by model-
ing spatially resolved polarimetric observations: Norris et al.
(2012) have detected 0.3 µm grains of Mg₂SiO₄ or MgSiO₃ for
W Hya, R Dor and R Leo at 2 stellar radii or less; Ohnaka
et al. (2016) detected 0.4–0.5 µm grains of Al₂O₃, Mg₂SiO₄ or
MgSiO₃ between ≈2–3 stellar radii in clumpy dust clouds for W
Hya; Adam & Ohnaka (2019) detected 0.1 µm grains of Al₂O₃,
Mg₂SiO₄, or MgSiO₃ in a clumpy shell around IK Tau between
2 and 5 stellar radii; Khouri et al. (2016) measured the asymmet-
ric photosphere of R Dor attributed to scattering by dust. Finally,
Khouri et al. (2018) observed clumpy structures in visible polar-
ized light around Mira at larger distance than our measurements
and within ≈8 photospheric radii. They used aluminum oxide
grains to model their data and concluded that aluminum could be
in dust phase if the grain size is less than 0.02 µm. All other exist-
ing observations, apart from our own and those of Khouri et al.
(2018), detect a few clouds of larger dust particles that account
for a good fraction of the mass loss. This is a discrepancy that is
to be resolved with more measurements. The greater grain size at
higher elevation satisfies the preferred conditions to drive mass
loss by radiation pressure on the grains (Höfner 2008), that is,
0.36–0.66 µm at 2–3 R<sub>*</sub>.

These findings give rise to a few more questions. The first
two regard the nature of the mechanism that could explain the

existence of such a clump and how long can it exist in the
atmosphere or above. A binary companion sufficiently close
to Mira A could reduce the local effective gravity enough to
enhance mass loss and it could also account for the symmetric
structures commonly observed in planetary nebulae. The known
companion Mira B is too far away to provide a gravity reduction
by more than ≈1%. A much closer companion would be needed,
even a massive planet, but there is no evidence for another very
close companion at this point. We find the dust clump to be on
the side opposite to Mira B as Chandler et al. (2007) did for
the excess of dust emission measured with ISI. Khouri et al.
(2018) did not find a systematic asymmetry at larger distance.
They found a dust trail connecting Mira A and Mira B; how-
ever, there is no clear evidence for an influence of Mira B on the
location of sites for dust production.

The size of the clump is larger than convective cells, such
as those simulated in 3D hydrodynamical codes. However, the
collective impact of complex convective structure could produce
larger scale irregularities. The coupling between pulsations and
convection could be invoked to explain such a clump of dust.
A pioneer study of this coupling was performed by Anand &
Michalitsanos (1976) with a mechanical analog of a giant or
a supergiant star. Their conclusion was that strong interactions
can be expected with asymmetrical fluctuations of the surface
brightness. Kiss et al. (2000) invoked this mechanism as a pos-
sible explanation for period and amplitude changes revealed
in photometric monitoring. A single pulsation with effective
convective coupling in the vicinity of a single or few large con-
vective cells could well enhance the efficiency of mass ejection
into the envelope, generating an over-density. This over-density
would presumably follow a ballistic trajectory, with much of the
material falling back toward the photosphere, where it would
encounter the next levitating shock. It is easy to suppose that
this pattern could preserve the over-density for some time, until
it is reduced by leakage up, especially into the extended molec-
ular layers, with a timeframe that is consistently similar to the
pulsational period.

Freytag & Höfner (2008) performed 3D numerical simula-
tions of an AGB star combining convection, pulsation, and dust
formation which confirm the Anand & Michalitsanos (1976) pre-
dictions. In their simulations, the bright convective cells provide
a strongly asymmetric brightness distribution at and above levels
they describe to be the top of the convective envelope which we
believe corresponds approximately to what we have referred to
as the core, or the depth to which near-IR imagery can penetrate.
Asymmetries in the density distribution appear most strongly
in the highest levels, which are detected in the near-IR through
their partial transparency and molecular absorption, which vary
through the pulsational period. Freytag & Höfner (2008) note a
typical lifetime for convective structure of years, compared with
a pulsational time of one year. Something that could prove par-
ticularly significant for the upper envelope opacity, Freytag &
Höfner (2008) predicted localized dust formation in regions of
enhanced density and depressed temperature which occur natu-
rally in the complex flows and shocks. This has been confirmed
by Höfner & Freytag (2019) with more recent 3D simulations
of atmospheric shock waves generated by large-scale convective
flows and pulsations. They find that the clumpy environment
of dust can form around M-type AGB stars. This is clearly
compatible with what we report here and with the result of
Ohnaka et al. (2016), Adam & Ohnaka (2019), Norris et al.
(2012). This is also compatible with the dimming of the visual
photometry at the time of our observations and during the few
cycles that came before.
## 7. Conclusion and outlook

The Mira image we obtained in *H* band at high spatial resolution offers supporting evidence for a possible outcome of pulsation-convection coupling and offers a possible direction for pursuing the evidence further.

Confirmation should be possible with a series of interferometric images obtained with spectral information and polarimetry, possibly combined with visible resolved polarimetry, as already shown in some works cited in this paper. The three requirements include: excellent image quality, a combination between spectral and spatial resolution, and a suitable time series with a few month time interval over a few periods, that is, over several years. The quality of the images can be addressed in two ways.

First, for the higher spatial frequency baselines, good (*u*,*v*) coverage is essential. This is possible at present only at the VLTI, where the large selection of baselines could be optimized for a target list, although the need for rather contemporaneous UV coverage may be a challenge. Second, as is clearly shown by the aperture masking and other results mentioned above, these stars have a lot of low spatial frequency power. It is essential to measure this correctly in order to properly interpret the high frequency measurements. This can be provided in the visible or near-infrared with extreme adaptive optics as shown by Ohnaka et al. (2016) with SPHERE or by aperture masking like the SAM mode of SPHERE on the VLT (Cheetham et al. 2016) and with MICADO on the ELT in the future (Lacour et al. 2014), or with short-exposure imaging with adaptive optics in a moderate speckle mode (Kervella et al. 2009). Mid-infrared wavelengths can be reached without adaptive optics (Poncelet et al. 2007) in this same mode. Near-IR polarimetric measurements are needed to simultaneously constrain scattering and absorption effects of dust. Spectral resolution is required to constrain temperatures and to separate the continuum (or pseudo-continuum) from the regions of the spectrum dominated by molecular opacity. Spectral resolutions as high as 4000 can be obtained with GRAVITY on the VLTI (GRAVITY Collaboration 2017) in the *K* band. With such data, we could envisage the detection of a correlation between the brightness distribution on the photosphere and the evolving opacity in the envelope over subsequent years. The *Freytag & Höfner* (2008) models specifically predict episodes of dust formation in the envelope, as close as 1–2 stellar radii from the surface. The combination of continuum and molecular opacities will help distinguish between dust and other continuous opacities. Some preliminary results on correlations between the visual light curve and the diameter measured in different band passes in the *K* band were recently published by *Wittkowski* et al. (2018) using GRAVITY. In addition, these data could be combined with mid-infrared four-telescope observations with the MATISSE instrument on VLTI (Lopez et al. 2014) to further characterize the properties of the dust.

*Acknowledgements.* This research has made use of NASA's Astrophysics Data System Bibliographic Services, the SIMBAD and AFOEV databases, operated at CDS, Strasbourg, France and of the AAVSO database. The work of STR is supported by NOIRLab, which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation. We thank the anonymous referee for suggestions which improved the content and presentation of the report.

## References

* Adam, C., & Ohnaka, K. 2019, A&A, 628, A132

* Anand, S. P. S., & Michalitsanos, A. G. 1976, Ap&SS, 45, 175

* Aronson, E., Bladh, S., & Höfner, S. 2017, A&A, 603, A116

* Baldwin, J. E., Beckett, M. G., Boysen, R. C., et al. 1996, A&A, 306, L13

* Berger, J.-P., Haguenauer, P., Kern, P. Y., et al. 2003, Proc. SPIE, 4838, 1099

* Bladh, S., & Höfner, S. 2012, A&A, 546, A76

* Bladh, S., Höfner, S., Aringer, B., & Eriksson, K. 2015, A&A, 575, A105

* Bobrowsky, M. 2007, Asymmetrical Planetary Nebulae IV (San Francisco: ASP), 22

* Bohren, C. F., & Huffman, D. R. 1983, Absorption and scattering of light by small particles (New York: John Wiley & Sons)

* Borde, P., Coudé du Foresto, V., Chagnon, G., & Perrin, G. 2002, VizieR Online Data Catalog: J/A+A/393/183

* Burns, D., Baldwin, J. E., Boysen, R. C., et al. 1998, MNRAS, 297, 462

* Chandler, A. A., Tatebe, K., Wishnow, E. H., Hale, D. D. S., & Townes, C. H. 2007, ApJ, 670, 1347

* Cheetham, A. C., Girard, J., Lacour, S., et al. 2016, Proc. SPIE, 9907, 99072T

* Coudé du Foresto, V., Ridgway, S., & Mariotti, J.-M. 1997, A&AS, 121, 379

* Creech-Eakman, M. J., & Thompson, R. R. 2009, ASP, 412, 149

* Freytag, B., & Höfner, S. 2008, A&A, 483, 571

* Freytag, B., Liljegren, S., & Höfner, S. 2017, A&A, 600, A137

* GRAVITY Collaboration (Abuter, R., et al.) 2017, A&A, 602, A94

* Hinkle, K. H., & Lebzelter, T. 2015, EAS Pub. Ser., 71, 249

* Hinkle, K. H., Scharlach, W. W. G., & Hall, D. N. B. 1984, ApJS, 56, 1

* Höfner, S. 2008, A&A, 491, L1

* Höfner, S., & Freytag, B. 2019, A&A, 623, A158

* Höfner, S., Bladh, S., Aringer, B., & Ahuja, R. 2016, A&A, 594, A108

* Ireland, M. J., Monnier, J. D., & Thureau, N. 2006, SPIE Conf. Ser., 6268, 62681T

* Jäger, C., Il'in, V. B., Henning, T., et al. 2003, J. Quant. Spectr. Rad. Transf., 79, 765

* Kamiński, T., Müller, H. S. P., Schmidt, M. R., et al. 2017, A&A, 599, A59

* Kervella, P., Verhoelst, T., Ridgway, S. T., et al. 2009, A&A, 504, 115

* Khouri, T., Maercker, M., Waters, L. B. F. M., et al. 2016, A&A, 591, A70

* Khouri, T., Vlemmings, W. H. T., Olofsson, H., et al. 2018, A&A, 620, A75

* Kiss, L. L., Szatmáry, K., Szabó, G., & Mattei, J. A. 2000, A&AS, 145, 283

* Knapp, G. R. 1985, ApJ, 293, 273

* Labeyrie, A., Koechlin, L., Bonneau, D., Blazit, A., & Foy, R. 1977, ApJ, 218, L75

* Lacour, S., Meimon, S., Thiébaut, E., et al. 2008, A&A, 485, 561

* Lacour, S., Thiébaut, E., Perrin, G., et al. 2009, ApJ, 707, 632

* Lacour, S., Baudoz, P., Gendron, E., et al. 2014, SPIE Conf. Ser., 9147, 91479F

* Le Bouquin, J.-B., Lacour, S., Renard, S., et al. 2009, A&A, 496, L1

* Liljegren, S., Höfner, S., Eriksson, K., & Nowotny, W. 2017, A&A, 606, A6

* Lopez, B., Danchi, W. C., Bester, M., et al. 1997, ApJ, 488, 807

* Lopez, B., Lagarde, S., Jaffe, W., et al. 2014, The Messenger, 157, 5

* Monnier, J. D. 2001, PASP, 113, 639

* Monnier, J. D., Berger, J.-P., Millan-Gabet, R., et al. 2006, ApJ, 647, 444

* Norris, B. R. M., Tuthill, P. G., Ireland, M. J., et al. 2012, Nature, 484, 220

* Ohnaka, K., Weigelt, G., & Hofmann, K.-H. 2016, A&A, 589, A91

* Ohnaka, K., Weigelt, G., & Hofmann, K. H. 2017, A&A, 597, A20

* Pedretti, E., Millan-Gabet, R., Monnier, J. D., et al. 2004, PASP, 116, 377

* Pedretti, E., Traub, W. A., Monnier, J. D., et al. 2005, Appl. Opt., 44, 5173

* Perrin, G. 2003, A&A, 398, 385

* Perrin, G., Ridgway, S. T., Mennesson, B., et al. 2004, A&A, 426, 279

* Perrin, G., Cotton, W. D., Millan-Gabet, R., & Mennesson, B. 2015, A&A, 576, A70

* Poncelet, A., Doucet, C., Perrin, G., Sol, H., & Lagage, P. O. 2007, A&A, 472, 823

* Ragland, S., Traub, W. A., Berger, J.-P., et al. 2006, ApJ, 652, 650

* Reid, M. J., & Menten, K. M. 1997, ApJ, 476, 327

* Renard, S., Thiébaut, É., & Malbet, F. 2011, A&A, 533, A64

* Ryde, N., Gustafsson, B., Eriksson, K., & Hinkle, K. H. 2000, ApJ, 545, 945

* Scholz, M. 2003, Proc. SPIE, 4838, 163

* Thiébaut, E. 2008, Proc. SPIE, 7013, 70131I

* Thiébaut, É., & Giovannelli, J.-F. 2010, IEEE Signal Process. Mag., 27, 97

* Thiébaut, É., & Young, J. 2017, J. Opt. Soc. Am. A, 34, 904

* Traub, W. A. 1998, Proc. SPIE, 3350, 848

* Traub, W. A., Ahearn, A., Carleton, N. P., et al. 2003, SPIE Conf. Ser., 4838, 45

* Tsuji, T. 2000, ApJ, 538, 801

* Whitelock, P., Marang, F., & Feast, M. 2000, MNRAS, 319, 728

* Wittkowski, M., Boboltz, D. A., Ohnaka, K., Driebe, T., & Scholz, M. 2007, A&A, 470, 191

* Wittkowski, M., Boboltz, D. A., Ireland, M., et al. 2011, A&A, 532, L7

* Wittkowski, M., Chiovassa, A., Freytag, B., et al. 2016, A&A, 587, A12

* Wittkowski, M., Rau, G., Chiovassa, A., et al. 2018, A&A, 613, L7

* Wong, K. T., Kamiński, T., Menten, K. M., & Wyrowski, F. 2016, A&A, 590, A127

* Wood, P. R. 2007, IAU Symp., 239, 343

* Woodruff, H. C., Tuthill, P. G., Monnier, J. D., et al. 2008, ApJ, 673, 418
## Appendix A: Image reconstruction tests

<figure><img src="image_10.png" /><figcaption>**Fig. A.1.** Images reconstructed with MACIM and MiRA. The name of the algorithm is given first followed by the type of regularization and the initial image. The ellipse at the bottom right is the FWHM of the central peak of the dirty beam.</figcaption></figure>

We tested two different softwares with various priors and initial images to reconstruct the image of Mira: MACIM (Ireland et al. 2006) and MiRA (Thiébaut 2008). We define the reconstructed image $\mathbf{x} \in \Omega \subset \mathbb{R}_+$, with $\Omega$ the set of non-negative and normalized n-pixel images, as the one that minimizes an objective function $\mathcal{J}(\mathbf{x})$. In this case, the objective function is the sum of three terms:

$$ \mathcal{J}(\mathbf{x}) = \mathcal{J}_{V^2}(\mathbf{x}) + \mathcal{J}_{\text{cp}}(\mathbf{x}) + \mu \mathcal{J}_{\text{rgl}}(\mathbf{x}), \quad (A.1) $$

where $\mathcal{J}_{V^2}(\mathbf{x})$ and $\mathcal{J}_{\text{cp}}(\mathbf{x})$ are two data fidelity terms measuring the distance between the model and the data (respectively for the squared visibilities and for the closures phases) while $\mathcal{J}_{\text{rgl}}(\mathbf{x})$ is a regularization term introduced to impose additional a priori constraints (so-called “priors” for short). With MiRA, the hyper-parameter, $\mu > 0$, is used to adjust the impact of the priors on the result. For MACIM, the hyper-parameter is set to $\mu = 1$. The data fidelity terms are weighted according to the variance of the measurements. The $\chi^2_{\text{data}}$ term in Sect. 4.1 is $\mathcal{J}_{V^2}(\mathbf{x}) + \mathcal{J}_{\text{cp}}(\mathbf{x})$. Dividing by the number of degrees of freedom yields the reduced $\chi^2$. Starting from a given initial image $\mathbf{x}_{\text{init}}$, MiRA (Thiébaut & Giovannelli 2010; Thiébaut & Young 2017) iteratively reduces $\mathcal{J}(\mathbf{x})$ to find a minimum. MACIM is a Monte-Carlo Markov chain algorithm who searches for the global minimum of $\mathcal{J}(\mathbf{x})$. MACIM uses a simulated annealing algorithm with the Metropolis sampler. For the kind of considered data, the objective function is not convex and the minimum

of $\mathcal{J}(\mathbf{x})$ over $\Omega$ may not be unique. The result produced by MiRA and MACIM thus depend on the initial image, $\mathbf{x}_{\text{init}}$, and on the regularization settings, $\mu$ and $\mathcal{J}_{\text{rgl}}(\mathbf{x})$.

For the reconstructions with MiRA, we considered three different regularizations. The first one, which measures the squared distance of the image to the initial image is given by:

$$ \mathcal{J}_{\text{SD2I}}(\mathbf{x}) = \sum_{j=1}^{n} \|x_j - x_{\text{init},j}\|^2. \quad (A.2) $$

The second one, favors compact images:

$$ \mathcal{J}_{\text{Compact}}(\mathbf{x}) = \sum_{j=1}^{n} \left(1 + (2 \|\theta_j\| / \gamma)^2\right) x_j^2, \quad (A.3) $$

with $\theta_j$ the 2-dimensional angular position of pixel j relative to the center of the field of view, $\gamma > 0$ the full width at half maximum (FWHM) of the assumed prior image ($\theta_j$ and $\gamma$ being in the same units). Attempting to minimize $\mathcal{J}_{\text{Compact}}(\mathbf{x})$ for $\mathbf{x} \in \Omega$ imposes that the brightest parts of the sought image be mostly concentrated around the center of the field of view. Hence this regularization corresponds to a “compactness” prior (Thiébaut & Young 2017). The third one is given by:

$$ \mathcal{J}_{\text{Edge}}(\mathbf{x}) = \sum_{j=1}^{n} \sqrt{\|\nabla_j \mathbf{x}\|^2 + \tau^2}, \quad (A.4) $$
<table><caption>**Table A.1.** Algorithms, initial images and regularizations used to reconstruct the 6 images of Fig. A.1.</caption><thead><tr><th>#</th><th>Algorithm</th><th>x<sub>init</sub></th><th>J<sub>rgl</sub>(x)</th><th>χ<sup>2</sup><sub>r,1</sub></th><th>χ<sup>2</sup><sub>r,2</sub></th></tr></thead><tbody><tr><td>1</td><td>MiRA</td><td>Circularly symmetric star+layer model <sup>(1)</sup></td><td>J<sub>SD2I</sub>(x), μ = 10<sup>8</sup></td><td>8.214</td><td>1.088</td></tr><tr><td>2</td><td>MACIM</td><td>Circular limb-darkened disk <sup>(2)</sup></td><td>J<sub>Dark Energy</sub>(x), μ = 1</td><td>11.40</td><td>1.348</td></tr><tr><td>3</td><td>MiRA</td><td>Image 2 reconstructed by MACIM</td><td>J<sub>Compact</sub>(x), μ = 2 × 10<sup>6</sup>, γ = 18 mas</td><td>6.734</td><td>0.8473</td></tr><tr><td>4</td><td>MiRA</td><td>Point-like source</td><td>J<sub>Compact</sub>(x), μ = 2 × 10<sup>6</sup>, γ = 18 mas</td><td>6.734</td><td>0.8473</td></tr><tr><td>5</td><td>MiRA</td><td>Image 2 reconstructed by MACIM</td><td>J<sub>Edge</sub>(x), μ = 10<sup>4</sup>, τ = 10<sup>-4</sup></td><td>6.504</td><td>0.8362</td></tr><tr><td>6</td><td>MiRA</td><td>Image 4 reconstructed by MiRA</td><td>J<sub>Edge</sub>(x), μ = 10<sup>4</sup>, τ = 10<sup>-4</sup></td><td>6.503</td><td>0.8354</td></tr></tbody></table>

Notes. <sup>(1)</sup>Circular star+layer model of Sect. 5 fitted to the $V^2$ data. <sup>(2)</sup>Circular limb-darkened disk model fitted to the $V^2$ data.

where $\nabla_j x$ denotes the 2D spatial gradient of $x$ at pixel $j$ and $\tau > 0$ is a selected threshold. This kind of regularization penalizes the image gradients and thus imposes a smooth image. The penalization is however less severe for the image gradients larger than $\tau$ in magnitude than for the image gradients smaller than $\tau$ in magnitude which helps preserving the edges in the image. Hence this regularization imposes an edge-preserving smoothness prior.

For MACIM, the regularizer $J_{\text{Dark Energy}}(x)$ is a dark interaction energy regularizer. This regularizer is the sum of all pixel boundaries with zero flux on either side of the pixel boundary. It gives higher weight to large regions of dark space in between regions of flux.

In addition to the image presented in Fig. 4 and reconstructed with the SD2I regularization and an initial image obtained by fitting a circularly symmetric star+ shell model to the $V^2$ data, we reconstructed four other images with MiRA and one with MACIM. These images are presented in Fig. A.1 together with the image of Fig. 4. The details for each image are listed in Table A.1. Images of 100×100 pixels were reconstructed with MiRA while 40×40 pixels were used for MACIM. The pixel size is 1 mas for MiRA and 1.44 mas for MACIM. As explained in Sect. 4, different $\mu$ values were tried for the reconstruction of image 1 with MiRA and the best result was obtained for $\mu = 10^8$. For images 3 to 6, also reconstructed with MiRA, the value of $\mu$ was adjusted to minimize the $\chi^2$. The parameter $\tau$ is in units of intensity per pixel and therefore depends on the pixel size and on the level of the flux normalization (Renard et al. 2011). For the image reconstructed with MACIM, $\mu$ was set to 1.

The reduced $\chi^2$ are given in Table A.1. The number of degrees of freedom is the number of independent parameters subtracted to the number of data (940, i.e. 705 $V^2$ and 235 closure phases). The number of independent parameters is not as easy to define as in classical $\chi^2$ fitting. The maximum number of independent parameters is the number of pixels, e.g. 10 000 for the images reconstructed with MiRA. This number is larger than the number of data but not all pixels are independent because (1) they oversample the angular resolution elements in the image and (2) the regularization induces correlations between pixels. What matters here are the relative values to compare the images, the absolute values of the reduced $\chi^2$ are therefore only indicative. We have chosen to define the reduced $\chi^2$ as: $\chi_r^2 = \chi^2/\#\text{data} = \chi^2/940$. Our definition of $\chi_r^2$ is thus an underestimation of the reduced $\chi^2$. The errors provided by the pipeline have been modified for image 2 reconstructed with MACIM: a 5% relative error was applied to the $V^2$ data with an additive default floor of 0.0001 and a minimum error of 0.5° was used for the closure phase data. Raw error bars were used for the images reconstructed with MiRA (1 and 4 to 6). We have computed the

reduced $\chi^2$ in two different ways: for $\chi_{r,1}^2$ we compared the reconstructed image with the $V^2$ and closure phase data using the raw error bars while for $\chi_{r,2}^2$ the modified error bars were used. Values are given in Table A.1. The lower values of the $\chi_{r,2}^2$ estimator are a good indicator that the error bars produced by the pipeline are underestimated as noted in Sect. 5.1. The reduced $\chi_{r,2}^2$ are all close to 1 meaning all images are reconstructed with the same quality with respect to the data.

In addition we have reconstructed the dirty beam with MiRA (Fig. A.2). All the visibilities of the sampled $(u, v)$ points were set to 1 for this image. The white ellipse on the image gives the FWHM of the central peak and gives the scale of the resolution in the reconstructed images.

In spite of the aforementioned issue of the solution not being unique, the reconstructions in Fig. A.1 clearly show that the result does not strictly depend on the given algorithm that was used. We see no major influence of the initial image on the final image, in particular, there is no influence coming from a circular or absence of symmetry. The Compact regularization with MiRA makes the brightest peak appear a bit larger compared to the other reconstructions but does not change the image substantially. The features in the image cannot be explained by the structures of the dirty beam. Overall, the reconstructed images are very similar and exhibit the same remarkable structures with a ≈20 mas north-south elongated central object surrounded by a fainter ≈40 mas circular environment and a dark patch in the environment that is west of the bright central peak. We conclude that the image reconstruction process is quite robust despite the limited $(u, v)$ coverage. Image 1 featured in Sect. 4 is fairly representative of all the reconstructed images.

<figure><img src="image_11.png" /><figcaption>**Fig. A.2.** Dirty beam reconstructed with MiRA. The ellipse at the bottom right is the FWHM of the central peak. It gives the scale of the resolution in the reconstructed images.</figcaption></figure>