Optical interferometry and Gaia measurement uncertainties reveal the physics of asymptotic giant branch stars

Context. Asymptotic giant branch stars are cool luminous evolved stars that are well observable across the Galaxy and populating Gaia data. They have complex stellar surface dynamics Aims. On the AGB star CL Lac, it has been shown that the convection-related variability accounts for a substantial part of the Gaia DR2 parallax error. We observed this star with the MIRC-X beam combiner installed at the CHARA interferometer to detect the presence of stellar surface inhomogeneities. Methods. We performed the reconstruction of aperture synthesis images from the interferometric observations at different wavelengths. Then, we used 3D radiative hydrodynamics simulations of stellar convection with CO5BOLD and the post-processing radiative transfer code Optim3D to compute intensity maps in the spectral channels of MIRC-X observations. Then, we determined the stellar radius and compared the 3D synthetic maps to the reconstructed ones focusing on matching the intensity contrast, the morphology of stellar surface structures, and the photocentre position at two different spectral channels, 1.52 and 1.70 micron, simultaneously. Results. We measured the apparent diameter of CL Lac at two wavelengths and recovered the radius using a Gaia parallax. In addition to this, the reconstructed images are characterised by the presence of a brighter area that largely affects the position of the photocentre. The comparison with 3D simulation shows good agreement with the observations both in terms of contrast and surface structure morphology, meaning that our model is adequate for explaining the observed inhomogenities. Conclusions. This work confirms the presence of convection-related surface structures on an AGB star of Gaia DR2. Our result will help us to take a step forward in exploiting Gaia measurement uncertainties to extract the fundamental properties of AGB stars using appropriate RHD simulations.


Introduction
The Gaia mission (Gaia Collaboration et al. 2016) is an astrometric, photometric, and spectroscopic space-borne mission, which in 2018 delivered high-precision astrometric parameters (i.e. positions, parallaxes, and proper motions) for over one billion sources (Gaia Collaboration et al. 2018). A considerable fraction of the intrinsically detected variable stars are long-period variables (LPVs, of which the survey contains more than 150 000, Holl et al. 2018) with large luminosity amplitudes and variability timescales, adequately sampled by Gaia (Gaia Collaboration et al. 2019). Among LPVs, there are cool luminous evolved stars with low to intermediate mass that reached a critical phase of their evolution at the end of the asymptotic giant branch (AGB). They are characterised by (i) large-amplitude variations in radius, brightness, and temperature of the star; (ii) a strong mass-loss rate (Ṁ = 10 −6 -10 −4 M ⊙ /yr, De Beck et al. 2010) driven by an interplay between pulsation, dust formation in the extended atmosphere, and radiation pressure on the dust (Höfner & Olofsson 2018); (iii) an inhomogenous visible surface (produced in the interior and shaped by the top of the convection zone: they travel outwards on timescales ranging from weeks to years (Freytag et al. 2017)); and (iv) strong molecular absorption bands in the optically thin region (e.g. Lançon & Wood 2000) and on the top of the convectionrelated surface structures.
The concordance of all these properties makes the measurement of accurate stellar parameters and mass-loss rate very com-A&A proofs: manuscript no. CL_Lac Notes. The first four columns display the G, Bp, Rp, and K photometric colours from Gaia and 2MASS; the two following columns indicate the luminosity (L ⋆ ) and effective temperature (T eff ); the next three show the parallax (p) with its error (σ ̟ ) in mas and in AU; the last column reports the pulsation period (P puls ). (a) Gaia DR2 (Gaia Collaboration et al. 2018) (b) 2MASS (Skrutskie et al. 2006) (c) Computed from K2MASS with a bolometric correction of 3.0 (d) The temperature is based on Gaia Bp/Rp photometric colours and was estimated for sources brighter than G=17 mag for 3000K<T eff <10000K.
Due to the degeneracy of the photometric temperature with interstellar extinction, an empirically trained machine learning algorithm was used for external spectroscopic datasets in low-extinction regions in order to derive temperatures (Andrae et al. 2018). However, due to the limited training set for cool stars, photometric temperatures have to be taken with extreme caution for these stars (g) Jura et al. (1993)  plex. In this context, a special emphasis has to be put on the problems with establishing their inter-dependence: a quantitative description of how the mass-loss rate depends on individual stellar parameters (and, consequently, how it changes as stars evolve) is essential for models of stellar and galactic chemical evolution (Höfner & Olofsson 2018). A noteworthy step forward in trying to find a way to quantitatively recover stellar parameters of AGB stars was proposed by Chiavassa et al. (2018) using Gaia DR2 data and stellar convection simulations (Freytag et al. 2017). The authors showed that the convection-related variability causes photocentre 1 displacements up to ≈ 11% of the stellar radius and accounts for a substantial part of the Gaia DR2 parallax error. They suggested that the fundamental properties of AGB stars could be measured directly from Gaia parallax errors. However, a final piece of information was missing, because the convection-related structures were quantified only indirectly (i.e. using photocentre displacement). Only with high angular resolution is it possible to detect the presence of stellar surface inhomogenities. This can be achieved using optical interferometry, which has already proven to be a powerful tool providing amazingly resolved aperture synthesis images of evolved stars (e.g. Chiavassa et al. 2010b;Montargès et al. 2014;Ohnaka et al. 2016;Montargès et al. 2016;Chiavassa et al. 2017;Montargès et al. 2017;Ohnaka et al. 2017;Wittkowski et al. 2017;Montargès et al. 2018;Paladini et al. 2018).
In this work, we show the presence of inhomogeneities on the stellar surface of the AGB star CL Lac, which is part of the sample of Gaia DR2 objects for which Chiavassa et al. (2018) showed that photocentre displacement explains the parallax error bars. The article is structured as follows: Section 2 introduces the interferometric observations and data reduction, while Section 3 evidences the interpretation of the data with 3D simulation of stellar convection. Section 4 summarises the conclusions and the impact on the determination of AGBs' stellar properties using Gaia parallax errors.

