BP3M: Bayesian Positions, Parallaxes, and Proper Motions derived from the Hubble Space Telescope and Gaia data

To demonstrate the performance of GaiaHub in sparse MW halo fields, we turn to the COSMOS field (Nayyeri et al., 2017; Muzzin et al., 2013) from the Cosmic Assembly Near-infrared Deep Extragalactic Legacy Survey (CANDELS; Grogin et al., 2011; Koekemoer et al., 2011, PIs: S. Faber, H. Ferguson). COSMOS is located at a high Galactic latitude and near the Galactic anti-center ($l=236.8~{}\mathrm{deg}$, $b=42.1~{}\mathrm{deg}$) and boasts HST imaging that covers a large area of sky: $\sim 2~{}\mathrm{deg}^{2}$ was imaged by ACS/WFC in 2003/2004, and 288 arcmin² of that area (i.e. $\sim 4\%$) was imaged again in 2011/2012. From simulated COSMOS-like data – presented in detail in Appendix B – PMs of halo and foreground disk stars in our magnitude range (i.e. $16<G<21.5~{}\mathrm{mag}$) can regularly be as large as $\sim 100~{}\mathrm{mas}~{}\mathrm{yr}^{-1}$ in size.

We analyse all COSMOS HST images within a $0.5~{}\mathrm{deg}$ radius of the field’s center, which is $\sim 39\%$ of the total COSMOS area. Figure 1 shows a histogram of the number of sources matched between the HST frames and Gaia for the 1619 HST images in our analysis. In these HST images, GaiaHub finds 2640 unique Gaia sources, with 4, 9, and 23 sources per image in the minimum, median, and maximum cases. To be clear, there are many repeat exposures at the same HST pointing as well as overlapping and multi-epoch regions of the COSMOS field, so all of the sources appear in more than one of the 1619 images.

A comparison of the Gaia- and GaiaHub-measured PMs for 1264 COSMOS sources is presented in Figure 2. The GaiaHub analysis was run using the --rewind_stars option, which iterates on the transformation fitting using Gaia (in the first iteration) and GaiaHub (in subsequent iterations) PMs to extract better solutions, especially when high-PM stars are present. However, the use of this option removes many sources that have poorly-measured Gaia PMs, resulting in more than half of the sources being rejected by GaiaHub in the COSMOS fields (i.e. 1659 used by GaiaHub from the 2640 possible unique sources). While more relaxed GaiaHub options produce PMs for all 2640 unique Gaia sources in the COSMOS images, these results are sensitive to small-number statistics and poorly-measured Gaia PMs that can shift the reference frame solution and yield extreme disagreements with Gaia.

In general, the GaiaHub and Gaia PMs follow the one-to-one line, and many PM measurements agree within their uncertainties. However, there are also some PMs that show significant disagreement between the two catalogs, as is evident in the rightmost panel. Using Mahalanobis distances between the PM measurements (i.e. accounting for their uncertainties) for sources with $<5\sigma$ disagreement between Gaia and GaiaHub PMs, we find that the GaiaHub PM uncertainties would need to be inflated by a factor of $\sim 1.7$ to explain the PM differences. For $\sim 10\%$ these PMs, the GaiaHub measurements are more than $5\sigma$ different from the corresponding Gaia values, whereas a sample size of $\sim 1000$ should expect to see no deviations this large.

Because GaiaHub is not tuned for sparse fields, it is not surprising to see these disagreements. Some of this disparity comes from the fact that GaiaHub is only able to propagate realistic uncertainties on the transformation parameters to the measured PMs when there are several repeat HST observations of the same field. For much of COSMOS, there are only $\sim 4$ HST observations per field, which is not enough samples to obtain statistically robust propagation using GaiaHub’s approach. While the transformation uncertainty is usually a small component of the total error budget when there are a large number of sources, the transformation uncertainties in the sparse images of COSMOS are significant, so many of the GaiaHub PM uncertainties are underestimated. Handling these sparse field cases therefore requires a new tool.

This subsection will be quite technical in detail, so readers who are less interested in the formal statistics are encouraged to skip to Section 3.2.

First, we define many of the key variables and parameters in Table 1. While we are ultimately concerned with the PMs in RA and Dec, it is convenient to work in a 2D plane projection (e.g. what a detector sees) when comparing images. Because Gaia’s measurements do not correspond to positions in a true image, we follow GaiaHub’s approach of transforming the Gaia coordinates into XY coordinates on a pseudo image using a tangent-plane projection

where the $(\alpha_{G,i}^{\prime},\delta_{G,i}^{\prime})$ are the Gaia-measured coordinates of source $i$ in radians, $\left(\alpha_{G,j,0},\delta_{G,j,0}\right)$ are the coordinates of the center of the Gaia pseudo image when considering HST image $j$, $\left(X_{G,j,0},Y_{G,j,0}\right)$ are the pixel coordinates of the center of the Gaia pseudo image, and $(X_{G,i,j},Y_{G,i,j})$ are the coordinates of source $i$ in the Gaia pseudo image in pixels. The $\left(\alpha_{G,j,0},\delta_{G,j,0}\right)$ and $\left(X_{G,j,0},Y_{G,j,0}\right)$ are inherited from the GaiaHub-measured values for each HST image, which come from an average of the source positions, and these values are not left as free parameters during the fitting process. The chosen pixel scale $\mathrm{PS}_{G,j}$ for each Gaia pseudo image is set to the nominal HST pixel scale that it is paired with (e.g. $50~{}\mathrm{mas}\cdot\mathrm{pixel}^{-1}$ for HST’s ACS/WFC). Then, changes in the RA and Dec can be transformed to changes in XY using the Jacobian matrix:

which is usually very close to

when $(\Delta\alpha*,\Delta\delta)$ are in mas and $(\Delta X_{G,i,j},\Delta Y_{G,i,j})$ are Gaia pseudo pixels. We note that this simplification only holds for small field-of-view detectors. In the cases of large PMs and large time baselines, the off-diagonal elements can start to become important, so we opt to use the more general version in our pipeline.

To set up the probability model, there are a few remaining key terms that we need to define:

which is the transformation matrix for the $j$-th HST frame to the $j$-th pseudo Gaia frame,

which is the GaiaHub-measured pixel position covariance matrix for the $i$-th source in the $j$-th HST frame,

which is the Gaia-measured position covariance matrix for the $i$-th source, and finally

