Resolving the extended stellar atmospheres of Asymptotic Giant Branch stars at (sub-)millimetre wavelengths

Mira A ($o$ Ceti) is the AGB star in the Mira AB system, and the archetypal Mira variable. Its companion is most likely a white-dwarf, but a K dwarf identification has also been suggested (e.g., Ireland et al., 2007). The two stars are physically separated by roughly 90 au. Mira B accretes material from the circumstellar envelope of Mira A, with observed interaction between the slow, dense outflow from Mira A and the faster, but thin outflow from Mira B (e.g. Ramstedt et al., 2014). The mass-loss rate of Mira A is estimated to be $\sim {10}^{-7}$ M${}_{\odot }$/year (Knapp et al., 1998; Ryde & Schöier, 2001), but the complex morphology make such determination rather uncertain. Despite its effect on the large-scale morphology of the circumstellar envelope, the companion is not expected to affect the close environment of the AGB star (Mohamed & Podsiadlowski, 2012). The distance to the Mira system derived from Hipparcos observations is $92\pm 10$ pc (van Leeuwen, 2007) and the distance from the period-luminosity relation is $110\pm 9$ pc (Haniff et al., 1995). Here we adopt $100$ pc. The visual magnitude of Mira A varies by roughly 8${}^{mag}$ over the pulsation cycle with a period of $\sim 332$ days (Watson et al., 2006).

Near-infrared observations between 1.0 and 4.0 $\mu$m (at $\phi \sim 0.7$ and $\phi \sim 0.0$) reveal uniform-disc diameters that vary as a function of wavelength from $\sim 30$ to $\sim 70$ mas (Woodruff et al., 2009). Monitoring the same wavelength range over eight pulsation cycles shows that the uniform-disc diameter varies in time at a given wavelength by at least 10 to 40% (Woodruff et al., 2008). A correction³³3In Vlemmings et al. (2015), a wrong interpretation of the uvmultifit results meant that the presented axis is the minor axis instead of the major axis. As a result, the position angle of the major axis should be ${90}^{\circ }$ rotated with respect to the angle presented in that paper. The derived brightness temperature $T_b$ is lower by (axis ratio)${}^2$. Hence, the results are fully consistent with Matthews et al. (2015). to the average uniform-ellipse sizes determined by Vlemmings et al. (2015) yield $(48.9\pm 0.7)\times (43.3\pm 0.1)$ mas at 229 GHz ($\phi \sim 0.5$) and of $(50.2\pm 0.6)\times (41.8\pm 0.4)$ mas at 94 GHz ($\phi \sim 0.45$). At 43 GHz, Reid & Menten (2007) measured an uniform-ellipse size of $54\pm 5\times 50\pm 5$ mas ($\phi =0.05$). We use $R_{\ast }=15$ mas, which is the smallest measurement of Woodruff et al. (2009) in the analysis. This corresponds to $\sim 1.5$ au.

The period of the semi-regular variable R Doradus appears to vary between two pulsation modes of 362 and 175 days on a time-scale of roughly 1000 days (Bedding et al., 1998). The non-unique period and the modest associated visual magnitude variation ($<2^{mag}$) makes it difficult to define the pulsation phase of R Dor at a given time. A distance of $\sim 60$ pc is obtained from observations using Hipparcos (59 pc, Knapp et al., 2003) and using aperture-masking in the near-infrared ($61\pm 7$ pc, Bedding et al., 1997). Models for observations of CO emission lines (Olofsson et al., 2002; Maercker et al., 2008; Van de Sande et al., 2018) imply gas mass-loss rates of $\sim {10}^{-7}$ M${}_{\odot }$ year${}^{-1}$. Based on oxygen isotopologue ratios, the initial stellar mass of R Dor is estimated to be between 1.3 and 1.6 $M_{\odot }$ (De Beck & Olofsson, 2018; Danilovich et al., 2017).

