Distances and parallax bias in Gaia DR2

TL;DR

This paper derives Bayesian distances for Gaia DR2 stars, identifies a parallax bias of about -0.054 mas, and emphasizes the importance of accounting for this bias in astrophysical analyses.

Contribution

It provides a direct validation of Gaia DR2 parallaxes using a Bayesian method across all stars in a sample, revealing a significant parallax bias and its dependence on measurement uncertainties.

Findings

Parallax bias delta p ≈ -0.054 mas in Gaia DR2.

Bias depends strongly on parallax uncertainty sigma_p.

Correcting bias influences astrophysical inferences from Gaia data.

Abstract

We derive Bayesian distances for all stars in the RV sample of Gaia DR2, and use the statistical method of Schoenrich, Binney & Asplund(2012) to validate the distances and test the Gaia parallaxes. In contrast to other methods, which rely on special sources, our method directly tests the distances to all stars in our sample. We find clear evidence for a near-linear trend of distance bias f with distance s, proving a parallax offset delta p. On average, we find delta p = -0.054 mas (parallaxes in Gaia DR2 need to be increased) when accounting for the parallax uncertainty under-estimate in the Gaia set (delta p = -0.048 mas on the raw parallax errors) with negligible formal error and a systematic uncertainty of about 0.006 mas. The value is in concordance with results from asteroseismic measurements, but differs from the much lower bias found on quasar samples. We further use our method to compile a comprehensive set of quality cuts in colour, apparent magnitude, and astrometric parameters. Last, we find that for this sample delta p appears to strongly depend on the parallax uncertainty sigma_p (when including the additional 0.043 mas) with a statistical confidence far in excess of 10\sigma and a proportionality factor close to 1, though the dependence varies somewhat with sigma_p. Correcting for the sigma_p dependence also resolves otherwise unexplained correlations of the offset with the number of observation periods n_{vis} and ecliptic latitude. Every study using Gaia DR2 parallaxes/distances should investigate the sensitivity of their results on the parallax biases described here and - for fainter samples - in the DR2 astrometry paper.

Tables

Table 1. Table 1: Parameters of the selection function (eq.  6 )

Parameter	value	unit
$a$	$97.9594$
$b$	$4.10228$	${\mathrm{kpc}}^{-1}$
$c$	$0.0159642$
$d$	$0.15$	${\mathrm{kpc}}^{-1}$
$l$	$2.22279$	${\mathrm{kpc}}^{-1}$
$l_2$	$2.97547$	$\mathrm{kpc}$
$k_2$	$2.0704$
$j$	$0.956364$	$\mathrm{kpc}$
$h$	$4.91193$
$z$	$0.024$	${\mathrm{kpc}}^{-1}$

Equations

$$f=(s^{\prime}-s)/s$,$$

$

$$\left(\begin{array}{c}U\\ V\\ W\end{array}\right)=\mathbf{M}(\mathbf{I}+f\mathbf{P}){\mathbf{M}}^t\left(\begin{array}{c}{U}_0\\ V_0\\ W_0\end{array}\right)$$

$$\left(\begin{array}{c}{T}_{\mathrm{uw}}\\ T_{\mathrm{vw}}\\ T_{\mathrm{ww}}\end{array}\right)=\left(\begin{array}{c}\cos l\sin b\cos b\\ \sin l\sin b\cos b\\ 1-\cos b^2\end{array}\right)$$

$$P(s)=N^{-1}s^{2}G(p,p_{0},{\sigma_{\rm p}})\rho(s({\rm p}),{l},{b})S(s({\rm p}))$,$$

$

$$N=\int ds{s}^{2}G(p,p_0,{\sigma }_{\mathrm{p}})\rho (s(\mathrm{p}),l,b)S(s(\mathrm{p}))$$

$$S(s)=aA(s)B(s)C(s)$$

$$A(s)$$

$$B(s)$$

$$C(s)$$

Full text

††pubyear: 2019††pagerange: Distances and parallax bias in Gaia DR2–References

1 Introduction

This paper presents the first statistical evaluation of parallax biases for a general stellar sample and the first derivation of unbiased stellar distances for the subsample of Gaia DR2 with line-of-sight velocity measurements (Gaia Collaboration et al., 2018a; Katz et al., 2019). Other than most other approaches, our method tests the distances of all stars in the selected sample and not just specific subgroups with different physical properties.

Since the first release of Gaia DR2 data, there has been clear evidence for offsets in the parallax measurements. In the data release papers, Lindegren et al. (2018) found a general offset of Gaia parallaxes by ${\delta_{\rm p}}=-0.029\,{\rm mas}$ (in the sense that the quoted values were too low) when studying known quasars in the Gaia catalogue. In addition, the median offset of the quasars is patchy on the sky (see also Fig. 15 in Arenou et al., 2018), and thus we also have to expect an error term ${\sigma_{\rm p}}$, which will in many aspects behave like a random error (as long as these patches are not resolved). The calibration involved in the data analysis for the astrometric solution of Gaia depends on both effective wavelength (colour) and magnitude (directly and through the ’window class’), so both offsets/biases and random errors cannot be expected to be homogeneous along the entire magnitude and colour range of the sample. The quasar catalogue, which is the best direct benchmark for Gaia parallaxes, is very faint (typically at apparent magnitudes ${\rm G}>17$), and the colour distribution does not match the stellar colour distribution very well. In contrast their internal checks between different parts of the analysis show significant changes in behaviour in particular at magnitudes ${\rm G}\sim 16,13,10$ (caused by their analysis windows, see Fig. 16 in Lindegren et al., 2018) and changes in their astrometric measurement accuracy (likely near ${\rm G}\sim 8$, see their Fig. 9). In addition, sources in several remote objects, in particular the LMC display a finely structured pattern in sky position of mean offsets related to the scanning law of Gaia.

In Gaia Collaboration et al. (2018d) the same offset can be seen when comparing the globular cluster mean parallaxes to the more accurate (at small parallax) values from the Harris (1996) catalogue. Gaia Collaboration et al. (2018d) also show that the Large and Small Magellanic clouds, and most of the satellite dwarf spheroidal galaxies they study, have negative average parallaxes. Consistent with these general offsets, but differing in the magnitude of the effect, Stassun & Torres (2018) found an offset of $(-0.082\pm 0.033)\,{\rm mas}$ for a sample of $89$ very bright ($G<12\,{\rm mag}$) binary systems, however with a large scatter of the single measurements. Yet, given the patchiness of the offset, and the very different magnitude range probed, these values could well be in agreement. Just when we were about to submit this paper, the latest study of Graczyk et al. (2019) found a value of ${\delta_{\rm p}}=-(0.054\pm 0.024)\,{\rm mas}$, which by an amusing coincidence is exactly our result, albeit with a much larger uncertainty. We also note that binarity was not worked into the astrometric solution, so it is to be expected that binary systems carry a different bias from the main stellar sample. Similarly, comparison with asteroseismic values points to significant parallax underestimates (Zinn et al., 2018), though again for rather specific subsets of stars in magnitude and colour/stellar evolutionary stage. However, their result of ${\delta_{\rm p}}=-0.05\,{\rm mas}$ bears high confidence, and again differs from the quasar result. Sahlholdt & Silva Aguirre (2018) looked at a smaller sample of dwarfs with asteroseismic data, and the offset they found was closer to that of the quasar sample, but with more significant uncertainty ($-0.035\pm 0.016\,{\rm mas}$). No test so far has measured what we are really concerned about: the bias of the full stellar sample and how it depends on the other properties of the observation.

