Heading Gaia to measure atmospheric dynamics in AGB stars

Context. Asymptotic Giant Branch stars are characterised by complex stellar surface dynamics that affect the measurements and amplify the uncertainties on stellar parameters. The resulting motion of the stellar photo-center could have adverse effects on the parallax determination with Gaia. Aims. We explore the impact of the convection-related surface structure in AGBs on the photocentric variability. We quantify these effects to characterise the observed parallax errors and estimate fundamental stellar parameters and dynamical properties. Methods. We use 3D radiative-hydrodynamics simulations of convection with CO5BOLD and the post-processing radiative transfer code Optim3D to compute intensity maps in the Gaia G band. From those maps, we calculate the intensity-weighted mean of all emitting points tiling the visible stellar surface and evaluate its motion as a function of time. We extract the parallax error from Gaia DR2 for a sample of semiregular variables in the solar neighbourhood and compare it to the synthetic predictions of photo-center displacements. Results. AGB stars show a complex surface morphology characterised by the presence of few large scale long-lived convective cells accompanied by short-lived and small scale structures. As a consequence, the position of the photo-center displays temporal excursions between 0.077 to 0.198 AU (5 to 11% of the corresponding stellar radius), depending on the simulation considered. We show that the convection-related variability accounts for a substantial part to the Gaia DR2 parallax error of our sample of semiregular variables. Finally, we put in evidence for a correlation between the mean photo-center displacement and the stellar fundamental parameters: surface gravity, radius, and pulsation. We denote that parallax variations could be exploited quantitatively using appropriate RHD simulations corresponding to the observed star.


Introduction
Low-to intermediate-mass stars evolve to red giant and asymptotic giant branch increasing the mass-loss during this evolution. They are characterised: (i) by large-amplitude variations in radius, brightness, and temperature of the star; and (ii) by a strong mass-loss rate driven by an interplay between pulsation, dust formation in the extended atmosphere, and radiation pressure on the dust (Höfner & Olofsson 2018). Their complex dynamics affect the measurements and amplify the uncertainties on stellar parameters.
Gaia (Gaia Collaboration 2016) is an astrometric, photometric, and spectroscopic space-borne mission. It performs a survey of a large part of the Milky Way. The second data release (Gaia DR2) in April 2018 (Gaia Collaboration 2018) brought highprecision astrometric parameters (i.e. positions, parallaxes, and proper motions) for over 1 billion sources brighter that G ≈ 20. Among all the objects that have been observed, the complicated atmospheric dynamics of asymptotic giant branch (AGB) stars affect the photocentric position and, in turn, their parallaxes (Chiavassa et al. 2011). The convection-related variability, in the context of Gaia astrometric measurements, can be considered as "noise" that must be quantified in order to better characterise any resulting error on the parallax determination. However, important information about stellar properties, such as the fundamental stellar parameters, may be hidden behind the Gaia measurement uncertainty.
In this work we explore the effect of convection-related surface structures on the photocentre to estimate its impact on the Gaia astrometric measurements.

