© ESO, 2011

1. Introduction

The Gaia mission (Perryman et al. 2001; Lindegren et al. 2008) will lead to the determination of the astrometric parameters of one billion objects with unprecedented precision. During its five years of satellite operations, every detected star will be observed on average 80 times, thus yielding its position (α₀, δ₀) at the reference epoch, its parallax (ϖ), and its proper motion (μ_α^∗ = μ_αcosδ, μ_δ). Additional parameters, such as the size of the orbit in the case of a binary, can also be adjusted. All these parameters derive from a χ² fit. To allow the comparison between models (with different number of parameters) the F2 statistics will be evaluated (Wilson & Hilferty 1931): $\mathit{F}\mathrm{2}\mathrm{=}\sqrt{\frac{\mathrm{9}\mathit{\nu }}{\mathrm{2}}}\left[\mathrm{3}\sqrt{\frac{{\mathit{\chi }}^{\mathrm{2}}}{\mathit{\nu }}}\mathrm{+}\frac{\mathrm{2}}{\mathrm{9}\mathit{\nu }}\mathrm{-}\mathrm{1}\right]\mathit{,}$(1)where ν is the number of degrees of freedom (the number of measurements minus the number of parameters of the fit). If χ² does follow a chi-square distribution with ν degrees of freedom, then F2 follows a N(0,1)-distribution irrespective of ν. A fit with a F2 larger than, say, 3 is considered poor.

Phenomena such as large convective cells or magnetic spots cause time-varying brightness asymmetries on the surface of stars. Their presence induces a time-dependent displacement of the photocentre with respect to the barycentre that possibly alters the astrometric parameters derived by Gaia (Bastian & Hefele 2005; Eriksson & Lindegren 2007). Svensson & Ludwig (2005) studied the dependence of the photocentric variability on the stellar parameters across the Hertzsprung-Russell diagram, from white dwarfs to red giants, using CO⁵BOLD hydrodynamical models. They found a linear relation between the standard deviation in the photocentre displacement (σ_P) and the pressure scale height at optical depth unity.

There is only one red-supergiant hydrodynamical 3D model available in literature, (Chiavassa et al. 2009), that predicts a standard deviation for the photocentre motion driven by convection of about 0.1 AU. This value, however, might be underestimated because of the approximations intrinsic to the model. Chiavassa et al. (2011) extended the analysis of Svensson & Ludwig (2005) with this red supergiant model. They argue that a radical change in the convective pattern is expected between the giant and supergiant stages but the large gap in T_eff/g between the two available models does not allow a law for σ_P as a function of T_eff/g to be unequivocally constrained. Hence, they proposed two different laws and performed a first analysis of the impact of this photocentric noise on the parallax, estimating the number of supergiants and giants for which it may be relevant.

The aim of this work is to investigate and provide an analytical prediction of the effects of these time-varying surface brightness asymmetries on the astrometric parameters derived by Gaia and to study their impact on the Gaia astrometric data processing. On the basis of the 3D supergiant model from Chiavassa et al. (2009), we simulated the photocentre motion due to convective cells and added it to Gaia simulated data (Sect. 2). We discuss the results of the astrometric fit of those data. In Sect. 3, we restrict our study to the single-star astrometric fit. Only five parameters (α₀, δ₀, ϖ, μ_α^∗, μ_δ) are derived in this solution. In Sect. 4, more complex, non single-star solutions are studied.

2. The model

Very large displacements of the photocentre are found in red supergiants, which are caused by the presence of extremely large convective cells. The flow of hot material towards the surface and the inflow of cold material around large convective cells produce huge surface brightness asymmetries. Evidence of this phenomenon are provided by interferometric observations of supergiants shch as α Ori (Buscher et al. 1990; Wilson et al. 1992; Tuthill et al. 1997; Young et al. 2000; Tatebe et al. 2007; Haubois et al. 2009; Ohnaka et al. 2011).

Fig. 1 Snapshot of the 3D supergiant model in the Gaia G band from Chiavassa et al. (2009). The colour scale denotes different temperatures. The 12 M⊙ star is modelled with an equidistant numerical mesh of 2353 grid points and a resolution of 8.6 R⊙ (or 0.040 AU). The luminosity average over spherical shells and over time (i.e., over 5 years) is L = 93   000 ± 1300 L⊙, Teff = 3490 ± 13 K, R = 832 ± 0.7 R⊙, and log g =  −0.337 ± 0.001.

A numerical 3D model of a red supergiant was developed by Chiavassa et al. (2009) using CO⁵BOLD (Freytag et al. 2002; Freytag & Höfner 2008). This model was used to predict spectra and intensity maps in the Gaia G band for the whole simulation time sequence, namely  ≈ 5 years with time steps  ≈ 23 days apart. A snapshot of the star is shown in Fig. 1.

