Overcoming Separation Between Counterparts Due to Unknown Proper Motions in Catalogue Cross-Matching

This paper is split into two parts. First, we detail the construction of a simple analytic model for the statistical distribution of potential proper motions of a source of a given magnitude and sky coordinates. In Section 2 we detail the constituent parts necessary to build the model of proper motions. We describe the process of building the distributions in Section 3, while Section 4 discusses the precision and accuracy of the model at various Galactic sightlines and brightnesses. Here we will return to Gaia, using its high-precision stellar proper motions across many Galactic sightlines to evaluate and corroborate our model distributions.

The second part of this work describes how to include this unknown proper motion distribution – or any distribution of proper motions, theoretical or poorly constrained yet detected – in the cross-matching process. We discuss the mathematical framework necessary for including proper motion drift in the evaluation of the separation between potential counterpart detections in Section 5. We also touch upon how to use these statistical distributions of proper motion drift as a discriminator between stars and galaxies. Concluding remarks are given in Section 7.

In Appendix A we outline the various coordinate systems and derive the transformation matrices used throughout this work, while in Appendix B we detail the inclusion of the proper motion PDFs within a probabilistic cross-match astrometric separation likelihood framework.

The first component in the Galactic motion is the unique velocity of the Sun, relative to the local standard of rest (LSR). The Sun’s motion through the Galaxy will induce a ‘secular’ parallax effect (i.e. distance-dependent, albeit non-periodic, as a trigonometric parallax would be) in the apparent movement of all other sources in the sky, with opposite sign. Hence we need to know the Sun’s motion, relative to this ‘zero point’ motion, the LSR, defined in the Heliocentric Cartesian coordinate frame, $(U,\,V,\,W)$ – velocities corresponding to the $(x,\,y,\,z)$ coordinate system. Here we use the values of Schönrich et al. (2010):

The main component of our model that will dictate the Galaxy-wide observed motions of sources is that of the Galactic rotation, and the stellar streaming. Descriptions of this motion go back to Oort (1927), with the Oort constants describing the motion of stars on closed orbits around the Galaxy. However, through modern kinematic surveys, obtaining the three-dimensional velocities and independent distances to a host of well-characterized objects, it is possible to directly measure the rotation curve of the Milky Way. Thus, we can derive the average tangential velocity of sources orbiting the Galactic center at a given radius, here following Model 3 of Mróz et al. (2019).

Obtaining the rotational velocity $\Theta$ at a given Galactocentric radius $R_{c}$, we can transform this Galactocentric Cylindrical azimuthal velocity into a Galactocentric Cartesian coordinate frame and subtract the Solar peculiar and LSR motion, obtaining $U_{1}$, $V_{1}$, and $W_{1}$ (Mróz et al., 2019, equations 5-7). To do so, we use the transformation

here $R_{\odot}$ is the Solar Galactocentric Cylindrical radius, $d_{\mathrm{ip}}$ is the in-plane distance from the Sun to the particular location, and $l$ is Galactic longitude – see Appendix A.1.1 for details. Deviating from Mróz et al. (2019), we set $(U_{s},\,V_{s},\,W_{s})$, the non-circular motion of the source, all to zero, and hence have

Once we have the Galactocentric Cartesian components of the rotational velocity, relative to the Sun, we can transpose into the in-plane Heliocentric radial and tangential velocities,

along with $v_{z}=W_{1}$, since the two axes are still in alignment. Finally, once we have appropriate Heliocentric Cylindrical velocities, we can construct our latitudinal velocity as

and relate the Heliocentric longitudinal and latitudinal velocities to their proper motions through

where $\pi$ is the parallax of the source. As our distances come from Galactic models, we assume they are not subject to any observational bias or uncertainty, and simply treat $\pi^{-1}=d$. The factor $k$ describes the translation from units of $\mathrm{km}\,\mathrm{s}^{-1}\,\mathrm{kpc}^{-1}$ to $\mathrm{mas}\,\mathrm{yr}^{-1}$, and is given by $k=0.2108\,\mathrm{mas}\,\mathrm{yr}^{-1}\,\mathrm{km}^{-1}\,\mathrm{s}\,\mathrm{kpc}$.