Unbiased distances are key to virtually every problem in modern astrophysics, and given the large sample sizes, we now need these distances at the $1\%$ level. For example the wave-pattern discovered by Schönrich & Dehnen (2018), which was confirmed by Huang et al. (2018) and Kawata et al. (2018), and which is likely related to the later findings of Antoja et al. (2018), has a total amplitude of below $1\,{\rm km}\,{\rm s}^{-1}$. To compare this: even at a perfect location towards the Galactic anticentre, the solar reflex motion will translate into a $\sim 0.1\,{\rm km}\,{\rm s}^{-1}$ bias for every $1\%$ in mean distance bias. Larger and more complex bias patterns will arise at other sky positions from cross-correlations between the velocity component measurements.

This work offers a solution to this problem. We will derive unbiased distances to the RV subsample (about $7$ million stars at magnitudes ${\rm G}\lesssim 15$) using the method proposed in Schönrich & Aumer (2017). By measuring the selection function directly from the sample, we derive a data-informed, and thus nearly unbiased prior, which suits a sample better than model-based priors (e.g. Astraatmadja & Bailer-Jones, 2016; Bailer-Jones et al., 2018) designed to fit pre-existing models of the entire Gaia sample, e.g. Bailer-Jones et al. (2018). As discussed in Schönrich & Aumer (2017), mismatch or neglect of the selection function will result in systematic bias of similar size to the measurement uncertainties for individual stars (i.e. of order $20\%$ for the common parallax quality cut of ${\rm p}/\sigma_{p}>5$). Our previous results have established that this method gives an unbiased translation from parallaxes to distances. Consequently, with these distances we can now directly measure and correct biases in the parallaxes using the statistical method of Schönrich et al. (2012).

Our paper is structured as follows: we start with a short description of the used data and coordinate system definitions in Section 2, followed by a description of our statistical distance estimator in Section 3. After this, we provide the formalism for deriving Bayesian distances in Section 4, including a derivation of the distance-dependent selection function $S(s)$. In Section 5 we quantify the different parallax biases in Gaia DR2, followed by a comparison to previous distance derivations and a comment on the distance to the Pleiades. Section 8 provides a summary of quality cuts necessary in Gaia DR2, followed by the Conclusions.

2 Data and definitions

2.1 Coordinate frame and definitions

Throughout this paper, we will use the standard definitions for Galactic coordinates and the Local Standard of Rest. We employ Galactic cylindrical coordinates $(R,z,\phi)$, where $R$ is the in-plane distance to the Galactic centre, $z$ is the altitude above or below the Galactic midplane, and $\phi$ is the Galactic azimuth, with the Sun placed at $\phi=0$. The distance of a star to the Galactic centre is termed $r=\sqrt{R^{2}+z^{2}}$, for the solar galactocentric distance, we use the value ${R_{0}}=8.27\,{\rm kpc}$ from Schönrich (2012), which is also in agreement with other determinations (Gillessen et al., 2009; McMillan, 2017), and only slightly in tension with the latest estimates from measurements of stellar orbits around Sgr A* from Gravity Collaboration et al. (2018). The vertical displacement of the Sun from the midplane, $\,{\rm z}_{\odot}=0.02\,{\rm kpc}$ is taken from Joshi (2007). We also tested that it does not have any significant impact on our results. The velocity vector in the heliocentric,111We use the somewhat negligent Galactic dynamicists’ term ”heliocentric”, while in truth, Gaia is measuring quantities in the solar system baryocentric frame. With relative motions between the two frames below $0.1\,{\rm km}\,{\rm s}^{-1}$, this difference does not matter. Cartesian frame is defined as $(U,V,W)$ with a right-handed set of components pointing radially inwards, in the direction of Galactic rotation, and upwards perpendicular to the plane. The velocity vector components in the Galactocentric cylindrical frame are analogously termed $(U_{g},V_{g},W_{g})$. To translate these velocity components, we use the motion of the Sun against the Local Standard of Rest as measured in Schönrich et al. (2010); Schönrich (2012): $({U_{\!\odot}},{{\boldsymbol{\upsilon_{\odot}}}},{W_{\!\odot}})=(11.1,250,7.24)\,{\rm km}\,{\rm s}^{-1}$, and if necessary, use the azimuthal velocity of the Sun against the Local Standard of Rest (${V_{\!\odot}}={{\boldsymbol{\upsilon_{\odot}}}}-V_{c}=12.24\,{\rm km}\,{\rm s}^{-1}$). For simplicity’s sake, we call $p$ the parallax of a star and ${\sigma_{\rm p}}$ the effective uncertainty of the parallax measurement assumed in that instance, which, depending on the examined set of assumptions, may contain the additional $\delta{\sigma_{\rm p}}=0.043\,{\rm mas}$ added in quadrature to the Gaia pipeline value ${\sigma_{\rm p,g}}$.

2.2 Data

Here we use the Gaia RV sample (Cropper et al., 2018; Gaia Collaboration et al., 2018c) from Gaia DR2 (Gaia Collaboration et al., 2018a) with more than $7$ million stars that have both astrometric and line-of-sight velocity measurements from the onboard spectrograph (Sartoretti et al., 2018) on the Gaia spacecraft (Gaia Collaboration et al., 2016). To ensure the quality of the data, we apply, if not stated otherwise for a specific task, a few quality cuts that were discussed in data release papers, e.g. Lindegren et al. (2018), namely: the number of visibility periods ${\rm n_{\rm vis}}\geq 5$ to ensure a full astrometric solution, a parallax quality cut of ${\rm p}/\sigma_{p}>5$ and lower limit ${\rm p}>0.1\,{\rm mas}$, which translates to approximately demanding a distance $s\lesssim 10\,{\rm kpc}$, an excess noise smaller than $1$, line-of-sight velocity limits of $|v_{\rm los}|<550\,{\rm km}\,{\rm s}^{-1}$ and $\sigma_{\rm los}<10\,{\rm km}\,{\rm s}^{-1}$. We usually remove the Galactic mid-plane from our sample, i.e. require $|{b}|>10\,{\rm deg}$. For our statistics the Galactic mid-plane carries no signal and we thus avoid problems with crowding and excessive reddening. We checked, though, that the measurement of the distance prior from low-$|b|$ data is similar to our higher latitude main sample, and provide distance estimates for these stars in the derived catalogue. In previous papers (see Schönrich & Aumer, 2017) we uncovered major problems with $v_{\rm los}$ measurements, in particular with LAMOST. Here, we just note that our tests of $v_{\rm los}$ accuracy and precision looked very decent on the Gaia sample, and we will concentrate on the more pressing issue of distances and Gaia parallaxes. We follow the convention of Gaia papers to call their apparent magnitudes in the three broad colour bands $({\rm G_{BP}},{\rm G},{\rm G_{RP}})$. Problems with capitalisation conventions do not arise, since we do not discuss absolute magnitudes through most of the paper.

3 Statistical Distance Estimation

3.1 General thought

