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

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

W. H. T. Vlemmings Department of Space, Earth and Environment, Chalmers University of Technology, Onsala Space Observatory, 439 92 Onsala, Sweden wouter.vlemmings@chalmers.se    T. Khouri Department of Space, Earth and Environment, Chalmers University of Technology, Onsala Space Observatory, 439 92 Onsala, Sweden    H. Olofsson Department of Space, Earth and Environment, Chalmers University of Technology, Onsala Space Observatory, 439 92 Onsala, Sweden
21-Feb-2019
Key Words.:
evolved stars, stars: AGB and post-AGB, stars: atmospheres, stars: activity, stars: imaging, Submillimetre: stars
Abstract

Context:The initial conditions for the mass loss during the Asymptotic Giant Branch (AGB) phase are set in their extended atmospheres, where, among others, convection and pulsation driven shocks determine the physical conditions.

Aims:High resolution observations of AGB stars at (sub-)millimetre wavelengths can now directly determine the morphology, activity, density, and temperature close to the stellar photosphere.

Methods:We use the Atacama Large Millimeter/submillimeter Array (ALMA) high angular resolution observations to resolve the extended atmospheres of four of the nearest AGB stars: W Hya, Mira A, R Dor and R Leo. We interpret the observations using a parameterized atmosphere model.

Results:We resolve all four AGB stars and determine the brightness temperature structure between and stellar radii. For W Hya and R Dor we confirm the existence of hotspots with brightness temperatures to  K. All four stars show deviations from spherical symmetry. We find variations on a time scale of days to weeks, and for R Leo we directly measure an outward motion of the millimetre wavelength surface with a velocity of at least  km s. For all objects but W Hya we find that the temperature-radius and size-frequency relations require the existence of a (likely inhomogeneous) layer of enhanced opacity.

Conclusions:The ALMA observations provide a unique probe of the structure of the extended AGB atmosphere. We find highly variable structures of hotspots and likely convective cells. In the future, these observations can be directly compared to multi-dimensional chromosphere and atmosphere models that determine the temperature, density, velocity and ionisation structure between the stellar photosphere and the dust formation region. However, our results show that for the best interpretation, both very accurate flux calibration and near-simultaneous observations are essential.

1 Introduction

Asymptotic Giant Branch (AGB) stars, with main-sequence masses between  M, lose significant amounts of mass in a stellar wind and are an important source of interstellar enrichment (Höfner & Olofsson, 2018). It is generally assumed that radiation pressure on dust, which formes at a few stellar radii, drives the AGB winds (e.g. Höfner, 2008). How the gas is transported from the stellar surface, through an extended stellar atmopshere, to the cooler regions where dust can form is less clear. Hence, our current understanding of the mass-loss processes around AGB stars depends critically on the description of convection, pulsations and shocks in the extended stellar atmosphere (e.g. Bowen, 1988; Woitke, 2006; Freytag & Höfner, 2008). Additionally, the chemistry before and during the dust formation is also strongly affected by shocks (e.g. Cherchneff, 2006).

Radio and (sub-)millimetre observations are uniquely suited to determining the properties of the extended AGB atmosphere. Continuum observations at centimeter wavelengths have shown emission in excess of stellar blackbody radiation. In Reid & Menten (1997, hereafter RM97), this emission was partly resolved and shown to arise from an extended region around the star where the opacity is dominated by electron-neutral free-free emission processes. Later radio continuum observations further resolved this region and showed that it encompasses most of the extended AGB atmosphere out to approximately two stellar radii (Reid & Menten, 2007; Matthews et al., 2015, 2018). These observations are inconsistent with the extended chromospheres that were produced in early AGB atmosphere models by the outward propagation of strong shocks (e.g. Bowen, 1988). Still, ultraviolet line and continuum emissions that appear to originate from an AGB chromosphere have been observed around the majority of AGB stars (e.g. Johnson & Luttermoser, 1987; Montez et al., 2017). The UV observations and radio observations can be reconciled if the extended atmospheres of these cool evolved stars are inhomogeneous, with the hot chromospheric plasma lying embedded in a cooler extended atmosphere and having a relatively low filling factor. This has previously been suggested for the supergiant Betelgeuse (Lim et al., 1998). Recent ALMA observations that resolved the submillimetre surface of the AGB star W Hya also support this hypothesis (Vlemmings et al., 2017).

Here we present further high-angular resolution ALMA observations at millimetre and sub-millimetre wavelengths to specifically probe the structure and physical properties of the inner extended atmosphere. These observations will produce direct benchmarks for the three-dimensional atmosphere models that currently serve as input for mass-loss prescription models (e.g. Freytag et al., 2017). We present the stars in our sample in § 2 and the observations and data analysis in § 3. The results of a uv-analysis are given in § 4. Subsequently, we present parameterized atmosphere models used to describe our observations in § 5.1, discuss the observed expansion of one of our sources in § 5.2, describe the observed hotspots and asymmetries in § 5.3, and analyse larger-scale excess emission in § 5.4. We end with a number of conclusions in § 6.