In practice, the TRILEGAL simulations do not provide either distance or parallax, but provide its distance modulus. Hence, to obtain a (three-dimensional) distance $d$ in kpc, we invert the absolute magnitude equation:

where $m-M$ is the distance modulus.

The above equations describe the average Galactic rotation velocity around the center of the Galaxy. However, objects have other sources of velocity that impact their observed proper motions, and this leads to a deviation away from the expected velocity, and thus proper motion. This is termed the asymmetric drift velocity, and essentially controls how much of the theoretical velocity a source should have is taken by other, random motions. Thus, we must include this component in the velocities.

We model three Galactic components (see Sections 2.4.1-2.4.3 for more details of their construction) in our simulations: the Galactic thin and thick discs, and a single Galactic (outer) halo. Each of these is given their own asymmetric drift, as a measure of the levels to which their motions are different from ‘pure’ streaming motion. We assume the thin disc of the Galaxy has an azimuthal asymmetric drift velocity of $10\,\mathrm{km}\,\mathrm{s}^{-1}$ (Robin et al., 2003). As Mróz et al. (2019) used Classical Cepheids in the derivation of their rotation curve, we also assume the rotation curve is valid for the thin disc and already folds in any drift velocity, and therefore just need to consider the relative drift velocities of the thick disc and the halo. We use a thick disc drift velocity of $49\,\mathrm{km}\,\mathrm{s}^{-1}$ (Pasetto et al., 2012a), and give the halo a drift velocity of $240\,\mathrm{km}\,\mathrm{s}^{-1}$ (Golubov et al., 2013), to essentially counteract the motion of the LSR, modelling the halo as stationary relative to the Galaxy. We therefore use $v_{a,\phi}=\{0\,,39\,,230\}\,\mathrm{km}\,\mathrm{s}^{-1}$ for the relative thin disc, thick disc, and halo drift velocities, respectively.

For a given location in the Galaxy, our decomposition of the drift velocity, from Galactocentric Cylindrical coordinates into Heliocentric Cylindrical coordinates, is given by the transformation

with

and

with $v_{a,\phi}$ taking on any one of the three given drift velocities above, depending on which component of the Galaxy is being considered. For details on the derivation of this transformation (rotation and mirror) matrix, see Appendix A.1.2.

The asymmetric drift velocity suggests that some of the motion that ought to be used by a given source in its rotation around the Galaxy is being used otherwise, in a random component. Hence, a collection of sources in a particular part of the Galaxy will have some spread of their velocities around some mean value. Thus, to be able to model our collection of sources in the Galaxy we require a description of the dispersion of the velocities.

Each component of the Galaxy modelled – thin and thick discs, and halo – have their own prescription of velocity dispersion. In addition, when simulating sources, we do not initially know to which component to assign a given simulated object. Therefore we simply generate a set of proper motion realisations for each of the three components, weighted according to their a priori density at that location. Hence, in the following sections we also describe the formulation of each component’s density profile, which are simply re-normalised by the sum of their densities to provide a prior probability, used as the weight for the proper motion distribution. We currently do not consider the Bulge, a common component of Galactic simulation models such as TRILEGAL, in our proper motion model. We chose to ignore this additional component for this initial, exploratory model, focussing on relatively bright sources, with detected Gaia proper motions to compare and verify our model against. This should mean they are sufficiently far from the Galactic center to not be influenced by the Bulge, avoiding the complexities that the inner region of the Galaxy impose on velocities of orbiting stars. However, with the upcoming LSST survey and the need to model much fainter, more distant objects in the next few years, we will investigate a robust Galactic Bulge/Bar model for inclusion within this proper motion framework. We simply conclude, for now, that it would be easy to add additional components, and we could model a simple Bulge component after e.g. Jackson et al. (2002).

Once a given Galactic component has its dispersion vector – its covariance matrix – in the Heliocentric Cylindrical coordinate system, then a realisation is drawn from a multivariate normal, given by