The projections P_x and P_y of the photocentre position on either axis are central to the present study. We assume that P_x and P_y are random variables (of zero mean) that are both Gaussian and Markovian¹. Therefore the process correlation function ${\mathit{C}}_{\mathit{x,y}}\mathrm{\left(}\mathrm{\Delta }\mathit{t}\mathrm{\right)}\mathrm{=}\mathit{E}\mathrm{\left[}\mathrm{\right(}{\mathit{P}}_{\mathit{x,y}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{-}\mathit{E}\mathrm{\left[}{\mathit{P}}_{\mathit{x,y}}\mathrm{\right]}\mathrm{\left)}\mathrm{\right(}{\mathit{P}}_{\mathit{x,y}}\mathrm{(}\mathit{t}\mathrm{+}\mathrm{\Delta }\mathit{t}\mathrm{)}\mathrm{-}\mathit{E}\mathrm{\left[}{\mathit{P}}_{\mathit{x,y}}\mathrm{\right]}{\mathrm{\left)}}^{\mathrm{\right]}}$(2)(where the operator E[] denotes the mean value over t) reduces to (Doob 1942) $\mathit{C}\mathrm{\left(}\mathrm{\Delta }\mathit{t,\tau }\mathrm{\right)}\mathrm{=}{\mathit{\sigma }}_{\mathit{P}}^{\mathrm{2}}\exp \mathrm{[}\mathrm{-}\mathrm{\Delta }\mathit{t}\mathit{/}\mathit{\tau }\mathrm{]}\mathit{.}$(3)The random variables P_x,y at time t_i are given by ${\mathit{P}}_{\mathit{x,y}}\mathrm{\left(}{\mathit{t}}_{\mathit{i}}\mathrm{\right)}\mathrm{=}\exp \mathrm{(}\mathrm{-}\mathrm{\Delta }{\mathit{t}}_{\mathit{i}}\mathit{/}\mathit{\tau }\mathrm{)}{\mathit{P}}_{\mathit{x,y}}\mathrm{\left(}{\mathit{t}}_{\mathit{i}\mathrm{-}\mathrm{1}}\mathrm{\right)}\mathrm{+}{\mathit{g}}_{\mathit{i}}$(4)where Δt_i = t_i − t_i − 1, g_i is extracted from a Gaussian distribution N(0,s) with ${\mathit{s}}^{\mathrm{2}}\mathrm{=}\left[\mathrm{1}\mathrm{-}\exp \mathrm{(}\mathrm{-}\mathrm{2}\mathrm{\Delta }\mathit{t}\mathit{/}\mathit{\tau }\mathrm{)}\right]{\mathit{\sigma }}_{\mathit{P}}^{\mathrm{2}}$. Thus, all statistical properties of such a Markov chain are determined by the three parameters of its mean $\mathit{E}{}^{\mathrm{[}}\mathit{P}_{\mathit{x,y}}{}^{\mathrm{]}}$ (which we here assume to be zero), its variance ${\mathit{\sigma }}_{\mathit{P}}^{\mathrm{2}}$, and its relaxation time τ, i.e. the time after which the correlation function becomes “negligibly small”, which in the present case is ${\mathit{\sigma }}_{\mathit{P}}^{\mathrm{2}}\mathit{/}\mathit{e}$.

Fig. 2 Upper panel: photocentre positions extracted from the snapshots over the five years of simulation. Lower panel: photocentre positions generated with a Markov chain with σP = 0.09 AU and τ = 4 × 106 s. Points are 23 days apart in both panels.

We computed the correlation function for the sequence of photocentre positions extracted from the snapshots of the 3D model in the Gaia G band shown in the top panel of Fig. 2. The result is shown as thick lines in Fig. 3.

Fig. 3 Thick lines: correlation function of the projections of the photocentre positions extracted from the 3D model. Thin lines: correlation function of simulated photocentre positions generated with an exponentially decreasing correlation function (see Eq. (3)) with σP ≃ 0.09 AU and τ ≃ 2 − 4 × 106 s sampled every 23 days.

Because of the time step of 23 days between two successive snapshots, we are only able to evaluate the correlation function on a discrete mesh with a 23-day time step. We see in Fig. 3 that the correlation function falls to zero for Δt ≈ 60 days. Assuming the correlation function to be exponentially decreasing as in Eq. (3), τ ≃ 46 days  ≃ 4 × 10⁶ s. The variance ${\mathit{\sigma }}_{\mathit{P}}^{\mathrm{2}}$ is taken as the arithmetic mean of ${\mathit{\sigma }}_{\mathit{P,x}}^{\mathrm{2}}$ and ${\mathit{\sigma }}_{\mathit{P,y}}^{\mathrm{2}}$, the values of the correlation function at Δt = 0: σ_P ≃ 0.09 AU. To check whether the assumed exponential correlation function is compatible with the one extracted from the 3D model data, we generated photocentre positions with σ_P = 0.09 AU and τ = 2 − 4 × 10⁶ s. An example of the generated positions is shown in the lower panel of Fig. 2. We then computed the correlation function with time-steps of 23 days. The results for some series of generated positions are shown as thin lines in Fig. 3 and are in good agreement with the correlation function for the 3D model points. Furthermore, it shows that we cannot distinguish between τ = 2 × 10⁶ s and τ = 4 × 10⁶ s.

We note that the two parameters σ_P and τ depend on the surface gravity and on the effective temperature of the star (Chiavassa et al. 2011). This means that the values we derived here are only valid for stars with physical parameters similar to the ones of the 3D model. In the analysis below, we also consider stars with different σ_P and τ in order to provide predictions for a wider class of objects.

In Fig. 4, we compare the distributions of P_x and P_y extracted from the snapshots of the 3D model and the ones generated with a Markov chain.

Fig. 4 Histograms of the Px (left panel) and Py (right panel) positions extracted from the snapshots of the 3D model (red striped histograms) and the ones generated with the Markov chain (blue solid histograms). The actual Px and Py are shown in Fig. 2. The black solid line represents the Gaussian distribution used in the Markov chain to generate the photocentre positions.

The solid black line is the Gaussian distribution N(0,σ_P) from which the generated positions are extracted. The 80 measurements considered here are not numerous enough to clearly reproduce the Gaussian distribution, though it does so if we simulate more photocentre positions with the Markov chain. We could not verify whether the distribution of the positions from the 3D model becomes Gaussian for longer times or with more measurements because no additional points are available.

The next step is to add these photocentre displacements generated with the Markov chain to simulated Gaia transit data.

The stellar population considered is a sample of  ≃ 5000 red supergiants distributed around the Galactic plane (0 ≤ l ≤ 180 and  − 20 ≤ b ≤ 20 in Galactic coordinates). These stars have absolute V magnitudes between –7 and –3.5 and apparent G magnitudes between 5.6 and 20. Almost no stars are found at less than 1 kpc (for example see Fig. 5). For such close distances, supergiants tend to appear brighter than the CCD saturation limit at G = 5.6 and will not be observed. The sample was generated with the Besançon model (Robin et al. 2004) and virtually observed with the Gaia Object Generator v7.0 (GOG; Isasi et al.2010). We then added to these “observations” photocentric shifts generated with a Markov chain. In this paper, we always add to all the stars in our sample a photocentric noise generated with Markov chains of the same σ_P and τ regardless of the physical parameters of the star.

The conclusions we draw in the analysis below should be interpreted in a statistical sense, as star-by-star effects are closely related to the particular Gaia scanning law and actual photocentre configuration.

3. Impact on the single-star parameters

As mentioned in the previous section, the goal of our work is to study the impact of surface brightness asymmetries on the stellar astrometric parameters derived by Gaia.

Fig. 5 Upper panel: relative variation in the single-star parallax () resulting from the addition of surface brightness asymmetries. Lower panel: same as upper panel but with the standard deviation predicted from Eq. (5) added (red points). The parameters of the Markov chain used to generate the data are σP = 0.09 AU and τ = 4 × 106 s.

The observations, as described in the previous section, are obtained by adding photocentre motions generated with a Markov chain to Gaia transit data. They were fitted with the five-parameter single-star model (${\mathit{\alpha }}_{\mathrm{0}}^{\mathrm{\ast }}$, δ₀, ϖ, μ_α^∗, μ_δ). In Fig. 5, we show some results concerning the parallax. The upper panel of Fig. 5 illustrates the difference between the parallaxes determined with () and without () the surface brightness asymmetries, in units of the statistical error in the parallax σ_ϖ. We note that σ_ϖ will be provided as output in the Gaia catalogue and is derived from the covariance matrix of the parameters obtained in the data fitting process. As shown in Appendix A of Chiavassa et al. (2011) and of this paper, the quantity σ_ϖ is independent of the presence of surface brightness asymmetries. The effect of the photocentre displacement can be quite significant as it can induce deviations in the derived parallaxes of up to 10σ for d ≈ 1 kpc. Nevertheless, we note that all effects become negligible for stars farther out than 5–6 kpc. At this distance, the epoch astrometric precision, which is about 0.03 mas for a star brighter than G = 12, is comparable to the angular dimension of the added photocentre displacement (0.09 AU on average). A smaller photocentre noise will clearly have an impact only for more nearby stars, and the opposite for larger σ_P. For example, at σ_P = 0.02 AU, only the most nearby stars (d < 2 kpc) have a greater than 3σ_ϖ. On the other hand, at σ_P = 0.2 AU, a deviation up to 30σ is reached.

The expected value of can be easily computed (details in Appendix A) assuming that the photocentre displacements are taken from a Gaussian distribution with variance ${\mathit{\sigma }}_{\mathit{P}}^{\mathrm{2}}$. We note that with this approach the positions of the photocentre are treated as uncorrelated, which, as we have seen, is not strictly true. The mean value of is zero because the photocentre positions, on average, coincide with the geometric centre of the star. This is also true for correlated photocentre positions if enough of them are considered. On the other hand, , the standard deviation in , indicates how much individual and values can differ: ${\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\varpi ̂}}\mathrm{=}\sqrt{\mathit{N}\frac{{\mathit{\sigma }}_{{\mathit{P}}_{\mathit{\alpha }}}^{\mathrm{2}}\mathrm{+}{\mathit{\sigma }}_{{\mathit{P}}_{\mathit{\delta }}}^{\mathrm{2}}}{\mathrm{2}}}\frac{{\mathit{\sigma }}_{\mathit{P}}}{\mathit{d}}\frac{{\mathit{\sigma }}_{\mathit{\varpi }}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathit{w}}^{\mathrm{2}}}\mathit{,}$(5)where d is the distance in kpc, N is the number of measurements, σ_w is the epoch astrometric precision, ${\mathit{\sigma }}_{{\mathit{P}}_{\mathit{\alpha }}}^{\mathrm{2}}$ and ${\mathit{\sigma }}_{{\mathit{P}}_{\mathit{\delta }}}^{\mathrm{2}}$ are the variances in the parallax factors (Green 1985, and Eqs. (A.19) and (A.20)). In the lower panel of Fig. 5, we have plotted the values estimated from Eq. (5) on top of the resulting from the single-star astrometric fit. The agreement between the two sets of points is remarkably good, meaning that can be used as a realistic estimate of the impact of the surface brightness asymmetries on the accuracy of .