which is the Gaia-measured PM covariance matrix for the $i$-th source. In these equations and the ones that will follow, our convention is to show matrices (and matrices of matrices) using a bold-faced typesetting.

The fundamental relationships between the prior Gaia measurements and the true parameter for source $i$ are given by:

which says that the Gaia-measured values are offset from the true values by noise dictated by their Gaia-measured uncertainties. While we have explicitly included a mean for the Gaia-measured position offset vector $\vec{\Delta\theta}_{i}^{\prime}$ to be as general as possible, in practice $\vec{\Delta\theta}_{i}^{\prime}=\vec{0}$ because Gaia has no expected offset from the positions it reports. In the form of the probability statements above, we have assumed that there is no correlation between, for example, parallax and $\mu_{\alpha*}$, which is not exactly correct because Gaia has measurements for these correlations. Our choices, however, make the following math easier, though it comes at the cost of some additional constraining power from the Gaia priors being ignored. The Gaia-measured correlations between these parameters are indeed included in future sections of this work when we compare the BP3M distributions to Gaia’s results. Future versions of the BP3M pipeline will work to incorporate these prior correlations, which will lead to even tighter posterior distributions on position, parallax, and PMs.

While it may seem nonphysical, Gaia observed parallax measurements can be negative (Lindegren et al., 2018), but it is important to remember that the observed parallax values define the mean of a distribution whose width/uncertainty usually places a large amount of probability in positive parallaxes; in this way, we must treat the Gaia astrometric measurements as distributions and not individual points. Luri et al. (2018), for example, explain how the definitions of motion on the sky (both by Gaia as well as in this work) technically allow for negative observed parallax values, which is especially likely to occur when the position uncertainty is relatively large compared to the size of the parallax motion (e.g. see their Section 3 and Figure 2); as a result, they remind the reader that the Gaia observed parallaxes should not be thought of as a direct measurement of distance, and instead, distances need to be estimated by proper statistical modelling of the information contained in the astrometric solution distributions.

Because some of the sources in an HST image have no Gaia-measured parallaxes or PMs, we find it useful to put a population/global prior on the PMs and parallaxes:

which says that there are some global distributions that the true parallaxes and PMs of the sources originate from. While we are free to play with the parameters of the population distributions, we note that they do need to be Gaussian in form so that we retain the necessary conjugacy. We choose to use diffuse hyperpriors, with the goal of minimally impacting the sources with Gaia-measured parallaxes and PMs while offering some guidance to the sources without. In the future, the parallax global prior could be made more constraining in cases where the distances are better known (e.g. clean populations of extra-Galactic stars) or using a magnitude-dependent parallax prior to better incorporate our understanding of the distribution of stars in the MW. Similar changes could also be made to the PM global prior. For the current version of the pipeline, which focuses on faint stars in the MW stellar halo and extra-Galactic sources, we choose the parallax prior to have a mean of 0.5 mas and width of 10.0 mas; when generating synthetic MW thick disk and halo stars (described in Appendix B), we find that these choices contain 99.9% of the synthetic parallaxes within $1\sigma$ and 100% of the synthetic parallaxes within $3\sigma$ (maximum parallax of $\sim 28$ mas). For the PM prior, we use the sources with Gaia-measured PMs to estimate a mean and covariance matrix, and then we multiply that covariance matrix by a factor of $10^{2}$ to guard against the possibility that some of the sources without Gaia PMs are significantly different in their PM from the other sources.

The information linking Gaia to HST image $j$ for source $i$ is then given by:

which says that the measured offset between the Gaia and HST positions (after applying the transformation) in the pseudo Gaia image is distributed around the offset implied by the sum of the motion from PM, parallax, and uncertainty in position. As mentioned back in the description of Equation 1, we are able to hold $(X_{0,j},Y_{0,j})$ fixed as constants while the $(W_{0,j},Z_{0,j})$ are left as free parameters during the fitting. We note that the relationship in Equation 6 ignores the impact that radial motion has on changing the distance to a star between observations (e.g. see see Section 2.3.1 of van de Ven et al., 2006) – and therefore the magnitude of the PM and parallax – though this is a safe assumption because the distance change is almost always extremely small (e.g. a star with LOS velocity of $1000~{}\mathrm{km}~{}\mathrm{s}^{-1}$ would only experience a distance change of 0.01 pc after 10 years).

The $\vec{\Delta\mathrm{plx}}_{i,j}$ term, which we refer to as the parallax offset vector, is a 2D vector that defines the direction and magnitude of the offset between two observations as a result of parallax motion for a source with a parallax of 1 mas; in this way, the parallax offset vector can be multiplied by any parallax value (i.e. $\mathrm{plx}_{i}\cdot\vec{\Delta\mathrm{plx}}_{i,j}$) to find the appropriate offset in $(\Delta\alpha*,\Delta\delta)$ between two observations as a result of parallax motion. To measure the $\vec{\Delta\mathrm{plx}}_{i,j}$, we use the HST and Gaia observation times, $(t_{H,j},t_{G})$, the position of the source in Gaia, $(\alpha_{i}^{\prime},\delta_{i}^{\prime})$, and built-in functions of astropy (Astropy Collaboration et al., 2013, 2018, 2022). Examples of the parallax motion for different positions on the sky are shown in Figure 3, where the Gaia observation time (i.e. J2016.0) is at the origin and the orbits trace out the parallax motion over the course of a year; the parallax offset vector at any time is simply the vector that connects the corresponding point on the ellipse to the origin.

With the above definitions in hand, we can use the distributions on the $\vec{\Delta d}_{G,i,j}$ vectors for all the sources in all images we are considering to construct the following likelihood distribution:

We note here that we have purposefully chosen to omit writing out the dependency on a few of the explanatory variables (e.g. $\Delta t_{j}$, the covariance matrices) for ease of reading the math.

From Bayes’ Law, we arrive at the following posterior:

where $p(a_{j},b_{j},c_{j},d_{j},W_{0,j},Z_{0,j})$ is the prior distribution on the transformation parameters and the remaining distributions have been defined above. In many of our initial analyses, we chose to use a flat prior on the transformation parameters, which yielded reasonable results. However, we soon realized from our analysis of a few hundred HST ACS/WFC images that there are relatively tight constraints for our expectations of the pixel scale ratio and skew terms. To that end, when using HST ACS/WFC data, we choose to employ a few relatively diffuse/not-overly-constraining priors on the skews, pixel scale ratio, and HST angle, as well as the $(W_{0,j},Z_{0,j})$ vector:

where we use the same prior on both the on- and off-axis skew terms. The prior means of each distribution come from the parameters measured by the previous iteration of transformation parameter fitting; in the case of the first iteration, we use the GaiaHub outputs. In this way, we tell the pipeline that the transformation solution is likely nearby the previous iteration’s solution, while the relatively large widths allows the values to change significantly if necessary. To transform from priors in angle, skews, and PSR, we need to account for a transformation Jacobian:

such that

It should be noted that these priors may not always apply to all sets of ACS/WFC data, so this may need to be adjusted for other fields/applications.

The complete forms of the posterior conditional distributions for each motion component or source $i$ are listed in Appendix A, defining the following important distributions:

where the ellipsis refers to the other variables that each distribution depends on. One key takeaway from the statistics that went into constructing the results in Appendix A is that the conjugacy between all distributions (i.e. all distributions are normally distributed/Gaussians) allows us to fairly easily combine multiple distributions to arrive at well-defined and closed-form posterior distributions.

With these posterior conditionals in hand, we can perform the following steps for each source for a given set of transformation parameters:

1. 1.

   Draw samples of $\mathrm{plx}_{i}$ from $p(\mathrm{plx}_{i}|\dots)$;

2. 2.

   Use those $\mathrm{plx}_{i}$ samples to draw samples of $\vec{\mu}_{i}$ from $p(\vec{\mu}_{i}|\mathrm{plx}_{i},\dots)$;

3. 3.

   Use those $\left(\vec{\mu}_{i},\mathrm{plx}_{i}\right)$ samples to draw $\vec{\Delta\theta_{i}}$ samples from $p(\vec{\Delta\theta_{i}}|\vec{\mu}_{i},\mathrm{plx}_{i},\dots)$.

Once we have samples of the PMs, parallaxes, and position offsets for each source, we can calculate the posterior probability of a set of transformation parameters. We do this by marginalizing over the individual samples of the PMs, parallaxes, and position offsets using Bayes’ Law such that

which is independent of the particular values of the PMs, parallaxes, and positions of each source. This relationship is the key that makes the simultaneous fitting of the 6 transformation parameters per image and 5 motion parameters per source feasible; because we can quickly draw the PM, parallax, and position samples directly from a known posterior distribution given a set of transformation parameters, we can efficiently calculate the posterior probability of the transformation parameters alone. Because we can calculate posterior probabilities for a given set of transformation parameters, we are able to sample those parameters from the posterior distribution using a Metropolis-Hastings (MH) MCMC algorithm. Then, for each sample of transformation parameters, we can sample from the posterior conditional distributions of parallaxes, PMs, and position offsets for each source. In this way, the uncertainty of the transformation fitting is propagated to the resulting PM, parallaxes, and positions.

In terms of implementing the statistics in the BP3M pipeline, the general steps are as follows:

1. 1.

   Read in the position, parallax, PM data from Gaia, where it exists, and the corresponding positions in the HST frames from the GaiaHub output catalogs³³3To be clear, GaiaHub analyses the HST images, and BP3M operates on the outputs of the HST image analysis from GaiaHub. BP3M does not interact with the raw HST images, though we will use “HST images” as a synonym for “GaiaHub-produced HST analysis catalog” throughout the text of this work.;

2. 2.

   Use the GaiaHub transformation parameter values as starting guesses;

3. 3.

   Using sources with Gaia priors on PM, estimate the global PM distribution;

4. 4.

   Draws samples for the transformation parameters using MH-MCMC, evaluating posterior probabilities using only sources with Gaia priors when measuring the transformation;

5. 5.

   Given a set of transformation parameters, draw samples of positions, parallaxes, and PMs for all sources;

6. 6.

   Identify outliers (e.g. bad cross-matches between Gaia and HST) as sources that have more than $2\sigma$ disagreement between the expected and observed positions in each HST image;

7. 7.

   Re-estimate the global PM distribution using the new posteriors;

8. 8.

   Repeat the fitting process using the non-outlier sources, including sources without Gaia priors;

9. 9.

   Compare the new list of outliers to a union of the new list and the most recent previous list of outliers. If the new list length is more than 10% different from the length of the union, repeat the fitting with the new list of outliers, and stop otherwise.

This process requires a minimum of two iterations to achieve good results, though it isn’t uncommon for the outlier list to change enough that a third iteration is required. In terms of processing time, running BP3M on a 2016 MacBook Pro for an HST image with $\sim 200$ sources takes approximate 15 minutes to complete its analysis, while an image with $\sim 10$ sources can be analysed in around 5 minutes. Increasing the number of HST images being analysed together effectively multiplies the computation time by the number of images; this is largely because we need to fit 6 additional transformation parameters for each additional image, and necessarily need to increase the number of MCMC walkers and number of MCMC iterations. Best practice is to analyse all images independently before combining them to save time on searching the transformation parameter space.

Because the HST images are mapped onto the global reference frame of Gaia, the transformation parameters that BP3M measures are also useful in constraining non-Gaia sources. Once the best set of transformation parameters between HST and Gaia have been measured, we can return to the list of all sources in the HST image as measured by GaiaHub, which can include sources much fainter than Gaia can see (i.e. $G>21.5~{}\mathrm{mag}$). When comparing multiple HST images at various epochs, these faint sources can be cross-matched with each other to recover their PMs and parallaxes to much fainter magnitudes. While the pipeline currently offers this feature, it is currently untested and optimized, but preliminary results suggest this approach will be fruitful. It may prove particularly useful in very sparse fields where the number of Gaia sources in each HST image is small, but the shared source list between the HST frames is large.

One key assumption that is built in to GaiaHub – and therefore BP3M – comes from how HST sources are cross-matched with Gaia. An HST source is matched with a Gaia source if they are nearest neighbours within some angular distance of each other (e.g. GaiaHub default of 5 HST pixels, $\sim 250~{}\mathrm{mas}$, though this can be changed), which ultimately sets upper limits on the sizes of PMs that can be measured. This works quite well in medium to low density regions where the stars are likely far enough from each other with respect to their motion between images. In high density regions or for fast moving stars, the closest pair of sources in successive images are not necessarily the correct matching. While future versions of our pipeline may adjust the cross-matching technique to reduce this confusion (e.g. rewind the Gaia sources using Gaia-measured PMs, where they exist, before cross-matching and defining matches probabilistically), we stress that our work focuses on the low to medium density regions where the cross-matching assumptions are valid.

