Estimating Distance from Parallaxes. IV. Distances to 1.33 Billion Stars in Gaia Data Release 2

Given the Gaussian likelihood in the parallax $\varpi$ with standard deviation ${\sigma }_{\varpi }$, plus the EDSD prior from Equation (1), the unnormalized posterior over the distance to a source is

The quantities $\varpi$, ${\sigma }_{\varpi }$, $l$, $b$ are all taken from the gaia_source table in GDR2 from the columns parallax, parallax_error, l, b respectively. ${\varpi }_{\mathrm{zp}}$ is the global parallax zeropoint, determined from Gaia's observations of quasars to be −0.029 mas (Lindegren et al. 2018). (The zeropoint actually varies as a function of sky position, magnitude, and color, but this is the best global estimate.) The length scale ${L}_{\mathrm{sph}}(l,b)$ has a Galactic longitude and latitude dependence as described in Section 2.3.

The properties of this posterior have been explored in some detail in Bailer-Jones (2015). For physical values of its three parameters—finite $\varpi$, positive ${\sigma }_{\varpi }$, positive ${L}_{\mathrm{sph}}$–it is always a proper (i.e., normalizable) density function. It can be either unimodal or bimodal, although for the vast majority of cases in GDR2 it is unimodal.

While the posterior in Equation (2) is the complete description of the distance to the source (fully specified by $\varpi$, ${\sigma }_{\varpi }$ in GDR2 and ${L}_{\mathrm{sph}}$ in our catalog), we sometimes want a point estimate along with some measure of the uncertainty. As known from previous work and seen from the examples in Section 3, this posterior is asymmetric, with a longer tail toward large distances. Being a proper (normalizable) posterior, all of the usual estimators—mean, median, mode, variance, quantiles, full-width at half-maximum—are defined. As the point estimator, ${r}_{\mathrm{est}}$, we prefer here the mode, ${r}_{\mathrm{mode}}$. This is found analytically by solving a cubic equation (Bailer-Jones 2015).

As a measure of uncertainty, we use the highest density interval (HDI) with probability p. The HDI is the span of the distance that encloses the region of highest posterior probability density, the integral of which is p. The span is defined by the lower and upper bounds, ${r}_{\mathrm{lo}}$ and ${r}_{\mathrm{hi}}$, respectively. Here, we set p = 0.6827, which is equal to the probability contained within $\pm 1\sigma$ of the mode for a Gaussian distribution. The distance posterior is asymmetric (sometimes significantly so) so ${r}_{\mathrm{hi}}$ − ${r}_{\mathrm{est}}$ and ${r}_{\mathrm{est}}$ − ${r}_{\mathrm{lo}}$ are unequal. Conceptually, the HDI can be found by lowering a horizontal line over the distribution until the area contained under the curve between its intercepts with the curve (i.e., ${r}_{\mathrm{lo}}$ and ${r}_{\mathrm{hi}}$) is equal to p (see Figure 5.10 of Bailer-Jones 2017). The HDI has the advantage over other confidence intervals in that it is guaranteed to contain the mode. Some other common measures are less flexible. The Fisher information approach, for example, assumes the posterior is locally Gaussian, the standard deviation of which is taken as the uncertainty (which is symmetric and therefore entirely inappropriate for positively constrained quantities like distances).

This HDI is unique if the posterior is unimodal. (A more general way of finding the HDI allows it to split into multiple intervals in the case of multimodality; see for example Hyndman 1996). For a small part of the parameter space defined by $(\varpi ,{\sigma }_{\varpi },L)$, the posterior is bimodal (see Bailer-Jones 2015). For these cases, we retain both the mode estimator and the HDI if possible (i.e., if the span does not include the minimum, thereby retaining the uniqueness; see below). Otherwise, we resort to using the median of the distribution as the point estimator, and report the 16th and 84th percentiles (i.e., $(1\pm p)/2$) as ${r}_{\mathrm{lo}}$ and ${r}_{\mathrm{hi}}$, respectively; together these two numbers form the equal-tailed interval (ETI), which has as much probability below the span as above, with p in between.