Observations with SPHERE/ZIMPOL revealed the size of the stellar disc diameter to vary strongly within the visible wavelength range, between $\sim 58$ mas and $\sim 72$ mas (Khouri et al., 2016a). The authors also reported a large variation in size (between $\sim 59$ mas and $\sim 72$ mas) within 48 days, which was accompanied by a striking change in morphology. Ireland et al. (2004) reported full-width-at-half-maximum (FWHM) sizes between 40 and 60 mas from fitting a Gaussian source to aperture-masked interferometry observations at visible wavelengths. At near-infrared wavelengths, a similar technique revealed a uniform-disc stellar radius of $\sim 27.5$ mas (Norris et al., 2012). Hence, we adopt $R_{\ast }=27.5$ mas, or $\sim 1.65$ au, in our analysis.

The distance to W Hydrae is found to be $\approx 100$ pc, from VLBI astrometry of maser spots (${98}_{-18}^{+30}$, Vlemmings et al., 2003) and Hipparcos observations (${104}_{-11}^{+14}$, van Leeuwen, 2007). Its mass-loss rate is estimated to be $\sim {10}^{-7}$ M${}_{\odot }$ yr${}^{-1}$, based on modelling of CO lines (e.g. Maercker et al., 2008; Khouri et al., 2014a), while its main-sequence mass has been constrained to be between 1.0 and 1.5 M${}_{\odot }$ (Khouri et al., 2014b; Danilovich et al., 2017), based on measured oxygen isotopic ratios. Over a period of $\sim 388$ days, the visual magnitude of W Hya varies between $\sim 6^{mag}$ and $9^{mag}$ (Zhao-Geisler et al., 2011). W Hya is sometimes classified as a Mira variable because of the relatively long period and large amplitude of its visible magnitude variations. However, a semi-regular classification is more appropriate based on studies of the gas-velocity variations in the atmosphere (Lebzelter et al., 2005).

The FWHM stellar size in the visible is found to reach values of $\sim 80$ mas at a pulsation phase of 0.44 (Ireland et al., 2004). Ohnaka et al. (2016) and Ohnaka et al. (2017) reported observations at maximum- and minimum-light phase, respectively, acquired in the near infrared with AMBER/VLTI and in the visible with SPHERE/ZIMPOL. Within one pulsation cycle, the authors measure a variation of 10% in the stellar size in the near infrared and a much larger variation of $\sim 50\%$ in visible wavelengths, which was accompanied by a strong change in morphology. As shown by Woodruff et al. (2008), its stellar radius also varies significantly in the near-infrared spectral range, by at least 10 to 50%. The smallest reported uniform-disc diameters of $\sim 30$ mas ($\phi \sim 0.79$) and $\sim 40$ mas ($\phi \sim 0.58$) were derived from observations in the near infrared by Woodruff et al. (2009), and the largest diameter (of $\sim 100$ mas) obtained in the mid-infrared (between $\phi \sim 0$ and $\sim 0.5$, Zhao-Geisler et al., 2011). ALMA images of the sub-mm continuum emission of W Hya ($\phi \sim 0.3$) revealed an asymmetric source with size $50.9\times 56.5$ mas and an unresolved hotspot on the surface (Vlemmings et al., 2017). At 43 GHz, Reid & Menten (2007) reported an uniform-ellipse size of $(69\pm 10)\times (46\pm 7)$ mas (at $\phi \sim 0.25$) and Matthews et al. (2018) found $(71\pm 4)\times (65\pm 3)$ mas ($\phi \sim 0.61$). At 22 GHz, RM97 found a uniform-disc diameter of $80\pm 15$ mas (average between measurements in two consecutive pulsation cycles at $\phi \sim 0.75$ and $\sim 0.8$). For consistency with (Vlemmings et al., 2017) we adopt $R_{\ast }=20$ mas $\approx 1.96$ au.

