Estimating Distances from Parallaxes. VI. A Method for Inferring Distances and Transverse Velocities from Parallaxes and Proper Motions Demonstrated on Gaia Data Release 3

C. A. L. Bailer-Jones

doi:10.3847/1538-3881/ad08bb

1. Introduction

Precise parallaxes may be simply inverted to get precise distance estimates of individual stars. This is not possible for low-precision parallaxes, however, because the inversion induces a disproportionately large random error and bias (Bailer-Jones 2015; Luri et al. 2018). Of the 1.47 billion stars with parallaxes in Gaia Data Release 3 (GDR3; Gaia Collaboration et al. 2023), 87% have a parallax S/Ns less than five, and 78% less than three. Even in the final Gaia data release the parallax S/Ns will only improve by a factor of two for most stars, leaving 75% still with parallax S/Ns less than five (60% less than three).

A more considered approach than parallax inversion is therefore required for distance estimation. This has been framed as a probabilistic inference (Bayesian) problem in the present series of papers. In the first paper (Bailer-Jones 2015), I looked at the general problem of inferring distances from parallaxes, and introduced the exponentially decreasing space density (EDSD) prior. The second paper (Astraatmadja & Bailer-Jones 2016a) introduced a more sophisticated distance prior based on a model of the Milky Way. Both this and simpler isotropic priors were applied in the third paper (Astraatmadja & Bailer-Jones 2016b) to estimate distances for all two million stars in GDR1 that had parallaxes. Paper IV (Bailer-Jones et al. 2018) built a direction-dependent EDSD prior for the whole Galaxy using a mock Gaia catalog and used this to infer geometric distances for all 1.33 billion stars in GDR2 with parallaxes. Bailer-Jones et al. (2021), the fifth in the series, did this again for the higher-precision parallaxes in GeDR3 using a more flexible direction-dependent distance prior (the generalized gamma distribution). This paper also extended the method to take advantage of the Gaia color and apparent magnitude of each star, which, together with a model for the color–absolute magnitude diagram (which also varied over the sky), gives information on distance (as already explored in Astraatmadja & Bailer-Jones 2016a). These distances, which we called photogeometric distances, were published for 1.35 billion stars, along with geometric distances for 1.47 billion stars.

Several other authors have published distance estimation methods or catalogs using the Gaia parallaxes. Some focus on specific objects, such as Wolf–Rayet stars (Rate & Crowther 2020) or Mira variables (Sanders 2023) or stellar clusters (Cantat-Gaudin et al. 2020; Olivares et al. 2020). Others have performed large-scale distance estimation of broad types of individual stars for which more general priors must be used. Both Fouesneau et al. (2022) and Anders et al. (2022) combined Gaia parallaxes with optical and infrared photometry from multiple surveys to infer intrinsic stellar properties, line-of-sight extinctions, and distances for over a hundred million stars. Andrae et al. (2023) used the Gaia BP/RP spectra and parallaxes to also estimate stellar parameters, extinctions, and distances, and published these as part of GDR3 for half a billion stars. Using the publicly released subset of the GDR3 BP/RP spectra, Zhang et al. (2023) also estimated stellar parameters (albeit parallax rather than distance, so the inversion problem remains). Although conceptually similar to the approach used to produce photogeometric distances in Bailer-Jones et al. (2021), these large-scale distance estimation projects differ in their methods, and by using the spectral energy distribution make stronger assumptions about the stars themselves.

Stellar proper motions can provide further distance information. A star moving with an angular velocity μ at distance r has a transverse velocity of

$\begin{eqnarray}&&\displaystyle \frac{v}{\mathrm{km}\,{{\rm{s}}}^{-1}}=k\displaystyle \frac{\mu }{\mathrm{mas}\,{\mathrm{yr}}^{-1}}\displaystyle \frac{r}{\mathrm{kpc}}{\rm{with}}\ k=4.740471.\end{eqnarray} \tag{ 1 }$

Assuming stars have a limited range of plausible velocities, a measurement of the proper motion puts a constraint on the plausible distance. A velocity prior quantifies what we mean by "plausible" and inference delivers a posterior over distance. Alone, this would provide a rather poor distance estimate, but when combined with a parallax measurement, we achieve a kinegeometric distance estimate that is potentially more accurate than a purely geometric one. This is of particular interest for distant stars in Gaia, as these tend to have very low-S/N parallaxes but still relatively accurate proper motions. I explore this approach in this paper.

The use of proper motions to aid distance estimation is not new (see Binney & Merrifield 1998, Section 2.2, for an overview). A simple case is to adopt a model of circular rotation in the Galaxy (a rather narrow prior), and to estimate the distance directly from the proper motion (e.g., Reid 2022). Proper motions have long been used to estimate the mean distance to a group of stars, often under the name of statistical parallaxes, by measuring the radial velocities as well as the proper motions (e.g., Strömberg 1936; Hemenway 1975). Schönrich et al. (2019) adopted a related approach to reduce the bias in parallax-based distance estimates, and they used this to estimate distances to seven million stars with RVs in GDR2. Zari et al. (2021) used a similar approach to the one in the present paper to estimate distances to OBA stars in the Galactic disk from GeDR3. As the target population was narrowly defined, the distance and velocity priors could be tailored to the expected distribution for young stars out to 3–4 kpc.

In the present paper, I use a direction and distant-dependent velocity prior derived from a Galaxy model to estimate kinegeometric distances to a random subset of stars from the entire GDR3 catalog (Gaia Collaboration 2022). I compare the results to geometric distances, and assess the performance of both types of estimate using a mock catalog. The kinegeometric method also provides estimates of the two transverse velocity components. Radial velocities are not used: When I refer to "velocity" in this paper, it means the total transverse velocity, not the 3D velocity. In contrast to the previous three papers in this series, I do not publish a catalog, because as we shall see, the improvements from using proper motions are substantial only in a limited part of parameter space.

2. Method

The prior is constructed for each HEALpixel (Górski et al. 2005) independently. As in Bailer-Jones et al. (2021), I use the equatorial, nested scheme at level 5, giving 12,288 HEALpixels with an area of 3.36 square degrees each. All sky plots in this paper are a Mollweide equal-area projection in Galactic coordinates, with north up, east increasing to the left, and (l, b) = (0°, 0°) in the center. They were plotted using healpy (Zonca et al. 2019) and show grid lines at 45° in latitude and longitude. All other plots, as well as the processing and analysis, were done in R (r-project.org).