where $i\in\{\mathrm{thin},\mathrm{thick},\mathrm{halo}\}$ and $\bm{\Sigma}{{}^{\prime}}_{\!i}$ is the Cylindrical frame rotated covariance matrix of the $i$th component.

The simulated proper motions are good across all sightlines and brightnesses; some examples are shown in Figures 2-4. We get good agreement in the mean proper motion, and shape of the distributions, of bulk source motions. For most sightline-brightness combinations the agreement is quantitative, while sometimes the shapes are merely broadly in agreement. Disagreement in modal proper motion drift is likely largely caused by our Galactic rotation model not capturing the fine detail of Galactic potentials or inaccuracies in our asymmetric drift velocity, while distribution width issues can mostly be explained by the extrapolation of the velocity dispersion vector. However, we stress that the model’s simplicity is one of its strengths in the context of inclusion within a larger cross-match process, and that these minor differences are more than acceptable for the purpose of improving Bayesian match likelihoods. These distributions are intended to reflect a wide range of potential positional shifts through time, rather than model any one specific proper motion. Hence so long as the rough widths – to within something like a factor two, which we achieve – and mean offsets – good to high accuracy using the Galactic rotation curve – are modelled to reasonable accuracy, our distributions are good enough for our work, and as intended.

Figure 5 shows some reduced statistics for the entire set of sightline-brightness combinations we tested in the Galactic plane – $G=\{12,\,15,\,18,\,19.5\}$, $l$ in the range from $0^{\circ}$ to $345^{\circ}$ in 15 degree intervals, and $b=\{-50^{\circ},\,-30^{\circ},\,0^{\circ},\,25^{\circ},\,40^{\circ}\}$, as well as the Galactic north and south poles at $-90^{\circ}\leq b\leq-80^{\circ}$, $-80^{\circ}\leq b\leq-70^{\circ}$, $70^{\circ}\leq b\leq 80^{\circ}$, and $80^{\circ}\leq b\leq 90^{\circ}$. Overall, we find that the widths of the Gaia data are approximately 80% that of our model (i.e. our model is too wide by 25%) across all positions and brightnesses. There is a roughly 10% spread in relative widths – middle column, Figure 5, cf. Figure 2, bottom right panel, where our red model has a slightly wider wing than the histogram of the black Gaia data. We also see evidence for overly narrow simulated proper motion distributions (Gaia-to-model ratios larger than one) along various sightlines, in approximately 8% of cases – but, again, get extremely good agreement along others. These slightly-too-wide distributions are likely related to our modelled radial and vertical dependencies of the thin disc dispersion vector, as the thin disc is the dominant term at the distances our Gaia data probe. Additionally, as can be seen in the top row of e.g. Figure 2, low-number statistics of brighter Gaia stars could be interpreted as lower standard deviations, as the ‘real’ distributions fail to probe the wings of the simulated distributions to high precision.

Relative to the standard deviation of the distributions, our biases – the mean motion offsets – are within $\sim 10$% two-thirds of the time, and almost always within $15-20$%, as shown by right-hand column, Figure 5. Absolutely, on a decade baseline, we find most of our mean motion offsets are within approximately 0.01 arcseconds (or 0.025 arcseconds on a 25-year baseline, as maybe be important for LSST), well within the ‘centroid’ precision of most photometric catalogues – left-hand column, Figure 5. These results – even where qualitative (e.g. Figure 4, bottom-left panel) as opposed to quantitatively (e.g. Figure 4, top-left panel) good fits – are satisfactory, and we therefore have chosen not to over-explore the residuals, as the subtleties of the spiral arm structure of the Milky Way are outside of the scope of this work.