The Mira-type variable R Leonis experiences visible magnitude variations of 7${}^{mag}$ over a period of 310 days (Watson et al., 2006). The distance to R Leo is not well determined, with published values ranging from $\sim 70$ to $110$ pc. The distance from the period-luminosity relation is $110\pm 9$ pc (Haniff et al., 1995). The recent geometric Gaia distance is ${71}_{-4}^{+5}$ pc, but the accuracy of this value is unclear since Gaia AGB distances are highly uncertain (see e.g. van Langevelde et al., 2019). Specifically for R Leo, the Gaia result has a goodness-of-fit value of $142$, a ${\chi }^2=53018$ and astrometric excess noise of $2.7$ mas. This implies a low-quality fit and the distance and especially the errors should thus be treated with extreme caution. Still, we use $71$ pc in our analysis. The mass-loss rate of R Leo has been determined to be $\sim {10}^{-7}$ M${}_{\odot }$/year (Danilovich et al., 2015; Ramstedt & Olofsson, 2014).

At visible wavelengths and $\phi =0.2$, the diameter was found to vary between $48.7\pm 2.3$ and $75.6\pm 3.7$ mas depending on the strength of TiO bands at the given wavelength (Hofmann et al., 2001). Norris et al. (2012) reported a uniform-disc diameter of $36.6\pm 0.6$ mas at 1.04 $\mu$m ($\phi =0.4$). Observations acquired over nine pulsation cycles showed that the uniform-disc diameter varies in time by at least 10 to 30% at a given wavelength in the near infrared (Woodruff et al., 2008). At $\phi =0.5$, the uniform-disc diameter between 1.0 and 4.0 $\mu$m was observed to vary between $\sim 27$ mas and $\sim 55$ mas (Woodruff et al., 2009). Observations by Tatebe et al. (2008) at 11.15 $\mu$m indicated a uniform-ellipse size $(31.81\pm 0.15)\times (30.67\pm 0.22)$ ($\phi =0.6$ to $0.8$). At 43 GHz, Reid & Menten (2007) measured the size of R Leo to be $(61\pm 10)\times (39\pm 6)$ ($\phi \sim 0.55$), while Matthews et al. (2018) reported $(54\pm 3)\times (45\pm 2)$ mas ($\phi \sim 0.23$). We adopt the smallest measurement of Woodruff et al. (2009), $R_{\ast }=13.5$ mas, corresponding to $\sim 0.96$ au, as the stellar radius.

In order to derive rough estimates for the physical conditions in the extended atmospheres probed by our ALMA observations, we construct a simple parameterized model atmosphere based on the radio-photosphere model of RM97. Since the density profile in the AGB atmosphere is significantly affected by pulsations and convective motions, we adopt a basic power-law description for the molecular hydrogen number density $\left(n_{H_2}\right)$:

Here $n_{H_2,0}$ is the number density (in $c{m}^{-3}$) at the stellar photosphere, $R_{\ast }$ is the photospheric radius, and $\alpha$ is the power law-index. For comparison, an atmosphere in hydrostatic equilibrium can, in the region of interest, be described with a power-law exponent $\alpha \approx 9$. In a dynamical atmosphere, the exact density vs. radius relation becomes complex and depends strongly on pulsation cycle with regions displaying slopes steeper than in the hydrostatic case, while others are significantly more shallow (e.g. Höfner et al., 2003). When convective motions are included, the average power-law exponent can also reach $\alpha \approx 6$ (Freytag et al., 2017). Density profiles with $\alpha <6$ are also typically used to describe the molecular atmosphere (e.g. Wong et al., 2016; Khouri et al., 2016b). We adopt $\alpha =6$ for our reference model.

For the temperature profile we either use the gray atmosphere profile also adopted in RM97, with:

or a power-law given by:

Here $T_{R_{\ast }}$ is the temperature (in K) at the photosphere and $\beta$ the power-law index. Finally, we approximate the ionisation profile also with a power-law in temperature:

Here $X_i=n_e/n_{H_2}$, where $n_e$ is the electron number density. We fix the ionisation at $T=3000$ K to be $X_{i,3000}=7\times {10}^{-6}$, which would correspond to single ionisation of K, Na, Ca and Al when assuming solar abundance. For the specific stellar models presented later, the power-law index $\gamma =8$, which fits the range provided by the detailed model calculations from RM97. The exact adopted value of $R_{\ast }$ for our sources (as discussed in § 2) does not affect any of our conclusions, and for different stellar radii the surface values of density and temperature can be scaled accordingly.