2.1. Distance and Transverse Velocity Posterior

Bailer-Jones et al. (2021) estimated geometric distances using a prior that varied only with distance (r) for each HEALpixel (p). The product of this prior and the likelihood is the unnormalized posterior probability distribution function (PDF)

$\begin{eqnarray}&&{P}_{{\rm{g}}}^{* }(r\,| \,\varpi ,{\sigma }_{\varpi },p)=\ P(\varpi \,| \,r,{\sigma }_{\varpi })\,P(r\,| \,p).\end{eqnarray} \tag{ 2 }$

The likelihood is a Gaussian in parallax (ϖ) with uncertainty σ_ϖ, as defined in Section 2.2 of Bailer-Jones et al. (2021). For comparison purposes, below I report geometric distances inferred from this posterior. These use the same input data, parallax zero point, and prior as in Bailer-Jones et al. (2021), so they are statistically identical.

Here, I introduce the kinegeometric distance, which uses the proper motions (μ_α*, μ_δ) in addition to the parallax. Proper motions are connected to the distance via Equation (1), so to exploit them we must introduce both a prior on the velocities (v_α*, v_δ) and a likelihood on the proper motions. Combining all this, we get the 3D kinegeometric posterior PDF over distance and velocity

$\begin{eqnarray}&&\begin{array}{l}{P}_{\mathrm{kg}}* (r,{v}_{\alpha * },{v}_{\delta }\,| \,\varpi ,{\mu }_{\alpha * },{\mu }_{\delta },{\rm{\Sigma }},p)\\ \quad =\,\mathop{\underbrace{P(\varpi ,{\mu }_{\alpha * },{\mu }_{\delta }\,| \,r,{v}_{\alpha * },{v}_{\delta },{\rm{\Sigma }})}}\limits_{\mathrm{likelihood}}\ \mathop{\underbrace{P(r,{v}_{\alpha * },{v}_{\delta }\,| \,p)}}\limits_{\mathrm{prior}},\end{array}\end{eqnarray} \tag{ 3 }$

with the * symbol indicating that the PDF is unnormalized. The likelihood is a 3D Gaussian with mean

$\begin{eqnarray}&&\left(\displaystyle \frac{1}{r},\displaystyle \frac{{v}_{\alpha * }}{{kr}},\displaystyle \frac{{v}_{\delta }}{{kr}}\right)\end{eqnarray} \tag{ 4 }$

and covariance Σ, which is the 3 × 3 variance–covariance matrix of the parallax and proper motions published in GDR3 for each star (Lindegren et al. 2021b; Gaia Collaboration et al. 2023). The form of the likelihood shows that the more precise the measured proper motion, the more precise the estimated value of v/r. As with the geometric distances, I use the parallax zero point of Lindegren et al. (2021a). I make no modifications to the proper motions. The prior will be introduced below. From the 3D posterior, we can estimate the distance and the transverse velocity components by marginalization. As I numerically sample the posterior (Section 2.3), this marginalization is trivial.

The main goal of this work is to assess the improvement in distance estimates by using the proper motions in addition to the parallaxes, but I will also investigate the velocities estimated by Equation (3). Transverse velocities can alternatively be obtained by multiplying the proper motions by the geometric distance. These I refer to as geometric velocities, and I will compute them for comparison purposes. I do not consider estimating velocities by dividing the proper motion by the parallax: this is very biased, and it fails completely when the parallax is not positive, which is the case for 24% of sources with parallaxes in GDR3.

2.2. Velocity Prior

The distance–velocity prior in Equation (3) can be decomposed without loss of generality as

$\begin{eqnarray}&&P(r,{v}_{\alpha * },{v}_{\delta }\,| \,p)=P({v}_{\alpha * },{v}_{\delta }\,| \,r,p)\ P(r\,| \,p).\end{eqnarray} \tag{ 5 }$

For the second term—the distance prior—I adopt the same one used for the geometric distances. This is described in Bailer-Jones et al. (2021, Section 2.3). It was built from the mock Gaia catalog of Rybizki et al. (2020) and is summarized in Figure 1. The first term in Equation (5), the velocity prior, I compute here from the same mock catalog, which is based on the Galaxy model of Robin et al. (2003) plus some later modifications. This model comprises several Galactic populations—bulge, thin and thick disks of a range of ages, stellar halo, dark matter halo, and interstellar medium—with velocities computed self-consistently with the mass distribution via age–velocity relations. The disk model includes a warp and flare and the outer bulge is strongly oblate (a bar). I exclude all star clusters.

**Figure 1.** Median of the distance prior in each HEALpixel, computed from the mock catalog of Rybizki et al. (2020).
Download figure:
Standard image High-resolution image

Stellar velocity distributions vary considerably with direction and distance, so a dependence on both is retained. The prior should be a continuous—and ideally, smooth—function of distance and direction. Here, I compute it independently for each HEALpixel (direction), so it is not smooth across HEALpixel boundaries, but it is smooth in distance. For a given HEALpixel, I investigated how the velocity distribution varied over distance shells, a shell being defined as the volume between two spheres of different radii centered on the Sun. With sufficient stars, the 2D distribution in each shell could be approximated well with a 2D Gaussian.

The task of building a prior then falls to determining how the mean and variance–covariance of this 2D Gaussian should vary with distance. Figure 2 supports the following explanation. I first construct a series of contiguous shells moving outward from the Sun, defined such that the volume of each shell is a fixed multiple λ of the volume of the previous shell. When the outer radius of the first shell is fixed (its inner radius is zero), this recurrence relation defines all shells. If stars were uniformly distributed in volume, then λ = 1 would produce shells with the same number of stars.