Interferometric observations with MIRC-X at CHARA
The main objective of this work is to detect the presence of stellar surface inhomogeneities for an object from Gaia's sample used in Chiavassa et al. (2018), and for which the convectionrelated variability accounts for a substantial part of the Gaia DR2 parallax error. We observed the LPV AGB star CL Lac with the MIRC-X beam combiner (Kraus et al. 2018;Anugu et al. 2018) and with the VEGA instrument (Mourard et al. 2009 On each night, we recorded three sets of MIRC-X data on CL Lac. Each set consisted of 10-15 min to record fringe data and another 10 min to record backgrounds and shutter sequences. In between these sets, we observed calibrator stars to calibrate the interferometric measurements. The calibrators and adopted angular diameters are listed in Table 2. In order to verify, we calibrated the calibrators against each other on each night. We saw no evidence for binarity (closure phases consistent with 0 • within their uncertainties). A fit to the calibrator visibilities provided angular diameters consistent with the adopted values. The data were reduced using the MIRC-X data reduction pipeline v1.2.0 2 to produce calibrated visibilities and closure phases. The MIRC-X detector is susceptible to bias in the bispectrum estimation as pointed out by Basden & Haniff (2004). The data pipeline corrects for this bispectrum bias with a method similar to the one outlined in Appendix C of that paper. Figure 1 displays the visibility and closure phase data for all the spectral channels as well as the UV-coverage. The data, and in particular A&A proofs: manuscript no. CL_Lac  Bourges et al. (2017) the closure phases, denote the presence of large departure points from the centrosymmetric case (i.e. closure phases not equal to 0 or ±180 • ). The calibrated OIFITS files are available on the Optical interferometry DataBase 3 . So as to obtain differential stellar radii in the optical and infrared, in order to retrieve information about the stellar dynamics and opacity run, we performed simultaneous observations with VEGA. We used the MIRC-X beam combiner as a fringe tracker. The longitudinal dispersion correctors were used to compensate for differences in dispersion between the visible and near-infrared. However, no fringes were spotted out in the optical with VEGA. This can be due by the combination of a large object (> 4.5 mas) and low squared visibilities (<0.1), very close to the limit of the instrumental capabilities. Nevertheless, we cannot definitively conclude on this point because the SNR was significantly lower than the confidence limit of three.