The Gaia end-of-mission astrometric accuracy σ_ϖ for a “perfect” single star can be written as (de Bruijne 2005) ${\mathit{\sigma }}_{\mathit{\varpi }}\mathrm{=}\mathit{m}\hspace{0.17em}{\mathit{g}}_{\mathit{\varpi }}\frac{{\mathit{\sigma }}_{\mathit{w}}}{\sqrt{\mathit{N}}}\mathrm{·}$(6)In this expression, m is an overall, end-of-mission contingency margin (m = 1.2 for a 20% margin), g_ϖ is a dimensionless geometrical factor relating to the scanning geometry, σ_w is the single measurement error (including all the error sources), and N is the number of transits. We can rewrite Eq. (5) for in terms of Eq. (6) ${\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\varpi ̂}}\mathrm{=}\frac{{\mathit{m}}^{\mathrm{2}}{\mathit{g}}_{\mathit{\varpi }}^{\mathrm{2}}}{\sqrt{\mathit{N}}}\frac{{\mathit{\sigma }}_{\mathit{P}}}{\mathit{d}}\sqrt{\frac{{\mathit{\sigma }}_{{\mathit{P}}_{\mathit{\alpha }}}^{\mathrm{2}}\mathrm{+}{\mathit{\sigma }}_{{\mathit{P}}_{\mathit{\delta }}}^{\mathrm{2}}}{\mathrm{2}}}\mathrm{\propto }\frac{{\mathit{\sigma }}_{\mathrm{\Theta }}}{\sqrt{\mathit{N}}}\mathrm{·}$(7)We see that is only related to the angular size of the photocentre displacement induced by the surface brightness asymmetries (${\mathit{\sigma }}_{\mathrm{\Theta }}\mathrm{=}\frac{{\mathit{\sigma }}_{\mathit{P}}}{\mathit{d}}$) and the scanning geometry. There is no dependence on the instrumental sensitivity σ_w because is intrinsic to the physics of the star. On the other hand, $\frac{{\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\varpi ̂}}}{{\mathit{\sigma }}_{\mathit{\varpi }}}$, the capability of Gaia to detect the noise induced by the surface brightness asymmetries, does depend on σ_w such that $\frac{{\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\varpi ̂}}}{{\mathit{\sigma }}_{\mathit{\varpi }}}\mathrm{\propto }\frac{{\mathit{\sigma }}_{\mathrm{\Theta }}}{{\mathit{\sigma }}_{\mathit{w}}}\mathrm{·}$(8)For red supergiants, σ_w is roughly constant because G < 12.6 in most of the cases.

Finally, if we consider $\frac{{\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\varpi ̂}}}{\mathit{\varpi }}$, which is the parallax inaccuracy (owing to the surface brightness asymmetries) relative to the parallax, we find that it is independent of the distance of the star $\frac{{\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\varpi ̂}}}{\mathit{\varpi }}\mathrm{\propto }{\mathit{\sigma }}_{\mathit{P}}\mathit{.}$(9)We can now compute $\frac{{\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\varpi ̂}}}{{\mathit{\sigma }}_{\mathit{\varpi }}}$ and $\frac{{\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\varpi ̂}}}{\mathit{\varpi }}$ for a red supergiant star using mean values for the scanning-law parameters g_ϖ = 1.91, m = 1.2, N = 78, σ_w = 30 μas (de Bruijne 2005; Lindegren 2010), and σ_P = 0.09 AU. For supergiants located at 1 kpc (σ_Θ = 0.09 mas), Eq. (8) predicts that $\frac{{\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\varpi ̂}}}{{\mathit{\sigma }}_{\mathit{\varpi }}}\mathit{<}\mathrm{4.3}$ for 68% of the sample and values larger than 13 for 0.3% of the sample. For supergiants at 5 kpc (σ_Θ = 0.018 mas), these numbers are 0.8 and 2.4, respectively. Similarly, Eq. (9) predicts that $\frac{{\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\varpi ̂}}}{\mathit{\varpi }}\mathit{<}\mathrm{3.4}\mathrm{\%}$ for 68% of the supergiants (independently of their distance) and values larger than 10.3% for 0.3% of the sample.

A similar behaviour is found for the four other parameters (${\mathit{\alpha }}_{\mathrm{0}}^{\mathrm{\ast }}$, δ₀, μ_α^∗, and μ_δ). Expressions equivalent to Eq. (5) are derived in Appendix A. We note that for position and proper motion, we find values that are up to 15σ off, while the photocentric noise was generated with the same Markov chain parameters as in Fig. 5.

We now analyse the accuracy (as opposed to ) of the “fitted single-star” parallax with respect to the true parallax.