For the determination of distance bias, we rely on the method of Schönrich et al. (2012) (hereafter SBA), which has been applied on various samples. The method relies on correlations between velocities, which depend on the position on the sky. The estimator is readily derived by writing down an estimate of stellar kinematics while allowing for a systematic distance bias $f=\langle s^{\prime}/s\rangle$, where $s^{\prime}$ denotes the estimated distance, $s$ the real distance to a star. This $f$ affects both tangential velocity components simultaneously, correlating them. To explain this with a simple example: imagine approaching a mountain horizontally. Knowing your velocity, your mind automatically has a clear estimate of the distance to the summit. This is because any incorrect estimate would translate your horizontal motion into a vertical component of motion of the mountain, i.e. the top of the mountain would have to be growing or shrinking (if you had over- or under-estimated the distance, respectively). Analogously, all parts of the mountain base below your level would appear to be moving downwards (upwards). And your brain knows this is not usually what mountains do. The SBA method allows us to extend this intuition and make it more robust against assumptions (i.e. we do not assume any mean motion or fixed velocity ellipsoid) to stellar samples.

The strength of the SBA method now lies in directly using the spatial dependence of this correlation of the heliocentric velocity components on galactic longitude (${l}$) and latitude (${b}$). As long as we have a sufficient sky coverage, we do not depend on classic assumptions of other methods: think of our mountain example. Observing many mountains around us, we gain a significant advantage in control of systematics over the use of just one single mountain. In a simplified picture, we just measure the pattern of apparent vertical velocities of all the mountain tops and mountain bases around us and try to find the distance correction that makes the angular dependence of this pattern disappear. Misjudging our own velocity, we would equally bias the motion of all summits and mountain roots, leaving the statistical distance evaluation unscathed. Similarly, we do not care about our horizontal velocity, since it can i) be measured and ii) would just affect our prediction for the strength of the effect around us; however, we just seek the distance factor at which the correlation of vertical velocity components with sky position disappears, making our own horizontal motion irrelevant. Translated to our real problem: Assumptions about the solar velocity do not matter for our method.

Similarly, let us assume that we are sitting in the middle of an orogeny (or reading Calvino’s Cosmicomics) and both mountain summits and roots are rapidly rising and sinking into the ground. Observing just one mountain in front of us, we would indeed infer that distances are overestimated, but behind us, the correlation term reverses sign, i.e. the apparent distance underestimate there cancels out the distance over-estimate inferred from the opposite direction. We learn from this that typically modes of the disc cancel out in a sample with large sky-area. Analogously, a wide halo stream passing through would cancel by the spatial terms. In short, galaxy physics can only affect our statistical measurements, if they vary across the sky in a way that correlates with the angle-terms of our method. In most cases (e.g. global breathing modes, streams) they will cancel out at first order.

Last, we note that other than our imaginary mountains, stars move horizontally, so in addition to the spatial correlations, we can benefit from two different horizontal velocity components with different dependence on sky position.

3.2 Formal argument and specific implementation

The formal derivation of our method (see SBA for a stringent treatment) is done by simply writing down what happens in the measurement. The vector of observed values $(s^{\prime}\mu_{l},s^{\prime}\mu_{b},v_{\rm los})$, where $s^{\prime}$ is the observationally inferred distance, $\mu_{l},\mu_{b}$ are the proper motions in Galactic longitude ${l}$ and latitude ${b}$. This vector is translated into the measured velocity components $(U,V,W)$ by a matrix ${\bf M}$ depending on ${l}$ and ${b}$. Since this is an orthogonal matrix, the inverse mapping (i.e. from the original velocity components) is done with the transverse ${\bf M}^{t}$. If we now assume that distances are changed by some relative bias

[TABLE]

where $s$ is the real distance, we can relate:

[TABLE]

where the index [math] indicates the real values, and $P$ is $\rm{diag}(1,1,0)$, which projects to the two proper motion components. Now, we see that the observed velocity components are correlated by $f$ via the matrix ${\bf T}={\bf M}{\bf P}{\bf M}^{t}$. Since the equations are linear, the average $f$ can thus be gained by a similar linear regression of any target velocity component $v_{i}$ onto the other velocity components $v_{j}$ multiplied with $T_{ij}$, the components of the matrix ${\bf T}$.

As discussed in SBA, using the in-plane velocity components $(U,V)$ mixes the statistics with a Galactic rotation estimate, and given the very large sample size here, we make the choice to avoid this possible source of systematic bias. We thus limit this study to using the correlation of both $U$ and $V$ velocity components with the vertical motion $W$. The relevant part of $T$ is thus

[TABLE]

Our method applies corrections for the following biases of this measurement:

• •

$v_{\rm los}$ determination errors, $\sigma_{\rm los}$, which would appear as distance underestimates (typically negligible due to the excellent precision and accuracy of the Gaia $v_{\rm los}$ estimates),

• •

proper motion determination errors, $\sigma_{\pm}$, which would appear as distance over-estimates, but are again negligible by more than an order of magnitude,

• •

the tilt of the velocity ellipsoid, which is of some importance for the statistics. This term is important, as the radially elongated velocity ellipsoid produces a locally changing correlation between the heliocentric velocity components, which can partially line-up with the $T_{\rm uw}$ and $T_{\rm vw}$ angle combinations.

We have to add the systematic uncertainties from these terms to our error budget, assuming that the uncertainties in these terms are statistically independent of each other. As already done in previous studies (Schönrich et al., 2012) we assume a systematic uncertainty of $10\%$ of the calculated correction value for the first two terms, and an uncertainty of $30\%$ for the turn of the velocity ellipsoid correction.

As stated above, due to the unprecedented precision of Gaia, the exact values of proper motion errors do not matter here as long as the order of magnitude of the uncertainty estimates in the Gaia pipelines is correct. Similarly, the error correlations are mostly inconsequential: two team members did independent tests on independently calculated mock samples, where we folded the mock measurements with the full error matrix between parallaxes and proper motions as given for each star in the sample. In these tests, the effect of error correlations on our statistics is more than one order of magnitude less than our systematic and formal error budget for whole-sky measurements. On pencil beams, like in tests of high $\beta$ stars, it contributes of order one tenth to the residual bias (see analysis below).

For the velocity ellipsoid correction, we assume that the velocity ellipsoid points to the Galactic centre at every position. Other than in previous applications of the Schönrich et al. (2012) method, the Gaia sample spans a large volume throughout the disc, and when we select by distance, stars in each sample will cover regions with vastly different values of the velocity dispersion. To optimise the estimate for this correction term, we directly measure the velocity dispersions in the Galactocentric spherical coordinate frame weighted by their impact on the distance estimator.

As described in Section 3.1 the method does not assume any velocity ellipsoid, and does not even require knowledge of the correlations between the velocity components (e.g. $U$ and $W$). The only important requirement is that Galactic structure does not infer a correlation between this velocity correlation and Galactic position. Realistic structure (e.g. a stream passing through the survey, or disc breathing modes) might produce a local velocity correlation (and so distance statistics on small patches of sky are uncertain, see the $\beta$ issue below), but cancels out to first order with large sky coverage. We have already tested and confirmed this on realistic simulations in the Appendix of Schönrich & Aumer (2017).

4 Distance Determinations and Selection Function for Gaia DR2

Before we can test distances/parallaxes with the SBA method, we have to derive distance expectation values for all stars, using the method of Schönrich & Aumer (2017). From parallax measurements we calculate the probability distribution in distance $P(s)$ for each star by

[TABLE]

where

[TABLE]