There is no analytic solution for the HDI for this posterior. We compute it with an iterative procedure involving a Taylor expansion. The basic idea is to take a small step in both directions away from the mode, compute the area under the curve covered by this step, and iterate this until the total area hits the desired limit.

We first normalize the posterior using Gaussian quadrature and denote this with $P(r)$. With $r$ set equal to the mode (where the first derivative is zero), and adopting a fixed negative ${\rm{\Delta }}P$, an initial step of size

is computed. Candidate lower and upper bounds are then defined as

The area under the curve between these bounds is computed using the trapezium rule. Assuming this area is less than p (${\rm{\Delta }}P$ is set small enough to ensure this is always the case), further steps are computed in both the negative and positive directions using the first-order term in the Taylor expansion, i.e.,

starting from i = 1, to give new candidate lower and upper bounds as

The additional area covered by these steps is estimated using the trapezium rule and added to the running total area. This process is iterated from Equation (5), with the result that both bounds move monotonically away from the mode. The procedure is stopped when the running total area computed exceeds p (or bimodality is detected; see below). For the distance catalog produced here, ${\rm{\Delta }}P$ is fixed to $-0.01P({r}_{\mathrm{mode}})$. For 98% of the sources in GDR2, between 39 and 73 iterations (the 1st and 99th percentiles) are required.

Note that the area computed by this method actually exceeds p, although generally by only a very small amount: for 99.9% of the sources, the area computed is still less than 0.6973. This is an acceptable degree of approximation. (Tests shows that the area calculation by the trapezium rule is generally very accurate: less than 0.06% of cases deviate from an accurately calculated area by more than 1 part in 1000.) In a few pathological cases, the posterior has a long, nearly flat tail extending to large distances. The numerical estimation of the HDI is then poor. In a very small fraction of cases, this can even lead to the computed area exceeding 0.9. If this happens, the mode and HDI are rejected as catalog entries and replaced with the median and equal-tailed quantiles, as explained below. (These solutions have result_flag=2 and modality_flag=1. A full description of the flags is given in the next section.)

Because the steps in the positive and negative directions are made independently, the values of the posterior at ${r}_{\mathrm{lo}}$ and ${r}_{\mathrm{hi}}$ are not identical: an imaginary line being lowered over the posterior does not remain exactly horizontal. Because P = 0 only for $r=0$ and $r=+\infty$ (Equation (2)), and because ${r}_{\mathrm{hi}}$ can never reach infinity, then in principle, ${r}_{\mathrm{lo}}$ should never reach exactly zero. With this numerical method, it can reach zero, although it only happens in 183,765 (0.014%) cases. The algorithm prevents steps to negative distances, so ${r}_{\mathrm{lo}}\geqslant 0$. (Note that this method of finding the HDI would not work without modification for distributions which drop to zero at finite values of r greater than zero.)

If the posterior is bimodal (which occurs in 0.09% of cases), then the search for the HDI is done starting from the highest mode. Provided the search does not encounter the minimum between the two modes, the HDI remains well defined, and this, together with the highest mode, are the posterior summary provided in the catalog. The result_flag has value 1 whenever the point estimate, ${r}_{\mathrm{est}}$, is the mode, and the lower and upper bounds of the confidence interval, ${r}_{\mathrm{lo}}$ and ${r}_{\mathrm{hi}}$, define the HDI.

The presence of bimodality (which is identified independently of the search, from the roots of ${dP}/{dr}=0$) is always indicated by the modality_flag, which is 1 for unimodal and 2 for bimodal. The minimum is encountered during the HDI search in less than 0.1% of all cases. If this happens, we stop the search and instead compute the median (the 50th percentile) and the ETI, which we report in the catalog. Such solutions are indicated with result_flag=2. The quantiles are found numerically by drawing 2 × 10⁴ samples from the posterior with a Markov Chain Monte Carlo (MCMC) method. Depending on the shape of the posterior, these estimates of the quantiles may not be very accurate, so should be treated with caution.

The HDI is a couple of orders of magnitude faster to compute than the ETI, because the former involves that many fewer calculations of the posterior density.