Aperture synthesis imaging
For the image reconstruction, we used the MIRA software package developed by Thiébaut (2008). This package is written in the scientific language yorick 4 . It reads interferometry data in the form of oifits files (Pauls et al. 2005), and the reconstructed image is compared to the visibility and closure phase data by means of Fourier transforms and a few additional steps. The pixel values are modulated to minimise a so-called cost function, which is the sum of a regularisation term plus data-related terms (χ 2 ). The data terms enforce agreement of the model image with the measured data (visibilities and closure phases). The regularisation term is a distance minimisation between the reconstructed image and an expected image, called the prior.
The data and the regularisation terms must be scaled to maximise the Bayesian evidence (the overall probability of the model integrated over all possible model parameter values). The scaling factor is called the 'hyperparameter' (µ).
We reconstructed images for the spectral channels individually (Fig. A.1) with: 0.2 mas/pixel, a 64x64-pixel field of view (FOV), a total variation regularisation (Renard et al. 2011), a circular Gaussian prior image with a FWHM equal to half the field of view, and a hyperparameter value of µ = 5 × 10 4 . This value of µ is confirmed by the L-curve plot shown in Fig. B.4, which shows the data term (χ 2 ) as a function of the hyperparameter (µ, all other image reconstruction parameters are random: pixel size, number of pixels, prior size, number of iterations): one can see that the χ 2 has a value close to 2 with µ ≤ 50000, and that the value raises by a large amount above. We show here the images convolved with a beam of FWHM λ/2B, meaning two times smaller as the theoretical angular resolution of the interferometer. This is to take advantage of the 'super-resolution' effect of image reconstruction with a sparse aperture reported in many works (e.g. Thiébaut 2008). For convenience, the image convolved with a beam of FWHM λ/B is also shown in Appendix B (Fig. B.5).
We also reconstructed images using 32×32 and 128×128 pixels with the same pixel size (see Fig. B.3), 64×64 images with different pixel sizes (from 0.1 to 0.4 mas, Fig. B.1), and hyperparameter values (from 5 to 5 × 10 6 , see Fig. B.2), to ensure that the images are not affected by reconstruction artifacts. We also checked that the global appearance of the image is similar when reconstructing a 'grey' image (i.e. computing an average across the spectral channels) and a chromatic image (Fig. B.3). Figure 2 displays two examples of aperture synthesis images at 1.52 and 1.70 µm. They were convolved with the resolution of the interferometer (Fig. 1, bottom panel). Both images clearly show a brighter area in the south-eastern part of the map that affects the position of the barycentre, which does not correspond to the geometrical one (dashed blue line). The morphology of CL Lac is complex and not centrosymmetric, largely affecting astrometric Gaia measurements . These stellar surface inhomogenities show a similar behaviour for all the spectral channels but with noticeable changes between 1.52 and 1.70 µm. Figure 2 (bottom panel) displays large differences (up to ∼ 20%), in particular close to the outer stellar disc borders.