Fig. 6 Histogram of the relative error in the single-star parallax (() resulting from the addition of surface brightness asymmetries for supergiant stars at different distances. The distribution for stars in the right panel (d > 5.5 kpc) closely reproduces the expected N(0,1). In the central panel (3.5 < d < 5.5 kpc), the distribution is wider with a best-fit N(0,2.5). Finally, for the left panel (d < 3.5 kpc) the σ of the distribution increases to four. Parameters of the Markov chain used here are σP = 0.09 AU and τ = 4 × 106 s. The x-axes in the three panels are at the same scale.

In the three panels of Fig. 6, we plot this difference, normalised by the estimated error in the parallax σ_ϖ, for different slices in distance. For stars far away (right panel), the distribution is clearly N(0,1) as expected. As we consider more nearby stars, the distribution becomes wider and wider. For stars in the range 3.5 < d < 5.5 kpc (central panel), the histogram is approximated well by a N(0,2.5) distribution, and for stars with d < 3.5 kpc (left panel) by N(0,4). If we increase σ_P, we obtain wider and wider distributions. For σ_P = 0.2 AU, only stars between 7.5 kpc and 8 kpc follow a N(0,1) distribution. On the other hand, lowering σ_P produces narrower distributions. For σ_P = 0.02 AU, the distribution becomes N(0,1) for stars at all distances indicating that at this level of photocentric noise, no effects are expected in Gaia data processing. We recall that actual stars have σ_P, which is dependent upon their absolute magnitude, but the exact relationship is not known so far (Chiavassa et al. 2011).

As already mentioned, we have used a five-parameter-model for our set of supergiants because we know we are dealing with genuine single stars. However, the F2 values for these solutions (Fig. 7, top panel) can exceed the threshold value three by more than one order of magnitude. An estimate of the increase in F2 induced by the presence of surface brightness asymmetries can be easily obtained by noticing that ${\mathit{\chi }}_{\mathit{a}}^{\mathrm{2}}$, which is the χ² in the presence of surface brightness asymmetries, can be written as (Chiavassa et al. 2011, Eqs. (11) and (12), and Appendix B of this paper) ${\mathit{\chi }}_{\mathit{a}}^{\mathrm{2}}\mathrm{\simeq }{\mathit{\chi }}_{\mathrm{na}}^{\mathrm{2}}\mathrm{+}\mathit{N}{\left(\frac{{\mathit{\sigma }}_{\mathit{P}}}{\mathit{d}{\mathit{\sigma }}_{\mathit{w}}}\right)}^{\mathrm{2}}\mathit{,}$(10)where ${\mathit{\chi }}_{\mathrm{na}}^{\mathrm{2}}$ is the χ² of the fit of the data without surface brightness asymmetries, N is the number of measurements, d is the distance in kpc, and σ_w is the Gaia single-measurement error. This expression can be easily converted into the analogous expression for F2 $\mathit{F}{\mathrm{2}}_{\mathit{a}}\mathrm{=}\mathit{F}{\mathrm{2}}_{\mathrm{na}}\mathrm{+}{\left(\frac{\mathrm{9}\mathit{\nu }}{\mathrm{2}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\left[{\left(\mathrm{1}\mathrm{+}\frac{{\mathit{\sigma }}_{\mathit{P}}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathit{w}}^{\mathrm{2}}}\right)}^{\mathrm{1}\mathit{/}\mathrm{3}}\mathrm{-}\mathrm{1}\right]\mathit{.}$(11)This quantity is plotted in the lower panel of Fig. 7 on top of the results obtained for the single-star solution.

Fig. 7 Upper panel: F2 for the single-star fit of the simulated data with the addition of surface brightness asymmetries. A F2 > 3 indicates a non-satisfactory fit. Lower panel: the same as the upper panel but with the theoretical prediction (Eq. (11)) added on top (see text).

There is on average remarkable agreement. As we previously pointed out, a case by case agreement is not to be expected. When we change the parameters of the Markov chain and set σ_P = 0.2 AU, F2_a increases to 100 for more nearby stars. No stars at d < 2.5 kpc have F2_a < 3, and we find values larger than three at all distances. In contrast if we look at what happens with σ_P = 0.02 AU, F2_a becomes larger than three only for a negligible fraction of the stars, and it never exceeds seven.

4. Non-single-star processing

We have so far treated the supergiants in our sample as single stars, as we know that they actually are. Nevertheless, when processing real Gaia data, we will not know in advance the kind of object we will be analysing. When we do not restrict the data processing flow to the five-parameter solution, stars with a single-star F2 larger than three are further processed. A chain of solutions of increasing complexity is applied until F2 becomes smaller than three. In the acceleration solution, a possible variation of the proper motion is added to better fit the data. In the orbital solution, the parameters of the orbit of a binary system are added. If neither the single-star solution nor the non-single-star solutions (acceleration or orbital) are acceptable on the basis of their F2, a stochastic solution is computed: some noise is added to the along-scan measurements in order to lower to zero the F2 of the single-star solution.

For our sample, an acceleration solution, an orbital solution, or a stochastic solution are found for some fraction of the stars, and this fraction depends on the parameters of the Markov chain. The distribution of solutions among these possible classes is plotted in Fig. 8 for three different sets of σ_P and τ.

Fig. 8 Distribution of the different types of solutions found as a function of the distance. In the three panels are shown the results for the same sample of stars to which different photocentric motions have been applied, as indicated by the corresponding σP and τ in the figure label.

The number of single-star solutions decreases drastically as we increase σ_P. Moreover, the impact on the astrometric parameters then increases, implying that the number of solutions in which they are unaltered decreases. We see as well that all the nearby (d ≲ 4 kpc) single-star solutions disappear. At the same time, we see an order-of-magnitude increase in the number of stochastic solutions. Finally, we point out the decrease of the fraction of orbital solutions as we decrease τ. This is a consequence of the loss of correlation between temporally close measurements that makes it harder to fit the points with an orbit.

That single-star solutions are rejected by the Gaia pipeline even though we are dealing with genuine single stars is unfortunate. Hence, we analysed these solutions to see whether it is possible to identify them and filter them out. We focused our attention on the orbital solutions as they can represent a significant fraction of the solutions for our sample.

As a first step, we analysed the significance of the orbital solution by looking at the ratio a/σ_a, where a is the semi-major axis of the orbit. If a is smaller or comparable to its error σ_a, the solution is not significant as this is an indication of null orbital parameters. In the left panel of Fig. 9, we show the ratio a/σ_a for a Markov chain with σ_P = 0.2 AU and τ = 4 × 10⁶ s and conclude that a large fraction of these solutions has a/σ_a < 1. They are thus spurious and may easily be flagged as such. For these solutions, the reduction in F2 under the acceptance value of 3 is caused by the χ² decreasing faster than the number of parameters increases.

Fig. 9 Left panel: distribution of the ratio a/σa, with a the semi-major axis of the orbit and σa its uncertainty. Solutions with a/σa < 2 are not significant and may be rejected. Right panel: ratio a/R, with R the stellar radius for stars with a/σa > 2. All orbits appear to have a semi-major axis smaller than the stellar radius and can thus be rejected. Markov chain parameters: σP = 0.2 AU and τ = 4 × 106 s.

As a second step, we studied the stars for which the orbital solution is significant (i.e. located in the peak at about a/σ_a = 6 in the left panel of Fig. 9) and we computed the ratio a/R, where R is the radius of the star. This quantity is plotted in the right panel of Fig. 9 and is always smaller than unity. We may thus reject all those solutions as they are unphysical, representing orbits with separations smaller than the radius of the stars under study. Moreover, we note that the size in AU of the semi-major axis found (Fig. 10) is always of the same order of magnitude as the photocentre motion added with the Markov chain, 0.2 AU being about the minimum value found.

Fig. 10 Semi-major axes a of the orbits derived with the orbital solution. The physical size of a in AU is comparable to the photocentre motion added to the simulated Gaia data. Markov chain parameters: σP = 0.2 AU and τ = 4 × 106 s.

Finally, we discuss the stochastic solutions. Stochastic noise is meant to represent any source of “noise” found in the data that cannot be modelled with any of the other models. We expect it to be correlated with the noise we added to the single-star Gaia data. A histogram of the stochastic noise (in AU) is shown in Fig. 11 for three different sets of parameters of the Markov chain. We show examples with extreme parameters in order to more clearly illustrate the occurring effects. In the left panel, we show the results for simulated photocentric motions with σ_P = 0.2 AU and τ = 4 × 10⁶ s. The blue horizontally striped histogram represents the distribution of the stochastic noise obtained by the data processing. It is clearly neither symmetric nor centred at the value of σ_P = 0.2 AU as expected (vertical dashed line). This is because in the fitting process, the position of the star is assumed to be at the mean photocentre position, which does not necessarily coincide with the geometrical centre of the star (example in the top panel of Fig. 2, where the mean photocentre position is in the first quadrant rather than at the origin). Once we quadratically add to the stochastic noise the difference between the mean photocentre position and the true position of the star, the resulting distribution (green diagonally striped histogram) becomes centred at 0.2 AU as expected. This positional shift and the asymmetry of the distribution is caused by the length of the Gaia mission (5 years) not being long enough for the distribution of photocentre positions to be Gaussian. Moreover, that measurements are taken at time intervals shorter than the relaxation time of about 40 days and the non-uniformity of the time sampling play a crucial role. The minimum time between two consecutive Gaia observations of the same target is about two hours. When we decrease τ to 4 × 10³ s (i.e. about one hour, meaning that all measured photocentre positions are almost uncorrelated), the histogram becomes almost symmetric and centred on 0.2 AU (central panel of Fig. 11). If we keep τ = 4 × 10⁶ s and decrease σ_P to 0.08 AU (right panel of Fig. 11), the stochastic noise distribution is still asymmetric and not peaked at 0.08 AU. However, the effects are mainly smaller than for the 0.2 AU case, merely because of the smaller size of the photocentric shifts generated.

Fig. 11 The blue horizontally striped histograms show the “stochastic” noise added to Gaia measurements in the computation of a stochastic solution. The green diagonally striped ones represent the stochastic noise plus the offset of the mean photocentre position with respect to the centre of the star. The three panels show results for three different sets of parameters of the Markov chain. The vertical dashed line corresponds to the standard deviation σP in the photocentre displacement.

From this discussion, we can conclude that the stochastic noise can be used as an estimate of the photocentric noise present on the surface of the star. This correspondence is at its closest for short correlation times and small photocentric motions. Under these circumstances, the distribution will be more symmetric and will peak at the correct value.

5. Conclusions

Surface brightness asymmetries are a common characteristic of stars and cause time-varying positional shifts between the photocentre and the barycentre. We have analysed their effects on the parameters derived from Gaia simulated data for supergiants. To achieve this goal, on the basis of a 3D model of a supergiant, we have simulated photocentric displacements to be added to along-scan measurements. We have then performed a five-parameter fit. For such a single-star solution, we have shown that the deduced parallaxes can be up to 13σ_ϖ off for a supergiant at 1 kpc with σ_P = 0.09 AU and $\frac{{\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\varpi ̂}}}{\mathit{\varpi }}$ up to 10%. We have shown as well that $\frac{{\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\varpi ̂}}}{{\mathit{\sigma }}_{\mathit{\varpi }}}$ decreases with distance as $\frac{{\mathit{\sigma }}_{\mathit{P}}}{\mathrm{d}{\mathit{\sigma }}_{\mathrm{w}}}$, but in contrast, $\frac{{\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\varpi ̂}}}{\mathit{\varpi }}$ does not vary with distance, being proportional to σ_P. At the same time, the quality of the fit can be very poor, with F2 parameters exceeding the threshold value of 3 by more than one order of magnitude. If more complex solutions are searched for, even though F2 becomes smaller than three, there is no improvement in the accuracy of the derived parameters and, at the same time, the computational effort is substantially increased. Spurious binary solutions can be rejected on the basis of significance tests and the comparison of the size of the semi-major axis with the radius of the star. When a stochastic solution is found, we observed a correlation between the stochastic noise and the Gaussian noise injected in the Markov chain to simulate the photocentric motion. This information could be exploited to obtain, at least, the order of magnitude of the convective noise present.