is the normalization, $s$ is the distance from the Sun, ${\rm p}$ denotes a parallax, $G(p,p_{0},\sigma_{p})$ is the (Gaussian) observational likelihood distribution in parallax, given the measurement $p_{0}$ and effective uncertainty $\sigma_{p}$, and $\rho(s(p),{l},{b})$ is the assumed density model. $S(s)$ denotes the selection function, i.e. number of stars detected in the sample divided by the number of stars actually there. As in previous papers, we use the simple density model from Schönrich & Bergemann (2014), which contains a thin disc, thick disc, and halo component (we can neglect the bulge because we only select stars with $|{b}|>10\,\,{\rm deg}$). All calculations are done with a self-adaptive trapezoid integration, where we start from the suspected maximum of the PDF and integrate to both sides in distance simultaneously, adapting the step-length upwards when the relative contribution of each segment to the integrated value falls below a threshold. We tested that our results are the same when lowering the initial step-length, and threshold by a factor $10$, or when using a simple step-wise integrator.

We tested by cutting the sample in Galactic latitude and longitude that this model is sufficiently close to the underlying distribution to ensure a good distance measurement. The most important quantity in the above equation, however, is the distance-dependent selection function $S(s(p))$, which arises from the magnitude-dependent cuts of the sample. Predicting $S(s)$ directly would require a full chemodynamical galaxy model including a three-dimensional reddening map. With this choice, we would be vulnerable to systematic uncertainties of the model choice, stellar evolution models, and the reddening map.

Here, we instead measure $S(s)$ from the data themselves. Nevertheless, it is useful to formulate a model-based expectation, which acts as a sanity check to our results. The top panel of Fig.  1 shows the selection function $S^{\prime}(s)$ calculated from simple population synthesis using the machinery of Schönrich & McMillan (2017) using B.A.S.T.I. stellar models (Pietrinferni et al., 2004, 2009) and a simple Salpeter (1955) Initial Mass Function (IMF). The shown $S^{\prime}(s)$ expresses the number of stars per solar mass (at birth) of a stellar population at a given distance $s$. Since the normalisation factor here is irrelevant for our purposes, we multiplied $S^{\prime}(s)$ with an arbitrary factor $80$ to facilitate direct comparison with the measured $S(s)$ in the bottom panel. For this figure, we suppose a simple magnitude-dependent selection, setting the selection probability constant below Johnson V-band $m_{V}<12.8\,{\rm mag}$ and then going linearly to zero for $m_{V}\geq 13.5\,{\rm mag}$.

$S(s)$ has three main regions: i) a quite steep roughly exponential decrease of $S^{\prime}(s)$ with distance $s$ in the near field, stemming from the magnitude limit moving up through the main sequence and turn-off region with increasing distance modulus. The exponential behaviour is expected, since the magnitude scale is logarithmic, the luminosity of stars is roughly the fourth power of the mass, and the IMF is close to a simple power law (as we can neglect the impact of stars below \sim$$0.5M_{\odot} where this is no longer true). This decrease is followed by ii) a weakly inclined plateau associated with the subgiant branches and red clump, followed by a downward knee and somewhat shallower decline at the largest distances due to iii) the red/asymptotic giant branch stars. There are some differences between different ages and metallicities: metal-poor stars are somewhat more luminous on average, shifting the features to the right, older populations have their subgiant branch at fainter magnitudes and less stars overall on the giant branch (the lengthening of stellar main sequence lifetimes outweighs the increase in the IMF). However, the common features let us predict the functional shape for the selection function rather well, when averaged over all stellar populations. We choose:

[TABLE]

where $a$ is a normalization constant, and the three multipliers are:

[TABLE]

The rationale behind this is to capture in $A(s)$ the general shape with two exponentials of scale-length $b^{-1}$ and $d^{-1}$ for the short (out to $\sim 1\,{\rm kpc}$) and long (beyond $3\,{\rm kpc}$), model the step in $S(s)$ with $\tan^{-1}$ in $B(s)$, as well as the (unimportant) drop-out of luminous or otherwise too close stars (potential loss due to proper motions) with $C(s)$.

We proceed as in Schönrich & Aumer (2017), iterating the distance calculation with the adapted prior. However, due to the good nominal quality of Gaia parallaxes, the general appearance of the measured $S(s)$ is already present from the first iteration using a flat prior. A fit of the function to data corrected for a mean parallax offset of ${\delta_{\rm p}}=-0.048\,{\rm mas}$ (see following below) after several iterations of the prior fitting is presented in Fig. 1, the values are provided in Table 1. The shape matches the prediction from the population synthesis. We also checked that the selection function does not vary strongly with galactic latitude, which signals that the spatial prior is sufficient and population differences throughout the galaxy will not strongly bias our distance estimates.

Apart from $d$, which sets the scalelength of the long exponential component, all parameters are well-constrained in the fit. The latter suffers from the fact that we cannot fit beyond the point where the parallax quality cut affects the sample. This cut must not enter $S(s)$, since the stars are culled after they have passed the magnitude limits. Inspection of Fig. 1 shows that given the relative drop-out rates of stars at different quality cuts, the fit must be close to the real shape. Extensive further tests of the far-distance end favoured an additional flattening of $S(s)$ by multiplying with $\exp(-(s^{\prime}/\,{\rm kpc}-4)/0.07)$, where $s^{\prime}={\rm min}(s,10\,{\rm kpc})$. The difference can be inspected between the top panel of Fig. 2, which does not contain this factor and the bottom panel, which does include it. It mostly serves to improve a slight kink in the distance statistics. The latter factor only has a minor effect on very remote stars beyond $s>3\,{\rm kpc}$, and we will apply it to all data shown from Fig. 3 onward. Both the model expectations and the derived distance statistics for very remote stars, which react strongly to the prior, support this choice.

5 Distance bias vs. distance: evidence for the parallax offset

5.1 Scanning vs. distance: Detection of the parallax offset

Fig. 2 measures the relative distance bias $1+f$ when using the Schönrich et al. (2012) method for the stellar sample when ordered and binned in the measured distance $s$. The top panel shows this scan of $1+f$ vs. $s$ for two different quality cuts in the relative parallax error, allowing a maximum ${\sigma_{\rm p}}/{\rm p}$ as given by Gaia DR2 of $10\%$ or $20\%$ respectively. It is apparent that $1+f$ increases nearly linearly with $s$, reaching a relative distance bias in excess of $5\%$ at $s\sim 1\,{\rm kpc}$ and errors in excess of $25\%$ near $s\sim 5\,{\rm kpc}$. Such a bias precludes a precision measurement of galactic kinematics even in the relative near field. The comparison of different ${\sigma_{\rm p}}/{\rm p}$ cuts proves that this cannot be an issue with the distance priors or other parts our measurement method. Otherwise, a tighter quality cut would drastically reduce the bias $1+f$ at a given distance. This leaves only one culprit: a bias in the Gaia parallaxes.

In fact, as we can see from the bottom panel in Fig. 2, this linear trend is a signature imprint of such a parallax bias. When we apply the correction ${\delta_{\rm p}}$ to all parallax measurements, we can minimise the trend for ${\delta_{\rm p}}=-0.048\,{\rm mas}$ with an uncertainty of about $0.006\,{\rm mas}$. The most distant bins beyond $s>3\,{\rm kpc}$ are in line with this estimate within the uncertainties. Small differences at these distances should be ascribed to uncertainties in $S(s)$ as discussed above. The offset is about double the amount found by Lindegren et al. (2018), but exactly in line with the asteroseismic evaluation of Gaia DR2 (Zinn et al., 2018). We note again that we do not believe the evaluations in Arenou et al. (2018) and Lindegren et al. (2018) on quasars to be applicable to our case, since those are in different apparent (magnitude) window classes with separate astrometric calibrations, and have all different kinematics/zero intrinsic parallax.