Three-dimensional simulations of AGB stars and synthetic images
We use the radiation-hydrodynamics (RHD) simulation of AGB stars (Freytag et al. 2017) computed with CO 5 BO1D (Freytag et al. 2012) code. The code solves the coupled nonlinear equations of compressible hydrodynamics and non-local radiative energy transfer in the presence of a fixed external spherically symmetric gravitational field on a 3D cartesian grid. The averaged stellar parameters of the simulation used in this work are reported in Table 3. The atmosphere in the simulation is characterised by global and small-scale shocks induced by large-amplitude, radial, and fundamental-mode pulsations (Freytag et al. 2017). The pulsation period is extracted with a fit of the Gaussian distribution in the power spectra of the simulations. The shocks contribute to the levitation of material and to the photocentre displacement (Fig. A.5 for the temporal behaviour, Chiavassa et al. 2018). A photocentre displacement of σ P = 0.131 AU (Table 3) corresponds to ≈ 27R ⊙ or ≈9% of the stellar radius.
To prepare the comparison with the interferometric data, we computed intensity maps in the same observed spectral channels (Fig. 2) for the 100 snapshots of the RHD simulation. For this purpose, we employed the code Optim3D (Chiavassa et al. 2009), which takes into account the Doppler shifts caused by the convective motions. The radiative transfer is computed in detail using pre-tabulated extinction coefficients per unit mass, like for MARCS (Gustafsson et al. 2008) as a function of temperature, density, and wavelength for the solar composition (Asplund et al. 2009). The temperature and density distributions are optimised to cover the values encountered in the outer layers of the RHD simulation. Notes. The first column shows the simulation name, then the next five columns the stellar parameters such as the mass (M ⋆ ), the average luminosity (L ⋆ ), radius (R ⋆ ), effective temperature (T eff ), and surface gravity (log g). The seventh column displays the stellar simulated time (t avg ) used to average all the quantities. A solar metallicity is assumed. The eighth and ninth columns report the pulsation period (P puls ) and the half of the distribution of the pulsation frequencies (σ puls ). Finally, the last two columns are the time-averaged value of the photocentre displacement P and its standard deviation (σ P ).  Table 4. CL Lac apparent diameters obtained with the best-fitting 1D-MARCS model, Θ 1D , and the one from the 3D temporally averaged intensity profile (Fig. 4, black curve), Θ 3D . The first step of our analysis is the determination of the stellar radius. In order to do this, we used (i) 1D, stationary, hydrostatic models; and (ii) the 3D RHD simulation of Table 3.

Stellar radius thanks to the Gaia parallax
For item (i), we selected several models from the MARCSgrid 5 (Gustafsson et al. 2008) with stellar parameters close to CL Lac (Table 1) For item (ii), we obtained the 100 azimuthally averaged intensity profiles from the synthetic maps using the method explained by Chiavassa et al. (2009). They were constructed using rings that are regularly spaced in µ (where µ = cos(θ), with θ being the angle between the line of sight and the radial direction). Figure 4 shows that the temporal fluctuation from one snapshot to another (green curves) is large, each snapshot being approximatively 23 days apart. Afterwards, we calculated synthetic visibility amplitudes from the different 1D-MARCS models (purple curve) and the 3D temporally averaged (black curve) intensity profiles using the Hankel transform. Finally, we fitted the visibility data from Fig. 1 separately for the two spectral channels at 1.52 and 1.70 µm to obtain the apparent diameter of CL Lac. The best-fitting MARCS model has T eff =2800K and log g=-0.5. Table 4 reports the 1D and 3D stellar diameters as well as the stellar radius at a distance of d = 870 ± 113 pc (i.e. inversion of the Gaia parallax, Table 1). The 1D-MARCS model returns apparent radii larger than the 3D, ∼ 1% and 5% at 1.52 and 1.70 µm, respectively. This is already noticeable from the intensity profiles in Fig. 4. Our result is in qualitative agreement with those obtained by several authors for dwarf and K giant stars To sum up, the RHD simulation used here is the closest available to CL Lac in terms of 1. stellar parameters with, in particular, a very similar radius (R ⋆ =306±18, Table 3) with respect to the value we found (Table 4) and pulsation period; 2. photocentre variability in the Gaia G band (σ P =0.131, Table 3) caused by convection-related structures that accounts, as demonstrated by Chiavassa et al. (2018), for a substantial part of the CL Lac parallax error (σ ̟ =0.130, Table 1).
It has to be noted that the T eff based on Gaia Bp/Rp photometric colours from Table 1 is almost 500K larger than the 1D MARCS and 3D RHD simulation, that, in turn, have overlapping values within error bars. However, due to the limited training set for cool stars in Gaia actual DR2, photometric temperatures have to be taken with extreme caution for these stars.