In a very small number of cases (about two in every million), the algorithm fails to give any summary of the posterior. In these cases, result_flag=0 and the three distances ${r}_{\mathrm{est}}$, ${r}_{\mathrm{lo}}$, and ${r}_{\mathrm{hi}}$ are set to NaN. This occurs when the parallax is negative yet highly statistically significant, ${\sigma }_{\varpi }/\varpi \ll -1$, which results in the unnormalized posterior density being numerically equal to zero everywhere on account of the finite machine precision. This precludes both the computation of the normalization constant (for the HDI) and an MCMC sampling (for the quantiles). Although we can still compute the mode (it is found algebraically), we choose not to report it. One could work around this problem by including scaling factors in the computation of the density, but the cases are very rare, plus the distance posterior is virtually identical to the prior, which is entirely specified by the length scale (always provided).

We publish only a single point estimate together with the bounds of a single confidence interval in order to keep the catalog manageable in size. Should users want other estimates, e.g., different size confidence intervals, it is quick and easy to recompute the posterior using $\varpi$ and ${\sigma }_{\varpi }$ in GDR2 and ${L}_{\mathrm{sph}}$ in our catalog. Processing all 1.33 billion sources took 75 CPU-core days, which on our multicore machine was just a few hours of wall-clock time. The code for the processing is available at https://github.com/ehalley/Gaia-DR2-distances.

We use a chemo-dynamical model to simulate the Galaxy (including extinction) as seen by Gaia, which here means all stars with apparent magnitude $G\leqslant 20.7$ mag. The resulting mock catalog and the Galaxy model on which it is based are described in Rybizki et al. (2018). From this catalog, we extract the distances to all stars in each of the 49152 equal-sized (0.84 square degrees) healpix cells (level 6) over the entire sky. We then fit the EDSD prior, Equation (1), to the stars in each cell. The maximum likelihood estimate for $L$ is just $\overline{r}/3$, where $\overline{r}$ is the mean of the distances in that cell. To avoid a few stars at very large distances dominating this estimate, we replace the mean with the median.

The resulting map, ${L}_{\mathrm{gal}}(l,b)$, is shown in Figure 1. The bulge, disk, and warp of the disk are quite prominent. Perhaps counter-intuitively, the median distance is generally smaller at high Galactic latitudes than at low latitudes. This is because of the very high density of stars far away in the Galactic disk which Gaia can see, despite the extinction. Although the maximum distance Gaia can observe to is larger at high latitudes, the density of star drops rapidly with height above the disk, so the median is determined primarily by nearer stars. The same map computed using $\overline{r}/3$ exhibits larger distances than this map at high latitudes.

Zoom In Zoom Out Reset image size

We do not use this map directly as the prior model, because it shows discontinuities across the healpix boundaries. This would lead to spatial discontinuities in the inferred distances. We therefore fit a spherical harmonic model of order n_max

where ${P}_{n}^{m}()$ is the associated Legendre polynomial of order (n,m), and $\{{s}_{n,m},{c}_{n,m}\}$ together with a are the free parameters to be determined. Note as ${P}_{n}^{m}(\sin b)\sin ({ml})=0$ when m = 0, the constants ${s}_{n,0}$ are unconstrained, so we set them to zero. The total number of free parameters is

We fit this model to the map of ${\mathrm{log}}_{10}L$ using linear least squares. After some experimentation, also with the healpix level of the input map, we settled on ${n}_{\max }=30$ as high enough to capture the essential spatial variations, but low enough to smooth out some of the rather sharp features in the model. The fitted map, which we denote ${L}_{\mathrm{sph}}(l,b)$, is shown in Figure 2. A histogram of the residuals of the fit is shown in Figure 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The distribution of the parallaxes is shown in Figure 4 (compare with the predictions in Figure 2 of Astraatmadja & Bailer-Jones 2016a). 24.9% of the sources have negative parallaxes. The distribution of the fractional parallax uncertainty, defined as ${\sigma }_{\varpi }/\varpi$, is shown in panel (b). Immediately striking is how many sources have very low S/N parallaxes: of those sources with positive parallaxes, only 15.7% have ${\sigma }_{\varpi }/\varpi \lt 0.2$, i.e., have an S/N greater than 5 (7.7% greater than 10, 0.98% greater than 50). These negative and/or low S/N parallaxes are due to the numerical dominance of relatively distant and faint sources. Cumulative distributions for ${\sigma }_{\varpi }/\varpi$ are shown in panels (c) and (d) for positive (black curve) and negative (red curve) parallaxes.