We further note that there is a significant positive bias $f\sim 5\%$ for the distances of the nearest stars or, equivalently, brightest parts of the sample (seen in the left-most green datapoint in the top panel Fig. 2). This bias is near impossible to explain with bad $v_{\rm los}$ measurements, which would feign a negative $f$. This bias is for these nearby stars orders of magnitudes larger than the previously discussed parallax offset, which for these stars is negligible. Some of this may be traced back to the magnitude dependent deviations (see below, Fig. 9), some to mis-identifications in the near field.

So far, we simply applied a parallax offset. However, from Lindegren et al. (2018) a position dependent variation in the parallax bias is expected. There has been discussion in the Gaia collaboration, if such a bias should be added to the error budget or not. Now, in our case, we are interested in the uncertainty for each single star. A priori, we do expect a random, but spatially correlated fluctuation of the parallax offset to enter single star parallax uncertainties just like an additional term that has to be added in quadrature to the formal uncertainty given by the pipeline, setting ${\sigma_{\rm p}}=\sqrt{{\sigma_{\rm p,g}}^{2}+(\delta{\sigma_{\rm p}})^{2}}$.

We thus quantify the correction ${\delta_{\rm p}}$ on two versions of the sample, once with and once without systematic $\delta{\sigma_{\rm p}}$, which we take to be $0.043\,{\rm mas}$ following the quasar analysis of Lindegren et al. (2018, table 4). One can argue for two ways of measuring the parallax bias: On one hand, we could demand that the average $f$ in the safe region $0.1<s/\,{\rm kpc}<3$ should be exactly [math]. While we have high confidence in the accuracy of the distance statistics in this area, one might still want to get rid of the need for a correct zero point. Indeed, since a parallax offset gives an almost linear $f(s)$ dependence, we can alternatively demand that the estimated slope $df/ds$ should be zero. However, within our systematic uncertainties, both methods, shown in Fig. 3 yield the same values for ${\delta_{\rm p}}$. Drawing a mean estimate from the shown results and further tests of varying quality cuts, we conclude that ${\delta_{\rm p}}=-0.048\,{\rm mas}$ with a negligible formal uncertainty and a systematic uncertainty of $\sim 0.006\,{\rm mas}$. The systematic uncertainty is a cautious estimate. The systematic uncertainties on velocity ellipsoid and measurement uncertainties are already priced into the formal errors, but we price in a correlated systematic error between evaluation bins, and performed a variation of fit parameters (region in $s$ to fit on), changes in quality cuts (see below), and comparison of the different methods. These would advise a mildly smaller number, and we added a budget for effects that we may have missed. Separate from these effects, accounting for $\delta{\sigma_{\rm p}}=0.043\,{\rm mas}$ further increases the $|{\delta_{\rm p}}|$ estimate by about $0.006\,{\rm mas}$ to ${\delta_{\rm p}}=-0.054\,{\rm mas}$. The small difference between the $f$ and the $df/ds$ estimation for ${\delta_{\rm p}}$ can be ascribed to bad luck as it is within $\sim 2\sigma$. We think, however, that it has contributions by one or several of the secondary problems identified below.

5.2 Quantifying parallax offset vs. uncertainty

Could we make different assumptions for parallax error visible with our distance method? In principle yes - the background is that the expectation value of a single star’s distance shifts systematically with the assumed error. Consequently, a misjudgement of uncertainties will typically show as a distance bias. However, as we saw in Fig. 3, this effect is subordinate to the other problems in the sample. Nevertheless, assessing the right value for the ${\rm p}/{\sigma_{\rm p}}$ quality cut requires a sample scan in this quantity. This is done in the top panel of Fig. 4. We remind that apart from stars at very small distances (corresponding to large ${\rm p}/{\sigma_{\rm p}}$) in this plot, the distance scan (see Fig. 2) shows no significant deviations after correcting for ${\delta_{\rm p}}$. Yet, the red line in the top panel of Fig. 4 reveals an impressive increase in $1+f$ towards lower ${\rm p}/{\sigma_{\rm p}}$. Note, however, that we have extended the sample to stars beyond our usual quality limit, admitting all stars with ${\rm p}/{\sigma_{\rm p}}>1$; below ${\rm p}/{\sigma_{\rm p}}<4$ there will definitely be a sizeable fraction of stars with catastrophic distance misestimates. Also, the large parallax uncertainty allows for major parts of the probability distribution $P(s)$ to cover highly uncertain regions of $S(s)$. However, extensive experimentation with $S(s)$ could not rectify the abnormalities in Fig. 4. One could now be tempted to argue that the suggested inclusion of $\delta{\sigma_{\rm p}}=0.043\,{\rm mas}$ removes the problem, but this is a deception caused by the left-shift of the graph, since we increased ${\sigma_{\rm p}}$. The $S(s)$ uncertainty can be resolved by limiting the ${\rm p}^{-1}$ range to contain the PDFs of even uncertain parallax measurements within the safe region $s<3.5\,{\rm kpc}$. The result of this is shown in the bottom panel of Fig. 4: there always remains a spike in $1+f$ on the left-hand side. Something else is going on here.

We now restore the parallax quality limit and plot $1+f$ against the pipeline parallax error ${\sigma_{\rm p}}$ itself in the top panel of Fig. 5. For ${\sigma_{\rm p}}>0.12\,{\rm mas}$ distance bias sky-rockets, even for the rather conservative cut ${\rm p}^{-1}<1.5\,{\rm kpc}$ applied here. This suddenly clarifies the tension in the plots of $1+f$ against ${\rm p}/{\sigma_{\rm p}}$, because we were just moving the region of large parallax errors via the distance/parallax cut. A further, unsettling observation is the uptrend of the distance bias with parallax error between ${\sigma_{\rm p}}=0.045\,{\rm mas}$ and ${\sigma_{\rm p}}=0.06\,{\rm mas}$, turning from a negative to a positive distance bias of order $2\%$.

The bottom panel of Fig. 5 attempts to qualify the nature of this failure. One could argue that stars with large parallax uncertainties noted by the pipeline should be binaries, affecting their proper motions, and thus our distance statistics. However, we see a clear increase of the distance bias with distance, putting the blame at a parallax offset exceeding by far the $\sim 0.05\,{\rm mas}$ of the entire sample, for which we already correct in the shown data. Given the displayed results, the only viable explanation is that $\delta{\rm p}$ is at least to some extent proportional to ${\sigma_{\rm p}}$.

With such an unusual finding, it is natural to point the search for an honest error at our own code. E.g. there could be a typo in our distance integral creating the dependence on ${\sigma_{\rm p}}$. Apart from controlling and testing our code, Fig.6 investigates this by checking different assumptions for the distance error. The purple line displays the original trend as found with a constant ${\delta_{\rm p}}=-0.048\,{\rm mas}$. The blue line corrects this to setting ${\delta_{\rm p}}={\sigma_{\rm p}}$ and the green line in addition adds $\delta{\sigma_{\rm p}}=0.043\,{\rm mas}$ in quadrature. If we had made a ${\sigma_{\rm p}}$ dependent error ourselves, the green line and the blue line should deviate in the same way (as the distance estimates depend very weakly on the assumed ${\sigma_{\rm p}}$). In contrast, we see that only with the full correction and assuming that ${\delta_{\rm p}}=-{\sigma_{\rm p}}$, we can rectify the trend in the sample.