Comparison with the 3D RHD simulation
In this section, we compare the synthetic maps to the aperture synthesis imaging. As a first step, we correct the limb-darkening (LD) effect dividing both the synthetic and reconstructed images. We proceeded as follows: (i) we created an LD model image (64×64 pixels and 0.2 mas/pixel resolution) with the same FOV as the reconstructed image using the 3D temporally averaged profile of Fig. 4; (ii) we rebinned the resolution of 0.2 mas/pixel the 3D synthetic maps so that they effectively possess the same FOV and resolution as the reconstructed images; and (iii) we convolved them to the resolution of the interferometric observations (Fig. 1, bottom panel) for both spectral channels at 1.52 and 1.70 µm.
To have strong constraints on the RHD simulation, we simultaneously inspected three different quantities for each map: [1] the intensity contrast in the maps; [2] the morphology of the surface structures; and [3] the position of the photocentre. We performed a minimisation procedure for each item to obtain the normalised χ 2 compiled = χ 2 contrast + χ 2 morphology + χ 2 photocentre . Starting from point [1], we divided pixel by pixel reconstructed images, and we convolved 3D synthetic maps by the LD one and used the definition of Tremblay et al. (2013) to estimate the intensity contrast as δI rms / I . The contrast depends on the radial cut chosen (i.e. the maximum radius adopted from the centre of the star from which outer pixels are not considered, Montargès et al. 2018). For CL Lac, we chose to cut at 0.80 stellar radii to avoid outer areas with rapid and steep increase. Figure 5 shows that several snapshots have very close values of contrast with respect to observations for both wavelength channels. The best matching 3D synthetic snapshot is the one with the closest intensity contrast with respect to the reconstructed images. Furthermore, we concentrated our analysis on the morphology of the surface structures (point [2]). To account for the unknown orientation of CL Lac with respect to the synthetic maps in the sky, we rotated the latter every 5 • around its centre, and we calculated the minimum difference with the reconstructed images. Several snapshots were found to be very similar to the observed images, but, as expected, none matched perfectly. Finally, we computed the position of the photocentre (point [3]) for each map and for the reconstructed images as the intensityweighted mean of the x − y positions of all emitting points tiling the visible stellar surface (Eq. 1 and 2 from Chiavassa et al. 2018). Then, we found the smallest difference between the photocentre in the observations and the synthetic maps. Among the different contributors to the χ 2 compiled , the morphology of the surface structures (χ 2 morphology ) is the most dominant. In addition to these three points, we also considered the difference between the reconstructed images at 1.52 and 1.70 µm (∆ reconstructed = image 1.52 −image 1.70 , bottom panel of Fig. 2) and for the simulation images (∆ 3D = image3D 1.52 − image3D 1.70 , bottom panel of Fig. 7). We aimed to have the closest difference (∆) between the two cases, and we defined χ 2 difference = ∆ reconstructed − ∆ 3D . This approach gives an additional and important wavelength dependence constraint for the simulation. In the end, the best matching snapshot is the one with the combination of the smallest χ 2 compiled versus χ 2 difference , with a higher weight for the latter. Figure 6 shows all the snapshots and all the different rotation angles considered, while Fig. 7 displays the best matching snapshot and rotation angle: the best combination returns χ 2 compiled =3.8 and χ 2 difference =0.8. Our 3D RHD simulation shows a good agreement with the observations both in terms of contrast an surface structures and wavelength dependence, meaning that it is adequate for explaining the inhomogenities observed in the reconstructed maps. This definitively confirms the presence of convectionrelated surface structures on the same Gaia DR2 object for which Chiavassa et al. (2018) proved that photocentre displacement explains the parallax error bars.