From the standpoint of the Gaia data processing, if a star could be identified as a supergiant independently of its parallax, a “supergiant flag” may then be used as a warning that the astrometric parameters may be far off. Classification methods independent of the parallax include the spectroscopic determination of the atmospheric parameters (T_eff,log g) from the Gaia spectrometer, or the classification as an irregular variable from the repeated photometric measurements. Irregular variables exist in two flavours, however: supergiants and stars at the base of the red giant branch (Jorissen et al. 1997). The two classes of irregular variables differ by several magnitudes in their respective absolute magnitudes, so that a 10% relative error in the parallax of the supergiant caused by the surface brightness asymmetries is not large enough to intermingle the two classes and introduce confusion.

To identify true binaries among supergiants despite the confusion introduced by the surface brightness asymmetries, an indirect way is provided by the detection of composite colour indices (in addition to the orbital separation probably now being significantly larger than the stellar radius). Such composite colour indices for systems involving supergiants were computed by Malkov et al. (2011), with the conclusion that these composite indices differ by at least 0.3 mag from the nearest single-star indices in the best two-colour diagram, thus allowing the identification of the binary systems.

This complementary non-astrometric information should thus allow for a more efficient and reliable evaluation of the significance of the astrometric solution in the presence of surface brightness asymmetries.

Acknowledgments

E.P. is supported by the ELSA (European Leadership in Space Astrometry) Research Training Network of the FP6 Programme and PRODEX. D.P. is FNRS Senior Research Associate. We thank L. Lindegren, B. Holl, A. Chiavassa, and C. Babusiaux for their support and useful discussions. Simulated data provided by the Simulation Unit (CU2) of the Gaia Data Processing Analysis Consortium (DPAC) have been used to complete this work. They are gratefully acknowledged for this contribution.

References

Appendix A: Estimation of the accuracy on the astrometric parameters

Gaia will provide N measurements (w) of the along-scan position of each star and hence the astrometric solution will be computed by solving an overdetermined system of N equations for β_na, the n unknowns (${\mathit{\alpha }}_{\mathrm{0}}^{\mathrm{\ast }}$, δ₀, ϖ, μ_α^∗, μ_δ in the case of a single-star solution discussed here): $\mathit{X}{{\beta }}_{\mathbf{na}}\mathrm{+}{ϵ}\mathrm{=}{w}\mathit{,}$(A.1)where ϵ is the vector of errors that we assume to be random and mutually uncorrelated, i.e. they have zero mean and their covariance matrix W is diagonal. We define X to be the matrix of known coefficients $\mathit{X}\mathrm{=}\left(\begin{array}{c}\\ \sin {\mathit{\theta }}_{\mathrm{1}}& \cos {\mathit{\theta }}_{\mathrm{1}}& {\mathit{f}}_{{\mathit{w}}_{\mathrm{1}}}& {\mathit{t}}_{\mathrm{1}}\sin {\mathit{\theta }}_{\mathrm{1}}& {\mathit{t}}_{\mathrm{1}}\cos {\mathit{\theta }}_{\mathrm{1}}\\ \sin {\mathit{\theta }}_{\mathrm{2}}& \cos {\mathit{\theta }}_{\mathrm{2}}& {\mathit{f}}_{{\mathit{w}}_{\mathrm{2}}}& {\mathit{t}}_{\mathrm{2}}\sin {\mathit{\theta }}_{\mathrm{2}}& {\mathit{t}}_{\mathrm{2}}\cos {\mathit{\theta }}_{\mathrm{2}}\\ \begin{array}{}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}& \begin{array}{}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}& \begin{array}{}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}& \begin{array}{}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}& \begin{array}{}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}\\ \sin {\mathit{\theta }}_{\mathit{N}}& \cos {\mathit{\theta }}_{\mathit{N}}& {\mathit{f}}_{{\mathit{w}}_{\mathit{N}}}& {\mathit{t}}_{\mathit{N}}\sin {\mathit{\theta }}_{\mathit{N}}& {\mathit{t}}_{\mathit{N}}\cos {\mathit{\theta }}_{\mathit{N}}\end{array}\right)\mathit{,}$(A.2)where θ is the position angle of the scan, t is the observing time, and f_w = r(P_αsinθ + P_δcosθ) is a linear combination of the parallax factors P_α and P_δ (Green 1985), and r is the radius vector of Gaia’s orbit at time t in AU ${\mathit{P}}_{\mathit{\alpha }}\mathrm{=}\cos \mathit{ϵ}\cos \mathit{\alpha }\sin \mathrm{\odot }\mathrm{-}\sin \mathit{\alpha }\cos \mathrm{\odot }\mathit{,}$(A.3)

${\mathit{P}}_{\mathit{\delta }}\mathrm{=}\mathrm{(}\sin \mathit{ϵ}\cos \mathit{\delta }\mathrm{-}\cos \mathit{ϵ}\sin \mathit{\alpha }\sin \mathit{\delta }\mathrm{)}\sin \mathrm{\odot }\mathrm{-}\cos \mathit{\alpha }\sin \mathit{\delta }\cos \mathrm{\odot }\mathit{,}$where  ⊙  and ϵ are, respectively, the longitude of the Sun and the obliquity of the ecliptic, both at time t.

