CROSS IDENTIFICATION OF STARS WITH UNKNOWN PROPER MOTIONS

The positional accuracy is characterized by a probability density function (PDF) on the celestial sphere. In a given model M, this $p(\mbox{\boldmath {$x$}}|\mbox{\boldmath {$r$}},M)$ function tells us where to expect $\mbox{\boldmath {$x$}}$ detections of an object that is at its true location $\mbox{\boldmath {$r$}}$. Throughout this paper, we use three-dimensional unit vectors for the positions on the sky, e.g., the aforementioned $\mbox{\boldmath {$x$}}$ and $\mbox{\boldmath {$r$}}$ quantities. Usually, the PDF is a very sharp peak and is assumed to be a normal distribution with some angular accuracy σ. The correct generalization to directional measurements is the Fisher (1953) distribution,

$$\begin{equation} F(\mbox{\boldmath {$x$}}|\mbox{\boldmath {$r$}},w) = \frac{w\,\delta (|\mbox{\boldmath {$x$}}|\!-\!1)}{4\pi \sinh w}\ \exp \left({w\,\mbox{\boldmath {$r$}}{}\mbox{\boldmath {$x$}}} \right) \end{equation} \tag{ 1 }$$

whose shape parameter w is essentially 1/σ² in the limit of large concentration; see details in Budavári & Szalay (2008).

The added complication comes from the fact that some objects are not stationary. If a given star is at location $\mbox{\boldmath {$r$}}$ now and has $\mbox{\boldmath {$\mu $}}$ proper motion, then Δt time later it would be at some other position $\mbox{\boldmath {$r$}}^{\prime }$,

$$\begin{equation} \mbox{\boldmath {$r$}}^{\prime } = \mbox{\boldmath {$r$}}^{\prime }(\Delta {}t; \mbox{\boldmath {$r$}}, \mbox{\boldmath {$\mu $}}), \end{equation} \tag{ 2 }$$

that is offset by a small displacement along a great circle. By substituting this position into our astrometric model, we create a new one M' with the added proper motion, $\mbox{\boldmath {$\mu $}}$, and time difference, Δt, parameters:

$$\begin{equation} p(\mbox{\boldmath {$x$}}|\Delta {}t,\mbox{\boldmath {$r$}},\mbox{\boldmath {$\mu $}},M^{\prime }) = F(\mbox{\boldmath {$x$}}|\mbox{\boldmath {$r$}}^{\prime }(\Delta {}t; \mbox{\boldmath {$r$}},\mbox{\boldmath {$\mu $}}),w). \end{equation} \tag{ 3 }$$

Naturally, there is nothing specific in this about the chosen characterization of the astrometry; one can use any appropriate PDF in place of the Fisher distribution instead.

At the heart of the probabilistic cross identification is the Bayes factor used for hypothesis testing. The question we are asking is whether our data D, a set of detected sources in separate catalogs with positions $\lbrace \mbox{\boldmath {$x$}}_i\rbrace$, are truly from the same object. For every catalog, we know its epoch and its astrometry characterized by a known p_i PDF. Let H denote the hypothesis that assumes that all measured positions are observations of the same object, and let K denote its complement, i.e., any one or more of the detections might belong to a separate object. By definition, the Bayes factor is the ratio of the likelihoods of the two hypotheses we wish to compare,

$$\begin{eqnarray} B(H,K|D) & = & \frac{p(D|H)}{p(D|K)} \end{eqnarray} \tag{ 4 }$$

that are calculated as the integrals over their entire parameter spaces.

If we assume that there is a single object behind the observations, we can integrate over its unknown proper motion and position to calculate

$$\begin{equation} p(D|H) = \int \!\!\!d\mbox{\boldmath {$r$}}\!\!\int \!\!\!d\mbox{\boldmath {$\mu $}}\,p(\mbox{\boldmath {$r$}},\mbox{\boldmath {$\mu $}}|H) \prod _{i=1}^{n} p_i(\mbox{\boldmath {$x$}}_i | \Delta {}t_i,\mbox{\boldmath {$r$}},\mbox{\boldmath {$\mu $}}, H), \end{equation} \tag{ 5 }$$

where the joint likelihood of H given the data is written as the product of the independent components and $p(\mbox{\boldmath {$r$}},\mbox{\boldmath {$\mu $}})$ is the prior on the parameters, which is the subject of the following section. The actual calculation of this likelihood depends on the prior and might only be accessible via numerical methods.

The complementary hypothesis is more complicated in the sense that the model has a set of independent objects with $\lbrace \mbox{\boldmath {$r$}}_i,\mbox{\boldmath {$\mu $}}_i\rbrace$ parameters, however, the result of the calculation turns out to be much simpler. Here, the integral separates into the product of