As a back-of-the-envelope test, the cross-matching technique is a good assumption when the average distance between stars in an HST image is twice as large as the average position change from PM:

where $\rho_{*}$ is the stellar number density in area on the sky. We can use this equation to determine different limits for choices of time baselines, densities, and average PM sizes. For example, with ACS/WFC’s $4096\times 4096$ pixel detector with pixel scale of $50~{}\mathrm{mas}\cdot\mathrm{pixel}^{-1}$, a time baseline of 15 years, and an average PM of $100~{}\mathrm{mas}~{}\mathrm{yr}^{-1}$, we find that there would need to be greater than $4660$ stars in the image for there to be a significant amount of confusion in the cross-matching; as these numbers are similar to HST images in the COSMOS field, this implies that sparse regions in the halo have very low risk of incorrect cross-matching. In denser regions, like nearby galaxies, a time-baseline of 15 years and average PM of $5~{}\mathrm{mas}~{}\mathrm{yr}^{-1}$ implies that a limiting stellar number density of $1/9~{}\mathrm{pixel}^{-2}=4.4\times 10^{-5}~{}\mathrm{mas}^{-2}$. In practice, we would like the factor of 2 to be even larger (e.g. 5 or 10) to be safe, but this sets the threshold. Future versions of the pipeline will likely explore improvements to the cross-matching between HST and Gaia using the posterior positions and motion measurements from BP3M.

While the comparison with real data will follow in Section 5, we begin by testing the pipeline on synthetic data to see how well we can recover the input transformation parameters and stellar motion (i.e. position, parallax, and PM) used to generate the synthetic HST image. Figure 4 shows the posterior distributions on the transformation parameters that BP3M measures for one synthetic HST image that has 200 sources and a 15 year time baseline from Gaia; the black points and histograms show the posterior draws, while the blue lines show the true values used to generate the image. The posterior distribution and the input values agree very strongly with each other. The chosen input transformation parameters are representative of real transformation solutions as measured by GaiaHub.

We repeat this process of generating synthetic images while changing the number of sources and the time baseline, and then measure how far the posterior distribution is from the truth. Specifically, we use the posterior samples of the transformation parameters (e.g. the black points in Figure 4) to define a 6D posterior mean vector, $\vec{v}$, and corresponding $6\times 6$ covariance matrix, $\boldsymbol{V}_{\vec{v}}$. The 6D distance between the posterior mean vector and the truth vector, $\vec{v}^{\prime}$, is then defined using the Mahalanobis distance metric

This is analogous to the 1D case of dividing the absolute difference between a mean value and the truth by the uncertainty. Likewise, the units of $D$ can be thought of as the number of $\sigma$ between the truth and mean, and we will use this definition (for varying dimensions of $\vec{v}$) when studying residuals throughout this work. When we analyse 300 synthetic HST images with time baselines between 5 and 15 years from Gaia and 5 to 200 sources per image, we get the blue histogram in Figure 5; the orange distribution shows the expected outcome – namely, a $\chi$ distribution with 6 degrees of freedom and a scale of 1.0 – which agrees well with the measured outputs. This result shows that the pipeline does a good job of recovering the input transformation parameters within the posterior uncertainties.

Next, we analyse how well the pipeline is able to recover the true motions of the sources in each synthetic image of MW halo stars. We also analyse the synthetic data using a GaiaHub-like approach, though we cannot pass the exact synthetic data to GaiaHub because it requires raw HST images that are not currently included in the simulations. Still, the GaiaHub analysis of synthetic data is similar to what GaiaHub would produce after iteration on the transformation parameter fitting (e.g. using the --rewind_stars option), which makes it the same as the analysis mode used to analyze the real COSMOS data in Figure 2.

An example comparison of the PMs for Gaia, GaiaHub, and BP3M is shown in Figure 6 for the same synthetic image that is considered in Figure 4, where in each panel we compare the observed PMs to the known truth. When showing GaiaHub and BP3M results (top middle and right panels), we color the sources by whether there are Gaia-measured parallaxes and PMs to be used as priors; green or blue points designate sources with (“Gaia Prior”) and magenta or orange points designate sources without (“No Gaia Prior”). To be precise, all of the sources have Gaia-measured priors on their position at J2016.0, but not all sources have Gaia priors on their parallaxes and PMs.

The first key takeaway from this figure is that the GaiaHub and BP3M PMs are clustered closer to the true values (i.e. the origin in the upper panels) when compared to Gaia alone: the black points in the Gaia panel represent the same sources as the green points in the GaiaHub panel and the blue points in the BP3M panel. In the GaiaHub panel, the green points near the origin have a systematic trend along the $y=x$ line that is not apparent in the BP3M panel. Investigating this trend using different synthetic data with different PM populations, it seems that large PMs are likely the cause. When focusing on smaller-magnitude PM populations (e.g. the extragalactic populations that GaiaHub was developed for), this systematic trend disappears. As a result, the PM uncertainties in the GaiaHub panel (and the corresponding points in the lower panels) need to be $\sim 2$ times larger on average to explain the difference between the measured PMs and the truth, whereas the BP3M PM uncertainties completely capture any differences with the truth. While GaiaHub and BP3M both show improvements over Gaia alone, BP3M handling of large PMs (e.g. as expected in the MW halo) and the accuracy of the PM uncertainties help motivate why BP3M’s creation was necessary.

Next, the bottom panels show that the size of the PM uncertainty has decreased significantly when HST information is combined with Gaia using BP3M; the left panel compares the PM uncertainties, while the right panel shows how much smaller the BP3M PM uncertainties are, compared to the Gaia PM uncertainties, as a function of magnitude. As mentioned in the previous paragraph, the GaiaHub uncertainties would need to be a factor of 2 larger than the values shown in this figure to explain the offsets of the measured PMs from the truth. Overall, this figure shows that combining the HST and Gaia data with BP3M not only increases PM precision, it also increases accuracy. To be clear, BP3M is not simply shrinking the Gaia PM distribution around the same mean as a result of an increased time baseline and number of images, it is truly improving the astrometric solution of each source.