All previous examples shown (Figures 2-4) were limited in Galactic latitude to $\lvert b\lvert\leq 50^{\circ}$, exploring primarily the proper motions of sources roughly in the Galactic plane. However, we must also verify that our model is good at high absolute Galactic latitudes, where we are viewing sources orbiting around the Galactic center ‘above’ us. As shown in Figure 6, we get good agreement for the Galactic longitudinal and latitudinal proper motions, after removing objects with parallax $\pi<0.05\,\mathrm{mas}$ ($d\gtrsim 20\,\mathrm{kpc}$) or $\pi/\sigma_{\pi}<2$. Here, close to $b=90^{\circ}$, our equations for the average rotational velocities (equations 3-7) simplify somewhat. First, looking straight up out of the Galactic plane, we have $d_{\mathrm{ip}}\approx 0$; we also have $R_{c}\approx R_{\odot}$, and hence $(R_{c}^{2}+R_{\odot}^{2}-d_{\mathrm{ip}}^{2})/(2\,R_{c}\,R_{\odot})\approx 1$, and can assume $\Theta(R_{c})=\Theta_{\odot}$, giving $U_{1}=-U_{\odot}$, $V_{1}=-V_{\odot}$. This further gives $v_{d}=-U_{\odot}\cos(l)-V_{\odot}\sin(l)$, $v_{l}=U_{\odot}\sin(l)-V_{\odot}\cos(l)$, and hence, along with a simplified $v_{b}=-v_{d}$,

This renders the orbital motions effectively just those of the Sun, with our dispersion vector giving good shape agreement to the Gaia data. Overall, we see good agreement in the shape of our model and the Gaia proper motions.

At this point it is worth briefly considering if, or at what brightnesses, this additional information is necessary. Figures 2-6 show a representative sample of Galactic sightlines and the various widths of the distributions of potential proper motions in each sightline-magnitude combination. Overall, the widths of these distributions are approximately $0.2-0.5$ arcsecond drifts over a 10 yr baseline ($20-50\,\mathrm{mas}\,\mathrm{yr}^{-1}$) at the bright end of our tests ($G=12$), reducing to $0.1-0.3$ arcsecond drifts in 10 years ($10-30\,\mathrm{mas}\,\mathrm{yr}^{-1}$) at $G=19.5$.

The first parameter we should compare the proper motions to is the precision on an individual position. If one or both of the positions in a given cross-match were highly uncertain, factors 10 or higher than the proper motion drift, this would dominate over the extra positional spread caused by the potential proper motion of the source. Its inclusion would then not contribute to the determination of potential counterparts. However, for a decade-long baseline, the spread of separations induced by unknown proper motion is at least a factor two or three higher than typical astrometric precisions, with even small motions over long enough baselines moving objects several astrometric precisions apart. For Gaia astrometric precisions are vastly higher; while 80% of its sources will also have incredibly high precision proper motions even the remaining sources will have coordinate positions significantly higher than the unknown proper motion distributions. A more typical ground-based survey, LSST should have at worse 0.07 arcsecond precision on each individual visit at $r=24$ (Ivezić et al., 2019). While this is a factor $\approx 3$ times smaller than the widths of our proper motion distributions on decade baselines, the real power of LSST lies in its repeated observations. Depending on whether the object is in the full ‘Wide-Fast-Deep’ (WFD) survey or in the Galactic Plane footprint, it will either be observed approximately 800 or 150 times across LSST’s full survey lifetime (Bianco et al., 2022). Hence the statistical precision on a co-added detection at $r=24$ is $0.003-0.006$ arcseconds depending on the exact number of visits. Even including $\approx 0.01\,\mathrm{arcsec}$ systematic precision (Ivezić et al., 2019) this is far below the widths of our models for proper motion drift. While those objects will likely also have proper motions after LSST DR3-4, $r=26.5$ coadded detections will have statistical astrometric precisions a factor $\sqrt{10}$ higher, $0.008-0.02$ arcseconds, still a factor 10 or more below our 10-year baseline drift spread. It will be therefore important to take these long-baseline drifts into account for faint LSST objects. For current-generation surveys such as SDSS, its very faintest sources have statistical positional uncertainties comparable to the tightest of our proper motion distribution widths ($\approx 0.2\,\mathrm{arcsec}$) so $r=24$ objects in SDSS may only see limited gains matching across a 10-year timespan. Of course, the drifts increase linearly with time, and so a 15-year baseline (2015-2030, for example, in the case of SDSS-LSST) increases the potential proper motion drifts to a larger impact than astrometric precision. Additionally, in the context of crowded field Bayesian cross-matching, even a ‘one-sigma’ positional movement will be enough to significantly disrupt match likelihoods.