We evaluate ${\beta ̂}{\mathbf{na}}_{}$, the best estimate of the parameters β_na in Eq. (A.1), by minimizing χ²${\mathit{\chi }}_{\mathrm{na}}^{\mathrm{2}}\mathrm{=}\mathrm{(}{w}\mathrm{-}\mathit{X}{{\beta }}_{\mathbf{na}}{\mathrm{)}}^{\mathit{T}}{\mathit{W}}^{-1}\mathrm{(}{w}\mathrm{-}\mathit{X}{{\beta }}_{\mathbf{na}}\mathrm{)}\mathit{,}$(A.4)which can be expressed by ${\beta ̂}\mathbf{na}\mathrm{=}{\mathrm{(}{\mathit{X}}^{\mathit{T}}{\mathit{W}}^{-1}{\mathit{X}}^{\mathrm{)}}}^{-1}\mathrm{\left(}{\mathit{X}}^{\mathit{T}}{\mathit{W}}^{-1}{w}\mathrm{\right)}\mathit{.}$(A.5)The difference between the true parameters and the estimated ones is easily obtained using Eqs. (A.5) and (A.1) ${\beta ̂}\mathbf{na}\mathrm{-}{{\beta }}_{\mathbf{na}}\mathrm{=}{\mathrm{(}{\mathit{X}}^{\mathit{T}}{\mathit{W}}^{-1}{\mathit{X}}^{\mathrm{)}}}^{-1}{\mathit{X}}^{\mathit{T}}{\mathit{W}}^{-1}{ϵ}\mathit{,}$(A.6)as well as Z, the covariance matrix of the parameters $\mathit{Z}\mathrm{=}\mathit{E}\mathrm{\left[}\left({\beta ̂}\mathbf{na}\mathrm{-}\mathrm{\left(}{{\beta }}_{\mathbf{na}}\mathrm{\right)}\right){\left({\beta ̂}\mathbf{na}\mathrm{-}\mathrm{\left(}{{\beta }}_{\mathbf{na}}\mathrm{\right)}\right)}^{\mathit{T}}\mathrm{\right]}\mathrm{=}{\mathrm{(}{\mathit{X}}^{\mathit{T}}{\mathit{W}}^{-1}{\mathit{X}}^{\mathrm{)}}}^{-1}\mathit{.}$(A.7)We note that Z is diagonal because all the off-diagonal terms contain sinθ, cosθ, or t factors that averages to zero. The mean parallax factors are also zero because of their dependence on  ⊙ .

We now consider a star with time-varying surface brightness asymmetries. Their presence causes a shift in the photocentre position and a consequent shift (Δw) in the observed along-scan position with respect to a star with a uniform disk.

This can be expressed as $\mathrm{\Delta }{\mathit{w}}_{\mathit{i}}\mathrm{=}\frac{\mathrm{1}}{\mathit{d}}\left[{\mathit{a}}_{\mathit{i}}\cos {\mathit{\theta }}_{\mathit{i}}\mathrm{+}{\mathit{b}}_{\mathit{i}}\sin {\mathit{\theta }}_{\mathit{i}}\right]\mathrm{=}\frac{\mathrm{1}}{\mathit{d}}{\mathit{D}}_{\mathit{i}}\mathit{,}$(A.8)where (a_i,b_i) are the coordinates of the photocentre in a reference system with its origin at the centre of the star (expressed in AU) and d is the distance of the star in kpc. We assume that a_i and b_i follow a Gaussian distribution with zero mean and standard deviation σ_P.

The astrometric problem to be solved now to get the parameters β_a is analogous to Eq. (A.1) $\mathit{X}{\beta }\mathit{a}\mathrm{+}{ϵ}\mathrm{=}{w}\mathrm{+}\mathbf{\Delta }{w}\mathit{.}$(A.9)The estimate of the parameters is again found by χ² minimization ${\mathit{\chi }}_{\mathit{a}}\mathrm{2}\mathrm{=}\mathrm{\left[}\mathrm{\right(}{w}\mathrm{+}\mathbf{\Delta }{w}\mathrm{)}\mathrm{-}\mathit{X}{\beta }\mathit{a}{\mathrm{]}}^{\mathit{T}}{\mathit{W}}^{-1}\mathrm{[}\mathrm{(}{w}\mathrm{+}\mathbf{\Delta }{w}\mathrm{)}\mathrm{-}\mathit{X}{\beta }\mathit{a}\mathrm{]}$(A.10)and the result, in analogy with Eq. (A.5), is ${\beta ̂}\mathit{a}\mathrm{=}{\mathrm{(}{\mathit{X}}^{\mathit{T}}{\mathit{W}}^{-1}{\mathit{X}}^{\mathrm{)}}}^{-1}\mathrm{\left(}{\mathit{X}}^{\mathit{T}}{\mathit{W}}^{-1}\mathrm{\right(}{w}\mathrm{+}\mathbf{\Delta }{w}{\mathrm{\left)}}^{\mathrm{\right)}}\mathit{.}$(A.11)With the use of Eqs. (A.5) and (A.11), we can produce an expression for the effect of the surface brightness asymmetries on the parameter estimates ${\beta ̂}\mathit{a}\mathrm{-}{\beta ̂}\mathbf{na}\mathrm{=}{\mathrm{(}{\mathit{X}}^{\mathit{T}}{\mathit{W}}^{-1}{\mathit{X}}^{\mathrm{)}}}^{-1}\mathrm{\left(}{\mathit{X}}^{\mathit{T}}{\mathit{W}}^{-1}\mathbf{\Delta }{w}\mathrm{\right)}\mathrm{=}\mathit{Z}\mathrm{\left(}{\mathit{X}}^{\mathit{T}}{\mathit{W}}^{-1}\mathbf{\Delta }{w}\mathrm{\right)}\mathit{.}$(A.12)For a non-variable star, the elements of Z and σ_wi, will be approximately all the same as they are a function of the apparent magnitude G and we can rewrite Eqs. (A.5), (A.7), and (A.12) for the derived parameters and the covariance matrix, as $\begin{array}{ccc}& & {\beta ̂}\mathbf{na}\mathrm{=}{\mathrm{(}{\mathit{X}}^{\mathit{T}}{\mathit{X}}^{\mathrm{)}}}^{-1}\mathrm{(}{\mathit{X}}^{\mathit{T}}{w}\mathrm{)}\mathit{,}\\ & & {\beta ̂}\mathit{a}\mathrm{=}{\mathrm{(}{\mathit{X}}^{\mathit{T}}{\mathit{X}}^{\mathrm{)}}}^{-1}\mathrm{(}{\mathit{X}}^{\mathit{T}}\mathrm{(}{w}\mathrm{+}\mathbf{\Delta }{w}{\mathrm{)}}^{\mathrm{)}}\mathit{,}\\ & & \mathit{Z}\mathrm{=}{\mathit{\sigma }}_{\mathit{w}}^{\mathrm{2}}{\mathrm{(}{\mathit{X}}^{\mathit{T}}{\mathit{X}}^{\mathrm{)}}}^{-1}\mathit{.}\end{array}$We note that the covariance matrix Z is independent of the presence of the brightness asymmetries, which means that the uncertainty in the deduced parameters will be the same. That is because what enters the covariance matrix are basically the measurement errors.