The particular pattern of PM uncertainties (i.e. a lower branch around $0.2~{}\mathrm{mas}~{}\mathrm{yr}^{-1}$ and an upper branch around $1.5~{}\mathrm{mas}~{}\mathrm{yr}^{-1}$) is a consequence of our choices in modelling the synthetic data, which is detailed in Appendix B. As a brief summary, the HST position uncertainties are based on real GaiaHub-analysed sources as a function of magnitude in COSMOS. At a given magnitude, some of the sources have well-measured positions, while others are less constrained, which leads to the PM uncertainty bifurcation shown in Figure 6.

When considering the uncertainty on a vector, there are a few different approaches one can take. In this work, we are concerned with 2D vectors (e.g. positions, PMs), 3D vectors (e.g. parallax plus 2D PM), and 5D vectors (i.e. 2D position, parallax, and 2D PM) of the different motion components. Because the different components of these vectors can have differing units, and because there can be substantial correlations in the covariance matrices we measure, we choose not to use the standard metric of the quadrature sum of the individual uncertainties of each component. Instead, when we compare vectors that have associated covariance matrices, we are interested in how much the size of that covariance matrix has changed, which includes the effect of correlation. We define the size of a vector’s uncertainty to be the determinant of the covariance matrix to the $1/(2d)$ power, where $d$ is the vector’s number of dimensions:

where $\boldsymbol{V}_{\vec{v}}$ is the covariance matrix of $\vec{v}$. With this definition, the area/volume created by the vector’s covariance matrix, i.e. the determinant $|\boldsymbol{V}_{\vec{v}}|$, is equal to the area/volume of a purely diagonal covariance matrix:

In this case, the resulting uncertainty size allows for the correlations between parameters to impact the certainty about the vector’s position, which serve to shrink the volume defined by that covariance matrix. For the remainder of this work, where we compare the uncertainty sizes of different vectors between Gaia and BP3M, we will be using the vector uncertainty size definition of Equation 14.

As an illustrative example, a highly-correlated measurement might have large uncertainties in all of the individual vector components, but the probability distribution implied by its covariance matrix covers only a small volume of parameter space owing to the high correlation between the components. Here, the quadrature sum of individual uncertainties would yield a large result, implying we know little about the true value of the vector, whereas the vector uncertainty size in Equation 14 would return a small value. This small value tells us that the volume of possible values that the vector can occupy is quite small because of the relationship between the dimensions. For a 2D vector $\vec{v}^{T}=(x,y)$, the vector uncertainty size is given by:

where $\sigma_{x}$ and $\sigma_{y}$ are the corresponding uncertainty in each of the vector components and $\rho_{x,y}$ is the correlation coefficient between $x$ and $y$.

To account for differences that individual realizations might have on the posterior BP3M PMs, we repeat the measurements of Figure 6; that is, we create 5 separate synthetic HST images that have 200 sources and a time baseline of 15 years. For the different realizations, the list of randomly chosen sources all come from the same synthetic catalog of COSMOS-like data, but each realization will have a slightly different set of true PMs, parallaxes, positions, magnitudes, and numbers of sources with/without Gaia priors. The differences of the posterior PMs from the truth for these $5\times 200=1000$ sources are shown in Figure 7, where the sources with Gaia-measured parallaxes and PMs are in blue and the converse are in orange. As expected for a well-behaved pipeline, the difference distribution clusters around the origin, and we see that the sources with Gaia priors on their motion have smaller posterior uncertainties than those without.

Of course, the PMs are only 2 of the 5 dimensions of the vector measured for each source, so we also compare the 5D vector (2D position, parallax, 2D PM) to the truth, again using the distance definition of Equation 13; this distribution for the 1000 stars in Figure 7 is given as the blue histogram in Figure 8, and it shows remarkable agreement with the expected $\chi$ distribution with 5 degrees of freedom and a scale of 1. These figures show that the pipeline recovers the true positions, parallaxes, and PMs of all sources within their posterior uncertainties.

Finally, we explore the impact that the number of sources in an image have on the PM improvement factor when comparing BP3M to Gaia. An example analysis of a synthetic HST image with 10 sources and a 15 year time baseline is given in Figure 9. While the BP3M pipeline again finds good agreement with the truth, the median PM improvement factor over Gaia alone is only 1.39 for the sources in this synthetic image, which is much smaller than the median factor of 8.3 seen for the sources in the $20.5<G<21~{}\mathrm{mag}$ range of Figure 6. This is largely because the smaller number of sources aren’t able to place as strong a constraint on the transformation parameters’ distribution, and the larger transformation parameter distribution propagates to the PM distributions of individual sources.

Comparing the GaiaHub and BP3M panels, there is good agreement between the PM means for this particular realization of synthetic data, which should not be taken for granted considering the large GaiaHub outliers that can appear in real COSMOS data (e.g. $\sim 10\%$ of the PMs in Figure 2). When we compare the PM uncertainties, however, we see significant differences between the two pipelines. In fact, the GaiaHub uncertainties would need to be increased by a factor of 8 to explain the differences with the true PMs, while the BP3M PMs agree with the truth within their uncertainties. Some of the PM uncertainty differences come from the fact that the GaiaHub analysis is only able to propagate realistic transformation solution uncertainties to the final PMs if there are several HST images of the same field, which is not the case when considering a single synthetic image. The impact of the transformation uncertainty is likely negligible for relatively large numbers of sources (e.g. $N>100$, as discussed in del Pino et al., 2022), but it becomes quite important in sparse fields and when there are only a handful of repeat HST observations of the same field (e.g. $N=4$ for many fields in COSMOS).

While it is fairly straightforward to estimate how much the GaiaHub uncertainties need to be inflated by using comparisons with the Gaia-measured PMs where they exist, it is much more difficult to assess and correct for the presence of systematics (e.g. as seen in the GaiaHub PMs of Figure 6) in sparse fields. Fortunately, these sparse halo conditions are what BP3M was developed for, and the synthetic data analysis provide evidence that the pipeline provides accurate PMs and uncertainties here.

We leave out example figures exploring the effect that a changing time baseline has on the posterior PMs because it follows expectations exactly; namely, for the same position uncertainty between two images, a smaller time baselines leads to larger PM uncertainties (e.g. Equation 2 of del Pino et al., 2022).

Our synthetic analysis also allows us to measure the improvement factors in the parallaxes and 2D position vectors of the sources. These results show, not unexpectedly, that the precision of these measurements is strongly tied (but not limited) to the magnitude of the source, the time baseline between images, the number of sources per image, and the effect of the position on the sky on the shape and orientation of the parallax motion. While we cannot change the time baselines of archival HST images, we can explore how changes in the observations times would have impacted the resulting BP3M parallaxes and positions. These lessons are particularly useful for planning of future observations using any telescope.