It is also useful to ask if proper motion precisions are ever comparable to the width of potential unknown proper motion. Again, for Gaia this is not the case due to its extremely high precision and repeated observations of all objects. At $G=20$ the median precision on its proper motions are of order $1\,\mathrm{mas}\,\mathrm{yr}^{-1}$ (Gaia Collaboration et al., 2021). For LSST, quoted proper motion uncertainties are also of order $1\,\mathrm{mas}\,\mathrm{yr}^{-1}$ (Ivezić et al., 2019) – but these assume $\approx 800$ visits. Hence the stellar proper motion precisions for Galactic Plane objects will be a factor $\approx 3$ higher due to the reduced number of observations within the same timeframe. However, even $5\,\mathrm{mas}\,\mathrm{yr}^{-1}$ is $0.05\,\mathrm{arcsec}$ over a 10-year timeframe and therefore a smaller, but sizeable, fraction than the unknown proper motion distribution widths. Right at the detection limit of proper motions with LSST it may be the case that it is preferable to not use the detected-but-unconstrained proper motions, though. SDSS has typical limiting proper motion precisions of $5\,\mathrm{mas}\,\mathrm{yr}^{-1}$ by $r\approx 20$. Hence, like LSST, where constrained individual SDSS proper motions brighter than about 20th magnitude are likely preferable to unknown proper motions, but below this limit and the single-visit detection limit of $r=22$ unknown proper motions may be the more precise constraint. In the IR, CatWISE (Eisenhardt et al., 2020) has proper motion uncertainties of $20\,\mathrm{mas}\,\mathrm{yr}^{-1}$ at $W1\approx 15$ (Marocco et al., 2021) and the VVV survey (Smith et al., 2018) cites uncertainties of $10\,\mathrm{mas}\,\mathrm{yr}^{-1}$ around $K_{s}\approx 16$; below these brightnesses the precisions on individual proper motions become comparable to or larger than the widths of typical unknown proper motion distributions.

As discussed in Section 2.4, we do not currently include a full prescription for the Galaxy. In particular, we do not model the Galactic Bulge (or Bar). This may have an effect at very low Galactic longitudes and latitudes. The Gaia data show a broader, almost flat distribution of longitudinal proper motions, where our simpler Galactic model, using the thin disc as the dominant term, has a bi-modal distribution of two narrower peaks, as shown in Figure 7, left hand panel. It can also be seen in the data (Figure 7, right hand panel) that there is a slightly too narrow distribution of latitudinal proper motions, as compared to the Gaia data. This likely comes back to the minor effects of radial dependencies of the $\sigma_{rr}^{2}$ term, either following a Gaussian- or Rayleigh-like distribution (as discussed in Section 2.4.1). However, it could also be the case that our radial and vertical dispersion scalings are failing at these smaller Galactic radii, as a significant fraction of Gaia sources ought to be sufficiently far removed from the Galactic center to be Bar or Bulge objects. Once again, we deem these minor issues beyond the scope of this preliminary work – the combined bi-modal longitudinal distribution almost entirely covers the distribution of Gaia motion drifts, to within better than a factor 1.5 or so, which is our goal. However, we highlight the issue that the very inner few degrees of the Galactic center may suffer systematic proper motion effects due to the nature of the Galactic Bulge and Bar complexities.

We also have not modelled the Magellanic Clouds, and indeed during testing found that several of our test fields are heavily ‘contaminated’ by sitting on the Small Magellanic Cloud (SMC) and Large Magellanic Cloud (LMC). These extra terms, as with the Bulge, would be easy to implement; provided a relative number density of sources, with some positional distribution, and bulk and dispersal proper motion – assuming the Magellanic Clouds are orbiting internally, and around the Milky way – the proper motions can be modelled in much the same way with the Galactic discs and halo. For now, we also simply urge the reader to take care when simulating sources centered on the Magellanic Clouds (SMC $l\sim 300^{\circ}$, $b\sim-45^{\circ}$; LMC $l\sim 280^{\circ}$, $b\sim-35^{\circ}$).