**Figure 2.** Velocity prior for HEALpixel 6200 at (l, b) = (2857, 348). The top left panel shows the positions of the midpoints of the shells and the number of sources in them (on a log scale). The two central panels show the distribution of the velocities (top v_α*, bottom v_δ) for some selected shells, with the colors denoting the shell as in the top left panel. The dashed lines show the fitted Gaussians. The two right panels show the values of the fitted means (upper) and fitted standard deviations (lower) for all shells, as red circles for v_α* and as blue triangles for v_δ. The lines show the spline fits, which extend to the outer radius of the last shell, here 7597 pc. The bottom left panel shows the fitted correlations of the Gaussians in each shell along with their spline fits, but these are not used (the correlation is set to zero in the model).
Download figure:
Standard image High-resolution image

fdg — **Figure 2.** Velocity prior for HEALpixel 6200 at (l, b) = (2857, 348). The top left panel shows the positions of the midpoints of the shells and the number of sources in them (on a log scale). The two central panels show the distribution of the velocities (top v_α*, bottom v_δ) for some selected shells, with the colors denoting the shell as in the top left panel. The dashed lines show the fitted Gaussians. The two right panels show the values of the fitted means (upper) and fitted standard deviations (lower) for all shells, as red circles for v_α* and as blue triangles for v_δ. The lines show the spline fits, which extend to the outer radius of the last shell, here 7597 pc. The bottom left panel shows the fitted correlations of the Gaussians in each shell along with their spline fits, but these are not used (the correlation is set to zero in the model).
Download figure:
Standard image High-resolution image

After experimenting to achieve a trade-off between having sufficiently small shells to capture the variation of the velocity distribution, yet sufficiently large shells to allow an accurate Gaussian fit, I selected λ = 1.3 and an initial outer shell radius of 200 pc. I generate a set of 45 shells such that the outer radius of the last shell is 15.3 kpc. Stars beyond this do not contribute to the prior. Gaia astrometry on more distant stars is poor, so larger distances are not worth accommodating. Some HEALpixels at high Galactic latitudes still have too few stars, so I merge successive shells (starting from the Sun) to attain at least 200 stars per shell. At higher latitudes, this leads to the outermost shell having an outer radius less than 15.3 kpc (as there are fewer than 200 stars beyond). The smallest outermost radius is 3.2 kpc with a median of 9.0 kpc. Among the 12,228 HEALpixels, the number of shells varies from 11 to 45 with a median value of 36. The number of stars per HEALpixel varies from 3200 to 3.3 million, with a median of 18 thousand. The top left panel of Figure 2 shows an example of the shell structure for one HEALpixel.

Once the shells have been fixed, I compute robust estimates¹ of the mean and standard deviation of the two velocities, as well as their correlation coefficient, for the stars in each shell. I then fit a smoothing spline to each of the parameters independently as a function of distance (right panels of Figure 2). The number of degrees of freedom of these fits is an increasing function of the number of shells, varying from three to ten with a median of eight. Beyond the outermost radius, the Gaussian parameters are kept constant at the value at the outermost radius. As the fit for the standard deviations could go to very small or even negative values, I force the fit of the two standard deviations never to drop below 2 km s⁻¹.

The results of this process for a relatively sparse field at high latitude, where some shells are merged, are shown in Figure 2. In this particular example, the variation of the correlation with distance is overfit; in other HEALpixels, there is no clear correlation. For this reason, I do not use these fits and I set the correlation to be zero for all distances for all HEALpixels; this is a conservative choice, making the prior less informative than it otherwise would be. Figure 3 summarizes the prior for a HEALpixel at a lower latitude toward the Galactic bulge, with many more stars where no shells are merged. One should note the large change in the mean velocities for the most distant shells.

**Figure 3.** As Figure 2, but now for HEALpixel 7593 at HEALpixel 7593 at (l, b) = (290, 77).
Download figure:
Standard image High-resolution image

fdg — **Figure 3.** As Figure 2, but now for HEALpixel 7593 at HEALpixel 7593 at (l, b) = (290, 77).
Download figure:
Standard image High-resolution image

When building this velocity prior, I excluded stars from the mock catalog that are fainter than the expected Gaia magnitude limit for that HEALpixel, as was done for the distance prior in Bailer-Jones et al. (2021; see Section 1 of that reference).

I experimented with making the velocity prior a function of magnitude or color. For magnitude, this accounted for little additional variance. Introducing a dependence on color was much harder, because some HEALpixels show gaps in their color distribution, which would necessitate additional modeling assumptions to get a continuous prior.

Figure 4 shows how the magnitude of the mean of the velocity prior varies over the sky for several distance slices. The velocities of nearby stars are small, but they increase at larger distances as we probe the faster orbits in the central regions of the Galaxy (one should note the change in color scale across the slices). The low-velocity regions near (l, b) = (64°, 16°) and (l, b) = (244°, − 16°) in the 0 kpc distance slice are the approximate solar apex and antapex directions (respectively), i.e., the direction in which the Sun is moving through the Galaxy.

**Figure 4.** Variation of the mean of the velocity prior over the Galactic sky at six different distances. The quantity shown is $\sqrt{{v}_{\alpha * }^{2}+{v}_{\delta }^{2}}$ , where v_α* and v_δ are the means of the two components of the velocity prior at the specified distance in each HEALpixel. In each distance slice, the color bar covers the full range of velocities plotted, and it is different for each slice.
Download figure:
Standard image High-resolution image

The inference (Equation (3)) uses the prior on the two components of the velocity, not the total transverse velocity. Plots like Figure 4 for the individual components are less astrophysically informative, however, because they show features of the (equatorial) coordinate system used. An artifact of this is nonetheless visible at the south equatorial pole (l, b) = (302 fdg 9, − 27 fdg 1) in some of the distance slices in Figure 4. This is a projection effect of the velocity vector onto the rapidly changing v_α* and v_δ basis vectors at the poles. In principle, we should not see it in the total transverse velocity, but because the prior is computed over the finite area of the HEALpixel, the average shows large jumps. The same effect occurs at the north equatorial pole, (l, b) = (122 fdg 9, + 27 fdg 1), but is barely visible. Using a smaller HEALpixel would mitigate this effect but create others, and this is a limitation of using a discrete prior.

Figure 5 shows how the square root of the sum of squares of the standard deviations of the velocity prior varies. This too is not directly used in the inference, but it gives an idea of the width of the velocity prior, which influences how precisely distance can be constrained for a given proper motion (discussed in Appendix A). The velocity constraint tends to be stronger at lower latitudes, as there is less dispersion about the mean rotational velocity of disk stars at a given distance and longitude. At higher latitudes, the dispersion increases significantly because the more isotropic orbits of halo stars dominate the tails of the samples. While not immediately visible, due to the variable color scale, the velocity constraint is generally weaker (larger standard deviation) at larger distances.