While we do not want to adhere to the notion that the parallax offset is perfectly equal to the parallax error, we do in Fig. 7 attempt to show a quantify this dependence. In both panels we use our usual quality cuts in colour, $g$-band magnitude and ${\rm p}/{\sigma_{\rm p}}>4$. In addition, we removed the parts of the sample with $\beta<-55\,{\rm deg}$, which, however, has very minor effect. The top panel demonstrates that if we would assume that ${\delta_{\rm p}}=-{\sigma_{\rm p}}$, both the average distance error and the trend of distance error with distance are within the systematic uncertainties, in line with assuming the usual additional parallax error $\delta{\sigma_{\rm p}}\sim 0.043\,{\rm mas}$ added in quadrature. The bottom panel tests different values of the proportionality constant $q$, when we set ${\delta_{\rm p}}=-q{\sigma_{\rm p}}$ after adding the additional term to ${\sigma_{\rm p}}$. Both statistics are in line with a value very close to $1$.

While the notion of ${\delta_{\rm p}}=-{\sigma_{\rm p}}$ is simple at first hand, it is rather unreasonable to believe that there should be such a perfect equality. To this end we tested a third option: While leaving $q$ free again, we added the offset found for the quasars to the parallax, i.e. set $-{\delta_{\rm p}}=0.029\,{\rm mas}+q{\sigma_{\rm p}}$. Since about half the offset in parallax needed is now again captured by a constant term, we can expect that a good solution will be found for $q\sim 0.5\,{\rm mas}$. This result is shown in Fig.  8, where we plot the distance bias against ${\sigma_{\rm p}}$ for both $q=1$ and $q=0.5$ plus the constant term. However, various experiments show that we cannot get rid of the trend of $f$ at small ${\sigma_{\rm p}}$, if $q$ is not close to 1 in this region. The equality is definitely not perfect, since we would require a $q$ closer to $0.5$ in all possible assumptions near ${\sigma_{\rm p}}\sim 0.09\,{\rm mas}$ and a larger $q>1$ for ${\sigma_{\rm p}}>0.012\,{\rm mas}$.

To summarize this: The best choice for the sample is to add $\delta{\sigma_{\rm p}}\sim 0.043\,{\rm mas}$ in quadrature, and to add $\sim 1$ times the parallax error ${\sigma_{\rm p}}$ to the parallax. When concerned about precision, we further advise to remove all stars with ${\sigma_{\rm p}}>0.08\,{\rm mas}$, and recommend to separately test stars with ${\sigma_{\rm p}}\lesssim 0.047\,{\rm mas}$ for anomalies.

5.3 Bias vs. colour and magnitude

After these more complicated considerations, the derivation of safe limits in colour and magnitude for the sample is quite straight-forward. We simply order our sample by ${\rm G}$ or ${\rm G_{RP}}$-band magnitude (for an overview of the photometric instrumentation and calibration in Gaia, see Riello et al., 2018; Evans et al., 2018), as well as by the given ${\rm G_{BP}}-{\rm G_{RP}}$ colour, and then iteratively tighten the quality cuts. The results after a first round of clean-ups is shown in Fig. 9. The top panel shows the distance bias vs. ${\rm G}$ magnitude after removal of all stars with ${\rm G}>14.5\,{\rm mag}$ and ${\rm G_{BP}}-{\rm G_{RP}}>1.5\,{\rm mag}$, while the bottom panel shows a scan in ${\rm G_{BP}}-{\rm G_{RP}}$ colour, after removal of all stars with ${\rm G}>14.5\,{\rm mag}$. If we had not censored the faintest stars in the top panel, they would be shown with $1+f\sim 0.8$, i.e. a very strong distance underestimate. Since we recall line-of-sight velocity measurement errors (even if unbiased) look like distance underestimates, the by far most likely explanation for the decline in $1+f$ is not a failure of Gaia parallaxes, but a much larger than indicated error from the $v_{\rm los}$ measurements. If would remain to be explored if the mild decline around ${\rm G}>13\,{\rm mag}$ may also be related to a change of magnitude window in the Gaia astrometry. Of course, the change of colour range implies a change in the selection function. We thus devised an automated measurement of the selection function $S(s)$, where we use about $40$ base points for $s<4\,{\rm kpc}$, on which we re-measure $S(s)$ and then calculate a grid of relative factors between we interpolate linearly on $\log{s/\,{\rm kpc}}$. As the reader can easily see from the difference between the red and blue points in the top panel of Fig. 9, this more appropriate but far more costly calculation does not significantly change the results. The largest change is around ${\rm G}\sim 8\,{\rm mag}$, where the re-calculation shortens all distances a little, since $S(s)$ drops more steeply with $s$, and thus exacerbates a little the negative bias in this area. We also note that the relatively sharp and borderline significant feature just below ${\rm G}\sim 11\,{\rm mag}$ resembles a lot what was shown in the re-calibration comparisons in Lindegren et al. (2018). However, the sample size and small amplitude of the effect precludes further investigation.

5.4 Bias vs. astrometric parameters

After discussing the main problems with the sample, it is worth looking at the sample quality vs. the different astrometric parameters. Fig. 10 shows the distance bias against the number of visibility periods ${\rm n_{\rm vis}}$ (top panel) and the ecliptic latitude $\beta$ (bottom panel). In both the top and bottom panel we compare the results obtained with the simpler assumption of a constant parallax offset ${\delta_{\rm p}}=-0.048\,{\rm mas}$ (red lines) vs. setting ${\delta_{\rm p}}=-{\sigma_{\rm p}}$ and limiting ${\sigma_{\rm p}}<0.12\,{\rm mas}$ (blue lines). The ${\sigma_{\rm p}}$ limit does not significantly alter the results. We were first rather astounded by the strong dependence of the distance bias on ${\rm n_{\rm vis}}$ in the naive evaluation (red lines). However, most of this effect is readily explained by the dependence of ${\delta_{\rm p}}$ on ${\sigma_{\rm p}}$. First, we note that due to the scanning law of Gaia, ${\rm n_{\rm vis}}$ very strongly correlates with $\beta$, so both trends have a common explanation. The culprit is easily found, when we correct for the apparent dependence of ${\delta_{\rm p}}$ on ${\sigma_{\rm p}}$ (blue lines), which diminishes the trend. Consistently with the top panel, most of the dependence of $1+f$ on $\beta$ in the bottom panel disappears when we apply this correction, with a mild failure of order $f\sim 5\%$ remaining around the ecliptic south pole. This is far less concerning than it looks since the ecliptic poles imply a pencil beam, where our method loses most of its advantages, and furthermore, the ecliptic poles are in a near-worst case location, since they are close to the Galactic plane almost exactly in the azimuthal direction: this implies weak statistics, and only the $V-W$ correlation term, which carries some mild signal from the galactic warp, no information from the $U-W$ term. We thus ran a few checks. The distance method did not reveal anything else than near-perfect distance bias when scanning along Galactic ${l}$ and ${b}$ (which would have argued against a failure here), but in turn, we did not see the apparent dip in the bottom panel of Fig. 10 when limiting the sample to the more robust stars with low $V_{g}$. We further compared the locus of the lower main sequence in absolute G band magnitude for stars with $\beta\lesssim-55\,{\rm deg}$ with equivalent fields in Galactic latitude and longitude, finding no significant difference. To summarize all this paragraph: Assuming the dependence of $-{\delta_{\rm p}}\sim{\sigma_{\rm p}}$ makes all problems in ${\rm n_{\rm vis}}$ and $\beta$ go away, with some possibility for a problem at $\beta\lesssim-55\,{\rm deg}$, so in sensitive studies, one might want to touch these stars with a long pole.