Our model for motions of objects in the plane of the sky assumes all sources are single stars, subject solely to the Galactic potential. However, half of objects are in some form of higher-order system (e.g. Raghavan et al., 2010) and should therefore be subject to additional on-sky motion. It is therefore reasonable to ask whether the non-inclusion of this effect, of unresolved binary objects, would have any impact on our derived proper motions.

If the binary were equal mass, any orbital motion of the two sources around their common barycentre would completely cancel by symmetry, and show no impact on the photocentre and proper motion of the blended sources. On the other hand, if the objects were very unequal in mass then both the barycentre and photocentre of the pair will be dominated by the larger, brighter main source, and effectively reduce to a singular object for our purposes. An object of approximately half the mass of the primary, however, contributes very little in luminosity but significantly in astrometric effects. Placing such an object on a worse-case orbit of approximately $10\,\mathrm{AU}$ would give an orbital period around 25 years, and a half-phase orbit on our key decade-long time interval between photometric catalogue ‘generations’. In that time, the primary object would travel halfway around the orbit, appearing to move a total of $\approx 6.5$AU, twice its orbital distance from the barycentre. At a typical distance of roughly $1\,\mathrm{kpc}$ this is $6.5\,\mathrm{mas}$ or $0.54\,\mathrm{mas}\,\mathrm{yr}^{-1}$. Such a perturbation is well below the of order $10\,\mathrm{mas}\,\mathrm{yr}^{-1}$ widths to the proper motion distributions observed for faint objects in our model. If the object were significantly closer – say $100\,\mathrm{pc}$ instead – the motion effects would be 10 times higher, and comparable to the model widths. In those cases the object would be much brighter, and likely have an individually measured proper motion or be known to be a multiple system through other means. We therefore believe the non-inclusion of higher-order systems is justifiable at the resolution we are aiming to achieve.

As noted in Section 3, the TRILEGAL simulations we use to construct our models of proper motions do not provide individual positions for simulated sources. To overcome this, we simply uniformly distributed sources within the rectangular area we sampled our Gaia proper motions in. This effect may explain some small disagreements between our simulated and Gaia proper motion distributions, as we are therefore not properly modelling any clustering, extinction effects, or other non-uniformity and correlations in the distances and positions of Galactic sources.

However, the effect on each individual proper motion should be relatively small, as the regions in question were mostly limited to several degrees in extent, and $\cos(x+5^{\circ})-\cos(x)\lesssim 0.08$, $\sin(x+5^{\circ})-\sin(x)\lesssim 0.08$ over the entire Galactic longitude. Thus our values for, e.g. the decomposition of $U_{\odot}$, or $\Theta$, within our proper motion equations, are of order 8% wrong at their most extreme, in the case of a simulated patch of sky five degrees wide. As an ensemble, however, this assumption should be a reasonable one, and the ‘incorrectness’ should average out, with uniformity of source distribution acceptable for small enough patches of sky, providing a statistical distribution of variations of velocity decomposition across the whole region.

To compare our ensemble proper motion distribution with the distribution of Gaia proper motions, we included the Gaia measurement uncertainty in our theoretical distribution of drifts. The Gaia data have individual uncertainties but to smooth the model with the uncertainty we had to select a single average value (the error bar included in the corners of sub-plots in e.g. Figure 2). A small part of the discrepancies between model and data in our analysis could therefore stem from this simplifying assumption, with no bearing on the model itself. If the Gaia data have a particularly broad distribution of measurement uncertainties, as they tend to at fainter magnitudes, our single uncertainty value would not produce an uncertainty-convolved motion drift distribution that reflected that of the Gaia data. Testing more complex treatments of Gaia uncertainty distributions in the comparison between model and data, we found that more fully describing the non-singular value of measurement precision did produce theoretical drift distributions that better matched the expected data, but not completely. We therefore still find a few sightlines with slight differences in central proper motion or distribution width, but perhaps 30% of these tensions are explainable by the different measurement precisions of faint Gaia proper motions. Ultimately, however, as mentioned in Section 3, we do not include this measurement uncertainty in the final model, just performing the convolution to compare to the Gaia distributions more accurately.