$$\begin{equation} p(D|K) = \prod _{i=1}^{n} \int \!\!\!d\mbox{\boldmath {$r$}}_i\!\!\int \!\!\!d\mbox{\boldmath {$\mu $}}_i\,p(\mbox{\boldmath {$r$}}_i,\mbox{\boldmath {$\mu $}}_i|K)\,p_i(\mbox{\boldmath {$x$}}_i | \Delta {}t_i,\mbox{\boldmath {$r$}}_i,\mbox{\boldmath {$\mu $}}_i, K). \end{equation} \tag{ 6 }$$

For each integral, we can select a reference time such that Δt_i = 0, hence the effect of the proper motion drops out, and we arrive at the same result as the stationary case discussed by Budavári & Szalay (2008).

To derive a more realistic $p(\mbox{\boldmath {$\mu $}}|\mbox{\boldmath {$r$}},H)$ prior, we choose to study the ensemble statistics of stars instead of approaching the problem with an analytic model. While the latter would have the advantage of providing a function at an arbitrary resolution, the formulae are difficult to derive and the analytic approximations might miss subtle details of the relation that could be relevant.

We study the properties of stars in the Sloan Digital Sky Survey catalog archive (York et al. 2000) that also contains accurate proper-motion measurements from the recalibrated United States Naval Observatory B1.0 Catalog (USNO-B; Munn et al. 2004). For this analysis, we pick stars from the Stripe 82 data set where multiple observations are available over a 300 deg² strip covering a narrow range in declination between ±125. These repeated observations were taken between June and December each year from 1998 to 2005 (Adelman-McCarthy et al. 2008). After rejecting saturated and faint sources (u should be in the range of 15–23.5 and g, r, i, and z in the range of 14.5–24), the number of stars is around 100,000. This size does not allow for a high-resolution determination of the prior, hence we analyze additional simulation data.

To extend the number of stars used in constructing the prior, we used the current state-of-the-art Besançon models (Robin et al. 2003) that match the SDSS distributions well. Assuming four different stellar populations in the Milky Way, using the Poisson equation and collisionless Boltzmann equation with a set of observational parameters (i.e., fitting parameters to the dynamical rotation curve), they compute the number of stars of a given age, type, effective temperature, and absolute magnitude, at any place in the Galaxy. The model has been successfully used for predictions of kinematics and comparison with observational data in studies, e.g., Bienaymé et al. (1992), Chareton et al. (1993), Ojha et al. (1999), Rapaport et al. (2001), and Soubiran et al. (2003). A total of 740, 000 stars were generated from the Besançon models using large-field equatorial coordinates, thus the prior is dominated by model data and not by SDSS measurements. The proper-motion distribution of SDSS data and that of the model are very consistent, with the model yielding somewhat wider distributions than the observations.

In preparation for binning the data, we separate the dependence of the prior on the different proper-motion components, $\mbox{\boldmath {$\mu $}}=(\mu _{\alpha },\mu _{\delta })$, and, omitting the explicit hypothesis, we write

$$\begin{equation} p(\mu _\alpha,\mu _\delta | \mbox{\boldmath {$r$}}) = p(\mu _\delta |\mbox{\boldmath {$r$}})\,p(\mu _\alpha |\mu _\delta,\mbox{\boldmath {$r$}}). \end{equation} \tag{ 10 }$$

Since the Stripe 82 data set in SDSS contains sources only in a narrow declination range between −126 and +126, we can safely neglect the dependence on declination in Equation (10) for the purpose of this study, thus

$$\begin{equation} p(\mu _\alpha,\mu _\delta |\mbox{\boldmath {$r$}}) \approx p(\mu _\delta |\alpha)\,p(\mu _\alpha |\mu _\delta,\alpha). \end{equation} \tag{ 11 }$$

We establish these relations one-by-one starting with the former using the basic property of conditional densities

$$\begin{equation} p(\mu _\delta |\alpha) = \frac{p(\mu _\delta,\alpha)}{\int \!p(\mu _\delta ^{\prime },\alpha)\,d\mu _\delta ^{\prime }}. \end{equation} \tag{ 12 }$$

To achieve a better signal-to-noise behavior across the entire parameter space, we do not use a uniform grid but vary bin sizes so that they follow the asinh() in both parameters. In this way, one can have higher resolution bins where more data are available (around the peak close to 0) and wider bins in the tail. We further improve the quality of the empirical prior by removing high-frequency noise with a convolution filter, whose characteristic width is approximately one pixel in size at any location. The integral in the denominator is evaluated by counting the stars in the appropriate bins using the widths of the bins to weight the counts.

Figure 1 shows the prior using the aforementioned nonlinear scale for μ_δ as a function of the position α. The distribution is centered on approximately −3 mas yr⁻¹, and the location of the mode is practically independent of the R.A. As one nears the direction of the Galactic Plane (α = 60° and α = 305° are the nearest regions in this stripe), the distribution becomes sharper. This is to be expected as, if we included a broader range in right ascension, the PDF would get even narrower as the velocity dispersion of stars decreases.