As before, we create synthetic HST images of COSMOS stars, but in this analysis, we do not let the HST position uncertainty of sources change as a function of magnitude based on GaiaHub examples (e.g. as discussed in Appendix B). Instead, we choose to set the HST position uncertainty to 0.01 pixels, which is the centroid uncertainty that previous studies have measured for well-behaved stars (e.g. Bellini et al., 2011). We adopt the single, well-behaved value to limit the effect that any particular choice/realization of position uncertainties may have on the conclusions we make about optimal observing strategies. Similarly, we choose to use 50 stars in each synthetic HST image to reduce the impact of transformation solution uncertainty as seen in low-source-count images.

Our initial tests suggest that there is minimal improvement on the Gaia-measured parallax or position vectors when using only one HST image (i.e. one HST in combination with Gaia). Additionally, multiple HST images taken very nearby in time (e.g. the multiple exposures taken to guard against cosmic rays) effectively act as a single epoch and are useful for improving the position uncertainty, but don’t offer much information to the parallax. When multiple HST images with common targets and substantial enough time baselines are analysed together, then the parallaxes and positions begin to show improvement. To better understand how parallax and position precision depend on observation times, we choose to simulate 3 epochs of COSMOS observations with HST, with time offsets from Gaia near 12, 8, and 4 years. Specifically, we make the observations occur at the following times:

1. 1.

   No Gaia Offset: All HST observations at year-multiples of the Gaia observation time,
   $\Delta t=\left[-12.0,-8.0,-4.0\right]~{}\mathrm{years}$;

2. 2.

   Half Gaia Offset: All HST observations at half-a-year offset from the Gaia date, $\Delta t=\left[-11.5,-7.5,-3.5\right]~{}\mathrm{years}$;

3. 3.

   Apocenter-then-Quarter: The earliest HST observations at the parallax orbit apocenter farthest from Gaia, then offset by a quarter year for the other epochs,
   $\Delta t=\left[-11.61,-7.36,-3.86\right]~{}\mathrm{years}$;

4. 4.

   Alternating Apocenters: The earliest HST observations at the parallax orbit apocenter farthest from Gaia, then alternating offsets by half a year for the other epochs,
   $\Delta t=\left[-11.61,-7.11,-3.61\right]~{}\mathrm{years}$.

The resulting position, parallax, PM uncertainties from analysing these HST observations concurrently is given in Figure 10, where the two bottom-most panels show the 3D (i.e. parallax, PM) and 5D (i.e. position, parallax, PM) vector uncertainties. This figure reveals many important lessons about the impact that relatively small time offsets of observations within a year have on the stellar measurements.

From the top three panels, we see that different observing strategies lead to significant differences in the resulting position and parallax precisions, though there are minimal differences between the PM uncertainties. This implies that improvements in parallax and position can be achieved with no cost to the PM precision. We emphasize, however, that this is only after analysing all 3 HST epochs together; when considering a single epoch at a time, sources without Gaia priors (i.e. $G>21~{}\mathrm{mag}$) do indeed have significantly larger PM uncertainties if there is an offset from the Gaia date as a result of the degeneracy between PM and parallax motions. Our results suggest that completed surveys with multiple HST epochs may be able to improve position and parallax measurements for free through careful design.

Of the observing strategies we consider (including many that are not represented in Figure 10), the best choice for improving positions and parallaxes was to use an alternating-parallax-orbit-apocenters approach, which is the 4th option in the list above and the red points in the figure. This strategy suggests that the parallax is most improved – and therefore the degeneracy between PM and parallax motions most disentangled – when we consider the average parallax orbit of the sources in a given field (e.g. the blue orbit in Figure 3 for COSMOS). The parallax ellipse can be used to identify the points in time that are most distant from each other (i.e. the extremes of the semi-major axis), which are therefore the easiest to tell apart for a given position uncertainty. Then, the best approach is to first observe at the time that is most offset from the Gaia observation time, and subsequent observations alternate by half a year such that they jump back and forth between the two apocenters. This approach also appears to have a large improvement on the position uncertainty, likely caused by the repeated sampling of the source’s position at the same time in the parallax orbit.

We analyse HST images using BP3M in nearby dSph galaxies to see how the posterior PMs and parallaxes compare to Gaia. This serves as a test with real data, as well as more proof that BP3M is improving the astrometric solutions for the sources it measures. We also show how analysing multiple HST images together impact the resulting PM and parallax precision.

A comparison of Gaia, GaiaHub, and BP3M measured PMs is presented in Figure 12 for the Fornax dwarf spheroidal. For the BP3M PMs, 6 HST images at the same epoch (time baseline from Gaia of 12.8 years) with a total of 198 unique sources are analysed together. The GaiaHub PMs are the same as in presented in del Pino et al. (2022) from analysis of the same 6 HST images, but in this case, a co-moving sample has been defined before transformation fitting to identify likely cluster members. While this figure is similar to Figure 6, the top panels are slightly different in that they show PM measurements instead of offsets from known PMs of synthetic data. Sources without Gaia-measured parallaxes and PMs are colored orange in the BP3M panel, which show larger uncertainties than the sources with Gaia priors in blue. The axes of the GaiaHub and BP3M panels have been set to cover the same range of values for easier visual comparison, and this range is $\sim 10$ times smaller than in the Gaia panel. As before, the lower panels compare the PM uncertainties. Figure 12 shows another PM comparison, this time using 10 HST images of the Draco dwarf spheroidal, taken at 3 different epochs (baselines of 11.2, 9.2, and 2.2 years from Gaia).

First, these figures show a significant tightening of the PMs when Gaia and HST are combined using either GaiaHub or BP3M. This suggests that the increase in PM precision over Gaia alone is not at the cost of PM accuracy. In comparing the PM uncertainties, we find that the GaiaHub values are slightly smaller than those of BP3M, but this is likely because BP3M propagates transformation parameter uncertainties to the final PMs and GaiaHub uses a single set of values for the transformation solution. We also see that the GaiaHub PMs form visually-tighter PM knots than the BP3M PMs, though we would expect them to be more similar based on the size of their PM uncertainties. These differences suggest that there are systematics in the BP3M PMs that are not present in the GaiaHub PMs. In testing possible explanations, we find that the BP3M results are able to look more like the GaiaHub ones when we restrict the samples to only the cluster members from del Pino et al. (2022). Theoretically, the BP3M statistics should not allow the presence of non-member sources to impact the PMs of member stars – in fact, their presence should help improve the transformation solution – but we find that this is not the case.