Methods
We used the radiation-hydrodynamics (RHD) simulations of AGB stars (Freytag et al. 2017) computed with CO 5 BOLD (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 in a three-dimensional (3D) cartesian grid. It is assumed that solar abundances are appropriate for M-type AGB stars. The basic stellar parameters of the RHD simulations are reported in Table 1. The configuration used is the "star-in-a-box", where the evolution of the outer convective envelope and the inner atmosphere of AGB stars are taken into account in the calculation. In the simulations, convection, waves, and shocks all contribute to the dynamical pressure and, therefore, to an increase of the stellar radius and to a levitation of material into layers where dust can form. No dust is included in any of the current models. The regularity of the pulsations A&A 617, L1 (2018)  Notes. The table shows the simulation name, the mass M , the average emitted luminosity L , the average approximate stellar radius R (we note that the radii vary by about 20% during one pulsation period), effective temperature T eff , surface gravity log g, the pulsation period P puls , the half of the distribution of the pulsation frequencies σ puls , and the stellar time t avg used for the averaging of the rest of the quantities. All these quantities are from Freytag et al. (2017). The last four columns are the time-averaged value of the photocentre displacement P = (P 2 x + P 2 y ) (1/2) , its standard deviation (σ P ), and the time-average of P x and P y .  Table 1. The number on the top indicates the stellar times at which the two snapshots were computed. decreases with decreasing gravity as the relative size of convection cells increases. The pulsation period is extracted with a fit of the Gaussian distribution in the power spectra of the simulations. The period of the dominant mode increase with the radius of the simulation ( Table 1 in Freytag et al. 2017).
We computed intensity maps based on snapshots from the RHD simulations integrating in the Gaia G photometric system (Evans et al. 2018). 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, as for MARCS (Gustafsson et al. 2008) as a function of temperature, density, and wavelength for the solar composition (Asplund et al. 2009). Micro-turbulence broadening was switched off. The temperature and density distributions are optimised to cover the values encountered in the outer layers of the RHD simulations.
The surface of the deep convection zone has large and small convective cells. The visible fluffy stellar surface is made of shock waves that are produced in the interior and are shaped by the top of the convection zone as they travel outward (Freytag et al. 2017). In addition to this, on the top of the convection-related surface structures, other structures appear. They result from the opacity effect and dynamics at Rosseland optical depths smaller than 1 (i.e. further up in the atmosphere with respect to the continuum-forming region). At the wavelengths in Gaia G-band, TiO molecules produce strong absorption. Both of these effects affect the position of the photocentre and cause temporal fluctuations during the Gaia mission, as already pointed out for red supergiant stars in Chiavassa et al. (2011). In the Gaia G photometric system (Fig. 1), the situation is analogous with a few large convective cells with sizes of a third of the stellar radii (i.e. ≈0.6 AU) that evolve on a scale of several months to a few years, as well as a few short-lived (weeks to month) convective cells at smaller scales ( 10% of the stellar radius) We calculated the position of the photocentre for each map (i.e. as a function of time) as the intensity-weighted mean of the x−y positions of all emitting points tiling the visible stellar surface according to where I (i, j) is the emerging intensity for the grid point (i, j) with coordinates x(i, j), y(i, j) of the simulation, and N = 281 is the total number of grid points in the simulated box. In the presence of surface brightness asymmetries, the photocentre position will not coincide with the barycentre of the star and its position will change as the surface pattern changes with time. This is displayed in the photocentre excursion plots (Figs. A.1-A.5) for each simulation in the Appendix together with the time-averaged photocentre position ( P ) and its standard deviations (σ P ) in Table 1. The value of σ P (the third column from the right in Table 1) is mostly fixed by the short time scales corresponding to the small atmospheric structures. However, the fact that P x and P y do not average to zero (last two columns of Table 1 and, e.g. Fig. A.6 in the Appendix), means that the photocentre tends not to be centred most of the time on the nominal centre of the star, because of the presence of a large convective cell. Depending on the simulation, σ P varies between 0.077 and 0.198 AU (≈5 to ≈11% of the corresponding stellar radius). This measure of the mean photocentre noise induced by the stellar dynamics in the simulations is compared in the following section to Gaia measurement uncertainty to extract information on stellar parameters from the astrometric measurements. It should be noted that the main information that is used to determine the astrometric characteristics of every star will be the along-scan measurement of Gaia. Chiavassa et al. (2011) showed that the projection of the star position along the scanning direction of the satellite with respect to a known reference point discloses similar, though slightly increased, values of σ P . At the current state of the DR2, it is not possible to perform this on the data and we assumed the conservative value of σ P directly extracted from the RHD simulations for the following comparisons.

Observations
Evolved late-type stars show convection-related variability that may be considered, in the context of Gaia astrometric measurement, as "noise". Chiavassa et al. (2011) demonstrated that RHD simulations can account for a substantial part of the supplementary "cosmic noise" that affects Hipparcos measurements for some prototypical red supergiant stars. As a consequence, the convection-related noise has to be quantified in order to better characterise any resulting error on the parallax.
We extracted the parallax error (σ ) from Gaia DR2 for a sample of semiregular variables (SRV) from Tabur et al. (2009), Glass & van Leeuwen (2007), and Jura et al. (1993) that match the theoretical luminosities of RHD simulations (Table 1). Moreover, we only considered stars with 4000 < L < 8000 in order to compare with our simulations. It has to be noted that σ may still vary in the following data releases because: (i) the mean number of measurements for each source amounts to 26 (Mowlavi et al. 2018) and this number will be 70-80 in total at the end of the nominal mission; (ii) and new solutions may be applied to adjust the imperfect chromaticity correction (Arenou et al. 2018). We cross-identified our sample stars with the Gaia DR2 as well as with the distance catalogue of Bailer-Jones et al. (2018) that uses a weak distance prior, varying as a function of Galactic longitude and latitude, to derive distances for the Gaia DR2 1 . The apparent K-band magnitudes were transformed to absolute K magnitudes using the Gaia distances of Bailer-Jones et al. (2018) and neglecting the interstellar absorption which should be very small in the local neighbourhood. The absolute K magnitudes were converted to bolometric magnitudes using the bolometric correction formula of Kerschbaum et al. (2010). The typical uncertainties in the bolometric correction (BC K ) are of the order of ±0.1 mag. Finally, the bolometric magnitudes were converted to luminosities assuming a solar M bol = 4.7 (Torres 2010). The error bars on luminosities were calculated using the upper and lower distance limits provided by Bailer-Jones et al. (2018) for each of our sources together with the error of ±0.1 mag for BC K . We do not correct for the variation of the K-band light curve as only amplitudes in the visual are available and the K-band amplitude of AGB stars is in general much smaller than in the visual.

Comparison and predictions
AGB stars are characterised by complex stellar surface dynamics that affect the measurements and the determined stellar parameters. The uncertainties in observed absolute magnitudes have been shown to originate mainly from uncertainties in the parallaxes. In this section we investigate if the parallax errors of our SRV sample can be explained by the resulting motion of the stellar photocentre seen in the RHD simulations. Figure 2 (left panel) displays the parallax errors against the luminosity and compares these results to the standard deviations of the photocentre displacement in the simulations. The latter show good agreement with the observations. There are two luminosity families in the 3D simulations. In general the more luminous models are larger, that is, they have lower surface gravity, which causes larger convection cells and a more "fluffy" atmosphere (Freytag et al. 2017). As a consequence, for the simulation with higher luminosity (≈7000 L ), the parallax error ratio, defined as σ /σ P (where σ is the average error for the stars with luminosity close enough to the corresponding simulation luminosity), is [0.75-0.98]. This attests that convection-related variability accounts for a substantial part of the parallax error in Gaia measurements. For lower luminosities (≈5000 L ), the ratio is [1.15-1.95], indicating that for those simulations the situation is less clear. However, the observed and simulated luminosities do not coincide exactly and the observed error bars are still very large. One limitation of the existing model grid is the restriction to 1 M . In the future, there will be models with other masses available. For a better comparison, we need simulations and observations with known luminosities, masses, and radii; neither is trivial.
An important piece of information is indeed hidden in the Gaia measurement uncertainty: using the corresponding stellar parameters of the RHD simulations of Table 1, we plotted (Fig. 2) σ P as a function of surface gravity (central panel) and pulsation (right panel). They display a correlation between the mean photocentre displacement and the stellar fundamental parameters. While effective temperature does not show particularly correlated points, simulations with lower surface gravity (i.e. more extended atmospheres) return larger  Table 1 (σ P , star symbol in red). Central panel: σ P against the surface gravity for the RHD simulations. Right panel: σ P against the logarithm of the period. excursions of the photocentre. This behaviour is explained by the correlation between the stellar atmospheric pressure scale height (H p ≈ T eff g ) and the photocentre displacement (Freytag 2001;Ludwig 2006;Chiavassa et al. 2011).
One property of AGB stars that is well constrained by observations is the period-luminosity (P-L) relation. The uncertainties in the determination of this relation are mainly based on the calculation of the distances and on the different P-L relations used. Figure 2 (right panel) reveals the correlation between the photocentre displacement and the logarithm of the pulsation: larger values of σ P correspond to longer pulsation periods. Global shocks induced by large-amplitude, radial, and fundamentalmode pulsations have an impact on the detailed stellar structure of the stellar atmosphere, together with small-scale shocks. Both contribute to the levitation of material and the detected photocentre displacement (Freytag et al. 2017). This result is likely associated with the P-L relation found by Freytag et al. (2017), who showed that the RHD simulations reproduce the correct period for a given luminosity compared to the observations of Whitelock et al. (2009). Given the fact that σ P explains Gaia measurement uncertainties on the parallaxes (left panel), we suggest that parallax variations from Gaia measurements could be exploited quantitatively using appropriate RHD simulations. However, the parameter space in our simulations is still limited (Table 1). In the future we aim to extend our RHD simulations' parameters to lower and higher luminosities (i.e. shorter and longer periods) which will enable a more quantitative comparison with respect to the upcoming Gaia data releases.

Summary and conclusions
We used the snapshots from RHD simulations of AGB stars to compute intensity maps in the Gaia G photometric system. The visible fluffy stellar surface is made of shock waves that are produced in the interior and are shaped by the top of the convection zone as they travel outward. The surface is characterised by the presence of few large and long-lived convective cells accompanied by short-lived and small-scale structures. As a consequence, the position of the photocentre is affected by temporal fluctuations.
We calculated the standard deviation of the photocentre excursion for each simulation and found that σ P varies between 0.077 and 0.198 AU (≈5 to ≈11% of the corresponding stellar radius) depending on the simulation. We compared the measurement of the mean photocentre noise induced by the stellar dynamics in the simulations (σ P ) to the measurement uncertainty on the parallax of a sample of AGB stars in the solar neighbourhood cross-matched with data from the Gaia DR2. We found good agreement with observations, suggesting that convection-related variability accounts for a substantial part of the parallax error. It should be noted that σ may still vary in the following data releases due to the increase of Gaia's measurements in number and further corrections to the parallax solution.
Here we present evidence for a correlation between the mean photocentre displacement and the stellar fundamental parameters: surface gravity and pulsation. Concerning the latter, we showed that larger values of σ P correspond to longer pulsation periods. This result, associated with the P-L relation found by Freytag et al. (2017), and the good agreement between simulations and observations (σ P vs. σ ), suggest that parallax variations from Gaia measurements could be exploited quantitatively using appropriate RHD simulations corresponding to the observed star.  Table 1 in the Gaia G band filter. The different snapshots are connected by the line segments; the total time covered is reported in the Table. The dashed lines intersect at the position of the geometrical centre of the images.  Table 1.  Table 1.  Table 1.  Table 1.  Table 1.  Table 1.