With these parameters, we can calculate the frequency dependent optical depth throughout the extended atmosphere by taking into account the relevant opacity sources. As found in RM97, of the possible processes that contribute to the opacity, the electron-neutral free-free emission dominates in the relevant temperature and density regime for the radio and (sub-)millimetre wavelengths. Hence we only considered this process. We did confirm that the electron-ion free-free emission contributes at least an order of magnitude less opacity even in the warmest, most ionized and densest regions of our standard model. We only consider interactions of the electrons with purely atomic hydrogen in our calculations. This is in contrast to the study in RM97 where an equilibrium state was assumed between atomic and molecular hydrogen. As was already noted in Bowen (1988) however, the formation rate of molecular hydrogen, in the shocked atmosphere of an AGB star, can be slow. RM97 argue that for high densities ($>{10}^{12}$ cm${}^{-3}$) in the extended atmosphere, the equilibrium time is $<250$ days. Since this is less than a typical pulsation period, RM97 assumed equilibrium conditions. However, for most of our sources, we find densities that are $<{10}^{12}$ cm${}^{-3}$ throughout a large part of the atmosphere. Also, recent models and observations (e.g. Khouri et al., 2016a; Freytag et al., 2017) have shown that the time-scales for significant changes due to convective motions and/or shocks can be measured in days or weeks. Hence we assume that equilibrium conditions are not reached. In case there is a significant molecular hydrogen component after all, the effective opacity decreases by up to a factor of five and our derived densities (or ionisation) would need to be increased. To calculate the opacity, we use the RM97 third-order polynomial fits to the absorption coefficients calculated by Dalgarno & Lane (1966). We note that another calculation and approximation is given in Harper et al. (2001), which is consistent within $\sim 10\%$ across the temperature range that we probe in the AGB atmosphere.

With the calculated opacities, we determine brightness temperature profiles using standard radiative transfer in the Rayleigh-Jeans limit. From these profiles we derive the frequency dependent stellar disc brightness temperature and size that can be directly compared with the observations. We assume a uniform stellar disc with a radius $R\left(\nu \right)$, which is the radius where the brightness temperature has dropped by $50\%$ in units of stellar radii $R_{\ast }$. The angular size, $\theta$, can be directly related to $R\left(\nu \right)$ by dividing it by the adopted angular stellar radius from Table 2. We focus in particular on the frequency-size and size-temperature relations as presented in Figs. 1, 2, 3, and 4 for the observations. As the opacity is proportional to the density (with exponent $2$) and ionisation (linearly) [in addition to the inverse quadratic dependence on frequency], it is clear that there are degeneracies between the various parameters used to describe the extended atmosphere. Using the above exponents and the power-law components from Eqs $1-5$, we see that the optical depth ${\tau }_{ff}\propto r^{-(\gamma \beta +2\alpha )}{\nu }^{-2}$. Since the temperature is best constrained by the observations at different wavelengths, the ionisation and number density are most degenerate. Although the power-law approximation is likely far from reality, with shocks and pulsations producing density and temperature profiles that cannot be described using simple formulae, a number of conclusions can still be drawn from a small parameter study. A coupling of a dynamic and/or convective atmosphere code with the full continuum radiative transfer is beyond the scope of this paper.