Fig. 11 shows the distance bias $1+f$ vs. the astrometric excess noise (top panel) and vs. the RUWE (a quantity based on the chi squared of the astrometric fit and stellar colour). To create these plots we use the usual limits in colour and magnitude, and waived the excess noise limit in the top panel. We note that only a very minor fraction of our stars have a non-zero excess noise value. These stars were censored from the sample. One could expect that at least some stars with larger excess noise should be binaries. In the discussion of Schönrich & Aumer (2017), the Gaia DR1 sample had a far larger temporal baseline for the astrometry compared to the $v_{\rm los}$ measurements, which raises the expectation of an apparent distance underestimate in the statistical distance estimator, since the $v_{\rm los}$ measurements would carry the additional velocity dispersion from the binary. Here, the case here is less clear-cut. However, we can still note: stars with a positive excess noise have a slight tendency to distance underestimates in our method (understandable, since the astrometric effect is not clear-cut and we still have some binary dispersion affecting the estimator). If we admit stars with very large ${\sigma_{\rm p}}$ (which correlates with the excess noise), we find that the sample shows clear signs of a break-down for excess noise values larger than $\sim 1$. Thus, we recommend applying a cut at this value. Further we note, the excess noise measurement is limited to stars with very good signal to noise in Gaia and thus leads to a concentrated sample in distance; neglecting this effect should cause slight distance over-estimates in our selection, so the true distance bias on the large–excessnoise stars is likely slightly underestimated.

The bottom panel of Fig. 11 provides our two suggestions for distance evaluation vs. the RUWE as defined in the additional release notes from Gaia. The usual argument is that stars with a large RUWE value are bad and should not be used. However, as we see from the figure, the only trend of $1+f$ that we can detect is in the evaluation where we hold the parallax offset fixed at ${\delta_{\rm p}}=-0.048\,{\rm mas}$. However, large values for RUWE (which is an expression for the quality of the astrometric fit) duly correlate with a larger ${\sigma_{\rm p}}$ given by the Gaia astrometric pipeline. If we correct for the trend in the Gaia parallax offset (green line), this slight bias vs. RUWE vanishes. In short, we cannot see any reason for applying a quality cut in RUWE - of course, stars with large RUWE have worse measurements, but the Gaia pipeline appears to be perfectly fine in not producing any bias vs. RUWE and mirrors the larger uncertainties in larger ${\sigma_{\rm p}}$ values.

Fig. 12 finally shows $1+f$ vs. the BP/RP flux excess factor $E_{\rm bprp}$. The quantity is often cited as an important quality measure for Gaia data. It compares the flux in the BP+RP bands to the total flux in the Gaia G-band and relies on their very similar coverage. Due to the Gaia passband definitions, $E_{\rm bprp}$ increases on average for red stars, but primarily expresses contamination by neighbouring stars and background, or misidentifications. We note that here we only measure the average distance error, i.e. we are not concerned with single outliers. The result is quite clear-cut: both very low values of $E_{\rm bprp}<1.172$ and large values of $E_{\rm bprp}>1.3$ signal compromised stars. We also note, however, that we needed to scan the sample with a very fine sample size mask of $12000$ stars moved in steps of $4000$ stars each, since the number of compromised stars is so low in the sample with ${b}>10\,{\rm deg}$. The error distribution within the bulk of the sample shows the usual number of outliers. Scanning with a larger sample size (100000) on the centre of the distribution shows that all sub-samples there have $|f|<2\%$, with a slightly suspicious region almost exactly at $E_{\rm bprp}\sim 1.2$.

6 Comparison to previous distance derivations

It is appropriate to ask how our distances differ from the two previous distance derivations by Bailer-Jones et al. (2018) and McMillan et al. (2018), so Fig. 13 shows a distance scan of $1+f$ vs. their distance $s$. In both cases the distances were estimated under the assumption that the Gaia DR2 parallax zero-point is $-0.029\,{\rm mas}$, based on the quasar results from Lindegren et al. (2018), which was the best available estimate at the times they were computed. We have shown that this is an underestimate for these stars. It is therefore no surprise that Fig. 13 shows that both studies overestimated stellar distances substantially, with the typical increase of $1+f$ with distance. In short, due to the parallax offset these distances are compromised and should not be used. The reason why the bias in the Bailer-Jones et al. (2018) distances is somewhat smaller than in the McMillan et al. (2018) distances is predominantly caused by the fact that we do not have proper expectation values for the Bailer-Jones et al. (2018). Instead, we have only the mode (maximum) of their posterior distance probability distribution, which, due to the skew probability distribution, is significantly smaller, for usual distance PDFs, than the expectation value. While this compensates for some of their intrinsic distance overestimate, the mode infers a difficult-to-predict and variable bias relative to the expectation value of a distribution; also one should not rely on two opposite biases to partly cancel. A longer discussion of this issue can be found in Schönrich & Aumer (2017).

7 The distance to the Pleiades

Since there has been so much discussion about the distance to the Pleiades, let us quickly analyse the effect that the offset has on their distance estimate. Historically, there has for a long time been a tension between a low astrometric distance estimate by Hipparcos (van Leeuwen, 2009), placing the distance of the Pleiades at $s=120\,{\rm pc}$ with a formal uncertainty of $1.5\%$, and results from isochrone fits to their photometry (Meynet et al., 1993; Stello & Nissen, 2001), as well as eclipsing binaries (Zwahlen et al., 2004; Southworth et al., 2005), which placed their distance in the range $\sim(130-137)\,{\rm pc}$ with similarly small uncertainties. The Gaia DR2 release (Gaia Collaboration et al., 2018b) estimates the Pleiades’ distance at $135.8\pm 0.1\,{\rm pc}$. Applying our parallax offset to this estimate brings the distance moderately down to $134.8\pm 0.2\,{\rm pc}$, and back towards the average of the stellar-physics-based estimates.

8 Summary of suggested quality cuts

Since the previous discussion was rather lengthy with many details, we provide here our suggestions for quality cuts on the Gaia sample, to ensure a minimization of the detected kinematic biases

• •

A colour cut ${\rm G_{BP}}-{\rm G_{RP}}<1.5\,{\rm mag}$. To be entirely safe, we suggest $0.5<{\rm G_{BP}}-{\rm G_{RP}}/\,{\rm mag}<1.4$.

• •

A magnitude cut for ${\rm G}<14.5\,{\rm mag}$, and ${\rm G_{BP}},{\rm G},{\rm G_{RP}}>0$. A safer limit is ${\rm G}<12.5\,{\rm mag}$ and ${\rm G_{RP}}<13.7\,{\rm mag}$.

• •

${\rm p}/{\sigma_{\rm p}}>4$, safer is ${\rm p}/{\sigma_{\rm p}}>10$.