**Figure 5.** As Figure 4, but now showing $\sqrt{{\sigma }_{\alpha * }^{2}+{\sigma }_{\delta }^{2}}$ , where σ_α* and σ_δ are the standard deviations of the two components of the velocity prior.
Download figure:
Standard image High-resolution image

2.3. Posterior Sampling

I sample the 3D posterior (Equation (3)) for each star with an MCMC method (Foreman-Mackey et al. 2013).² The distance is generally strongly correlated with the velocity components, because for a given proper motion, an increase in the distance can be compensated for by a decrease in the velocity. Such correlations did not pose a problem for the sampler, and estimating instead the true angular velocities—which would not be correlated—did not confer particular advantages. After some experimentation, I found that 60 burn-in iterations followed by 200 sampling iterations with 20 walkers gave good convergence. The acceptance rate was around 0.6, and the sample chains were thinned before parameters were estimated. For each parameter, I estimate the 16th, 50th, and 84th percentiles, the median providing the parameter estimate and the other two the lower and upper 1σ-like uncertainty estimates.

3. Performance on the Mock GDR3 Catalog

To estimate the accuracy of the distance and velocity estimates, I apply the method to a noisy version of the mock Gaia catalog. This is achieved by adding Gaussian random noise to the parallaxes and proper motions at the expected noise level in GDR3 (Rybizki et al. 2020). The noise-free version of this same catalog was used to build the priors, so this assessment shows what the method can achieve before being affected by any mismatch between the prior and reality.

I analyze the results over the whole sky using two samples. The constant fraction sample selects 0.8% of the stars at random from each HEALpixel to give a total of 10.3 million sources. Like the full Gaia sample, this is dominated by faint and low-accuracy parallax stars, many of which lie in the Galactic plane and bulge. While this is useful for examining trends in performance with distance and magnitude, it is not as useful for the examining the variation over the sky. For that purpose, I use the constant number sample, which comprises approximately 900 sources selected at random per HEALpixel, for a total of about 11.5 million sources. Both samples include only stars that are brighter than the prior magnitude limit in each HEALpixel (Section 2.2).

3.1. Overall Performance

Figure 6 compares the estimated median distances and velocities against their true values for the constant fraction sample. The top left panel shows the geometric distance estimated from Equation (2). The other panels show the kinegeometric distance and two velocity components (v_α*, v_δ) from Equation (3). The variations and differences are discussed in more detail below. Notable for now, though, is the smaller spread for kinegeometric distances compared to geometric ones at larger distances. From now on, I will discuss and show results for the total kinegeometric velocity $\sqrt{{v}_{\alpha * }^{2}+{v}_{\delta }^{2}}$ rather than its individual components. This will be compared with the total velocity derived from the geometric method, ${{kr}}_{\mathrm{geo}}\sqrt{{\mu }_{\alpha * }^{2}+{\mu }_{\delta }^{2}}$ .

Performance statistics for the constant fraction sample are summarized in the top block of Table 1. For the distances, I use the fractional residuals, (estimated-true)/true, because distance accuracy is a strong function of distance. For the total velocities, I use the actual residuals, (estimated-true). We see that the residuals for the two methods are broadly similar, with a slightly heavier tail to negative residuals. For the distances, the fractional bias—median of (estimated-true)/true—is −0.0128 for the geometric distance and −0.0156 for the kinegeometric distance. For the velocities, the bias—median of (estimated-true)—is −0.41 km s⁻¹ for the geometric method and −0.88 km s⁻¹ for the kinegeometric one. Thus, the biases are quite small on average.

Table 1. Median Performance for the Two Distance and Two Velocity Estimates according to Four Different Statistics, Computed on the GDR3 Mock Catalog Using the Two Different Samples (Section 3) in the Top and Middle Blocks

		Geo Distance	Kinegeo Distance	Geo Velocity	Kinegeo Velocity
constant fraction sample
Bias	est-true	−20.9 pc	−29.9 pc	−0.410 km s⁻¹	−0.875 km s⁻¹
Scatter	∣est-true ∣	977 pc	779 pc	19.8 km s⁻¹	16.0 km s⁻¹
Fractional bias	(est-true)/true	−0.0128	−0.0156	−0.0098	−0.0178
Fractional scatter	∣(est-true)/true ∣	0.228	0.185	0.238	0.192

constant number sample
Bias	est-true	−0.788 pc	−2.53 pc	−0.030 km s⁻¹	−0.192 km s⁻¹
Scatter	∣est-true ∣	199 pc	177 pc	6.14 km s⁻¹	5.68 km s⁻¹
Fractional bias	(est-true)/true	−0.0018	−0.0047	−0.0009	−0.0055
Fractional scatter	∣(est-true)/true ∣	0.125	0.112	0.128	0.114

constant fraction sample		for expected GDR5 astrometric precisions
Bias	est-true	−10.1 pc	−18.1 pc	−0.229 km s⁻¹	−0.461 km s⁻¹
Scatter	∣est-true ∣	722 pc	595 pc	14.7 km s⁻¹	11.6 km s⁻¹
Fractional bias	(est-true)/true	−0.0072	−0.0105	−0.0057	−0.0106
Fractional scatter	∣(est-true)/true ∣	0.168	0.139	0.178	0.139

Note. The bottom block shows the results for simulations at the higher precisions expected in GDR5, as discussed in Section 5.

Download table as: ASCII Typeset image

The biases, which we can think of as the systematic deviation about the identity line in Figure 6, are part of the estimation error. The total estimation error I characterize by the median of the absolute values of the residuals, more specifically as median(∣estimated-true)∣/true) for the distance and median(∣estimated-true∣) for the velocities. These robust versions of the error I refer to as the scatter. One should note that this includes the bias. For the distances, the fractional scatters are 0.228 for geometric distances and 0.185 for kinegeometric distances. For the velocities, the scatters are 19.8 km s⁻¹ for the geometric method and 16.0 km s⁻¹ for the kinegeometric one. These are quite large, but it must be remembered that most of the stars in the sample are distant and faint.