In Fig. 9, we show the size vs. frequency and brightness temperature vs. size relations for models with different $\alpha ,\beta ,$ and $\gamma$. The fixed model parameters for the plotted models are given in Table 4. In the figures we also indicate the observed range of slopes $\left(\eta \right)$ of the size vs. frequency relation derived from the observations. From this is it clear that the lowest of the derived values of $\eta$, observed for Mira A, R Dor and R Leo, cannot be reproduced for any but the steepest density profiles ($\alpha >9$), and for none of the temperature or ionisation profiles. One of the most straightforward ways to match the shallow size vs. frequency relations observed is to include a transition region with a sudden increase in opacity. Such a transition could naturally be related to a shock traveling through the extended atmosphere and mimics the density behaviour in dynamical atmopshere models (e.g. Bowen, 1988; Höfner et al., 2003), with density variations up to a factor of $\sim 10$. We thus produce models, which we describe as a ’layer’ model, with a $10$ times increase in ${\tau }_{ff}$ at a single ’shock’ radius ($R_s$) to investigate if this could reproduce the observations. Including such an opacity increase specifically changes the stellar profile by introducing a region that is sufficiently optically thick so that there is no noticeable change of size with increasing frequency (within the observed frequency range). It can also produce limb brightening effects or ring-like structures in the residual maps once the layer has moved further out. As a result, the size determined using a uniform stellar disc will increase, thereby producing a sharp drop in the temperature vs. size relation and a leveling off of the size vs. frequency relation. Ring-like structures related to this effect have possibly been seen around Betelguese (O’Gorman et al., 2017) and IRC+10216 (Matthews et al., 2018).

Here we discuss the modeling results for the four stars separately. We adopted the standard ionisation description for all stars. We note that we treat the different epochs presented in this paper for all stars except R Leo as sufficiently simultaneous in this analysis. The results are shown in Fig. 10.

ALMA observed R Leo within two weeks on four and two occasions in Band 4 and Band 6, respectively, at fairly similar frequencies in each given band. We are thus able to investigate its time variable behavior in detail. Looking at the flux density variations, we find differences of on average $5\%$. These can mostly be attributed to absolute flux calibration uncertainties, even though the same flux calibrator was used. However, we also note a clear trend of increasing radius and ellipticity with time. In Fig. 11 we show a comparison of the uv-data and model fits for the first and last observing epoch in Band 4. The fit indicates that the size change is robustly determined. In Fig. 12 we display the measured size against the time after the first observations for both Band 4 and Band 6. In this plot we adjust all measurements to 140 and 225 GHz for the Band 4 and Band 6 data, respectively, using the individually fitted slopes to the frequency-size relations. We can fit a linear expansion of $0.17\pm 0.02$ mas day${}^{-1}$ to the Band 4 data and find that this is consistent with the size measurements during the two epochs at 225 GHz. These measurements imply that the continuum surface seen in the millimetre wavelength observations is expanding at $10.6\pm 1.4$ km s${}^{-1}$ when assuming a distance to R Leo of 71 pc. The fact that the 140 and 225 GHz surfaces undergo a similar expansion implies a bulk motion of material that encompasses almost $0.2R_{\ast }$. A motion of the continuum surface can occur because of changes in density, ionisation and/or temperature. The temperature in the region that we are probing does not seem to change by more than $\pm 5\%$ (most likely due to absolute flux calibration uncertainties). Such a variation would produce a change in ionisation of $\sim 50\%$ and a similar change in opacity. According to our models for R Leo, this would result in a change of measured radius, at 140 GHz, of $\sim 0.05R_{\ast }$. This assumes the temperature changes throughout the entire envelope which requires radiative heating, which is unlikely. The radius change would be significantly less if the temperature change is only local. Hence, the motion that we are seeing is likely due to density or ionisation changes (by factors $>2$ or $>4$, respectively).

The measured velocity matches the velocity dispersion seen in the gravitationally bound molecular gas around a number of stars (e.g. Hinkle et al., 1997; Wong et al., 2016; Vlemmings et al., 2017, 2018). Since the star also increases its ellipticity, the motion is not spherically symmetric. The expansion velocity taking the major axis values alone would correspond to $12\pm 1$ km s${}^{-1}$ and that of the minor axis to $9\pm 1$ km s${}^{-1}$. Note that if the movement occurs mostly in only one direction and not equally in opposite directions (which is unlikely if it is a result of convective motions), the inferred velocities are up to a factor of two higher. In any case, the velocity is less than the escape velocity, which is $\sim 34$ km s${}^{-1}$ at $1.6R_{\ast }$ for a star with $M=1$ M${}_{\odot }$. The velocity exceeds the local sound speed in the atmosphere, which is $\sim 5$ km s${}^{-1}$. It is also larger than the velocities of $\lesssim 8$ km s${}^{-1}$ found at $R\sim 1.6R_{\ast }$ in current convection models (Freytag et al., 2017).