Throughout this work we have used the TRILEGAL simulations to provide a set of theoretical distances for sources of a particular Galactic sightline and magnitude. We could, of course, use any model of the Milky Way to achieve this, such as the Besançon model (Robin et al., 2003, 2012, 2014, 2017; Czekaj et al., 2014; Bienaymé et al., 2015). With the Besançon models, however, we receive simulated proper motions for the objects returned in our query, unlike with TRILEGAL. We can use these simulated proper motions to further verify the robustness of our proper motion model; but this then raises the question of why we simply do not use these simulated proper motions for use in our cross-matches. We will address both of these issues in the next two sections.

Our model for proper motions assumes the source in question is a star – objects orbiting the Galactic center in some fashion. However, for an all-sky catalogue cross-match we will also, at high Galactic latitudes, be matching a considerable number of galaxies. We therefore need to model the two cases. First, that the sources being matched are stars, and hence have the statistically modelled unknown proper motion distribution, with which we wish to ‘blur’ out our potential match separations. Second, they are galaxies, which have zero proper motion, being altogether too distant to have visibly moved anywhere. We therefore have a slightly different probability of match (sources being ‘counterparts’, under hypothesis $c$) given separation $d$, now also conditioned on the ‘type of source’ hypothesis, which we will denote as $p(c|d,\mathcal{S})$ and $p(c|d,\mathcal{G})$ for a ‘star’ and ‘galaxy’ pairing respectively.

When matching, we are generally only concerned with the overall probability of the two sources having a given sky separation under the hypothesis of their being matched, $p(d|c)$ – this term is denoted $G$ by Wilson & Naylor (2018a). Note that this differs from $p(c|d)$, the probability of the two sources being counterparts given their sky separation, Wilson & Naylor (2018a)’s $g$. $p(d|c)$ we can obtain by the marginalisation over the two hypotheses:

where $P(\mathcal{S}|c)$ is the prior probability that these counterparts (with given sky positions, brightnesses, etc.) are stars (or galaxies, in the opposite case). We work under the assumption there is no third type of object – crudely labelling objects as ‘in the Milky Way’ or ‘outside the Milky Way’ – and thus $P(\mathcal{S}|c)+P(\mathcal{G}|c)=1$.

However, we can also ask a related but separate question: ‘what is the probability that these two detections are of a star, given that they are counterparts with a given separation?’, which looks like

Here we have the likelihood of the separation given the hypothesis that the sources are counterparts and stars, multiplied by the prior chance of the sources being stars given they are counterparts, normalised by the overall chance of either a galaxy or star pair having this particular detection offset. To calculate both $P(\mathcal{S}|d,c)$ and $p(d|c)$ we therefore need both prior and likelihood terms.

Calculating the likelihood terms $p(d|c,\mathcal{S})$ and $p(d|c,\mathcal{G})$ is relatively straightforward, simply being the convolution of the respective AUFs (containing all relevant AUF components for the two catalogues) of the sources in question. For $p(d|c,\mathcal{G})$ this does not include any proper motion terms, as the ‘proper motion model’ for galaxies is a static one – mathematically, this is equivalent to the convolution of $G$ and a delta function at zero proper motion, with $f*\delta=f$ – and hence $p(d|c,\mathcal{G})=G$. For $p(d|c,\mathcal{S})$, however, we wish to include the motion of Galactic sources, and hence subsequently convolve by the $h{{}^{\prime}}_{\!\mathrm{pm}}$ PDF, describing the potential additional on-sky movement due to the epoch difference between the two sets of observations; $p(d|c,\mathcal{S})=G{{}^{\prime}}\equiv G*h{{}^{\prime}}_{\!\mathrm{pm}}$.