Zoom In Zoom Out Reset image size

Figure 5 shows examples of several posteriors and their corresponding priors. (More examples can be found in Figure 12 of Bailer-Jones 2015.) Recall that the mode of the posterior is at $2{L}_{\mathrm{sph}}$. We do not plot the likelihood because it is an improper distribution in $r$, so can only be plotted with an arbitrary scaling. The exact shape of the posterior depends, in a somewhat complex manner, on the value of all three of its variables: $\varpi$, ${\sigma }_{\varpi }/\varpi$, and ${L}_{\mathrm{sph}}$. These examples have been chosen to show a range of behaviors; they do not represent the frequency of such posteriors among the results. Panels (a)–(c) are for small positive values of ${\sigma }_{\varpi }/\varpi$: these posteriors are dominated by the likelihood (as opposed to the prior) and are therefore approximately Gaussian with a mode close to the inverse parallax. Panels (d) and (e) are cases where the likelihood and prior have similar size influence on the posterior. As ${\sigma }_{\varpi }/\varpi$ increases, the positive distance tail extends and the posterior becomes increasingly asymmetric. Panel (d) is a relatively rare example of a posterior with a high, narrow peak accompanied by a long tail.

Zoom In Zoom Out Reset image size

Broadly speaking, the larger ${\sigma }_{\varpi }/\varpi$ is, the more the prior dominates the posterior. This is seen in panels (f) and (h). However, it is not simply the case that the larger the value of ${\sigma }_{\varpi }/\varpi$, the closer the posterior is to the prior. We see this in panel (g), where ${\sigma }_{\varpi }/\varpi$ = 2.64 (larger than in panel f), yet the posterior mode is at more than twice the value of the prior mode. The posterior is always the (normalized) product of prior and likelihood, so if these peak at very different values, as is the case in panel (g) ($2{L}_{\mathrm{sph}}=1.0$ kpc versus $1/\varpi =10.4$ kpc), the likelihood can still play a significant role, even when the data are poor (i.e., ${\sigma }_{\varpi }/\varpi$ is large). In other words, prior dominance does not increase monotonically with ${\sigma }_{\varpi }/\varpi$. The actual values of the parallax and the prior length scale are also important.

Panel (i) is an example of a negative parallax. Panel (j) shows the relatively rare case of a bimodal posterior. In this particular case, the HDI cannot be uniquely defined, so our algorithm returns the median as the point estimate, and the equal-tailed quantiles to define the asymmetric confidence interval.

Figure 6 summarizes the properties of the distance estimates for unimodal posteriors with mode and HDI estimates. We discuss this figure for most of the rest of this section.

Zoom In Zoom Out Reset image size

Panel (a) shows the distribution of the distance estimate ${r}_{\mathrm{est}}$, with the distribution of the length scale, ${L}_{\mathrm{sph}}$, in panel (d) beneath it for comparison. A characterization of the fractional distance uncertainties is shown by the histogram in panel (b), where $({r}_{\mathrm{hi}}-{r}_{\mathrm{lo}})/2$ is taken as a single symmetrized measure of the uncertainty. This can be compared with the distribution of the fractional parallax uncertainties shown in Figure 4(b). The fractional distance uncertainties are generally smaller because the prior prevents the posterior distribution becoming arbitrarily wide.

Panel (c) indicates the skewnesses of the posteriors. The upper part of the confidence interval (${r}_{\mathrm{hi}}$ − ${r}_{\mathrm{est}}$) is always larger than the lower part (${r}_{\mathrm{est}}$ − ${r}_{\mathrm{lo}}$). In 2%, 0.6%, and 0.15% of the cases, it exceeds it by factors of more than 5, 10, and 20, respectively. As mentioned earlier, if the asymmetry is so large that the positive tail is nearly flat, then the HDI may not have been estimated very accurately.