Previous observations of AGB stars at radio and sub-mm/mm wavelengths have revealed a variety of asymmetries and hotspots in the extended atmospheres. Radio observations at different epochs have shown significant ellipticity, brightness asymmetries and non-uniform surfaces (e.g., RM97 Reid & Menten, 2007; Matthews et al., 2018). ALMA observations of Mira A and W Hya also revealed a mottled surface structure at sub-mm/mm wavelengths (Vlemmings et al., 2015; Matthews et al., 2015; Vlemmings et al., 2017). The most prominent of these was the detection of a hotspot with $T_b>50.000$ K on the surface of W Hya. As can be seen in the residuals from a uniform stellar disc model and the uv-model fits for the sources presented here (Figs 5, 6, 7, and 8), all stars display asymmetries and/or hotspots to various degree. In some cases, negative residuals are seen that are similar in strength to the positive features. This is a natural result when determining the residuals after subtracting a best-fit uniform-disc model, which minimises the intergrated flux across the disc.

Mira A: Our observations of Mira A indicate that, at the current epoch, the star is relatively circular with a marginally significant ellipticity of a few percent seen in the Band 6 data. This is in contrast with previous ALMA observations that indicate an ellipticity between $10-20\%$ (Vlemmings et al., 2015; Matthews et al., 2015). While the previous observations had significant residuals that could indicate a strong hotspot, the current data only have a weak spotted structure, predominantly in the limb of the disc. This could be related to a (non-uniform) enhanced opacity layer, such as included in the atmosphere models, with changes in density or ionisation by a factor of $\sim 3$ or $\sim 10$, respectively.

R Dor: R Dor shows strong residual hotspots, especially in the Band 7 observations. At those frequencies we can determine a lower limit to $T_b$ for both brightest hotspot of $>3300$ K above the average stellar disc temperature. The actual temperature could be significantly higher if the hotspot sizes are smaller than the formal limit of $<6\times 3$ mas obtained from a fit to the residual image and uv-data. The brightness variation across the stellar surface appears to correlate with structure seen at visible wavelengths (Khouri, 2018, and in prep.). In a comparison between the Band 6 and 7 observations, there also appears to be a correspondence between the structure seen at both wavelengths. The two frequencies probe a slightly different region of $\sim 0.03$ R${}_{\ast }$ in depth. If the structures are indeed the same, they would thus have a vertical extent of at least $8\times {10}^{11}$ cm ($\sim 0.05$ au). Since R Dor appears to be the star with the smallest (sub-)millimetre size compared to its stellar radius as determined from infra-red observations we are potentially probing very close to a chromosphere if it is present. Notably, the ellipticity of R Dor increases towards the shorter wavelengths in Band 7, with an ellipticity up to $\sim 6\pm 1\%$ at a position angle of $44\pm 2^{\circ }$. There is no correspondence between this position angle and the apparent rotation axis (Vlemmings et al., 2018).