Thus, the likelihood for our new question is easy to calculate; we are therefore left with the derivation of the prior, $P(\mathcal{S}|c)$. The ‘conditioned on the fact that the sources are counterparts’ aspect of the prior is tricky to implement in practice. We wish to know, analogous to Wilson & Naylor (2018a)’s derivation of photometric likelihoods, the distribution of stars and galaxies as a function of the two bandpasses in question – e.g. $r$ and $J$, for a match between optical and infrared data. Thus, $P(\mathcal{S}|c)$ is really ‘what is the probability that these two sources are stars given that they are counterparts with magnitude limits (or dynamic ranges) in their respective bandpasses?’, $P(\mathcal{S}|c,m_{\mathrm{lim,r}},m_{\mathrm{lim,J}})$. Due to the nature of the simulated objects – being derived from one-sided distributions, a function of just a single magnitude in one bandpass in one of the two catalogues – we are unable to create two-dimensional relationships between stars and galaxies in the construction of these priors. In fact, our likelihoods, $p(d|c,\mathcal{S})$, should implicitly assume counterparts for sources, but are built from the full distribution of sources of just a single magnitude. Here we could have, for example, a case where sources of $J=17$ either have optical brightnesses $r=18$ or $r=25$ (being two classes of objects at differing distances, say); this distance distribution is blurred into a bimodal proper motion distribution in the IR, but one class of object is rejected if we consider the dynamic range of the optical data for an example $m_{\mathrm{lim},r}=20$.

At present, the explicit dependency on the two-sided, magnitude-magnitude relationship between sources in our two catalogues is beyond the scope of this work, due to the nature of the outputs available from most Galactic simulations being limited to a particular set of bandpasses for a specific catalogue. We therefore simply note here that for now, the construction of these models is one-sided – in contrast to the cross-matching algorithms of Wilson & Naylor (2018a), taking into account both catalogues symmetrically, in both AUF-based astrometry and photometry. We thus sidestep this dependency by constructing our priors on star and galaxy counts on single magnitude source counts, effectively creating $P(\mathcal{S})$ and $P(\mathcal{G})$, removing the dependency on $c$ within the priors. We can still, however, account for the dynamic range of each bandpass within its given catalogue on a per-filter basis, and hence implicitly use $P(\mathcal{S}|m_{\mathrm{lim}})$ in a practical implementation.

With this minor practical dependency removed, we conclude that with the inclusion of a distribution of unknown proper motions for Galaxy-based stars, it is possible to discriminate between stars and galaxies in photometric catalogues. The equations

and

allow for the drift of Galactic sources with time, recovering them as non-static sources. This is the most certain question that can be answered; stars, as shown in e.g. Figure 3, can have a very high probability of small proper motions in certain sightlines in the Galaxy. Thus, zero proper motion does not necessarily mean galaxy; but a combination of delta-function likelihood for Galactic proper motion and imbalanced priors at high Galactic latitudes mean that zero proper motion objects will bias towards being extragalactic. On the other hand, if a source has a proper motion distribution which is significantly non-zero, as is the case for Galactic longitude proper motions at $l=270^{\circ}$, $b=0^{\circ}$ (Figure 8), then we should see a breaking of this degeneracy. The offset between the two sources being considered as potential counterparts should now be able to tell whether the sources are further apart than their respective AUFs would suggest – at which point they are very likely detections of a star – or if they have an offset compatible with their AUFs – at which point they are very likely a galaxy.

We also note that we should also consider the proper motions within the context of non-matches, but it is easy to see that this results in a trivial case, effectively ignoring the proper motions. For the counter hypothesis of ‘these sources are unrelated to one another, and separate detections of two physical sky objects’, each source can have its own proper motion, based on its own statistical distribution of potential motions. In these cases, we need to compare to the hypothesis that these sources are not related to one another given the separation between them. This involves the double, but separate, marginalisation over all possible unknown locations and proper motions, for both objects. Ultimately, as the integrals are separate the proper motions do not affect the end result – see Appendix B.2. This is obvious intuitively: the distribution of separations of unrelated, randomly placed objects is independent of the unknown motion history of those objects.