Zoom In Zoom Out Reset image size

The second term of Equation (11) is constructed similarly using the same adaptive binning and smoothing but in even higher dimensions. It is difficult to visualize a four-dimensional PDF, hence, in Figure 2, we plot slices of the prior at various α values. Both axes are shown in the transformed scale. The values α = 575 and α = 3075 represent the two edges of Stripe 82, which are closest parts to the Galactic Plane. The same effect can be seen on these panels as in Figure 1: looking out of the Plane the PDF gets more disperse. The boxy (squared) shape of the contour lines arises from the asinh() transformation; on a linear system, the contours would appear to be more circular.

Zoom In Zoom Out Reset image size

As mentioned earlier, Stripe 82 was observed repeatedly from 1998 to 2005, between June and December of each year. Thus, we can obtain multi-epoch observations to test our method. We choose a range of stars with different proper motions observed at different epochs. To be sure that for our tests the observed stars are the same in each epoch, we select bright stars (magnitude r < 17) with tolerances in all magnitudes (typically 0.7 in u, 0.4 in g, r, i, and 0.5 in z). The query gives us on the average 20 epochs per star, from which many were observed with small time intervals while the biggest time interval between the epochs is approximately 6.5 years. We divide the time interval into three approximately equal parts and thus get four observations of each star with much the same time intervals between them. According to Equations (5) and (6), only R.A. and decl. coordinates are used for calculating the Bayes factors; USNO measurements of proper motion on the forthcoming figures and tables are only shown as a reference. We randomly select a dozen stars for the following tests.

We calculate the integrals of the Bayes factor numerically. Our Monte Carlo implementation generates independent random positions $\lbrace \mbox{\boldmath {$r$}}_n\rbrace$ (three-dimensional unit vectors) and two random components of the velocity that yield the $\mbox{\boldmath {$\mu $}}_{n}$ vector in the tangent plane of each $\mbox{\boldmath {$r$}}_n$, i.e.,

$$\begin{equation} \mbox{\boldmath {$r$}}_n^{\prime } = (\mbox{\boldmath {$r$}}_n + \mbox{\boldmath {$\mu $}}_n\,\Delta {}t) / | \mbox{\boldmath {$r$}}_n + \mbox{\boldmath {$\mu $}}_n\,\Delta {}t |. \end{equation} \tag{ 13 }$$

In theory, one has to integrate the position over the whole celestial sphere, and the proper motion out to infinity, but the integrands always drop sharply in practice, hence one can bound all the relevant parameters easily to reduce the computational need and to use the above approximation to the motion. For more efficient implementations, one can utilize more sophisticated Markov chain Monte Carlo (MCMC) methods.

The uncertainty estimates include two separate sources of errors. The numerical imprecision is tuned by the number of generated random parameters and can be estimated in the process of the integration. In our calculations, this error term is kept at a low level and contributes 10⁻² order of magnitude to the value of the weight of evidence.

Another source of error comes from the uncertainty in position measurements. While this is small for the SDSS detections, σ = 01, the short time differences could yield large relative errors in the proper motion. To get the order of this error, we generate 100 random realizations of the position of every star from the appropriate Gaussian distribution and recalculate the Bayes factor to derive the rms error in the weight of evidence. The figures later in this section contain error bars that represent these 1σ deviations.

After the static case, the simplest is a uniform prior as introduced in Equation (9). This analytic formula may appear at first not to favor any particular proper motion, yet the Bayes factor has some non-trivial scaling properties that are worth considering.

As the displacement of the source is a product of the time difference and the proper motion, associations at the same distances but with varying time intervals will indeed have different Bayes factors. In the case of a longer Δt, only smaller proper motions will contribute to the integral in Equation (5) shrinking the integral domain. This yields a scaling by a factor of Δt⁻², as seen in the left panel of Figure 3. This means that associations will be assigned lower qualities if they are farther in time even if their angular separations are identical.

Zoom In Zoom Out Reset image size

Another interesting aspect is the selection of the limiting μ_max value. Our choice of 600 mas yr⁻¹ is admittedly somewhat arbitrary and was selected to cover the stars in our sample. If one decreased its value, then stars moving at faster speeds would quickly get lower Bayes factors and clearly not be associated. Increasing limits make the value of the constant prior drop, which in turn will lower the quality of the associations. For illustrating this effect, we compute the Bayes factors for a star with μ = 360 mas yr⁻¹ and Δt = 2 years as a function of μ_max, see the right panel of Figure 3. As the prior is proportional to μ⁻²_max, the curve follows the same trend.