W Hya: As was found previously (Vlemmings et al., 2017), W Hya shows strong deviations from a uniform disc. Hotspots as well as negative patches are seen in both the Bands 4 and 7 observations. As discussed above, in the event of strong hotspots, negative bowls also naturally occur when a uniform disc with spots is fitted by a disc alone. Negative patches towards the star can also occur due to an increase of opacity along the line of sight. In the limb they can be the result of a convolution of the beam with the residuals that remain when a top-hat stellar disc is subtracted from a more shallow stellar profile. In the Band 4 observations, we find that the prominent hotspot can be constrained to $<5.3\times 5.3$ mas, which implies $T_b>3750$ K. In Band 7, the two spots are very compact ($<2.6\times 2.6$ mas) with $T_b>10.000$ K. This is consistent with the previously observed hotspot. We find no correspondence in the positions of the spots seen in Bands 4 and 7, nor with the position of the spot seen in previous observations. Since the region probed at the shorter wavelengths (at $\sim 1.2R_{\ast }$) is still optically thick at longer wavelengths (where we find a radius of $\sim 1.33R_{\ast }$), the Band 7 spots are likely obscured at Band 4. Additionally, the current observations were taken almost a month apart, while the previous observations were taken two years earlier. Since significant changes in the optical are seen on time-scales of weeks (Khouri et al., in prep.), and convective motions and associated shocks operate on similar time-scales. it is likely that the lifetime of the hotspots is of the same order.

R Leo: Although the multi-epoch observations of R Leo show a clear signature of an expanding layer, there is no sign of strong hotspots. Only in the Band 6 observations do we see an extension in the north-east. This structure appears marginally resolved ($\sim 10\times 5$ mas implying $T_b\sim 1000$ K) along the limb of the star, which means it could be related to the layer invoked to explain the observed sizes and temperatures. The ellipticity also increases from nearly circular to $\sim 7\%$ ellipticity with a position angle due north. The fact that the residuals remain towards the north-east implies that the true shape of the stellar disc is not quite elliptical and that an expansion of material with a size of up to an entire quarter of the stellar disc is occurring. Such a large size would not be inconsistent with that of a large convective cell.

As noted previously, in particular the shortest wavelength Mira data show a clear deviation from a simple stellar disc model at the largest angular scales. We show this residual, for one of the spws, in Fig. 13. The same emission is seen in all spws, with an almost linear increase of integrated flux density from $12$ mJy at $330$ GHz to $31$ mJy at $344$ GHz. It is immediately apparent that the extended emission peaks towards the east of the star and extends over almost $5R_{\ast }$. The extent and morphology bears resemblance to the emission seen in SO, CO and TiO${}_2$ lines shown in Khouri et al. (2018) and extends out to the dust ridges seen with SPHERE/ZIMPOL in the same paper.

We suggest three possible origins of the extended emission. If the emission originates from circumstellar dust emission close to the star it would represent a dust mass of $\sim 2\times {10}^{-7}$ M${}_{\odot }$ at $T=1000$ K adopting circumstellar dust properties from Knapp et al. (1993). The large flux density difference between the spws would argue against the emission being due to dust. However, since the emission only corresponds to a few percent of the subtracted stellar emission, the flux density measurements might be quite uncertain. If the emission is due to dust, it is however unclear why the VLT/SPHERE observations (Khouri et al., 2018) show such a clear lack of dust scattering emission in the region.

A second possibility would be that the emission is due to many weak residual spectral lines remaining in the data. Although a significant effort was made to remove the obvious lines, the line density at the shortest wavelengths around Mira A is significant. If we are seeing stacked lines this could explain the correspondence with the aforementioned stronger lines arising in the same region. It would hint at an asymmetric mass loss process that is not surprising considering the complexity of the Mira system. It is not clear why the flux density of the emission changes almost linearly with frequency.

Finally, the emission could represent the electron-neutral free-free emission from the optically thin extended atmosphere. Although we find that the brightness profile is generally well represented by a stellar disc model, the true profile could have weaker emission extending out to several stellar radii. The Mira 46 GHz observations indicate an extent of almost $4R_{\ast }$. Furthermore, although we were unable to produce a similar high fidelity map at Band 6, also at those frequencies a short-baseline excess was found. This hypothesis does not naturally explain the strong flux density increase with frequency, but potentially rapid changes of opacity could be involved. If this is indeed the cause for the emission, it would also indicate quite a non-symmetric envelope.

At the moment we cannot rule out any of the possible origins with certainty, although dust emission appears unlikely. To distinguish between these possibilities, sensitive observations at high frequencies and improved models will be needed.