Obs. date ALMA project Band Ref. freq. Beam size, P.A. Rms Flux Cal. Gain cal.
[GHz] [mas], [] [Jy
beam
Mira A
2017-09-21 2016.1.00004.S 4 , J0006-0624;  , J0228-0337
2017-09-22 2016.1.00004.S 6 , J0006-0624;  , J0228-0337
2017-11-09 2017.1.00191.S 7 , J0208-0047;  , J0217+0144
R Dor
2017-09-26 2016.1.00004.S 4 , J0519-4546;  , J0428-6438
2017-09-25 2016.1.00004.S 6 , J0519-4546;  , J0428-6438
2017-11-09 2017.1.00191.S 7 , J0522-3627;  , J0516-6207
W Hya
2017-10-28 2017.1.00075.S 4 , J1427-4206;  , J1351-2912
2017-11-25 2016.A.00029.S 7 , J1337-1257;  , J1351-2912
R Leo
2017-09-30 2016.1.00004.S 4 , J0854+2006;  , J1002+1216
2017-10-03 2017.1.00862.S 6 , J0854+2006;  , J1002+1216
2017-10-05 2017.1.00862.S 4 , J0854+2006;  , J1002+1216
2017-10-09 2017.1.00862.S 4 , J0854+2006;  , J1002+1216
2017-10-11 2017.1.00075.S 6 , J0854+2006;  , J1002+1216
2017-10-13 2017.1.00075.S 4 , J0854+2006;  , J1002+1216
111 $a$$a$footnotetext: Briggs weighted image of a single spw at the reference frequency. $b$$b$footnotetext: Calibrator name, flux densities (Jy at the reference frequency) and spectral index.
Table 1: Observational details

2 The Sources

Star Distance Radius Period
[pc] [mas] [days]
Mira A 100 15 332
R Dor 60 27.5 362/175
W Hya 98 20 388
R Leo 71 13.5 310
222 $a$$a$footnotetext: References listed in the text
Table 2: Adopted parameters of the sample

We observed four nearby AGB stars with ALMA at high angular resolution in order to describe their extended atmosphere at (sub-)millimetre wavelengths. Here we briefly present the four stars and in particular provide the distances used in our analysis. The observed stars pulsate with periods of the order of one year and, over each pulsation cycle, experience variations of orders of magnitude in visual light. The visual-light variation is caused by the strongly time-variable molecular opacity in this wavelength range (Reid & Goldston, 2002). The opacity (and, consequently, the stellar size) vary strongly with wavelength, as does the dominating opacity source (from atoms to molecules to dust grains to free electrons). Moreover, the stellar size at a given wavelength is also observed to vary significantly in time. This is a consequence of changes in density and temperature in the inner regions that are caused by stellar pulsations and motion of large convective cells.

2.1 Mira A

Mira A ( 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  M/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  pc (van Leeuwen, 2007) and the distance from the period-luminosity relation is  pc (Haniff et al., 1995). Here we adopt  pc. The visual magnitude of Mira A varies by roughly 8 over the pulsation cycle with a period of  days (Watson et al., 2006).

Near-infrared observations between 1.0 and 4.0 m (at and ) reveal uniform-disc diameters that vary as a function of wavelength from to  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 correction333In 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 rotated with respect to the angle presented in that paper. The derived brightness temperature is lower by (axis ratio). Hence, the results are fully consistent with Matthews et al. (2015). to the average uniform-ellipse sizes determined by Vlemmings et al. (2015) yield  mas at 229 GHz () and of  mas at 94 GHz (). At 43 GHz, Reid & Menten (2007) measured an uniform-ellipse size of  mas (). We use  mas, which is the smallest measurement of Woodruff et al. (2009) in the analysis. This corresponds to  au.

2.2 R Dor

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 () makes it difficult to define the pulsation phase of R Dor at a given time. A distance of  pc is obtained from observations using Hipparcos (59 pc, Knapp et al., 2003) and using aperture-masking in the near-infrared ( 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  M year. Based on oxygen isotopologue ratios, the initial stellar mass of R Dor is estimated to be between 1.3 and 1.6  (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  mas and  mas (Khouri et al., 2016a). The authors also reported a large variation in size (between  mas and  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  mas (Norris et al., 2012). Hence, we adopt  mas, or  au, in our analysis.

2.3 W Hya

The distance to W Hydrae is found to be  pc, from VLBI astrometry of maser spots (, Vlemmings et al., 2003) and Hipparcos observations (, van Leeuwen, 2007). Its mass-loss rate is estimated to be  M yr, 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 (Khouri et al., 2014b; Danilovich et al., 2017), based on measured oxygen isotopic ratios. Over a period of  days, the visual magnitude of W Hya varies between and (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  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 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  mas () and  mas () were derived from observations in the near infrared by Woodruff et al. (2009), and the largest diameter (of  mas) obtained in the mid-infrared (between and , Zhao-Geisler et al., 2011). ALMA images of the sub-mm continuum emission of W Hya () revealed an asymmetric source with size  mas and an unresolved hotspot on the surface (Vlemmings et al., 2017). At 43 GHz, Reid & Menten (2007) reported an uniform-ellipse size of  mas (at ) and Matthews et al. (2018) found  mas (). At 22 GHz, RM97 found a uniform-disc diameter of  mas (average between measurements in two consecutive pulsation cycles at and ). For consistency with (Vlemmings et al., 2017) we adopt  mas  au.

2.4 R Leo

The Mira-type variable R Leonis experiences visible magnitude variations of 7 over a period of 310 days (Watson et al., 2006). The distance to R Leo is not well determined, with published values ranging from to  pc. The distance from the period-luminosity relation is  pc (Haniff et al., 1995). The recent geometric Gaia distance is  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 , a and astrometric excess noise of  mas. This implies a low-quality fit and the distance and especially the errors should thus be treated with extreme caution. Still, we use  pc in our analysis. The mass-loss rate of R Leo has been determined to be  M/year (Danilovich et al., 2015; Ramstedt & Olofsson, 2014).

At visible wavelengths and , the diameter was found to vary between and  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  mas at 1.04 m (). 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 , the uniform-disc diameter between 1.0 and 4.0 m was observed to vary between  mas and  mas (Woodruff et al., 2009). Observations by Tatebe et al. (2008) at 11.15 m indicated a uniform-ellipse size ( to ). At 43 GHz, Reid & Menten (2007) measured the size of R Leo to be (), while Matthews et al. (2018) reported  mas (). We adopt the smallest measurement of Woodruff et al. (2009),  mas, corresponding to  au, as the stellar radius.

Frequency Obs. date Axis ratio P.A. Spectral
[mJy] [mas] [] [K] index
Mira A
2017-09-21
2017-09-22
2017-11-09
R Dor
2017-09-26
2017-09-25
2017-11-09
W Hya
2017-10-28
2017-11-25
R Leo
2017-09-30
2017-10-03
2017-10-05
2017-10-09
2017-10-11
2017-10-13
444 $a$$a$footnotetext: Spectral index fit combining all epochs/bands. $b$$b$footnotetext: Fit restricted to baseline lengths  M (Band 4) and  M (Bands 6 and 7). $c$$c$footnotetext: Combined spectral index fit with previous epoch.
Table 3: Flux density, size, temperature and spectral index results.

3 Observations and data reduction

The observations presented in this paper were taken between 21 September and 25 November 2017 in some of the longest baseline ALMA configurations with maximum baselines between 10 and 16.2 km. The data were obtained as part of several different projects. Innitially we analysed data from the projects; 2016.1.00004.S, 2016.A.00029.S, 2017.1.00075.S (PI: Vlemmings), and 2017.1.00191.S (PI: Khouri). When we noted an apparent increase in stellar size for R Leo in our three observing epochs (see § 5.2), we retrieved archive data of R Leo taken between our observations as part of project 2017.1.00862.S (PI: Kaminski). All observations were taken in spectral line mode, with four spectral windows (spws) with a bandwidth of  GHz each. The initial calibration (of the bandpass, amplitude, absolute flux-density scale and phase) for several of the datasets was done using the ’Common Astronomy Software Applications’ (CASA) package pipeline (McMullin et al., 2007), while others were manually calibrated by observatory staff using different version of CASA. Observation details, including the flux and gain calibrators used, are given in Table. 1.

We inspected the results of the initial calibrations, in particular the absolute flux calibration and flagging. We identified a few cases where extra flagging was required. Additionally, because of the multi-epoch coverage of R Leo in Band 4, we were able to determine systematic flux density differences between calibrators for the flux and gain that appeared to affect the measurements of the star. Hence, the flux densities in these observations were re-calibrated using a renormalized flux density for the gain calibrator (J1002+1216) of  mJy at  GHz. This reduced the flux densities by in the epochs taken on 30 September and 9 October. Based on the remaining scatter in the R Leo observations, we estimate a remaining absolute flux calibration uncertainty in Band 4 of for R Leo. For the other sources only a single epoch was available at each observing band and we assume a uncertainty for all Bands.

After the initial calibration, we performed self-calibration. For the data sets taken as part of projects 2016.1.00004.S and 2017.1.00075.S, one of the spws contained a strong SiO maser that was used for the first two steps of phase self-calibration. These solutions were then applied to all other spws. We subsequently carefully removed obvious spectral lines and performed spectral (20 channels) and time smoothing (10 sec) to reduce the data size. We then performed one or two, depending on signal-to-noise (SNR) improvements, further self-calibration steps for phase solutions on the continuum. These solutions were also applied to all the data. No amplitude self-calibration was performed and for the data-sets without SiO masers, only continuum self-calibration was done. The self-calibration generally improved the SNR from to or better. In Table 1 we give the rms noise and beam size for the lowest frequency spw in the data using Briggs robust weighting (Briggs, 1995). The maximum recoverable scale depends on the shortest baselines in the antenna configuration and ranges from at the longest wavelengths to at the shortest wavelengths. Results for a selection of molecular lines that were included in the observations have previously been presented in Vlemmings et al. (2018) and Khouri et al. (2019).

4 Results

Figure 1: Top left: Flux density measurements determined by uv-fitting a uniform elliptical stellar disc to Mira A. Only the fitted flux density uncertainties are included. The absolute flux density uncertainty between the different observing bands are discussed in the text. The symbol size is generally larger than the fitted flux density uncertainties. The solid and dotted lines represent a spectral index fit for two observing epochs. Top right: Derived temperature vs. radius profiles for Mira A. The solid and dotted lines correspond to the epochs in the top left figure. For illustration, the long-dashed line corresponds to a gray atmosphere temperature profile (RM97) with a temperature at of  K. Bottom left: The size vs. frequency relation at both observing epochs. The size denotes the average of the major and minor axis. Bottom right: The radially averaged real and imaginary visibilities and the best fit elliptical stellar disc model of Mira A for the three observed bands.
Figure 2: As Fig. 1 for R Dor. In the top right panel, the long-dashed line indicates the gray atmosphere temperature profile with a temperature at of  K.
Figure 3: As Fig. 1 for W Hya. In the top right panel, the long-dashed line indicates the gray atmosphere temperature profile with a temperature at of  K.
Figure 4: As Fig. 1 for R Leo. The six different epochs are denoted by the symbols and colors as indicated in the top left panel. No fits are indicated for the individual epochs in the temperature vs. radius (top-right) panel. Here, the long-dashed line indicates the gray atmosphere temperature profile with a temperature at of  K. In the bottom left, the lines indicate the maximum and minimum power-law fit to the size vs. frequency relation when correcting each epoch for the observed expansion (§ 5.2). In the uv-plot only two epochs, one for each observing band, are indicated for illustration.

Most of our analysis is performed using the uv-fitting code uvmultifit (Martí-Vidal et al., 2014) in a similar manner as described in Vlemmings et al. (2017). We fit a uniform elliptical stellar disc for each of the spws in our observations. This provides us with the flux density, major axis, axis ratio and position angle as a function of frequency and time. For Mira A, we performed the fits after its companion Mira B was subtracted from the data. The results of our fits are listed in Table 3 along with the observing date and stellar phase . In Matthews et al. (2015), an additional uncertainty was added to the size measurements to account for residual phase errors on the long baselines that arise from tropospheric phase variations. This uncertainty depends on the baseline length and rms phase fluctuations for the longest baselines. When using self-calibration, the phase solutions are directly determined for each integration interval on source. Hence, the rms phase fluctuations are significantly reduced. Adopting the approach described in Matthews et al. (2015) we conservatively find that a relative error term, with respect to the measured size, of on average in Band 4, in Band 6 and in Band 7 should be added. The exact fractions depend on the measured size and baseline lengths. The errors in Table 3 reflect this added uncertainty.

The spectral indices provided in Table 3 are determined in two ways, depending on whether different epochs were combined or not. When epochs are combined for a fit, we include a absolute flux density uncertainty for the different bands (except for R Leo where we adopt a uncertainty in Band 4). In that case, all spws in the band are assigned the same fractional error and the spectral index and its associated error are determined by a Monte Carlo scheme. For several epochs we also present a spectral index fit using only the four (or in some cases three) spectral windows. In that case, no additional flux density errors are assumed besides the measurement uncertainties. In neither case did we include the uncertainty in the assumed spectral index of the flux calibrator. This error is of the order of for our calibrators.

Although a uniform disc model provides good fits to all our observations, residual structure clearly remains. This is discussed in more detail in § 5.3. While in most cases this concerns small scale structure, we found that for a number of spws, particularly those in Band 7 of the Mira observations, excess flux at short baselines affected the fits. This excess can be seen in the lower right panel of Fig. 1 for baseline lengths below  M. For the affected spws we restricted our fits to uv-distances , or  M for Band 4 and for Bands 6 and 7, respectively. The possible origin of the excess, particularly for Mira A, is discussed in § 5.4.

The measured sizes and flux densities are used to obtain the source brightness temperature using:

(1)

Here is the stellar flux density in mJy, the major axis in mas, the axis ratio (), the frequency in GHz, and the light-speed in m s. The results for the flux density vs. frequency, size vs. frequency, and temperature vs. size relations are shown in Figs  123, and 4. These figures also show the quality of the uv-fits.

For illustration, we also produced images for all of the observations of Mira A, R Dor, and W Hya, as well as three of the epochs for R Leo. These are shown in Figs 567, and 8. The images were produced using uniform weighting in order to reach the highest resolution. Hence, the rms in the images is significantly higher than the rms reached when using Briggs weighting with a robust parameter of (presented in Table 1). The figures also show the residuals after subtracting the best fit stellar disc model from the data in each spw.

Figure 5: Top row: Uniformly weighted continuum maps of Mira A in the three observed bands. The images serve as an illustration since the analysis is performed on the visibilities by uv-fitting. The images are produced for a single spw, using the highest frequency spw in each observing band. Contours are drawn at and , where and  mJy beam from left to right. The beam sizes, from left to right, are , , and  mas. Bottom row: Residual images after subtracting the best fit uniform disc model presented in Table 3. The models are fitted for the spectral windows individually and the subsequent residual image is produced using all spectral windows. The beam sizes, from left to right, are , , and  mas. Contours are drawn at and , where and  mJy beam from left to right. The black ellipse indicates the size of the stellar disc.
Figure 6: As Fig.5 for R Dor. Here and  mJy beam from left to right in the top panels, and and  mJy beam from left to right in the bottom panels. The beam sizes are , , and  mas on the top row and , , and  mas on the bottom row.
Figure 7: As Fig.5 for the two observations of W Hya, with and  mJy beam from left to right in the top panels and and  mJy beam from left to right in the bottom panels. The beam sizes are , and  mas on the top row and , and  mas on the bottom row.
Figure 8: As Fig.5 for R Leo. Only three of the epochs are shown for illustration. The Band 6 epoch from 11 October was restored with a circular beam in order to highlight the northwest extension. In the top row, and  mJy beam from left to right. In the bottom row, and  mJy beam from left to right. The beam sizes are , , and  mas on the top row, and , , and  mas on the bottom row.

5 Discussion

5.1 Parameterized model atmospheres

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 :

(2)

Here is the number density (in ) at the stellar photosphere, is the photospheric radius, and 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 . 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 (Freytag et al., 2017). Density profiles with are also typically used to describe the molecular atmosphere (e.g. Wong et al., 2016; Khouri et al., 2016b). We adopt for our reference model.

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

(3)

or a power-law given by:

(4)

Here is the temperature (in K) at the photosphere and the power-law index. Finally, we approximate the ionisation profile also with a power-law in temperature:

(5)

Here , where is the electron number density. We fix the ionisation at  K to be , 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 , which fits the range provided by the detailed model calculations from RM97. The exact adopted value of 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 ( cm) in the extended atmosphere, the equilibrium time is  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  cm 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 across the temperature range that we probe in the AGB atmosphere.

 cm
 K
gray atmosphere
Table 4: Model atmosphere: fixed parameters

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 , which is the radius where the brightness temperature has dropped by in units of stellar radii . The angular size, , can be directly related to 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. 123, and 4 for the observations. As the opacity is proportional to the density (with exponent ) 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 , we see that the optical depth . 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.

Figure 9: The brightness temperature vs. size (left column) and size vs. frequency (right column) relations for model atmospheres (see text) with different power-law slopes for the density profile , top), temperature profile (, middle), and ionisation (, bottom). The two dashed line indicate the range of the observed slopes of the size vs. frequency relation (). Standard parameters are given in Table 4. We vary from (blue) to (red) in steps of . We vary from (blue) to (red) in steps of and include the gray atmosphere model (dark blue). We vary from (blue) to (red) in steps of .

In Fig. 9, we show the size vs. frequency and brightness temperature vs. size relations for models with different and . The fixed model parameters for the plotted models are given in Table 4. In the figures we also indicate the observed range of slopes of the size vs. frequency relation derived from the observations. From this is it clear that the lowest of the derived values of , observed for Mira A, R Dor and R Leo, cannot be reproduced for any but the steepest density profiles (), 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 . We thus produce models, which we describe as a ’layer’ model, with a times increase in at a single ’shock’ radius () 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.

Figure 10: The extended atmosphere models for Mira A (top left), R Dor (top right), W Hya (bottom left) and R Leo (bottom right) that can reproduce our observations. For each star the top panel shows the temperature vs. radius relation and the bottom panel shows the size vs. frequency relation. The solid circles are the measurements presented in this paper, and are taken as an average of all spws per observations. The open circles are 46 GHz radio observations (Matthews et al., 2015, 2018) shown for comparison. The radio observations were taken several years before the ALMA observations. In the panels for W Hya, the open square indicates the ALMA observations taken approximately 2 years earlier (Vlemmings et al., 2017). The thick solid lines are the models described in the text. In the top panels for each star, the red dashed line represents the assumed temperature profile. In the bottom panels for each star, the red dashed lines represent the size vs. frequency power-law relation fitted to the observations. For Mira A and R Dor, the dotted line represents a model without a high-opacity layer. For W Hya the dotted line represents the gray atmosphere model with that comes closest to describing the radio observations. For R Leo the dotted line represents the model for which the layer location has shifted (see text).

5.1.1 Mira A

The slope of the size vs. frequency relation for Mira A is shallow, with determined from the Band 4 and 6 observations, and between the Band 7 observations taken  month later and the Band 4 observations. The observations can be reproduced by a typical density profile with  cm and . For the temperature, a gray atmosphere with  K is consistent with the measured brightness temperatures. Additionally, the low and apparent abrupt drop of between the Band 6 and 7 observations requires, in our model framework, a sudden increase (moving inwards) in opacity by at least a factor of 10 at . We do not cover a sufficient range of radii to determine the thickness of the region with elevated opacity. Considering the relation between opacity and the modeling parameters, the opacity change can be produced by an increase of the density by at least a factor of three or a change of ionisation by a factor of or a combination of these. A significant change in temperature would have affected the brightness temperature vs. size relation. To indicate how our model compares with radio-wavelength observations at 46 GHz (Matthews et al., 2015), we also present the average of the major and minor axis of these observations in the figure despite that the radio epoch was three years before the ALMA observations. We find that our model underestimates the brightness temperature measured at 46 GHz by almost while it underestimates the radius by . Considering that the temperature uncertainty of the radio observations is and the ellipticity is almost , the model could be consistent. The size measurements of the previous ALMA observations at 96 and 230 GHz are also consistent but the derived brightness temperatures are higher (Vlemmings et al., 2015; Matthews et al., 2015).

The layer of enhanced opacity we derive for Mira A is likely related to shocks originating from stellar pulsations and/or convection. Since we do not have sufficient observing epochs spread over a few weeks time-scale taken at similar frequencies, we are unable to determine a motion of the layer. Alternatively, it might not be an actual spherical layer, as we observe asymmetries such as a mild elliptical shape and compact structures in our residual maps. Our models are not yet able to determine the effect of such structures. The derived density and temperature relations should thus be seen as rough indications only. Still, deviations from the densities we find by a factor of three or more would, in our framework, produce significantly different stellar sizes and are thus unlikely unless compensated by differences in ionisation or another opacity source.

5.1.2 R Dor

The sizes of R Dor observed at the different frequencies are, in stellar radii, the smallest of our sample. Similar to Mira A, the change of size with frequency is low ( between Band 4 and Band 6, and in Band 7). A combination of a fast drop in temperature and the low also implies the existence of a layer with increased opacity around R Dor. Since the measured size is small, we can reproduce the profiles with a lower density of  cm and . The temperature is well reproduced by a gray atmosphere with  K. The opacity enhancement (by a factor of ) around R Dor occurs at . For R Dor, the ellipticity increases markedly towards the higher frequencies in Band 7 and also several clear residual emissions occur both toward the star as well as in the limbs. Although the similar caveats should be considered as for Mira A, the observations imply that the layer is related to the asymmetric structures that are likely (convective) shocks. Alternatively or additionally, the increased opacity so close to the photosphere could arise due to the presence of a higher ionisation region related to a possible chromosphere. No long wavelength observations of R Dor exist for comparison because of its low declination.

5.1.3 W Hya

W Hya is the only star in our sample that displays a slope that can be reproduced without adding a region with enhanced opacity ( in Band 4, and in Band 7). For this, we need a density profile with  cm and . A temperature profile with  K and produces a slightly better correspondence with the data than a gray body atmosphere model. However, a gray body atmosphere with  cm and is better able to reproduce the radio observations as was already noted for W Hya by RM97. This likely reflects the complex nature of the AGB extended atmosphere where a simple single power-law is not sufficient to characterize the density and temperature. W Hya also shows significant residuals from a simple stellar disc model, both in the observations presented here and the previous ALMA observations (Vlemmings et al., 2017).

5.1.4 R Leo

R Leo presents a complex case. As will be discussed in § 5.2, a clear expansion is seen. This points to an outward motion of a shock wave. When adjusting each epoch for the fitted outward motion of  km s and fitting a single power-law to the adjusted size vs. frequency relation, we find . We can reproduce the range of observed sizes and temperatures by introducing an increased opacity layer as shown in Fig. 10. The layer moves from to and the underlying density profile has  cm and . The temperature is described by a gray atmosphere model with  K. Since it appears that at least part of the expansion of R Leo is asymmetric, with significant residuals from a uniform stellar disc towards the north-east, the apparent layer might in reality be a directional shock-wave propagating through the extended atmosphere.

5.2 Expansion of R Leo

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 . 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  mas day 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  km s 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 . 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 (most likely due to absolute flux calibration uncertainties). Such a variation would produce a change in ionisation of 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 . 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 or , 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  km s and that of the minor axis to  km s. 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  km s at for a star with  M. The velocity exceeds the local sound speed in the atmosphere, which is  km s. It is also larger than the velocities of  km s found at in current convection models (Freytag et al., 2017).

Figure 11: A uv-plot comparison of the real part of the visibilities between R Leo Band 4 observations taken on 2017-09-30 (black solid line and circles) and 2017-10-13 (red dashed line and squares). In most cases the error bars are less than the symbol size. The lines represent the uniform disc model, where the long-dashed blue line corresponds to the best fit model of the 2017-10-13 observations scaled to the flux density of the 2017-09-30 observations. Although the difference in flux density is most likely predominantly caused by the uncertainty in absolute flux calibration, the long-dashed blue line shows that the change in radius of the model is robustly determined.
Figure 12: The average diameter vs. observing time after the first epoch (2017-09-30) for the star R Leo. The black circles indicate the radius at 140 GHz in Band 4 and the red squares are the fitted radii at 225 GHz in Band 6. The black solid line is the error weighted fitted expansion of the optically thick surface. The expansion velocity corresponds to  km s (assuming a distance of  pc). The red dashed line is the fit to the  GHz radii shifted by a time-lag of days.

5.3 Asymmetries and hotspots

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  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 567, 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 (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 or , 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 for both brightest hotspot of  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  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  R in depth. If the structures are indeed the same, they would thus have a vertical extent of at least  cm ( 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 at a position angle of . 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  mas, which implies  K. In Band 7, the two spots are very compact ( mas) with  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 ) is still optically thick at longer wavelengths (where we find a radius of ), 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 ( mas implying  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 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.

5.4 Extended Emission Surrounding Mira A

Figure 13: The residual map of the extended emission around Mira A at 332 GHz. We have subtracted the best fit uniform disc model after fitting to the data with uv-distance  M. The residual image is made using Briggs weighting and restricting the baselines to  M, which results in a synthesized beam of  mas at a position angle of (indicated in green in the bottom-right corner). The black contours are the corresponding map of Mira A using all data drawn at and times the peak intensity. The beam size for all data is shown in red in the bottom-right corner. The white contour is drawn at , with  mJy beam. The residual emission has a peak of  mJy beam which corresponds to of the stellar flux density.

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  mJy at  GHz to  mJy at  GHz. It is immediately apparent that the extended emission peaks towards the east of the star and extends over almost . The extent and morphology bears resemblance to the emission seen in SO, CO and TiO 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  M at  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 . 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.

6 Conclusions

We have shown that multi-epoch and multi-frequency high angular-resolution observations with ALMA can be used to probe the extended atmospheres of AGB stars. This has allowed us to constrain the density, ionisation and temperature between . The main source of uncertainty on the temperature determination has been shown to be the absolute flux calibration accuracy of ALMA. Specifically, multi-epoch observations of the same flux- and gain-calibrator pair over four epochs of ALMA Band 4 observations indicate that an uncertainty of can only be reached with sufficient observations within a short time-scale. For single epoch observations, the uncertainty, even at low frequency, is closer to . Still, size measurements as well as surface maps can be made with great accuracy. We find that:

  • Parameterized atmosphere models can describe the observations for typical surface hydrogen number densities between  cm and photospheric temperatures of  K. These models assume the electron-neutral free-free process as the dominant opacity source as previously found in RM97. However, for R Dor, Mira A, and R Leo a layer of increased continuum opacity is needed to reproduce the shallow slope of the size vs. frequency relation. This layer, of unknown thickness, could be due to a sudden increase of neutral hydrogen density or effective ionisation or a combination of this. Alternatively an unknown additional opacity source is needed. It is likely that the layer is non-uniform and it might be related to convective shocks.

  • A mottled structure of spots can be seen most clearly on the surface of R Dor and W Hya. The brightness temperature of the brightest spot of R Dor is  K while the two brightest spots on the surface of W Hya have  K, consistent with previous observations (Vlemmings et al., 2017). The spots likely represent the shocks that cause the formation of a low-filling factor chromosphere. The time-scale of these spots is of the order of a few weeks at most.

  • Our six epochs of observations of R Leo show a linear increase in measured size with an expansion velocity of  km s. We suggest that the expansion is related to the motion of the aforementioned high opacity layer. The velocity is consistent with the typical line-widths seen in the stellar atmosphere but does not exceed the escape velocity.

  • In particular for Mira A, we find significant excess flux at the shortest baseline lengths compared to a uniform stellar disc model. The most likely explanation is that we are seeing either an asymmetric distribution of optically thin neutral-ion free-free emission out to almost  R, or the effect of stacking weak emission lines in what were thought to be line-free regions of the spectrum. However, although unlikely, we cannot yet fully rule out that the emission comes from warm circumstellar dust.

Acknowledgements.
This work was supported by ERC consolidator grant 614264. WV, TK, and HO also acknowledge support from the Swedish Research Council. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2016.1.00004.S, ADS/JAO.ALMA#2016.A.00029.S, ADS/JAO.ALMA#2017.1.00075.S, ADS/JAO.ALMA#2017.1.00191.S and ADS/JAO.ALMA#2017.1.00862.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), NSC and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. We also acknowledge support from the Nordic ALMA Regional Centre (ARC) node based at Onsala Space Observatory. The Nordic ARC node is funded through Swedish Research Council grant No 2017-00648.

References

  • Bedding et al. (1998) Bedding, T. R., Zijlstra, A. A., Jones, A., & Foster, G. 1998, MNRAS, 301, 1073
  • Bedding et al. (1997) Bedding, T. R., Zijlstra, A. A., von der Luhe, O., et al. 1997, MNRAS, 286, 957
  • Bowen (1988) Bowen, G. H. 1988, ApJ, 329, 299
  • Briggs (1995) Briggs, D. S. 1995, in Bulletin of the American Astronomical Society, Vol. 27, American Astronomical Society Meeting Abstracts, 1444
  • Cherchneff (2006) Cherchneff, I. 2006, A&A, 456, 1001
  • Dalgarno & Lane (1966) Dalgarno, A. & Lane, N. F. 1966, ApJ, 145, 623
  • Danilovich et al. (2017) Danilovich, T., Lombaert, R., Decin, L., et al. 2017, A&A, 602, A14
  • Danilovich et al. (2015) Danilovich, T., Teyssier, D., Justtanont, K., et al. 2015, A&A, 581, A60
  • De Beck & Olofsson (2018) De Beck, E. & Olofsson, H. 2018, A&A, 615, A8
  • Freytag & Höfner (2008) Freytag, B. & Höfner, S. 2008, A&A, 483, 571
  • Freytag et al. (2017) Freytag, B., Liljegren, S., & Höfner, S. 2017, A&A, 600, A137
  • Haniff et al. (1995) Haniff, C. A., Scholz, M., & Tuthill, P. G. 1995, MNRAS, 276, 640
  • Harper et al. (2001) Harper, G. M., Brown, A., & Lim, J. 2001, ApJ, 551, 1073
  • Hinkle et al. (1997) Hinkle, K. H., Lebzelter, T., & Scharlach, W. W. G. 1997, AJ, 114, 2686
  • Hofmann et al. (2001) Hofmann, K. H., Balega, Y., Scholz, M., & Weigelt, G. 2001, A&A, 376, 518
  • Höfner (2008) Höfner, S. 2008, A&A, 491, L1
  • Höfner et al. (2003) Höfner, S., Gautschy-Loidl, R., Aringer, B., & Jørgensen, U. G. 2003, A&A, 399, 589
  • Höfner & Olofsson (2018) Höfner, S. & Olofsson, H. 2018, A&A Rev., 26, 1
  • Ireland et al. (2007) Ireland, M. J., Monnier, J. D., Tuthill, P. G., et al. 2007, ApJ, 662, 651
  • Ireland et al. (2004) Ireland, M. J., Tuthill, P. G., Bedding, T. R., Robertson, J. G., & Jacob, A. P. 2004, MNRAS, 350, 365
  • Johnson & Luttermoser (1987) Johnson, H. R. & Luttermoser, D. G. 1987, ApJ, 314, 329
  • Khouri (2018) Khouri, T. 2018, in Imaging of Stellar Surfaces, 17
  • Khouri et al. (2014a) Khouri, T., de Koter, A., Decin, L., et al. 2014a, A&A, 561, A5
  • Khouri et al. (2014b) Khouri, T., de Koter, A., Decin, L., et al. 2014b, A&A, 570, A67
  • Khouri et al. (2016a) Khouri, T., Maercker, M., Waters, L. B. F. M., et al. 2016a, A&A, 591, A70
  • Khouri et al. (2019) Khouri, T., Velilla-Prieto, L., De Beck, E., et al. 2019, A&A, 623, L1
  • Khouri et al. (2018) Khouri, T., Vlemmings, W. H. T., Olofsson, H., et al. 2018, A&A, 620, A75
  • Khouri et al. (2016b) Khouri, T., Vlemmings, W. H. T., Ramstedt, S., et al. 2016b, MNRAS, 463, L74
  • Knapp et al. (2003) Knapp, G. R., Pourbaix, D., Platais, I., & Jorissen, A. 2003, A&A, 403, 993
  • Knapp et al. (1993) Knapp, G. R., Sandell, G., & Robson, E. I. 1993, ApJS, 88, 173
  • Knapp et al. (1998) Knapp, G. R., Young, K., Lee, E., & Jorissen, A. 1998, ApJS, 117, 209
  • Lebzelter et al. (2005) Lebzelter, T., Hinkle, K. H., Wood, P. R., Joyce, R. R., & Fekel, F. C. 2005, A&A, 431, 623
  • Lim et al. (1998) Lim, J., Carilli, C. L., White, S. M., Beasley, A. J., & Marson, R. G. 1998, Nature, 392, 575
  • Maercker et al. (2008) Maercker, M., Schöier, F. L., Olofsson, H., Bergman, P., & Ramstedt, S. 2008, A&A, 479, 779
  • Martí-Vidal et al. (2014) Martí-Vidal, I., Vlemmings, W. H. T., Muller, S., & Casey, S. 2014, A&A, 563, A136
  • Matthews et al. (2015) Matthews, L. D., Reid, M. J., & Menten, K. M. 2015, ApJ, 808, 36
  • Matthews et al. (2018) Matthews, L. D., Reid, M. J., Menten, K. M., & Akiyama, K. 2018, AJ, 156, 15
  • McMullin et al. (2007) McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in Astronomical Society of the Pacific Conference Series, Vol. 376, Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell, 127
  • Mohamed & Podsiadlowski (2012) Mohamed, S. & Podsiadlowski, P. 2012, Baltic Astronomy, 21, 88
  • Montez et al. (2017) Montez, Jr., R., Ramstedt, S., Kastner, J. H., Vlemmings, W., & Sanchez, E. 2017, ApJ, 841, 33
  • Norris et al. (2012) Norris, B. R. M., Tuthill, P. G., Ireland, M. J., et al. 2012, Nature, 484, 220
  • O’Gorman et al. (2017) O’Gorman, E., Kervella, P., Harper, G. M., et al. 2017, A&A, 602, L10
  • Ohnaka et al. (2016) Ohnaka, K., Weigelt, G., & Hofmann, K.-H. 2016, A&A, 589, A91
  • Ohnaka et al. (2017) Ohnaka, K., Weigelt, G., & Hofmann, K.-H. 2017, A&A, 597, A20
  • Olofsson et al. (2002) Olofsson, H., González Delgado, D., Kerschbaum, F., & Schöier, F. L. 2002, A&A, 391, 1053
  • Ramstedt et al. (2014) Ramstedt, S., Mohamed, S., Vlemmings, W. H. T., et al. 2014, A&A, 570, L14
  • Ramstedt & Olofsson (2014) Ramstedt, S. & Olofsson, H. 2014, A&A, 566, A145
  • Reid & Goldston (2002) Reid, M. J. & Goldston, J. E. 2002, ApJ, 568, 931
  • Reid & Menten (1997) Reid, M. J. & Menten, K. M. 1997, ApJ, 476, 327
  • Reid & Menten (2007) Reid, M. J. & Menten, K. M. 2007, ApJ, 671, 2068
  • Ryde & Schöier (2001) Ryde, N. & Schöier, F. L. 2001, ApJ, 547, 384
  • Tatebe et al. (2008) Tatebe, K., Wishnow, E. H., Ryan, C. S., et al. 2008, ApJ, 689, 1289
  • Van de Sande et al. (2018) Van de Sande, M., Decin, L., Lombaert, R., et al. 2018, A&A, 609, A63
  • van Langevelde et al. (2019) van Langevelde, H., Quiroga-Nuñez, L. H., Vlemmings, W., et al. 2019, arXiv e-prints
  • van Leeuwen (2007) van Leeuwen, F. 2007, A&A, 474, 653
  • Vlemmings et al. (2017) Vlemmings, W., Khouri, T., O’Gorman, E., et al. 2017, Nature Astronomy, 1, 848
  • Vlemmings et al. (2018) Vlemmings, W. H. T., Khouri, T., Beck, E. D., et al. 2018, A&A, 613, L4
  • Vlemmings et al. (2015) Vlemmings, W. H. T., Ramstedt, S., O’Gorman, E., et al. 2015, A&A, 577, L4
  • Vlemmings et al. (2003) Vlemmings, W. H. T., van Langevelde, H. J., Diamond, P. J., Habing, H. J., & Schilizzi, R. T. 2003, A&A, 407, 213
  • Watson et al. (2006) Watson, C. L., Henden, A. A., & Price, A. 2006, Society for Astronomical Sciences Annual Symposium, 25, 47
  • Woitke (2006) Woitke, P. 2006, A&A, 460, L9
  • Wong et al. (2016) Wong, K. T., Kamiński, T., Menten, K. M., & Wyrowski, F. 2016, A&A, 590, A127
  • Woodruff et al. (2009) Woodruff, H. C., Ireland, M. J., Tuthill, P. G., et al. 2009, ApJ, 691, 1328
  • Woodruff et al. (2008) Woodruff, H. C., Tuthill, P. G., Monnier, J. D., et al. 2008, ApJ, 673, 418
  • Zhao-Geisler et al. (2011) Zhao-Geisler, R., Quirrenbach, A., Köhler, R., Lopez, B., & Leinert, C. 2011, A&A, 530, A120
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
""
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
   
Add comment
Cancel
Loading ...
354907
This is a comment super asjknd jkasnjk adsnkj
Upvote
Downvote
""
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters
Submit
Cancel

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test
Test description