The middle block of Table 1 shows the performance statistics on the constant number sample, where we see much smaller bias and scatter in all cases. This is because this sample has more stars at higher latitudes that are generally nearer. This is not as representative, so for global statistics I report the worse results on the constant number sample elsewhere in this paper. The final block is for GDR5, discussed in Section 5.

The goal of the present paper is to introduce the kinegeometric estimates. The above statistics show that kinegeometric distances are overall more accurate than the geometric ones by a factor of only 1.25 (the ratio of the scatters). Within this, they have a bias 1.2 times larger than the geometric estimates, although the average bias is only −1.5%. In the constant fraction sample, 58% of the stars have a kinegeometric error (scatter) that is smaller than the geometric error. Twenty-seven percent of stars have a kinegeometric error that is smaller than the geometric error by a factor of two or more (compared to 17% the other way around). The is no obvious, narrow part of parameter space where one estimate is always better than the other, but there are broad variations with distance and position, which we now explore.

3.2. Performance as a Function of Distance

Figure 7 shows how the bias and scatter of the kinegeometric distance estimate varies as a function of distance as a density plot, with the orange line indicating the median at each distance. We see a complex dependence on distance, in part because it is an average over sources over the whole sky. Compared to the geometric method (blue line), the kinegeometric bias and scatter are slightly smaller at most distances. The reason the geometric distance nonetheless showed a slightly smaller bias overall (Table 1) is because a large number of stars lie in the range where the geometric bias is smaller. Beyond about 7 kpc, the kinegeometric distances are more accurate. This is the regime where the kinegeometric distances are most useful.

Figure 8 shows the performance as a function of parallax S/ N (defined as the parallax divided by the parallax uncertainty). The performances of the two methods are similar in terms of bias, but the kinegeometric distances have lower scatter for low parallax S/N. This is precisely the region where the proper motions are helping, because even though the proper motion S/ Ns are generally lower too, they are still high enough to provide additional information on distance. At higher parallax S/N, the differences are small because the parallax dominates, and the priors are less important.

Figure 9 compares the bias and scatter of the two total velocity estimates as a function of distance. Again, we see a large range of these statistics at a given distance, as evidenced by the dotted lines in these plots. The velocity estimates have similar performances (in both bias and scatter) for stars with true distances less than about 7 kpc. At larger distances, the kinegeometric velocity estimates are more accurate, due to the assistance from the proper motions and velocity prior.

3.3. Performance as a Function of Direction

We now turn to look at how the performance varies over the sky using the constant number sample.

Figure 10 shows the fractional bias and scatter, by computing the statistic for each star and then taking the median over each HEALpixel. Looking first at the upper panels (the bias), we see that both methods tend to underestimate distances in lines of sight in the disk and bulge. This is because these directions are dominated by distant stars (see Figure 1), which generally have low parallax S/Ns, which in turn introduces a bias (explained in Appendix B). The kinegeometric estimates show larger negative biases by up to an average of 13% in two regions above and below the Galactic center. This is related to the large changes in the velocity prior beyond a few kpc (see Figure 4). The geometric distance estimates are less affected.

The scatter plots (bottom row of Figure 10) show a larger scatter in both distance estimates at low latitudes. This is a consequence of most stars being at large true distances, so they have low-S/N parallaxes. The kinegeometric distances show a smaller scatter, though, because the proper motions improve the distance accuracy primarily at low parallax S/N (as already discussed in the context of Figure 8).

Figure 11 shows the differences between the kinegeometric and geometric distance estimates over the sky (computing the median per HEALpixel of the differences for each star). At higher Galactic latitudes, the estimates are very similar on average. In the central regions of the Galaxy, the kinegeometric distances are smaller, by up to 11% on average within a HEALpixel. In the rest of the Galactic disk, the kinegeometric distances are larger by up to 8% on average.

Figure 12 shows how the performances of the two velocity estimates vary over the sky. The bias patterns are broadly similar to what we saw with the distances, especially for the kinegeometric estimates, e.g., the blue blobs at l ≃ 0° mentioned above. The bias map for the geometric velocities, however, shows two regions in the disk with a small but significant overestimate, at around l = − 155° and l = 150° (the "red cheeks" in the top left panel of Figure 12). They are slightly displaced from the Galactic plane, so they are presumably related to the warp of the disk in the Galaxy model on which the mock catalog is based. They do not appear in the bias map for kinegeometric distances (top right panel of Figure 12). Some structures in these regions are seen in the nearer slices of the velocity prior (Figure 4), but not in the distance prior (Figure 1), suggesting that not accounting for these in the geometric prior is problematic. Naturally, if these features do not exist in the real Galaxy, the (unknown) biases in the inference in GDR3 could be different.

The lower panels of Figure 12 show the scatter in the two velocity estimates, which are similar in pattern, but slightly larger for geometric velocities. Figure 13 shows the average differences per HEALpixel between the two velocity estimates.

**Figure 13.** Median differences between the kinegeometric and geometric velocities per HEALpixel (computed on the individual stars), for the constant number sample in the mock catalog. The color bar has white at zero and diverges on different linear scales to the maximum negative and positive differences.
Download figure:
Standard image High-resolution image

3.4. Uncertainty Estimates

From the posterior PDF for each star, I compute the 16th and 84th percentiles for each parameter, then take their difference from the median to obtain lower and upper 1σ-like uncertainty estimates. Half their difference gives a single uncertainty estimate for the parameter. If the residuals (median minus true) have an unbiased Gaussian distribution, and if the uncertainty is a statistically correct estimate of the residual, then the ratio of the residual to the uncertainty—the normalized residuals—should have a unit Gaussian distribution. These are shown in Figure 14. All four distributions have a mean close to zero and a standard deviation close to unity, although the distributions for both distance measures show a slight negative skew. This is not surprising, because we know there is a negative distance bias for distant, low parallax S/N stars (see Appendix B). But overall, this shows that the uncertainty estimates are good statistical estimates of the true error.

**Figure 14.** Normalized residuals for the geometric posterior distance (top left) and for the three parameters of the kinegeometric posterior (other panels), computed for all sources in mock HEALpixel 6200. The smooth red lines show the unit Gaussian for comparison.
Download figure:
Standard image High-resolution image