Panel (e) shows how the fractional distance uncertainty varies with distance. The complex shape is a consequence of the true distribution in the Galaxy, the variation of ${L}_{\mathrm{sph}}(l,b)$, and the properties of the posterior.

The relation between the fractional parallax uncertainty and distance uncertainty is shown in panel (f). For small (positive) values of both parameters, they are highly correlated. For large positive values of ${\sigma }_{\varpi }/\varpi$, the fractional distance uncertainty can never get very large, due to the prior. For the prior itself, $({r}_{\mathrm{hi}}-{r}_{\mathrm{lo}})/(2{r}_{\mathrm{mode}})\simeq 0.74$ (independent of $L$), and this is the value that the upper envelope of the distribution asymptotes to as ${\sigma }_{\varpi }/\varpi \to \infty$. For some intermediate values of ${\sigma }_{\varpi }/\varpi$, the fractional distance uncertainties are quite large. These are due to posteriors with long positive tails like that shown in Figure 5(d): ${r}_{\mathrm{hi}}\gg {r}_{\mathrm{est}}$ which results in $({r}_{\mathrm{hi}}-{r}_{\mathrm{lo}})/2{r}_{\mathrm{est}}$ being large. If we used the median or mean as a distance estimate instead, this measure of the fractional distance uncertainty would not extend to such large values (although those estimators would introduce other problems).

Panel (g) plots the ratio of the mode of the posterior (${r}_{\mathrm{est}}$) to the mode of the prior ($2{L}_{\mathrm{sph}}$) as a function of the fractional parallax uncertainty, ${\sigma }_{\varpi }/\varpi$, to give some sense of how dominant the prior is. The posterior is a (renormalized) product of the prior and likelihood, so the location of the posterior mode depends not only on the modes of the prior and likelihood—the latter is at $1/\varpi$ when considered as a function of distance—but also on their widths, i.e., how peaked/flat the distributions are. This generates a rather complex behavior in panel (g). For small positive values of ${\sigma }_{\varpi }/\varpi$ (below 0.1–0.2), we see a wide range of values of ${r}_{\mathrm{est}}/2{L}_{\mathrm{sph}}$, as we would expect: there is little constraint from the prior. As ${\sigma }_{\varpi }/\varpi$ increases, the prior plays more of a role and so the posterior mode deviates from it by less. We see a sharp lower boundary in the scatter plot, showing that ${r}_{\mathrm{est}}$ cannot drop below some finite multiple of the prior mode. This is due to the shape of the prior density function, which always drops to zero at zero distance: for a likelihood of finite width (i.e., ${\sigma }_{\varpi }/\varpi \gt 0$), its product with the prior always produces a mode above zero. The upper envelope, in contrast, is blurred because there is no hard upper limit on the location of the prior mode.

A comparison of the distances to the naive inverse parallax estimates is shown in panel (h) by plotting their ratio (i.e., ${r}_{\mathrm{est}}\varpi$), and zoomed in for the higher-precision parallaxes in panel (i). Once ${\sigma }_{\varpi }/\varpi$ exceeds around 0.1–0.3, the inverse parallax deviates considerably from ${r}_{\mathrm{est}}$. For large fractional parallax uncertainty, ${r}_{\mathrm{est}}$ tends to be less than $1/\varpi$. This is because a noisy inverse parallax can extend to arbitrarily large distances, something that the prior constrains. The inverse parallax is known to be a biased and very noisy estimator (e.g., Bailer-Jones 2015; Luri et al. 2018), which tends to overestimates the true distance.

The left panel of Figure 7 shows the distribution of estimated distances in Galactic coordinates, for those stars where the posterior is unimodal and this mode is the distance estimator. Compare this to the right panel, which shows the distribution of the modes of the corresponding priors. The overall features of the prior remain evident in the left panel, partly because such features correspond to real, previously known features in the Galaxy, but also because most stars have large fractional parallax uncertainties. Yet we clearly see features that are not present in the prior map, such as the Large and Small Magellanic Clouds. Our mean distances to these satellite galaxies are of course underestimated, because the stars have very large fractional parallax uncertainties, so our estimates are prior-dominated. The largest length scales in the prior model are only about 2.6 kpc—corresponding to prior modes of 5.2 kpc—whereas these galaxies are in fact about 50–60 kpc away. The same problem applies to distant giants, because the prior will normally be dominated by the nearer dwarfs in the model.