In reality, the accuracy (to be distinguished from the precision) of the parameters can be quite different in the two cases. Indeed, the parameters derived in presence of brightness asymmetries will be different from the ones found for a uniform-disk star, their difference being given by: $\mathbf{\Delta }{\beta ̂}\mathrm{=}{\beta ̂}\mathit{a}\mathrm{-}{\beta ̂}\mathbf{na}\mathrm{=}\mathit{Z}\mathrm{\left(}{\mathit{X}}^{\mathit{T}}{\mathit{W}}^{-1}\mathbf{\Delta }{w}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{{\mathit{\sigma }}_{\mathit{w}}^{\mathrm{2}}}\frac{\mathrm{1}}{\mathit{d}}\left(\begin{array}{c}\\ {\mathit{\sigma }}_{{\mathit{\beta }}_{\mathrm{1}}}^{\mathrm{2}}{\sum }_{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{N}}{\mathit{X}}_{\mathit{i}\mathrm{1}}{\mathit{D}}_{\mathit{i}}\\ \begin{array}{}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}\\ {\mathit{\sigma }}_{{\mathit{\beta }}_{\mathrm{5}}}^{\mathrm{2}}{\sum }_{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{N}}{\mathit{X}}_{\mathit{i}\mathrm{5}}{\mathit{D}}_{\mathit{i}}\end{array}\right)\mathit{,}$(A.16)where ${\mathit{\sigma }}_{{\mathit{\beta }}_{\mathit{j}}}^{\mathrm{2}}$ are the diagonal elements of the covariance matrix Z, X is the matrix of the known coefficients described in Eq. (A.2), and D describes the photocentre motion (see Eq. (A.8)). In this equation, all the terms are known apart from ${\sum }_{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{N}}{\mathit{X}}_{\mathit{ij}}{\mathit{D}}_{\mathit{i}}$. These terms are on average zero because the factors are all uncorrelated and  ⟨ D_i ⟩  = 0. As we expect to see some effects on the deduced parameters because of the presence of surface brightness asymmetries, it is interesting to compute ${\mathit{\sigma }}_{\mathbf{\Delta }{\beta ̂}}^{\mathrm{2}}$, the standard deviation of $\widehat{{{\beta }}_{{a}}}\mathrm{-}{\beta ̂}{}_{\mathbf{na}}$. This quantity gives an estimate of the errors in the value of the derived parameters caused by the presence of the brightness asymmetries. For the parameter β_j, this means that we need to estimate the variance in N ⟨ X_ijD_i ⟩ _i. We can easily compute those factors for all the parameters where  ⟨ a ⟩  =  ⟨ b ⟩  = 0, ${\mathit{\sigma }}_{\mathit{a}}^{\mathrm{2}}\mathrm{=}{\mathit{\sigma }}_{\mathit{b}}^{\mathrm{2}}\mathrm{=}{\mathit{\sigma }}_{\mathit{P}}^{\mathrm{2}}$ and  ⟨ t ⟩  = 0.

For the first two parameters ${\mathit{\alpha }}_{\mathrm{0}}^{\mathrm{\ast }}$ and δ, we have to compute the variance of  ⟨ X_ijD_i ⟩ _i = (sinθ)D and  ⟨ X_ijD_i ⟩ _i = (cosθ)D, which leads to ${\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\alpha ̂}\begin{array}{}\mathrm{\ast }\\ \mathrm{0}\end{array}}\mathrm{=}{\mathit{\sigma }}_{\mathit{\alpha ̂}\begin{array}{}\mathrm{\ast }\\ \mathrm{0}\mathrm{-}\mathit{a}\end{array}\mathrm{-}\mathit{\alpha ̂}\begin{array}{}\mathrm{\ast }\\ \mathrm{0}\mathrm{-}\mathit{na}\end{array}}\mathrm{=}\sqrt{\frac{\mathit{N}}{\mathrm{2}}}\frac{{\mathit{\sigma }}_{{\mathit{\alpha }}_{\mathrm{0}}^{\mathrm{\ast }}}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathit{w}}^{\mathrm{2}}}\frac{{\mathit{\sigma }}_{\mathit{P}}}{\mathit{d}}$(A.17)and ${\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\delta ̂}\mathrm{0}}\mathrm{=}{\mathit{\sigma }}_{\mathit{\delta ̂}\mathrm{0}\mathrm{-}\mathit{a}\mathrm{-}\mathit{\delta ̂}\mathrm{0}\mathrm{-}\mathit{na}}\mathrm{=}\sqrt{\frac{\mathit{N}}{\mathrm{2}}}\frac{{\mathit{\sigma }}_{{\mathit{\delta }}_{\mathrm{0}}}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathit{w}}^{\mathrm{2}}}\frac{{\mathit{\sigma }}_{\mathit{P}}}{\mathit{d}}\mathrm{·}$(A.18)For the parallax  ⟨ X_ijD_i ⟩ _i = f_wD, we need the variance in the parallax factors, which is easily obtained by noting that the only term that is not constant in time is  ⊙  and neglecting the variations in r: $\begin{array}{ccc}& & {\mathit{\sigma }}_{{\mathit{P}}_{\mathit{\alpha }}}^{\mathrm{2}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}{\cos }^{\mathrm{2}}\mathit{ϵ}{\cos }^{\mathrm{2}}\mathit{\alpha }\mathrm{+}{\sin }^{\mathrm{2}}{\mathit{\alpha }}^{\mathrm{)}}\mathit{,}\\ & & {\mathit{\sigma }}_{{\mathit{P}}_{\mathit{\delta }}}^{\mathrm{2}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{[}{\left(\sin \mathit{ϵ}\cos \mathit{\delta }\mathrm{-}\cos \mathit{ϵ}\sin \mathit{\alpha }\sin \mathit{\delta }\right)}^{\mathrm{2}}\mathrm{+}{\cos }^{\mathrm{2}}\mathit{\alpha }{\sin }^{\mathrm{2}}{\mathit{\delta }}^{\mathrm{]}}\mathit{.} \end{array}$Hence ${\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\varpi ̂}}\mathrm{=}{\mathit{\sigma }}_{\mathit{\varpi ̂}\mathit{a}\mathrm{-}\mathit{\varpi ̂}\mathrm{na}}\mathrm{=}\sqrt{{\mathit{\sigma }}_{{\mathit{P}}_{\mathit{\alpha }}}^{\mathrm{2}}\mathrm{+}{\mathit{\sigma }}_{{\mathit{P}}_{\mathit{\delta }}}^{\mathrm{2}}}\sqrt{\frac{\mathit{N}}{\mathrm{2}}}\frac{{\mathit{\sigma }}_{\mathit{\varpi }}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathit{w}}^{\mathrm{2}}}\frac{{\mathit{\sigma }}_{\mathit{P}}}{\mathit{d}}\mathrm{·}$(A.21)Finally for the two components of the proper motion μ_α^∗ and μ_δ, we have to compute the variance in  ⟨ X_ijD_i ⟩ _i = t(sinθ)D and  ⟨ X_ijD_i ⟩ _i = t(cosθ)D, which leads to ${\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\mu ̂}{\mathit{\alpha }}^{\mathrm{\ast }}}\mathrm{=}{\mathit{\sigma }}_{\mathit{\mu ̂}{\mathit{\alpha }}^{\mathrm{\ast }}\mathrm{-}\mathit{a}\mathrm{-}\mathit{\mu ̂}{\mathit{\alpha }}^{\mathrm{\ast }}\mathrm{-}\mathit{na}}\mathrm{=}\sqrt{{\mathit{\sigma }}_{\mathit{t}}^{\mathrm{2}}}\sqrt{\frac{\mathit{N}}{\mathrm{2}}}\frac{{\mathit{\sigma }}_{{\mathit{\mu }}_{{\mathit{\alpha }}^{\mathrm{\ast }}}}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathit{w}}^{\mathrm{2}}}\frac{{\mathit{\sigma }}_{\mathit{P}}}{\mathit{d}}$(A.22)and ${\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\mu ̂}\mathit{\delta }}\mathrm{=}{\mathit{\sigma }}_{\mathit{\mu ̂}\mathit{\delta }\mathrm{-}\mathit{a}\mathrm{-}\mathit{\mu ̂}\mathit{\delta }\mathrm{-}\mathit{na}}\mathrm{=}\sqrt{{\mathit{\sigma }}_{\mathit{t}}^{\mathrm{2}}}\sqrt{\frac{\mathit{N}}{\mathrm{2}}}\frac{{\mathit{\sigma }}_{{\mathit{\mu }}_{\mathit{\delta }}}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathit{w}}^{\mathrm{2}}}\frac{{\mathit{\sigma }}_{\mathit{P}}}{\mathit{d}}\mathrm{·}$(A.23)