Some further testing suggests the likely culprit comes from BP3M’s estimate of the source position uncertainty in the HST images. As discussed previously, BP3M uses the same HST position uncertainties as GaiaHub, and this is calculated as a fraction of a PSF quality of fit statistic. When analysing a kinematically-cleaned, co-moving population of stars (i.e. the case that GaiaHub was designed for), this estimate of uncertainty performs well because the sample is restricted to faint sources that follow a linear relation between the quality of PSF fitting and HST position uncertainty. The validity of this approximation is evident from the resulting GaiaHub PMs and the results when BP3M uses the same subsample. However, this uncertainty approximation can occasionally underestimate the true uncertainty – especially for bright-but-not-saturated sources – leading to cases where the transformation solution becomes torqued away from the truth by a few unfairly-high-weighted sources. As more HST images are analysed together (e.g. as we present here with the dSph data), the amount that HST information contributes to the posterior measurements begins to outweigh the Gaia-measured priors, and it becomes more important for the HST uncertainties to be more accurate. It will be important for future versions of BP3M to identify other methods of measuring the HST position uncertainties to be as statistically robust as possible.

For all the BP3M analyses presented in this work, we emphasize that we use all Gaia sources found in each HST image to test the general nature of the pipeline. Our goal is to retain as many sources as possible for transformation fitting, which is especially important in sparse field cases. For BP3M, foreground bright stars – which often have well-measured Gaia positions, parallaxes, and PMs – that do not belong to the distant populations of interest are often useful anchors for the transformation fitting.

We next repeat the BP3M analysis while varying the number of HST images considered together to explore the effect this has on the posterior PMs and parallaxes. The resulting improvement factors as a function of magnitude are given in Figures 13 and 14 for Fornax and Draco respectively, with the legend displaying the time baselines from Gaia for the different sets of HST images. In the case of Fornax, the most significant improvement in PM occurs when analysing a single HST image with Gaia (i.e. a median PM precision increase by a factor of $\sim 8.6$ for $20.5<G<21~{}\mathrm{mag}$), though there is still an additional improvement when the remaining five HST images are brought in ($\sim 14\%$ improvement for for $20.5<G<21~{}\mathrm{mag}$). The parallax uncertainties, however, do not improve very significantly for any of the images, though there is a slight, $\sim 1\%$, precision increase at the faintest magnitudes. This is likely because the HST images are taken at a single epoch and the time baseline is quite close to a one year multiple of Gaia’s observation date (i.e. only offset by 20% of a year); this means that the parallax orbit is being sampled repeatedly in a small area at the same approximate time, which does not lead to significant increases in parallax precision.

For Draco, increasing the number of HST images does substantially improve the PMs for each image added; considering 10 HST images together with Gaia improves the PM precision by a median factor of $\sim 13$ for $20.5<G<21~{}\mathrm{mag}$ compared to a median factor of $\sim 7.8$ when only one HST image is used. The parallaxes also see a noticeable change in precision as a function of HST images analysed. While the time offsets from Gaia are again only $\sim 20\%$ of a year, the HST images are spread over three epochs. In the case of Fornax, we effectively have measurements of positions for each source at two different times (i.e. a single epoch of HST images plus Gaia) while the Draco data have measurements of positions at four different times; for the Draco stars then, there is smaller set of parallax and PM combinations that can explain those four position measurements per source than if there were only two position measurements. It is likely that the improvement in the parallax precision also contributes to the increase in PM precision because the degeneracy between the different types of apparent motion become less entangled. However, the time offsets of the HST images from Gaia and from each other are not optimally spaced in the parallax orbit (i.e. near the apocenters), resulting in a median parallax precision increase by a factor of $\sim 1.18$ for the faintest magnitudes ($G>20.5~{}\mathrm{mag}$).

As an additional method of validating the BP3M posterior PMs and their uncertainties, we cross-match the $\sim 2000$ unique COSMOS sources mentioned in Section 2 with the MILLIQUAS QSO catalog⁴⁴4https://heasarc.gsfc.nasa.gov/W3Browse/all/milliquas.html (Flesch, 2023). There are 52 sources in the QSO list that have COSMOS sources within $10~{}\mathrm{mas}$, so we run the HST images that contain these sources through the BP3M pipeline, with the resulting posterior PMs shown in Figure 15. As expected for extremely distant sources, the PMs are closely distributed around the origin.

Next, we use the PM uncertainties to measure the 2D distances of the PMs from stationary, and then compare these measurements to the expected distribution; the result of this process is shown in Figure 16. The expected distribution (orange histogram) for an equivalent sample size agrees quite well with the measured one (blue histogram). This is additional and external evidence that the BP3M posterior PM distributions are reasonable and trustworthy.

Like in Section 5.2, we identify sources in our COSMOS sample that are shared with an external survey; in this case, we cross-match with the COSMOS sample of the HALO7D survey (Cunningham et al., 2019a, b; McKinnon et al., 2023). The HALO7D PMs are described in Cunningham et al. (2019b); to summarize, the PMs are measured using multiple epochs of HST imaging with the absolute reference frame being defined by registration with background galaxies. These HST+HST PMs are therefore a result of traditional PM measurement techniques.

We analyse the HST images in the multi-epoch region of COSMOS that have Gaia sources in common with the HALO7D survey, finding 31 cross-matches. In Figure 17, we see that the BP3M PMs have strong agreement with their HST-only counterparts. As in the other comparisons we have presented, the distribution of 2D differences between the PM measurements agree statistically with our expectations, providing the final piece of validation for the BP3M outputs.

Using the HST images of the COSMOS field discussed in Section 2.1, we run BP3M on all images within $0.5^{\circ}$ of the field’s center. This includes all regions of COSMOS with multi-epoch HST imaging (near the center), as well as single-epoch-only HST regions. The posterior PMs we measure for the 2640 unique sources that are found in common between the HST images and Gaia are shown in Figure 18, where the sources are colored by whether or not they have Gaia-measured PMs and parallaxes. As a quick sanity check, this PM distribution visually agrees with our expectations based on the synthetic COSMOS data created using previously-measured velocity distributions for the MW thick disk and stellar halo⁵⁵5See Figure 23 in Appendix B, but note that our “visual comparison” here means looking at the range and spread of the PMs, not comparing the different colored populations in each figure.. There may be, however, more BP3M PMs clustered at the origin than originally expected (i.e. the orange, “No Gaia Prior” points), but many of these sources may truly be background stationary sources, which were not included as a component of the synthetic data.