4. Results on GDR3

There are 1.46 billion stars in GDR3 that have parallaxes and proper motions. As with the mock catalog, I apply the geometric and kinegeometric inference methods to two randomly selected subsets: a constant number sample with 0.8% of all sources (10.3 million sources), and a constant fraction sample with 900 sources per HEALpixel (10.8 million sources), in both cases restricted to sources brighter than the prior magnitude limit in each HEALpixel (Bailer-Jones et al. 2021, Section 1), to enable comparison to the estimates on the mock catalog in the previous section.

4.1. Distances

Figure 15 shows the median distance per HEALpixel. The median distance increases to lower Galactic latitudes in both methods, because of the higher density of stars at larger distances within the disk. The median distance is smaller in the Galactic plane, however, because Gaia's depth is limited by strong interstellar dust extinction. Many of the patterns we see close to the plane are due to dust density variations. At high Galactic latitudes, the depth through the disk is small, and although distant stars in the halo are visible, they are rare, so the median remains relatively small. The bulge is quite prominent in both distance estimates as a larger, approximately circular region of more distant stars. Close to the Galactic center, we see dust-free lines of sight that extend the median distance to over 7 kpc. Neither the mock catalog as we use it nor the prior contain star clusters, so cluster distances will not be correct. This is most noticeable for stars in the Large and Small Magellanic Clouds, which we see have severely underestimated distances.

The differences between the two distance estimates are shown in the bottom panel of Figure 15 (median of differences per star). The largest differences are in parts of the Galactic plane, parts of the bulge, and in several star clusters. Most prominent is that the kinegeometric distance is smaller in the bulge above and below the Galactic center. This is a result of the measured proper motions combined with the velocities of the Galaxy prior. In parts of the Galactic plane, the kinegeometric distances are larger; in others, the differences are negligible; and in the northern part of the disk near the anticenter, the kinegeometric distances are smaller. These patterns are presumably related to the Galactic bar and warp, which unsurprisingly differ between GDR3 and the Galactic model used to make the prior.

The larger kinegeometric distances in the upper and lower parts of the bulge were also seen in the results on the noisy mock data in Figure 11, although they were less prominent there and also extended through the Galactic disk. We should not overinterpret differences between Figures 11 and 15, however: they can occur for many reasons, including differences in the underlying stellar spatial and velocity distributions or in the extinctions between the mock catalog and real Galaxy. The larger color bar scale on the difference plot for the GDR3 results shows that kinegeometric and geometric distances differ much more in the real inference than they do in the mock inference, no doubt a result of these differences.

Figure 16 compares the two distance estimates as a function of geometric distance. The density plot and the quantiles (dashed lines) show that there is a large spread, but the median trend (solid orange line) is for the kinegeometric distances to be on average smaller than the geometric ones beyond 6 kpc, by 20% at 10 kpc.³ This is not true in all lines of sight, however, such as the Galactic center and parts of the disk, as we saw in Figure 15. When we consider just the low parallax S/N stars (<3), we see a large difference in the distance estimates on average also for nearby stars (the dashed pink line in Figure 16).

We saw in Section 3.4 (Figure 14) that the width of the posterior (the precision) was a good statistical estimate of the actual error. We can therefore compare the posterior width (the difference between the 84th and 16th percentiles) for the two types of distance estimate in order to predict which is more accurate. Using the constant fraction sample, we find that, for 62% of the stars, the kinegeometric posterior is narrower than the geometric one. Fifteen percent of stars have a kinegeometric posterior that is less than half the width of the geometric one, compared to just 4% the other way around. This suggests kinegeometric distances will be more accurate for a significant number of stars. Figure 17 plots the ratio of these posterior widths over the sky, taking the median over all stars in each HEALpixel. We see that in almost all parts of the sky—but in particular in the bulge and disk—the kinegeometric distances are more precise (but we must remember that each HEALpixel is an average over a large range of distances). We also find that kinegeometric distances tend to be more precise for lower parallax S/Ns, as we would expect.

**Figure 17.** Median per HEALpixel of the ratio (kinegeo/geo) of the posterior widths for the constant fraction sample in GDR3. Numbers less than 1 (in orange) show where the width of the kinegeometric distance posterior is narrower than the width of the geometric distance posterior, and these are therefore statistically likely to be more accurate. We see that this is the case in most parts of the sky, but particularly so in the disk and bulge. The bilinear color bar with white at unity and diverges on different linear scales to the minimum and maximum values.
Download figure:
Standard image High-resolution image

4.2. Velocities

Figure 18 shows the median velocity per HEALpixel inferred by both the geometric method (geometric distance times proper motion) and the kinegeometric method. These are transverse velocities relative to the Solar System barycentre. The largest velocities occur within about 90° in longitude of the Galactic center, in particular at low latitudes. This is where we see mostly disk stars on smaller radius orbits than the Sun's. The exception is in the Galactic plane, where dust extinction means Gaia only sees nearby stars that are orbiting the Galaxy with low velocities relative to the Sun. At larger longitudes, we see the outer part of the disk—where, due to the flattening of the Galactic rotation curve, stars are orbiting at speeds more similar to the Sun's.

The differences between the two velocity estimates are shown in the bottom panel of Figure 18. This is broadly similar to the difference plot obtained from the mock catalog results in Figure 13, although the scales are different. The main difference (other than the absence of star clusters in the mock results) is around the Galactic center: whereas for mock results the geometric velocities were several km s⁻¹ larger on average, in GDR3 the kinegeometric velocities are slightly larger.

5. Expected Performance with the Final Gaia Data Release (GDR5) and the Limits of Proper Motions

The distance accuracy that this method can attain is determined by the accuracies of the parallaxes and proper motions, as well as the width of the distance and velocity priors. The accuracy of the data will improve in later Gaia data releases. If t is the time span of observations, then whereas the parallax accuracy improves as t^1/2 on account of the increase in the number of observations, the proper motion accuracy improves as t^3/2, due to the additional increase in the observational baseline. This assumes systematic errors do not dominate, which is the case for most stars. GDR3 is based on 34 months of data, whereas the final release, GDR5, will be based on about 10.5 yr (126 months) of data. Using these factors, the parallax and proper motion accuracies should therefore increase by respective factors of 1.9 and 7.1 between GDR3 and GDR5. This suggests that, while both the kinegeometric and geometric distance estimates should be more accurate with GDR5 data, the improvement for the kinegeometric ones might be larger.