• •

${\sigma_{\rm p}}<0.1\,{\rm mas}$ with ${\sigma_{\rm p}}$ as given by the Gaia pipeline, safer is ${\sigma_{\rm p}}<0.07\,{\rm mas}$.

• •

${\rm n_{\rm vis}}>5$ as pointed out in Lindegren et al. (2018) and excess noise $<1$.

• •

For the BP-RP excess flux factor, use $1.172<E_{\rm bprp}<1.3$. Tighter cuts might apply if the number of outliers is important.

• •

$s>80\,{\rm pc}$ for studies that need assurance of distance systematics $<4\%$.

Some notes: The ${\sigma_{\rm p}}$ dependence of the parallax offset is comparably well controlled as long as we choose a safe limit on ${\sigma_{\rm p}}$. If one needs to use a looser limit on ${\sigma_{\rm p}}$, we strongly advise to control for the dependence of ${\delta_{\rm p}}$ on ${\sigma_{\rm p}}$. If needed, the upper limit on the excess flux factor $E_{\rm bprp}$ can be relaxed a bit with proper caution, see Fig. 12. As detailed in Section 5.4, we find no reason to test or cut in RUWE, when applying our set of quality cuts on this sample. Note further that the faint magnitude limit, and likely the colour limit are an imprint of the $v_{\rm los}$ measurement quality in the Gaia RV subset, i.e. these should be validated separately for different samples.

9 Conclusions

We have used the methods of Schönrich et al. (2012) and Schönrich & Aumer (2017) to derive bias-free expectation values for stellar distances in the RV subset of Gaia DR2.222Please find the datasets with the derived distances and simple estimates of stellar velocities and positions in Galactic coordinates either in the MNRAS online materials or at https://zenodo.org/record/2557803.. While Gaia parallaxes have been extensively tested in Lindegren et al. (2018), their quasar sample has no overlap with our study in magnitude, thus the quasar sample is inadequate for testing distances in the important Gaia RV subset for two reasons: quasars have all ${\rm p}=0$, i.e. differences in the astrometric pipeline for ${\rm p}>0$ cannot be tested, and basically none of their quasars shares the same evaluation cohort as the stellar sample.

We derived Bayesian distances for all stars in the subset and validated the distance expectation values for stars with ${\rm p}/{\sigma_{\rm p}}>4$, the minimum justifiable quality requirement for relative parallax error. All distances and derived kinematics will be made available with this work.

Our study provides clear proof for an average parallax offset ${\delta_{\rm p}}=-0.054\,{\rm mas}$ (Gaia parallaxes are too small) with negligible formal uncertainty and a systematic uncertainty of $\sim 0.006\,{\rm mas}$. The parallax offset is clearly identifiable as such as it results in an almost perfectly linear uptrend of distance bias with distance $s$ that reaches values in excess of $f>30\%$ for $s>3\,{\rm kpc}$. This offset is comparable to the findings of Zinn et al. (2018) using asteroseismic data. It is significantly larger than the value of $-{\delta_{\rm p}}=0.029\,{\rm mas}$ found by Lindegren et al. (2018) for quasars. This also advocates a re-analysis of cluster distances. Even the very nearby Pleiades are pulled back by this offset by $\sim 1\,{\rm pc}$. Similarly, even a comparably benign bias of $10\%$ creates larger deviations of mean velocities, than found e.g. for the warp/wave pattern in the local Galactic disc. Every study using Gaia DR2 parallaxes/distances should investigate the sensitivity of their results on the parallax biases described here and - for fainter samples - in the DR2 astrometry paper.

We evaluated different assumptions for the parallax error in the Gaia pipeline and found that our estimate for ${\delta_{\rm p}}$ is nearly unaffected by changing the parallax error. Not adding the additional error $\delta{\sigma_{\rm p}}=0.043$ to the astrometric pipeline value ${\sigma_{\rm p}}$ in quadrature, decreases our estimate for $-{\delta_{\rm p}}$ by about $0.006\,{\rm mas}$.

As we used a self-informed prior for the distance-dependent selection function $S(s)$, the method provides a good approximation for $S(s)$, which we provide in equation (6). Assuringly, $S(s)$ displays the behaviour expected from population synthesis and the magnitude limits of the RV sub-sample of Gaia: an almost perfectly exponential decrease for $s<1\,{\rm kpc}$ related to the main sequence, a knee at intermediate distances, when the magnitude cut passes the level of the subgiant branch, and a slower decrease of $S(s)$ towards large $s$.

After correction for a constant parallax offset we still found a highly significant correlation of the distance bias $f$ with ${\sigma_{\rm p}}$. While this could point to a problem with the assumptions in our Bayesian distances, this explanation is unlikely since $S(s)$ can be measured to high confidence from the data. To the contrary, our results suggest that $-{\delta_{\rm p}}$ is roughly proportional to ${\sigma_{\rm p}}$ (best-fit value $q=1.05$) after adding the $0.043\,{\rm mas}$ additional error to the Gaia parallaxes, and further show that stars with ${\sigma_{\rm p}}\gtrsim 0.1\,{\rm mas}$ should be discarded from analysis. It is unreasonable to think that the parallax uncertainty is added to the parallax value in a simple way, and so it is no surprise that the required factor is not constant in parallax: $q$ depends significantly on $p$ and around $p\sim 0.09\,{\rm mas}$, it is closer to $q=0.5$.

Resolving this dependency of ${\rm p}$ on ${\sigma_{\rm p}}$ also removes a highly significant trend of our measured distance bias with the number of visibility periods ${\rm n_{\rm vis}}$, and consequently with ecliptic latitude $\beta$. We note that our method would still flag a distance problem at the ecliptic south pole. However, when used on a very narrow area on the sky, we lose most of our statistical corrections and do not trust the evaluation. In fact consistent with this expectation, when limiting the sample to the more robust stars with low azimuthal velocity, this dependency was not confirmed, consistent with an evaluation of the derived HR diagram.

We further used the method to evaluate safe limits to be imposed both on apparent magnitude and stellar colour, finding that red stars at ${\rm G_{BP}}-{\rm G_{RP}}>1.5$ are compromised as well as stars with ${\rm G}\gtrsim 14\,{\rm mag}$ are flagged for distance underestimates. The most likely explanation is a decline in quality of the otherwise very well determined $v_{\rm los}$.

We further tested for astrometric parameters, finding no biases related to $RUWE$ (after removal of the aforementioned problems), and no strong correlation of $f$ with astrometric excess noise values smaller than $1$. At least for getting distance expectation values in the Gaia RV sample, this strongly argues for not using $RUWE$ as a quality indicator. A mild decrease in distance estimates could point to stellar binaries.

A summary of all quality cuts is provided in Section 8.

Acknowledgements

We thank our referee, U. Bastian, for very thorough and insightful comments to the paper. It is a pleasure to thank Lennart Lindegren, J. Magorrian, J. Binney, F. van Leeuwen, and A. Mora for helpful comments. RS is supported by a Royal Society University Research Fellowship. This work was performed using the Cambridge Service for Data Driven Discovery (CSD3), part of which is operated by the University of Cambridge Research Computing on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The DiRAC component of CSD3 was funded by BEIS capital funding via STFC capital grants ST/P002307/1 and ST/R002452/1 and STFC operations grant ST/R00689X/1. DiRAC is part of the National e-Infrastructure. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.