While each source may appear in multiple HST frames, we only show the results from analysing each HST image separately, and the displayed PM is from the posterior distribution with the smallest uncertainty. Currently, it is very computationally expensive to analyse multiple HST concurrently, making the complete simultaneous analysis of COSMOS HST images intractable. Fortunately, as we determined from our dSph test data, analyses of single HST images do not show noticeable systematics in BP3M PMs as a result of the HST position uncertainty estimates, so the PMs we present for the COSMOS sample are robust. After increasing computation efficiency and HST position uncertainty estimates in planned updates to the pipeline, we will be able to present a complete COSMOS analysis in future work.

We also compare the BP3M PMs and uncertainties to the corresponding Gaia PMs in Figure 19, where we again see excellent agreement between the two samples (top panels). Because the sources in this figure include the same sources shown in Figure 2, we can also compare the sparse-field PM performance between the GaiaHub and BP3M pipelines. In general, BP3M returns much closer agreement with the Gaia PMs, as is expected because it uses those measurements as prior constraints. In addition, BP3M is able to measure consistent PMs for all 2640 Gaia sources in the COSMOS images, which is a large increase from the 1659 PMs catalog of GaiaHub. These results suggest that BP3M’s approach has successfully addressed the issues producing the outliers of Figure 19, indicating it is a valid method of extracting useful and trustworthy PMs in sparse fields.

The middle and bottom panels of Figure 19 show the improvement on the PM uncertainty that BP3M measures by combining HST and Gaia, especially for the $\sim 25\%$ of sources that had no Gaia-measured PMs; the median PM uncertainty for sources without Gaia priors is found to be $\sim 1.16~{}\mathrm{mas}~{}\mathrm{yr}^{-1}$. While many of the brightest magnitude sources – those most likely to be at or near saturation in the HST images – show no improvement on Gaia’s PMs, the faintest magnitude sources ($20.25<G<20.75~{}\mathrm{mag}$) have a median PM improvement factor of 2.6. As mentioned for the small-source-count case in Figure 9, the less extreme improvements in PM we see in sparse fields are likely the result of larger uncertainty in the transformation parameters, which propagates to the PM uncertainties of all sources in an image.

The BP3M transformation parameters could likely be pinned down much more precisely by identifying stationary sources in each image. For example, by finding the centroids in each HST image that are clearly associated with quasars and background galaxies, then updating the priors in parallax and PM to a narrow distribution around $\vec{0}$, the transformation solution can be improved by these anchor points. This would reduce the posterior width of the transformation parameter distributions for each image, especially for particularly sparse images, which would also reduce the posterior widths in the PMs for all of the sources in an HST image. Exploring the impact that known stationary targets can have on the BP3M PMs will be the focus of future work.

The median PM uncertainties in the $20<G<21.5~{}\mathrm{mag}$ range are summarized for the COSMOS sample in Table 2. First, the BP3M sample that has Gaia parallax and PM priors has a significant increase in precision compared to Gaia alone. Second, the Gaia COSMOS sample ($20<G<20.75~{}\mathrm{mag}$) has a similar median PM uncertainty to the BP3M sample without Gaia priors ($20.75<G<21.5~{}\mathrm{mag}$), but the latter contains significantly more stars. Effectively, BP3M substantially increases the number of stars in a single field with Gaia-quality PMs at the faintest magnitudes.

A typical MW stellar halo isochrone (i.e. ${\mathrm{age}}=12~{}\mathrm{Gyr}$ and $[\mathrm{Fe}/\mathrm{H}]=-1.2~{}\mathrm{dex}$, see Appendix B) has a red giant branch (RGB) base around an absolute magnitude of $M_{G}\approx 3~{}\mathrm{mag}$. For stars at the base of the RGB with apparent magnitudes in the $20<G<20.75~{}\mathrm{mag}$ range, this translates to a distance range of $\sim 25-35~{}\mathrm{kpc}$. Using the median PM uncertainties in Table 2, this would imply an uncertainty in velocity on the sky in the range of $\sim 123-173~{}\mathrm{km}~{}\mathrm{s}^{-1}$ for the Gaia-alone sample versus a $\sim 52-74~{}\mathrm{km}~{}\mathrm{s}^{-1}$ range for the BP3M sample with Gaia priors. For the BP3M sample without Gaia priors, the apparent magnitude range of $20.75<G<21.5~{}\mathrm{mag}$ translates to a distance range of $\sim 35-50~{}\mathrm{kpc}$, and results in a sky velocity uncertainty range of $\sim 190-268~{}\mathrm{km}~{}\mathrm{s}^{-1}$ when using the median PM uncertainty for this subsample.

Like with the dwarf spheroidal data, many of these COSMOS sources appear in multiple HST images. As a result, the BP3M PM uncertainties presented in this section are significantly larger than the possible PM uncertainty that could be achieved by analysing the images concurrently, which will be the focus on future work. This may also likely yield improvements on the parallax precision, especially for sources in the multi-epoch HST regions of the field. Finally, there are many HST sources that are too faint to have Gaia counterparts, but they nonetheless appear in multiple HST frames. For these sources, future versions of BP3M can use the HST-to-Gaia transformation parameters for each image to identify the positions of these sources in a global reference frame as a function of time, and thus determine their best fit PMs and parallaxes.

In tests of a handful of COSMOS images where we analyse all overlapping HST concurrently, we measure a median PM uncertainty of $0.064~{}\mathrm{mas}~{}\mathrm{yr}^{-1}$ for $20<G<20.75$ mag (sources with Gaia-measured PMs) and a median PM uncertainty of $0.25~{}\mathrm{mas}~{}\mathrm{yr}^{-1}$ for $20.75<G<21.5$ mag (sources with Gaia-measured PMs). In future work, we plan to use the extreme precision PMs to investigate substructure in the MW stellar halo at large distances ($D>40$ kpc). Other future applications including measuring internal kinematics of perturbed stellar streams – whose stars are often in sparse fields that are analogous to COSMOS, and where the 3D velocity error budget is dominated by PM terms (e.g., GD-1 Price-Whelan & Bonaca, 2018; Malhan & Ibata, 2019) – to constrain the dark matter subhalo mass function.