I tested this expectation using the mock catalog with similar simulations as before (Section 3), but now scaling down the parallax and proper motion uncertainties by the above factors (1.9 and 7.1, respectively). The results for the constant fraction sample are shown in the bottom block of Table 1. Looking first at the geometric distances—which depend on the parallaxes but not the proper motions—we see that the fractional bias is decreased from GDR3 to GDR5 by a factor of 1.8, and the fractional scatter is decreased by a factor of 1.35. These improvements are not as large as the improvement in the parallaxes. This is expected, because the prior also plays a role: The distance accuracies for the numerous stars with very low parallax S/Ns (≲1) are only weakly dependent on their parallax S/Ns, because their distances are prior-dominated. For such stars, doubling the parallax S/N hardly changes their distance estimates—and therefore hardly improves their accuracies.⁴

For the kinegeometric distances, the fractional bias and scatter are decreased by factors of 1.55 and 1.35, respectively. Contrary to possible expectations, the kinegeometric distances do not improve more than the geometric distances: On average, at least, the larger improvement in the proper motions over the parallaxes has not helped. The reason is that the proper motion only provides a distance estimate when combined with an assumed velocity distribution (the prior) via Equation (1). Once the proper motions are sufficiently precise, such that the width of this velocity distribution is the limiting factor in determining the distance precision, then improving the proper motions further will not help. This is explained in Appendix A. It appears that, on average at least, this limit has already been achieved in GDR3.

The variation of the distance accuracy (bias and scatter) with distance for GDR5 is shown as the dashed lines in Figure 19. The profiles are broadly similar to GDR3 (solid lines), but we see that the improvements in the accuracy vary with distance. Below 2 kpc, the scatter in both the geometric and the kinegeometric distances is reduced by a factor of nearly two, whereas at around 5 kpc, there is little improvement.

For the velocities, the situation is similar (Table 1). The scatter in the geometric velocities is expected to decrease by a factor of 1.35 times between GDR3 and GDR5, and by a factor of 1.4 for kinegeometric velocities. We might have expected much more improvement in both, because the proper motions improve by a factor of 7.1. But the mapping from proper motion to velocity depends also on the estimated distance, and the improvement in this, which is coming mostly from the parallaxes, is also a limiting factor in the improvement in the velocities.

I conclude that, in GDR5, kinegeometric distances and velocities should be about 1.35 times more accurate than they are in GDR3, when averaged over the entire Gaia sample in GDR3. These will still be only about 1.25 times more accurate than geometric distances and velocities. The limited improvement is because the velocity priors are relatively broad.

While it may be worth using proper motions in GDR5 to improve the distance estimates in some regions, this comes at the cost of building in assumptions on velocities. If we are happy with this, then a better accuracy should be obtained by combining parallaxes and proper motions also with colors, magnitudes, and even spectra.

6. Summary and Conclusions

The kinegeometric method introduced in this paper estimates distances and transverse velocities together in a 3D posterior PDF using the parallaxes and proper motions. It uses a direction-dependence distance prior and a direction- and distance-dependent velocity prior. Geometric distances, in contrast, are obtained from just the parallaxes and the distance prior.

The performance was assessed on a random subset of the mock GDR3 catalog. This showed that on average, kinegeometric distances are 1.25 times more accurate than geometric ones (Table 1). The median fractional distance residual (scatter) of the kinegeometric distances is 19%, averaged over all sources; averaged over the sky, it is 11%. The main improvement brought by the proper motions is where the parallax S/N is low (≲4; Figure 8). There is a considerable spread in performance with position in the Galaxy, however. Close to the Sun the kinegeometric distance errors are very small, then increase steadily to 20% at 2 kpc and beyond (Figure 7). Kinegeometric distances have a low bias at high Galactic latitudes, but up to several percent in the bulge and disk, with complex patterns emerging at low latitudes due to the stellar density and prior variations (Figure 10). The biggest advantage of kinegeometric distances over geometric ones is for stars beyond about 7 kpc from the Sun.

The kinegeometric method also provides transverse velocity estimates. These can be compared to "geometric" velocities obtained by multiplying the geometric distances by proper motion. Both approaches have similar overall performance for stars within 7 kpc, but beyond that the kinegeometric velocities are considerably more accurate, including a lower bias (Figure 9).

When applied to real GDR3 data, we find that differences in the two distance estimates depend on direction and distance. Beyond several kpc, kinegeometric distances tend to be smaller on average (Figure 16), which in particular is reflected by two regions in the bulge (Figure 15). In the disk, the picture is more complicated, a reflection of how the disk kinematics in the prior combine with the measured proper motions to inform the distance estimates. The precisions (widths) of the posteriors are a statistically good estimate of the accuracies (Figure 14), so we can use them to investigate the expected accuracy of the distance estimators in the absence of a ground truth. The kinegeometric distances are more precise than the geometric ones for 62% of stars, and (on average) over most of the sky (Figure 17).

Overall, we see only a modest average improvement in distance accuracy when including proper motions. This conclusion is based on the 34 months of Gaia data in GDR3. The S/N of the parallaxes and proper motions are expected to improve in future Gaia data releases, by respective factors of 1.9 and 7.1 in the final release (GDR5), assuming this is based on 10.5 yr of data. I find that this will improve the accuracy of the kinegeometric distances and velocities by a factor of 1.35. The geometric distances improve by a similar factor, meaning kinegeometric distances are still only about 1.25 times more accurate than geometric ones in GDR5. The reason for this limited improvement is that proper motions only constrain distances when combined with an expected velocity, which we introduce as a distribution of velocities. If this prior is broad, then even very precise proper motions do not constrain the distances much (see Appendix A). This is very different from the situation with parallaxes, in which an increased accuracy of the parallax translates directly into an increased accuracy in the distance (provided the distance is not prior-dominated). Hence, for many stars in Gaia, the velocity distributions in the Galaxy are too broad for even high-precision proper motions to provide strong constraints on distance.

Acknowledgments