Appendix B: Expression for χ² in the presence of surface brightness asymmetries

We estimate the increase in the χ² caused by the presence of surface brightness asymmetries. We can rewrite Eq. (A.10) for ${\mathit{\chi }}_{\mathit{a}}^{\mathrm{2}}$ in a slightly different form ${\mathit{\chi }}_{\mathit{a}}\mathrm{2}\mathrm{=}\frac{\mathrm{1}}{{\mathit{\sigma }}_{\mathit{w}}^{\mathrm{2}}}\mathrm{[}{w}\mathrm{-}\mathit{X}{{\beta }}_{\mathbf{na}}\mathrm{+}\mathbf{\Delta }{w}\mathrm{-}\mathit{X}\mathbf{\Delta }{\beta }\mathrm{)}{\mathrm{]}}^{\mathit{T}}\mathrm{[}{w}\mathrm{-}\mathit{X}{{\beta }}_{\mathbf{na}}\mathrm{+}\mathbf{\Delta }{w}\mathrm{-}\mathit{X}\mathbf{\Delta }{\beta }\mathrm{)}\mathrm{]}$(B.1)\newpagewhere Δβ is β_a − β_na. The minimum ${\mathit{\chi }}_{\mathit{a}}^{\mathrm{2}}$ will be, by definition, the one for the deduced parameters, hence with (Eq. (A.16)) and with (Eq. (A.5)). Using those values, we obtain

${\mathit{\chi }}_{\mathit{a}}\mathrm{2}\mathrm{=}\frac{\mathrm{1}}{{\mathit{\sigma }}_{\mathit{w}}^{\mathrm{2}}}\mathrm{[}{\left({w}\mathrm{-}\mathit{X}{\beta ̂}\mathbf{na}\right)}^{\mathit{T}}\left({w}\mathrm{-}\mathit{X}{\beta ̂}\mathbf{na}\right)\mathrm{+}\mathbf{\Delta }{w}\mathit{T}\mathbf{\Delta }{w}\mathrm{+}{\left(\mathit{X}\mathbf{\Delta }{\beta ̂}\right)}^{\mathit{T}}\left(\mathit{X}\mathbf{\Delta }{\beta ̂}\right){\mathrm{+}}^{}$$\mathrm{-}\mathrm{2}\mathbf{\Delta }{w}\mathit{T}\left(\mathit{X}\mathbf{\Delta }{\beta ̂}\right)\mathrm{+}\mathrm{2}{\left({w}\mathrm{-}\mathit{X}{\beta ̂}\mathbf{na}\right)}^{\mathit{T}}\mathbf{\Delta }{w}\mathrm{-}\mathrm{2}{\left({w}\mathrm{-}\mathit{X}{\beta ̂}\mathbf{na}\right)}^{\mathit{T}}{\left(\mathit{X}\mathbf{\Delta }{\beta ̂}\right)}^{\mathrm{]}}\mathit{.}$(B.2)The last four terms are negligible with respect to the first two, leading to ${\mathit{\chi }}_{\mathit{a}}^{\mathrm{2}}\mathrm{\simeq }{\mathit{\chi }}_{\mathbf{na}}^{\mathrm{2}}\mathrm{+}\mathit{N}\frac{\mathrm{1}}{{\mathit{\sigma }}_{\mathit{w}}^{\mathrm{2}}}{\left(\frac{{\mathit{\sigma }}_{\mathit{P}}}{\mathit{d}}\right)}^{\mathrm{2}}\mathrm{·}$(B.3)

All Figures

Fig. 1 Snapshot of the 3D supergiant model in the Gaia G band from Chiavassa et al. (2009). The colour scale denotes different temperatures. The 12 M⊙ star is modelled with an equidistant numerical mesh of 2353 grid points and a resolution of 8.6 R⊙ (or 0.040 AU). The luminosity average over spherical shells and over time (i.e., over 5 years) is L = 93   000 ± 1300 L⊙, Teff = 3490 ± 13 K, R = 832 ± 0.7 R⊙, and log g =  −0.337 ± 0.001.

In the text

	Fig. 2 Upper panel: photocentre positions extracted from the snapshots over the five years of simulation. Lower panel: photocentre positions generated with a Markov chain with σP = 0.09 AU and τ = 4 × 106 s. Points are 23 days apart in both panels.
In the text

	Fig. 3 Thick lines: correlation function of the projections of the photocentre positions extracted from the 3D model. Thin lines: correlation function of simulated photocentre positions generated with an exponentially decreasing correlation function (see Eq. (3)) with σP ≃ 0.09 AU and τ ≃ 2 − 4 × 106 s sampled every 23 days.
In the text

	Fig. 4 Histograms of the Px (left panel) and Py (right panel) positions extracted from the snapshots of the 3D model (red striped histograms) and the ones generated with the Markov chain (blue solid histograms). The actual Px and Py are shown in Fig. 2. The black solid line represents the Gaussian distribution used in the Markov chain to generate the photocentre positions.
In the text

	Fig. 5 Upper panel: relative variation in the single-star parallax () resulting from the addition of surface brightness asymmetries. Lower panel: same as upper panel but with the standard deviation predicted from Eq. (5) added (red points). The parameters of the Markov chain used to generate the data are σP = 0.09 AU and τ = 4 × 106 s.
In the text

Fig. 6 Histogram of the relative error in the single-star parallax (() resulting from the addition of surface brightness asymmetries for supergiant stars at different distances. The distribution for stars in the right panel (d > 5.5 kpc) closely reproduces the expected N(0,1). In the central panel (3.5 < d < 5.5 kpc), the distribution is wider with a best-fit N(0,2.5). Finally, for the left panel (d < 3.5 kpc) the σ of the distribution increases to four. Parameters of the Markov chain used here are σP = 0.09 AU and τ = 4 × 106 s. The x-axes in the three panels are at the same scale.

In the text

	Fig. 7 Upper panel: F2 for the single-star fit of the simulated data with the addition of surface brightness asymmetries. A F2 > 3 indicates a non-satisfactory fit. Lower panel: the same as the upper panel but with the theoretical prediction (Eq. (11)) added on top (see text).
In the text

	Fig. 8 Distribution of the different types of solutions found as a function of the distance. In the three panels are shown the results for the same sample of stars to which different photocentric motions have been applied, as indicated by the corresponding σP and τ in the figure label.
In the text

Fig. 9 Left panel: distribution of the ratio a/σa, with a the semi-major axis of the orbit and σa its uncertainty. Solutions with a/σa < 2 are not significant and may be rejected. Right panel: ratio a/R, with R the stellar radius for stars with a/σa > 2. All orbits appear to have a semi-major axis smaller than the stellar radius and can thus be rejected. Markov chain parameters: σP = 0.2 AU and τ = 4 × 106 s.

In the text

	Fig. 10 Semi-major axes a of the orbits derived with the orbital solution. The physical size of a in AU is comparable to the photocentre motion added to the simulated Gaia data. Markov chain parameters: σP = 0.2 AU and τ = 4 × 106 s.
In the text

Fig. 11 The blue horizontally striped histograms show the “stochastic” noise added to Gaia measurements in the computation of a stochastic solution. The green diagonally striped ones represent the stochastic noise plus the offset of the mean photocentre position with respect to the centre of the star. The three panels show results for three different sets of parameters of the Markov chain. The vertical dashed line corresponds to the standard deviation σP in the photocentre displacement.

In the text