Zoom In Zoom Out Reset image size

Figure 8 compares the projected distribution of stars in the Galactic plane from our distance estimates (left) with those obtained from a naive parallax inversion (right). This shows the tendency of inverse parallax to extend to implausibly large distances. The radial lines from the Sun in both panels are due to nearby dust clouds diminishing the number of stars observed along certain sight-lines.

Zoom In Zoom Out Reset image size

A simple validation of our distance estimates can be obtained using stellar clusters, on the assumption that the members of a physical (gravitationally bound) cluster span just a small range of distances.

Given a list of clusters with known sky positions, we select cluster members based on their proper motions. That is, we assume all members follow a common space motion (with some dispersion) that is different from the bulk of the field stars in the same region. We model the proper motion plane using Gaussian mixture models for both the cluster region and a field region. This allows us to assess the membership probability to the cluster on a star-by-star basis. For present purposes, we focused on purity by selecting only the most probable members. This method gives very a similar selection of members to that found and reported in Gaia Collaboration (2018b).

Figure 9 shows some results of validating our distance estimates using clusters. Each panel shows our distances to individual cluster members as a function of the fractional parallax uncertainty. (In all cases the estimates turn out to have result_flag=1 and modality_flag=1.) These distances can be compared to the inverse of the mean parallax of the cluster members, shown as the blue line. This mean is computed using all members regardless of ${\sigma }_{\varpi }/\varpi$ (including negative parallaxes if present). Although the parallax uncertainties are correlated on the angular scale of most stellar clusters ($\lesssim 1^\circ$), these correlations need not be considered when forming a simple mean. They would need to be accommodated if we estimated the formal uncertainty in this mean. Yet due to the large number of cluster members, the uncertainty in the mean is instead limited by the systematic parallax error, which is generally below 0.1 mas (Lindegren et al. 2018). Thus, for clusters within one kpc or so, the inverse mean parallax is a reasonable approximation of the cluster distance. We are not suggesting that this is the best way to estimate true cluster distances. We are just using it to give some baseline against which to compare our distance estimates (and uncertainties therein).

Zoom In Zoom Out Reset image size

Figure 9 shows that for nearby clusters (first five panels), our distance estimates for the individual stars agree with the cluster mean. This is as expected, because for nearby clusters, all of the members have small ${\sigma }_{\varpi }/\varpi$, and for small ${\sigma }_{\varpi }/\varpi$, the posterior mode is similar to inverse parallax. In all of these cases, the prior has negligible impact (the modes of the priors are much larger than the estimated distances). These panels also show that the estimated uncertainties (the error bars plotted) are reasonable, considering that clusters have a finite depth.

As we move to more distant clusters, we see some informative behavior. First, we see that some cluster members have large ${\sigma }_{\varpi }/\varpi$. As discussed earlier, the posteriors for stars with large ${\sigma }_{\varpi }/\varpi$ are often dominated by the prior (e.g., Figure 5(f)). This is seen nicely in the panel for NGC 2360, where the estimates move up toward the prior for large ${\sigma }_{\varpi }/\varpi$. The error bars also get larger. The next cluster, NGC0884, does not show such a clear trend, and in fact at large ${\sigma }_{\varpi }/\varpi$, the distance estimates are systematically larger than the prior. This is an example where even for poor parallaxes the posterior is not entirely dominated by the prior, because the likelihood, even though broad, suggests quite a different distance than the prior does. This was also seen in Figure 5(g). The final two clusters are dominated by stars with large ${\sigma }_{\varpi }/\varpi$. Here, the posteriors tend to be closer to the priors. Both clusters are so distant that many of their members have negative parallaxes. Yet our inferred distances to these are broadly consistent with the mean cluster parallax. Large uncertainties notwithstanding, this demonstrates that useful information can be extracted from negative parallaxes. Of course for these more distant clusters, the mean cluster parallax has a significant systematic error (order 0.1 mas), as do the individual parallaxes used for our distance estimates.