First, we analyze the quality of the associations as a function of time difference between the observations. The top panels of Figure 4 show the logarithm of the Bayes factor, aka the weight of evidence for all stars in our test sample as a function of their proper motion. Open circles represent the results from the static model that can be obtained analytically as in Budavári & Szalay (2008), and crosses show the new measurements from the numerical integration of the improved model using the empirical prior introduced in this study. Triangles signal the value for a simple model of constant proper-motion prior with μ_max = 600 mas yr^-1. If we correct for small relative differences in time intervals between the epochs taking 2, 4.5, and 6.5 years as a reference, respectively, this prior yields practically constant weights of evidence. For reference, the W = 0 threshold is plotted as the dashed horizontal line. This is the theoretical dividing line above which the observations support the hypothesis of the match. All panels contain the same objects but the calculations are based on different detections that are farther apart in time as we go from left to right. What we see immediately is that as the time difference increases, the models provide increasingly different results: the static model starts rejecting stars with larger proper motions much faster than models that accommodate the possibility of the sources moving.

Zoom In Zoom Out Reset image size

While the only objective measure of the quality for the match is the Bayes factor, its interpretation for the uninitiated is admittedly not as obvious at first as a probability value would be, where one has a good sense of the meaning of the values. From the Bayes factor, we can calculate the posterior P(H|D), if we have a prior P(H) via the equation

$$\begin{equation} P(H|D)=\left[ 1 + \frac{1-P(H)}{B(H,K|D)\,P(H)} \right]^{-1}. \end{equation} \tag{ 14 }$$

Assuming a constant prior over the sky with the value of P(H) = 1/N, where N = 10⁹ is the estimated total number of stars on the sky as computed from the average density in SDSS, we can plot the matching probabilities for comparison. Note that large posterior probabilities are not sensitive to small modulations in the density; it changes linearly only for small values of the Bayes factor. The bottom panels of Figure 4 use the same symbols as the top ones to illustrate the derived posteriors using the above constant prior. The difference between the models is possibly even more striking here. While the left panel has very similar estimates from the different models, with time the separations grow large enough to quickly zero out the probabilities for stars with proper motions larger than 100 mas yr⁻¹, whereas the new models keep the probabilities significantly larger. Table 1 contains the measurements for all stars as a function of the time difference. The first column of the table is the identifier of the star, the ObjID in SDSS Data Release 6. The reference proper-motion values are taken from ProperMotions table of the SDSS Catalog Science Archive, which combine the SDSS and the recalibrated USNO-B astrometry for a precise and reliable determination.

We see two important features of the associations using the new proper-motion priors. The constant prior yields lower and lower Bayes factors as the elapsed time increases and the probability drops regardless of the proper motion. Even when the separations are small enough for the static model to perfectly recover the object, the constant proper-motion prior yields a lower 90% probability. The empirical prior always outperforms the static model but, in this two-epoch observation case, the probabilities of the fast stars fall below the constant case or any reasonable probability threshold.

Next, we turn our attention to the potential improvements from including additional epochs to the data sets. For this comparison, we keep the first and last observations in time, hence the baseline is the same for all cases. We add to these two observations additional measurements whose epochs are between them in time. These two-, three-, and four-way associations are shown in the left, middle, and right panels of Figure 5, respectively. It is apparent that adding new detections significantly improves the proper-motion models: the reasonably good associations of only two detections are promoted to essentially certain matches by including intermediate detections. In contrast, the static model continues to reject the associations of all high proper-motion stars. We see that one of the stars with μ ∼ 100 mas yr^-1 actually gets a high probability even in the static model, when the angular separation for only a few years between the epochs is small enough to recover the star.

Zoom In Zoom Out Reset image size

Table 2 shows the measurements as a function of the number of epochs used in the calculations. We see that the empirical prior of the improved model assigns 100% probability of all stars when considering all four epochs, and even the three-epoch computations would yield close to that with the μ = 300 mas yr^-1 star getting a lower 97%. The exception from this is the fastest star at μ = 555 mas yr^-1 in the case of the empirical prior, whose probability is essentially zero in all panels. The reason for this is that this star is one of the highest proper-motion stars in Stripe 82 and even with the generated model stars, which appear in the prior, we have very few (roughly 40) high proper-motion stars.

It is worth reiterating the reason for and the consequence of these results. Associations of more than two detections benefit dramatically more from the proper-motion prior because two points can always be connected with a straight line unlike three or more. In other words, the prior probability of two detections being on a great circle is 100% but for three or more it is small, hence such combinations will get boosted by the alignment. Having seen the convincingly large probabilities for the three-way cases and assuming the same maximum time difference between observations, one can conclude that, for the time intervals we consider here, surveying strategies with three epochs are superior to those with only two, but adding more would not noticeably improve our ability to correctly cross identify the detections.