I thank Morgan Fouesneau, Rene Andrae, and Sara Jamal for discussions on this project, and the referee for several useful suggestions. This work was funded in part by the DLR (German space agency) via grant 50 QG 2102 and has made use of data from the European Space Agency (ESA) mission Gaia (http://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC; http://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.

Facility: Gaia - .

Appendix A: How the Proper Motion and Velocity Prior Constrain Distance

The kinegeometric posterior of Equation (3) is a 3D distribution coupling distance and velocity. We can understand how the proper motion and velocity prior combine to produce distance estimates by considering a simpler version with just one velocity component and no parallax. In this case, the posterior over distance is

$\begin{eqnarray}&&\begin{array}{l}P(r\,| \,\mu )=\displaystyle \int P(r,v\,| \,\mu )\,{dv}\\ \quad =\ \displaystyle \int \displaystyle \frac{P(\mu \,| \,r,v)P(r,v)}{P(\mu )}\,{dv}\propto P(r)\displaystyle \int \\ P(\mu \,| \,r,v)P(v\,| \,r)\,{dv}.\end{array}\end{eqnarray} \tag{ A1 }$

We see that the proper motion likelihood constrains the distance only through the integral of its product with the velocity prior. This is quite different from the parallax, where P(r∣ϖ) ∝ P(r)P(ϖ∣r) without any integral.⁵ Consequently, even if we measure a very precise proper motion, such that the likelihood P(μ∣r, v) for some v in Equation (A1) is narrow, if the corresponding velocity prior P(v∣r) is broader, the contribution to the posterior will also be broad.

This combination is demonstrated in Figure 20. I use the v_δ prior from HEALpixel 7593 (Figure 3), except that I take the absolute value of the mean of the velocity prior, to more conveniently deal with positive values. The prior (blue curves) is shown for 10 and 30 km s⁻¹. Both velocities show a preference for smaller distances. The two orange curves are the likelihoods corresponding to these same two velocities assuming a proper motion measurement of 1 ± 0.25 mas yr⁻¹. The values of 1 mas yr⁻¹ and 10 km s⁻¹ correspond to a distance of 2.1 kpc (and 30 km s⁻¹ to 6.3 kpc; Equation (1)), but we see that this low-S/N measurement does not constrain the distance much for either velocity. Each gray curve is the product of the prior and likelihood, and therefore a term under the integral in Equation (A1) for a given velocity.

If we repeat this for a continuous range of velocities, then the sum of all the resulting gray curves is the integral in Equation (A1), which is the unnormalized distance posterior (but omitting a distance prior). This is illustrated in the top left panel of Figure 21, where only a few of the gray curves with different velocities are shown, but the resulting integral (dashed black line) has been computed using a broad, dense grid of velocities. We see from this plot that, even though each contribution to the integral is relatively narrow, their sum—and therefore the posterior—is broad. Consequently, even if we have higher-S/N proper motions—and so narrower likelihoods and thus narrower gray distributions—the integral is more or less unchanged. This is demonstrated in the bottom left panel of Figure 21, where the S/N is increased by a factor of four. This is the main point: improving the proper motion S/N does not necessarily improve the distance accuracy. This is why we expect little improvement in kinegeometric distances from GDR3 to GDR5 (Section 5). If the prior is narrower, then we can get a narrower posterior for given data. This is illustrated in the right columns of Figure 21, where the standard deviation of the velocity prior has been halved at all distances. The proper motion likelihoods are unchanged (compared to the panel on the left), but some are now more suppressed by the prior.

**Figure 21.** Demonstration of how the likelihood—prior product, P(μ∣r, v)P(v∣r), shown as the gray curves, sum to make the unnormalized distance posterior (black dashed line) via Equation (A1). The gray curves are shown just for 10 velocities from 5 to 50 km s⁻¹ in steps of 5 km s⁻¹, but some are so flat as to be invisible. In all panels, the measured proper motion is 1 mas yr⁻¹. The top row is for a proper motion S/N of 4, and the bottom for 16. The right column is for a narrower prior with half the standard deviation of that in the left column. None of the functions shown are probability density functions in r. The gray ones have been normalized to unit area over the range plotted. The black dashed curve has been scaled to have the same maximum as the largest gray peak.
Download figure:
Standard image High-resolution image

Appendix B: Random Errors Cause Bias

One might think that, because the priors are built from the mock catalog, there should be no bias in the results on the mock catalog in Section 3 (for either distance estimate). This is incorrect for at least three reasons. First, the mock catalog used in the performance assessment is noisy. Although the noise is symmetric (Gaussian) in the parallax, it is not symmetric in the distance, due to their nonlinear relation (r ∼ 1/ϖ; see Bailer-Jones 2015). Second, the distance and velocity priors are fits, so they do not represent all the variances in the mock data. We may in fact not want a prior tuned exactly to all the details and specific choices of the mock catalog, because this would be a rather strong prior. Third, in reflecting the true Galaxy, the distance prior has the bulk of its probability mass within a few kpc. Stars beyond these distances tend to have low parallax S/Ns, i.e., the data are relatively uninformative, so the prior will tend to pull the inferred distance below the true distance for these stars. This is why we see larger biases at larger distances in Figure 7. Nearby low parallax S/N stars will tend to have their distances pulled up by the prior for the same reason, but this is much less common because there are fewer nearby stars and nearby stars tend to have higher parallax S/Ns on account of their larger parallaxes and brighter apparent magnitudes.

"Not using a prior" is not a solution to this problem, because a prior is invariably present implicitly. Indeed, if we just invert a parallax to estimate the distance, the biases are much larger (Bailer-Jones 2015). A prior uniform in distance or velocity is a poor choice, and much more biased than the present choice: A uniform distance prior is equivalent to assuming that the stellar number density drops as 1/r² from the Sun, which is demonstrably wrong. A uniform velocity distribution relative to the Sun contradicts the existence of the different components of the Galaxy. Unless one is prepared to use highly tuned priors, low-S/N data give low-S/N inferences. The only other solution is to acquire better data.

Estimating Distances from Parallaxes. VI. A Method for Inferring Distances and Transverse Velocities from Parallaxes and Proper Motions Demonstrated on Gaia Data Release 3

Article metrics

Author affiliations

Dates

Abstract

1. Introduction

2. Method