Discussion and conclusions
Stellar dynamics induce an intrinsic noise that increases the measurement uncertainty of the parallax of AGB stars in Gaia DR2 , among which there is the star CL Lac. However, the effect of convection-related surface structures was estimated only indirectly using the displacement of the photocentre (σ P ). In addition to this, the authors showed also a clear correlation between stellar parameters (surface gravity and pulsation period) and the standard deviation of the photocentre movements: 3D simulations with a lower surface gravity (i.e. more extended atmospheres) return larger excursions of the photocentre.
Thanks to optical interferometry, we managed, in this work, to discover the presence of stellar inhomogenities on the AGB CL Lac. Its morphology is complex, not centrosymmetric, and likely affects astrometric Gaia measurements. Using stateof-the-art 3D global simulations of stellar convection with CO 5 BOLD, we were able to interpret them as a consequence of small-and global-scale shocks that contribute to the levitation of material. This definitively confirms the presence of convectionrelated surface structures on the same Gaia DR2 object for which Chiavassa et al. (2018) discussed the photocentre displacement.
Our result will help us to make a step forward in exploiting Gaia measurement uncertainties using appropriate RHD simulations. The long-term aim is to uniquely extract the fundamental properties of AGB stars directly from Gaia errors and make it possible to systematically study the properties of convection in stars other than the Sun. These properties are paramount to understanding the physics of these stars: for example, the surface gravity controls the size of granules that, in turn, contribute to the stellar photometric variations. Similarly, the pulsation period provides important information about stellar (mean) interior as well as global shocks induced by large-amplitude, radial, and fundamental-mode pulsations. To achieve our goal, a series of steps are required. We list them here: 1. An additional study of the photometric colours from Bp and Rp photometric systems; Chiavassa et al. (2011) showed that red supergiant stars 6 have magnitude excursions up to 0.28 and 0.13 magnitudes in Bp and Rp, respectively. This impacts the determination of their parameters. We expected similar (or higher) values for AGBs. 2. One limitation of the existing model grid (Freytag et al. 2017) is the restriction to 1 M ⊙ . A significantly larger RHD simulation grid with a larger range of stellar masses is necessary to explore the large number of AGBs in Gaia (Holl et al. 2018). 3. With Gaia DR2, the mean number of measurements for each source amounts to ≈ 26 (Mowlavi et al. 2018). Finally, when Gaia DR4 is available in 2022, the number of measurements will increase to 70-80, possibly changing the parallax error. Additionally, DR4 will provide time-dependant measurements allowing a direct and more precise comparison with the time-dependent RHD simulations. We present here a subset of all the images that were reconstructed on CL Lac with the MIRA software to illustrate the robustness of the features discussed in this work. Image reconstruction with optical interferometry is known to have flaws due to the sparse sampling of the measurements, and also due to the intrinsically non-linear properties of the available interferometric observables (i.e. squared visibilities and closure phases). Therefore, we tested several of the input parameters of MIRA code, namely the pixel size (Fig. B.1), the value of the hyperparameter µ (Fig. B.2), and the FOV (i.e. changing the number of pixels, Fig. B.3). Figure B.1 shows that the image behaviour does not change drastically as a function of the pixel size: the disc of the star is the same size, and the bright area in the south-east is always located in the same place. We chose a pixel size of 0.2 mas. Then, Fig. B.2 displays that too low µ values (≤ 10000) return inhomogenous images, while too large values (≥ 100000) tend to flatten the surface structures. We selected a median value of µ = 50000. Finally, Fig. B.2 indicates that changing the image FOV does not affect its appearance.
Moreover, it has to be noted that the resolution of images reconstructed with MiRA is not necessarily λ/B, where B is the largest baseline of the interferometric beam. Depending on the signal-to-noise and UV-plane coverage of the data, a degree of super-resolution (up to a factor of 2-3) is possible. In Fig. B.5, we show two examples at 1.52 µm with other convolutions, namely an interferometric beam of FWHM λ/(B) and λ/(4B). These images have to be compared with those in Fig. 2. A. Chiavassa et al.: Optical Interferometry and Gaia measurement uncertainties      The intensity is normalised between [0, 1] and the contour lines are the same as in Fig. 2.