# Supplement: An isolated mass gap black hole or neutron star detected with
astrometric microlensing

Casey Y. Lam Casey Y. Lam<EMAIL_ADDRESS>University of California,
Berkeley, Department of Astronomy, Berkeley, CA 94720  Jessica R. Lu
University of California, Berkeley, Department of Astronomy, Berkeley, CA
94720  Andrzej Udalski OGLE collaboration  Astronomical Observatory,
University of Warsaw, Al. Ujazdowskie 4,00-478 Warszawa, Poland  Ian Bond  MOA
collaboration  School of Mathematical and Computational Sciences, Massey
University, Private Bag 102-904 North Shore Mail Centre, Auckland 0745, New
Zealand  David P. Bennett MOA collaboration  Code 667, NASA Goddard Space
Flight Center, Greenbelt, MD 20771, USA  Department of Astronomy, University
of Maryland, College Park, MD 20742, USA  Jan Skowron OGLE collaboration
Astronomical Observatory, University of Warsaw, Al. Ujazdowskie 4,00-478
Warszawa, Poland  Przemek Mróz OGLE collaboration  Astronomical Observatory,
University of Warsaw, Al. Ujazdowskie 4,00-478 Warszawa, Poland  Radek Poleski
OGLE collaboration  Astronomical Observatory, University of Warsaw, Al.
Ujazdowskie 4,00-478 Warszawa, Poland  Takahiro Sumi  MOA collaboration
Department of Earth and Space Science, Graduate School of Science, Osaka
University, Toyonaka, Osaka 560-0043, Japan  Michał K. Szymański OGLE
collaboration  Astronomical Observatory, University of Warsaw, Al. Ujazdowskie
4,00-478 Warszawa, Poland  Szymon Kozłowski OGLE collaboration  Astronomical
Observatory, University of Warsaw, Al. Ujazdowskie 4,00-478 Warszawa, Poland
Paweł Pietrukowicz OGLE collaboration  Astronomical Observatory, University of
Warsaw, Al. Ujazdowskie 4,00-478 Warszawa, Poland  Igor Soszyński OGLE
collaboration  Astronomical Observatory, University of Warsaw, Al. Ujazdowskie
4,00-478 Warszawa, Poland  Krzysztof Ulaczyk OGLE collaboration  Astronomical
Observatory, University of Warsaw, Al. Ujazdowskie 4,00-478 Warszawa, Poland
Department of Physics, University of Warwick, Gibbet Hill Road,Coventry, CV4
7AL, UK  Łukasz Wyrzykowski OGLE collaboration  Astronomical Observatory,
University of Warsaw, Al. Ujazdowskie 4,00-478 Warszawa, Poland  Shota
Miyazaki  MOA collaboration  Department of Earth and Space Science, Graduate
School of Science, Osaka University, Toyonaka, Osaka 560-0043, Japan  Daisuke
Suzuki  MOA collaboration  Department of Earth and Space Science, Graduate
School of Science, Osaka University, Toyonaka, Osaka 560-0043, Japan  Naoki
Koshimoto  MOA collaboration  Code 667, NASA Goddard Space Flight Center,
Greenbelt, MD 20771, USA  Department of Astronomy, University of Maryland,
College Park, MD 20742, USA  Department of Astronomy, Graduate School of
Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033,
Japan  Nicholas J. Rattenbury MOA collaboration  Department of Physics,
University of Auckland, Private Bag 92019, Auckland, New Zealand  Matthew W.
Hosek, Jr. University of California, Los Angeles, Department of Astronomy, Los
Angeles, CA 90095, USA  Fumio Abe  MOA collaboration  Institute for Space-
Earth Environmental Research, Nagoya University, Nagoya 464-8601, Japan
Richard Barry  MOA collaboration  Code 667, NASA Goddard Space Flight Center,
Greenbelt, MD 20771, USA  Aparna Bhattacharya  MOA collaboration  Code 667,
NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA  Department of
Astronomy, University of Maryland, College Park, MD 20742, USA  Akihiko Fukui
MOA collaboration  Department of Earth and Planetary Science, Graduate School
of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033,
Japan  Instituto de Astrofísica de Canarias, Vía Láctea s/n, E-38205 La
Laguna, Tenerife, Spain  Hirosane Fujii  MOA collaboration  Institute for
Space-Earth Environmental Research, Nagoya University, Nagoya 464-8601, Japan
Yuki Hirao  MOA collaboration  Department of Earth and Space Science, Graduate
School of Science, Osaka University, Toyonaka, Osaka 560-0043, Japan
Yoshitaka Itow  MOA collaboration  Institute for Space-Earth Environmental
Research, Nagoya University, Nagoya 464-8601, Japan  Rintaro Kirikawa  MOA
collaboration  Department of Earth and Space Science, Graduate School of
Science, Osaka University, Toyonaka, Osaka 560-0043, Japan  Iona Kondo  MOA
collaboration  Department of Earth and Space Science, Graduate School of
Science, Osaka University, Toyonaka, Osaka 560-0043, Japan  Yutaka Matsubara
MOA collaboration  Institute for Space-Earth Environmental Research, Nagoya
University, Nagoya 464-8601, Japan  Sho Matsumoto  MOA collaboration
Department of Earth and Space Science, Graduate School of Science, Osaka
University, Toyonaka, Osaka 560-0043, Japan  Yasushi Muraki  MOA collaboration
Institute for Space-Earth Environmental Research, Nagoya University, Nagoya
464-8601, Japan  Greg Olmschenk  MOA collaboration  Code 667, NASA Goddard
Space Flight Center, Greenbelt, MD 20771, USA  Clément Ranc MOA collaboration
Zentrum für Astronomie der Universität Heidelberg, Astronomiches Rechen-
Institut, Mönchhofstr. 12-14, 69120 Heidelberg, Germany  Arisa Okamura
Department of Earth and Space Science, Graduate School of Science, Osaka
University, Toyonaka, Osaka 560-0043, Japan  Yuki Satoh  Department of Earth
and Space Science, Graduate School of Science, Osaka University, Toyonaka,
Osaka 560-0043, Japan  Stela Ishitani Silva  MOA collaboration  Department of
Physics, The Catholic University of America, Washington, DC 20064, USA  Code
667, NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA  Taiga Toda
MOA collaboration  Department of Earth and Space Science, Graduate School of
Science, Osaka University, Toyonaka, Osaka 560-0043, Japan  Paul J. Tristram
MOA collaboration  University of Canterbury Mt. John Observatory, P.O. Box 56,
Lake Tekapo 8770, New Zealand  Aikaterini Vandorou  MOA collaboration  Code
667, NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA  Department of
Astronomy, University of Maryland, College Park, MD 20742, USA  Hibiki Yama
MOA collaboration  Department of Earth and Space Science, Graduate School of
Science, Osaka University, Toyonaka, Osaka 560-0043, Japan  Natasha S. Abrams
University of California, Berkeley, Department of Astronomy, Berkeley, CA
94720  Shrihan Agarwal University of California, Berkeley, Department of
Astronomy, Berkeley, CA 94720  Sam Rose University of California, Berkeley,
Department of Astronomy, Berkeley, CA 94720  Sean K. Terry University of
California, Berkeley, Department of Astronomy, Berkeley, CA 94720

###### Abstract

This Supplement provides supporting material for Lam et al. (2022). We briefly
summarize past gravitational microlensing searches for black holes and present
details of the observations, analysis, and modeling of five black hole (BH)
candidates observed with both ground-based photometric microlensing surveys
and Hubble Space Telescope astrometry and photometry. We present detailed
results for four of the five candidates that show no or low probability for
the lens to be a black hole. In these cases, the lens masses are $<2M_{\odot}$
and two of the four are likely white dwarfs or neutron stars. We also present
detailed methods for comparing the full sample of five candidates to
theoretical expectations of the number of BHs in the Milky Way ($\sim
10^{8}$).

††software: Galaxia (Sharma et al., 2011), astropy (Astropy Collaboration et
al., 2013, 2018), Matplotlib (Hunter, 2007), NumPy (van der Walt et al.,
2011), SciPy (Virtanen et al., 2019), SPISEA (Hosek et al., 2020), PopSyCLE
(Lam et al., 2020), dynesty (Speagle, 2020), PyMultiNest (Buchner et al.,
2014; Feroz et al., 2009), dustmaps (Green, 2018), hst1pass, xym2mat, xym2bar
(Anderson & King, 2006), ks2 (Anderson et al., 2008; Bellini et al., 2018),
celerite (Foreman-Mackey et al., 2017)

##  1 Supplemental Introduction

This paper is a supplement to Lam et al. (2022), in which we present results
from a search for stellar mass black holes in the Milky Way using
gravitational microlensing. Lam et al. (2022) and this supplement describe our
analysis of five Milky Way black hole (BH) candidates, which constitute a
sufficient sample to place preliminary constraints on the number of isolated
black holes in the Milky Way.

This supplement is organized as follows. In §2, we present a detailed review
of past microlensing searches for BHs. In §3, the datasets and reduction
processes for our BH search are described, and in §4 the photometric and
astrometric analysis of the HST data is explained. In §5, the microlensing
modeling and fitting procedure are described in detail, while §6 describes how
to combine the high-resolution imaging and results of the microlensing
modeling to constrain the lens’ luminosity. §7 presents the results of the
modeling for the five candidates. Note that results for the single mass-gap
black hole or neutron star candidate are presented in the main paper Lam et
al. (2022). Details on constraints on the Milky Way BH population and future
BH microlensing searches are discussed in §8.

##  2 Past Microlensing Searches for Black Holes

The advent of ground-based microlensing surveys provided a new avenue to
search for isolated black holes (BHs). Notable BH candidates identified with
photometric microlensing include MACHO-96-BLG-5, MACHO-98-BLG-6 (Bennett et
al., 2002), and MACHO-99-BLG-22/OGLE-1999-BUL-32 (Mao et al., 2002).
Poindexter et al. (2005) found that MACHO-99-BLG-22 is likely a BH,
MACHO-96-BLG-5 is possibly a BH, and MACHO-98-BLG-6 is most likely not a BH,
with BH lens probabilities of 78%, 37%, and 2%, respectively. On the other
hand, Abdurrahman et al. (2021) found that MACHO-96-BLG-5 and MACHO-98-BLG-6
are still good BH candidates, ruling out non-BH lenses for source-lens
relative proper motions larger than 2.5 mas/yr. However, mass estimates for
these lensing events cannot be made without invoking a Galactic model. These
candidates only had photometric microlensing observations, which alone cannot
constrain the mass of the lens, unless rare higher-order effects such as
finite-source effects are detected.

As mentioned briefly at the end of Paczynski (1986), microlensing also has an
astrometric signature, in which the centroid of the image is displaced from
the source’s true position (Hog et al., 1995; Miyamoto & Yoshii, 1995; Walker,
1995). In contrast to the now-routine measurements of photometric
microlensing, detections of astrometric microlensing are still at the
forefront of our technical capabilities. Typical astrometric shifts toward the
Bulge are $O(0.01-1)$ mas, and few existing facilities are currently capable
of the astrometric precision to perform this measurement. Only a handful of
astrometric measurements of the gravitational deflection of light have ever
been made, all for nearby ($<10$ pc) lenses that were astrometrically
anticipated (Eddington, 1919; Sahu et al., 2017; Zurlo et al., 2018). However,
a combination of photometric and astrometric microlensing together can
determine the mass of the lensing object, making detection of astrometric
microlensing important for BH searches.

To date, there have been two endeavors to measure lens masses by combining
photometric and astrometric microlensing. Lu et al. (2016) attempted a
measurement with Keck laser guide star adaptive optics (LGS AO), but no
detections of astrometric microlensing were made. Kains et al. (2017) reported
a detection of astrometric microlensing made with the Hubble Space Telescope
(HST), but the signal was very weak and no lens masses were well constrained.
The currently operating Gaia mission is also anticipated to make measurements
of astrometric microlensing, which can be searched for once per-epoch
astrometry is released (McGill et al. (2020) and references therein). A
handful of these deflections should be due to BHs (Rybicki et al., 2018).

##  3 Observations in Detail

###  3.1 Targets and Selection Criteria

Table 1: Target Summary

Short Name | OGLE Alert Name | MOA Alert Name | RA (J2000.0) | Dec (J2000.0)
---|---|---|---|---
MB09260 | – | MOA-2009-BLG-260 | 17:58:28.561 | -26:50:20.88
MB10364 | – | MOA-2010-BLG-364 | 17:57:05.401 | -34:27:05.01
OB110037 | OGLE-2011-BLG-0037 | MOA-2011-BLG-039 | 17:55:55.83 | -30:33:39.7
OB110310 | OGLE-2011-BLG-0310 | MOA-2011-BLG-332 | 17:51:25.39 | -30:24:35.0
OB110462 | OGLE-2011-BLG-0462 | MOA-2011-BLG-191 | 17:51:40.19 |

Five BH candidate microlensing events111Three additional other targets were
initially observed with HST, but dropped from the target list after a year
(Sahu, 2012). were selected from the OGLE Early Warning
System222http://ogle.astrouw.edu.pl/ogle4/ews/ews.html (Udalski et al., 2015)
and MOA Alerts333https://www.massey.ac.nz/$\sim$iabond/moa/alerts/ to be
imaged with HST. These targets were selected to have long ($t_{E}>200$ days)
duration, no light contribution from the lens, and high magnifications to
allow detection of parallax signals, making them good isolated BH candidates
(Sahu, 2009). Preliminary results by Sahu et al. (2017) reported that all five
candidates were low-mass ($<0.5M_{\odot}$) stars, with no BH detections in the
sample. We chose to re-analyze these targets in order to use both BH
detections and non-detections in the sample to constrain the total number of
BHs in the Milky Way.

Three of the targets, OGLE-2011-BLG-0037/MOA-2011-BLG-039,
OGLE-2011-0310/MOA-2011-BLG-332, and OGLE-2011-BLG-0462/MOA-2011-BLG-191
(hereafter OB110037, OB110310, and OB110462), were alerted by both OGLE and
MOA. The other two targets, MOA-2009-BLG-260 and MOA-2010-BLG-364 (hereafter
MB09260 and MB10364), were only alerted by MOA. Table 3.1 lists their
coordinates. Figures 3 \- 8 show the lightcurves of the targets.

###  3.2 MOA

The MOA-II survey is carried out with a 1.8-m telescope at Mt. John University
Observatory in New Zealand (Hearnshaw et al., 2006; Sumi, 2008). The seeing at
the site ranges from $\sim 1.8^{\prime\prime}-3.5^{\prime\prime}$, with the
median seeing being $\sim 2.5^{\prime\prime}$. The telescope has a 2.2 deg2
FOV, with a 10-chip CCD camera with a plate scale $0.57^{\prime\prime}$/pixel.
The main observations are taken using the MOA-Red (630-1000 nm) filter (Bond
et al., 2001).

MOA data for MB09260 and MB10364 were reduced as described in Bond et al.
(2017).444OB110037, OB110310, and OB110462 also have MOA lightcurves. For
simplicity we only present the OGLE lightcurve fits for those events, since
the seeing at OGLE is better. The MOA lightcurves are photometrically
calibrated to the OGLE-III I-band.

###  3.3 OGLE

The OGLE-IV survey is carried out at the 1.3-m Warsaw telescope at Las
Campanas Observatory in Chile (Udalski et al., 2015). The seeing at the site
ranges from $\sim 1.0^{\prime\prime}-2.0^{\prime\prime}$, with the median
seeing being $\sim 1.3^{\prime\prime}.$ The telescope has a 1.4 deg2 FOV, with
a 32-chip CCD camera with a plate scale of $0.26^{\prime\prime}$/pixel. The
main observations are taken using the OGLE-I filter, which is similar to
Cousins I-band. The data was reduced using the Difference Image Analysis
technique as implemented by Wozniak (2000). In addition, we rescaled the
photometric uncertainties of OB110462 according to the method described in
Skowron et al. (2016).

OGLE data is only available for OB110037, OB110310, and OB110462. The
magnification of MB09260 was not observed by OGLE as it occurred during the
OGLE-III to OGLE-IV upgrade. MB10364 is located in a gap in the detectors of
the OGLE camera.

###  3.4 HST

Table 2: HST WFC3-UVIS Observations

Note. — Asterisk (∗) denotes observations excluded from analysis.

Table 3: HST WFC3-UVIS Observations

Note. — Asterisk (∗) denotes observations excluded from analysis.

HST observations come from a multi-year campaign following up these five
targets (GO-11707, GO-12322, GO-12670, GO-12986, GO-13458, GO-14783; PI: K. C.
Sahu). Observations were taken with the UVIS channel on the Wide Field Camera
3 (WFC3) in two different wide-band filters, F606W ($V$-band) and F814W
($I$-band). Table 3.4 summarizes the HST observations.

The WFC3 UVIS channel is composed of two 2k$\times$4k CCDs and has a
$162^{\prime\prime}\times 162^{\prime\prime}$ field of view with a plate scale
$0.04^{\prime\prime}$/pixel. WFC3 UVIS supports sub-arraying, in which only a
portion of the entire detector is read out, which can reduce data volume or
exposure time and increase observational efficiency. All observations prior to
2011-07-22 were taken with the UVIS1-2K4-SUB subarray mode. Beginning HST
Cycle 18, more subarray sizes were made available, and observations after
2011-07-22 were taken with the UVIS2-2K2C-SUB subarray mode, a 2k$\times$2k
subarray.

Additional observations of OB110462 commenced in Cycle 29 (GO-16760; PI: C.
Lam). These were taken in as similar a configuration as possible to the later
epochs of the archival program, using WFC3 UVIS in UVIS2-2K2C-SUB subarraying
mode, with observations in F606W and F814W filters. The first set of
observations from this program was taken October 2021 and is presented here;
an additional set of observations is anticipated to be taken Fall 2022 (Lam &
Lu, 2021).

Figure 1:  Central $3^{\prime\prime}\times 3^{\prime\prime}$ of HST WFC3-UVIS
F814W combined images of the observed fields, centered on the target
(circled). These images are of the target at or near baseline, i.e.
unmagnified. The color stretch is logarithmic. Note that the color scale is
not the same across panels. Figure 2:  CMDs for each field. The target at
baseline magnitude and color is marked as a red star.

###  3.5 Gaia

Gaia is an all-sky scanning astrometric space mission (Gaia Collaboration et
al., 2016). All of the targets, with the exception of OB110462, are found in
Gaia Early Data Release 3 (EDR3, Table 3.5). Gaia EDR3 covers the period from
25 July 2014 to 28 May 2017 (Gaia Collaboration et al., 2021). MB10364,
OB110037, and OB110310 have proper motions and parallaxes, while MB09260 does
not. OB110462 was not in Gaia as it is too faint. We note that there is a Gaia
source located $\sim 0.35^{\prime\prime}$ away from OB110462 (Gaia EDR3 Source
ID 4056442477683080960), which coincides with the bright star directly west of
the target seen in the HST images (Figure 1).

Table 4: Gaia EDR3 Values

Note. — Magnitude uncertainties are estimated from the Gaia reported flux
errors. Zero-point (ZP) correction comes from Lindegren et al. (2021a), ∗
denotes values that are extrapolations. For full descriptions we refer the
reader to Gaia EDR3 documentation (van Leeuwen et al., 2021), Section 13.1.1
gaia_source. OB110462 is not in Gaia.

Figure 3: _Top panel_ : Detrended MB09260 lightcurve, as seen by MOA and HST
F814W and F606W. The maximum likelihood model (MLE, described in §7) is
plotted over the data. _Second from top panel_ : The residuals to the MLE
model. The Gaussian Process (GP) model is plotted on top of the residual. We
emphasize the residual is not independently fit by the GP, but is
simultaneously fit with the model parameters; this is purely to visualize the
data (also see Golovich et al. (2022)). See §5 for more details about the
fitting procedure. _Second from bottom panel_ : Same as top panel, but zoomed
into the three most magnified years. _Bottom panel_ : Same as second from top
panel, but zoomed into the three most magnified years. Figure 4: Same as
Figure 3, but for MB10364. Note we do not include GP in the MB10364 fit.
Figure 5: Same as Figure 3, but for OB110037. Instead of MOA data, the blue
data are OGLE data. Figure 6: Same as Figure 5, but for OB110310. Figure 7:
Same as Figure 5, but for OB110462. This maximum likelihood model was
calculated using the default weighted likelihood described in §5.3. Figure 8:
Same as Figure 5, but for OB110462. This maximum likelihood model was
calculated using the equal weight likelihood described in §5.3.

##  4 HST Data Analysis

###  4.1 Reduction

The HST archival data was accessed from the Mikulski Archive for Space
Telescopes555https://archive.stsci.edu/hst/ (MAST) in June 2021. The Cycle 29
data was accessed October 2021. For the following analysis, we employed the
calibrated, flat-fielded, individual exposures corrected for charge transfer
efficiency (CTE, HST files with suffix _flc666See Gennaro (2018) for a full
description of the different file name suffixes.). The archival data was
processed with Version 3.6.0 (Dec-31-2020) of the calwf3 pipeline, using
Version 2.0 of the CTE correction algorithm.

CTE can alter astrometry at the milliarcsecond level, hence it is important to
use _flc files. However, even the _flc files do not necessarily fix all
problems associated with CTE (Kuhn & Anderson, 2021). Exploration of other
methods of CTE correction will be explored in future work. At the present we
mitigate CTE effects via other methods (§4.2.2) and validate our astrometry to
ensure it is not distorted by CTE.

Images were converted into calibrated star lists via the following steps.

  1. 1.

_Star list extraction from individual frames._ Star lists were extracted from
the individual _flc exposures by modeling the PSFs of sources with hst1pass,
an updated version of the software described in Anderson & King (2006).
Empirical filter-dependent PSF models as described in Anderson (2016) and
geometric distortion solutions as described in Bellini et al. (2011) were used
when performing source extraction with hst1pass.

  2. 2.

_Combined star list for one epoch._ Within a single epoch and filter, multiple
star lists were aligned to a common coordinate system in an iterative manner
using xym2mat and xym2bar (Anderson & King, 2006), which include the
distortion solution for the WFC3 camera and filters (Bellini et al., 2011) to
produce a single matched star list.

  3. 3.

_Photometric calibration._ Lastly, a zero-point is applied to the star lists
to convert from instrumental to Vega magnitudes. Star lists were calibrated
against photometrically calibrated star lists on the Hubble Legacy Archive,
Data Release 10 (HLA DR10). A magnitude offset is applied later during the
astrometric alignment (§4.3) to obtain more precise relative photometry.

Note that data taken in F606W and F814W filters are treated as independent
measurements. That is, observations taken on the same date are treated as
distinct epochs, and are not combined into a single star list, as the
importance of filter dependence in astrometry is not well established. See
§4.2.5 and Appendix H for further details.

Certain epochs were excluded from the analysis; these are marked with an
asterisk in Table 3.4. The reason for their exclusion is detailed as follows.

  * •

_Epochs with only a single frame._ Observations with only a single frame per
filter cannot produce any useful photometric or astrometric constraints using
hst1pass. This is the case for the MB09260 2009-10-01 F814W and 2009-10-19,
2010-03-22, and 2010-06-14 F606W epochs. In the MB09260 2010-10-20 F606W
epoch, a cosmic ray in one of the exposures interfered with the extraction of
the target, effectively leaving only a single usable frame.

  * •

_Multiple exposure times._ Although in principle mixing multiple exposure
times in a single epoch is possible, in practice most of the data was obtained
with several long exposures and only a single frame with a shorter exposure.
Rather than analyze the few short frames with different detection thresholds,
PSF reference stars, and astrometric reference stars, which can lead to
systematic errors, we choose to only use frames with the same exposure times
within an epoch. For this reason, some frames from the MB09260 2009-10-01
F606W, MB10364 2010-09-13 F814W, OB110037 2011-08-15 F814W, OB110310
2013-10-21 F814W, and OB110462 2011-08-08 F606W and F814W epochs were not
used.

  * •

_Saturation of target._ No useful astrometric limits can be placed when the
target is saturated. The target is saturated in both filters in the MB10364
2011-04-13 and 2011-07-22 epochs.

  * •

_Telescope pointing issues._ The observations of OB110462 on 2017-08-11
suffered a telescope drifting issue, resulting in streaked images.

  * •

_Astrometric alignment systematics._ Although there are no standalone issues
with the observations of OB110462 on 2013-05-13, astrometric systematics are
apparent in the reference stars when this epoch is astrometrically aligned
along with the other epochs using the methodology described in §4.2. This is
due to the difference in position angle of the observations taken, as the
2013-05-13 epoch was taken at PA = 99.9 deg, while all the other epochs were
taken with a PA different by $\sim 180$ deg, with PA = 255.2 - 276.1 deg.
Thus, the 2013-05-13 epoch is left out of the analysis.

The other targets (MB09260, OB110037, OB110310) with $\sim 180$ deg
differences in PA across observations do not suffer this same problem, as
there are multiple observations at each PA. This allows the systematics due to
the $\sim 180$ deg PA flip to be calibrated out during the astrometric
alignment.

###  4.2 HST Astrometric Analysis

The positional measurements extracted from the different epochs of HST data
(§4.1) must be transformed into a common reference frame in order to derive
the motion of the target. This is an iterative process, with multiple “passes”
at refining the reference frame, allowing for the best relative astrometry
possible to be extracted. We follow a similar procedure as described in §4.2
of Lu et al. (2016). The final positions and magnitudes of the targets
resulting from the analysis in this section is presented in Table 4.3. The
photometry is shown in Figures 3 \- 8 and the astrometry is shown in Figures
14 \- 19.

####  4.2.1 Alignment procedure

Following standard image processing techniques, a 2-D polynomial
transformation of the form

$\displaystyle x^{\prime}$
$\displaystyle=a_{0}+a_{1}x+a_{2}y+a_{3}x^{2}+a_{4}xy+a_{5}y^{2}+...$ (1)
$\displaystyle y^{\prime}$
$\displaystyle=b_{0}+b_{1}x+b_{2}y+b_{3}x^{2}+b_{4}xy+b_{5}y^{2}+...$ (2)

is applied to the images in order to match them to a reference image. A first
order 2-D polynomial transformation ($x^{\prime}=a_{0}+a_{1}x+a_{2}y$;
$y^{\prime}=b_{0}+b_{1}x+b_{2}y$) is an affine transformation777An affine
transformation maps points to points, lines to lines, planes to planes, and so
on. Affine tranformations preserve collinearity and ratios of distances.
Parallel lines also remain parallel after an affine transformation., which can
be used to model translation, rotation, scaling, and shearing introduced by
the camera. A higher-order polynomial can correct for additional distortions,
but going beyond second order generally does not improve results as the number
of free parameters quickly increases and results in overfitting.

In the first pass, the HST images are aligned to the absolute reference frame
of Gaia with a first order 2-D polynomial transformation to roughly establish
the transformation. The Gaia EDR3 catalog is matched to the HST catalog using
the pattern matching algorithm of Groth (1986).

In subsequent passes, the HST images are aligned to themselves, using a 2-D
polynomial transformation going up to second order. It is empirically
determined that 3-4 passes gives optimal results. In each successive pass, the
HST images are aligned to the reference frame derived in the previous pass,
which continually refines the reference frame and derived proper motions.

To calculate the optimal transformation, a set of reference stars
$r^{R}_{ref}(t_{0})$ are selected from the stars in reference frame $R(t_{0})$
observed at time $t_{0}$. The reference stars $r^{R}_{ref}(t_{0})$ are matched
to corresponding stars $u^{U}_{ref}(t)$ in the untransformed $U(t)$ frame
observed at time $t$. The transformation $T:U(t)\rightarrow R(t)$ is found by
least-squares minimization of the $x$ and $y$ position residuals from the
alignment

$\displaystyle x_{res}$
$\displaystyle=\sum_{i}w_{x,i}(x^{R}_{r,i}(t)-T(x^{U}_{u,i}(t)))^{2}$ (3)
$\displaystyle y_{res}$
$\displaystyle=\sum_{i}w_{y,i}(y^{R}_{r,i}(t)-T(y^{U}_{u,i}(t)))^{2}$ (4)

where $w_{(x,y),i}$ is the weight for the $i$-th reference star,
$(x^{R}(t)_{r,i},y^{R}(t)_{r,i})$ and $(x^{U}_{u,i}(t),y^{U}_{u,i}(t))$ are
the positions of the stars in the reference $R(t)$ and untransformed $U(t)$
frames at time $t$. The positions of the reference stars $r^{R}_{ref}(t)$ in
the reference frame $R$ at time $t$ are propagated from time $t_{0}$ using the
proper motions

$\displaystyle x^{R}_{r,i}(t)$
$\displaystyle=x^{R}_{r,i}(t_{0})+v^{R}_{x,r,i}(t-t_{0})$ (5) $\displaystyle
y^{R}_{r,i}(t)$ $\displaystyle=y^{R}_{r,i}(t_{0})+v^{R}_{y,r,i}(t-t_{0})$ (6)

where $v^{R}_{(x,y),r,i}$ are the proper motions of the $i$-th reference stars
in reference frame $R$. After the transformation $T$ is derived, it is applied
to all the stars $u^{U}$ in the $U$ frame to obtain a transformed star list
$u^{R}(t)=T(u^{U}(t))$, where the stars $u^{R}(t)$ are now in the frame $R$.
This yields starlists for all $N$ epochs
$u_{1}^{R}(t_{1}),\,...\,,u_{N}^{R}(t_{N})$, where the positions of all the
stars are now in the same reference frame $R$. For each star $j$, a proper
motion is derived by finding the best-fit straight line via non-linear least
squares through the $n\leq N$ observations.888Some stars are not detected in
all epochs, which is why it is possible to have $n<N$ observations.

Lu et al. (2016) examined several different weighting schemes and showed that
the resulting astrometry is not affected. We choose to use weights
$w_{(x,y),i}=1/\sigma_{(x,y),i}$, where $\sigma_{(x,y),i}$ are the positional
uncertainties of the stars in the untransformed frame. For the positional
uncertainties, instead of using the root mean square values $\sigma_{RMS}$
directly returned by hst1pass, we follow Hosek et al. (2015) and use the error
on the mean $\sigma_{RMS}/\sqrt{N}$ where $N$ is the number of frames the
source is detected in, with an additional empirical additive error. The
uncertainties, as well as the procedure used to determine them, are detailed
in Appendix A. The additive error term dominates over the
$\sigma_{RMS}/\sqrt{N}$ term for bright stars, which makes the positional
errors more uniform across epochs and magnitude as would be expected for
systematic errors.

As we are interested in the astrometry of the target, the target itself is not
used to establish the transformation into a common reference frame (§4.2.2)
nor to judge the quality of the final transformation (Figure 9).

####  4.2.2 Reference star selection

Reference stars are stars assumed to have linear proper motions, which are
used to derive the reference frame transformation and the motions of the other
stars. The selection of reference stars depends on multiple considerations,
such as the stellar density, amount of geometric distortion, instrumental
systematics, and number and brightness of targets of interest, to name a few.
The goal is to balance having enough stars to establish the reference frame,
while excluding stars which would produce a non-stable reference frame. The
criteria for reference star selection for each target is summarized in Table
4.2.3. We choose reference stars with brightness similar to the target,
relatively large radial separations from our target of interest, and exclude
likely foreground stars. The target itself is also excluded from being a
reference star. We detail the reasoning for these choices below.

_Brightness range_ : Due to the nature of CTE, there are strong magnitude-
dependent astrometric residuals, even when using the latest CTE-corrected _flc
images. However, this is not unexpected (Kuhn & Anderson, 2021). For this
reason, stars are chosen to be in a brightness range similar to the one
spanned by the target as narrow as possible. For OB110310 and OB110462, which
are relatively faint and where there are many stars of similar brightness, all
stars falling within $\pm 0.1$ mag of the target’s brightest and faintest in
the HST data are used. MB09260, on the other hand, is brighter, with less
stars of similar brightness, so the range is larger, with all stars falling
within $\pm 1$ mag of the target’s brightest and faintest in the HST data are
used. MB10364 is so bright that many stars of comparable brightness are
saturated in the longer exposures. Because of this, only bright (F814W
$<18.0$, F606W $<19.2$) and unsaturated stars were selected; the bright limit
on the magnitude range differs between epochs because of the different
exposure times.

_Spatial separation_ : Only reference stars within $30^{\prime\prime}$, or
$20^{\prime\prime}$ for the denser field around OB110462, are used as
reference stars. This minimizes the impact of geometric distortion residuals
and spatially dependent PSF variations.

_Foreground stars_ : A key assumption in the astrometric alignment process is
that reference stars have linear proper motions, and parallax effects are
ignored. For a typical bulge star 8 kpc from Earth, this is a reasonable
assumption as the parallax will be 1/8000 arcsec = 0.125 mas, below our
achievable astrometric precision. However, for nearby stars, ignoring parallax
is an issue when trying to derive an accurate transformation. As all the
target fields are toward the highly extincted Galactic bulge, bright blue
stars as identified on a CMD (Figure 2) are likely to be nearby and have a
non-negligible parallax, and are excluded from the set of reference stars. The
color-magnitude exclusion criteria are listed in the last column of Table
4.2.3. For MB10364, no bright blue stars were removed as reference stars, as
all the observations came from within 6 weeks of the same time of year. Hence,
any type of yearly parallax signal would be negligible within this time span.

_Number of detections_ : We require reference stars to be detected in most, if
not all, epochs. If there are $N_{\textrm{ep}}$ total epochs, we require
reference stars to be detected in $N_{\textrm{ep,detect}}=N_{\textrm{ep}}-2$
epochs.

Lastly, as the motion of the target is the quantity we are interested in, we
do not use it as a reference star.

####  4.2.3 Derived stellar proper motions

To evaluate the goodness of the fits of the derived stellar proper motions, we
consider the $\chi^{2}$ distributions of the position residuals

$\displaystyle\chi^{2}_{x}$
$\displaystyle=\sum_{t}\Big{(}\frac{x_{t}-x_{t,fit}}{\sigma_{x_{t}}}\Big{)}^{2}$
(7) $\displaystyle\chi^{2}_{y}$
$\displaystyle=\sum_{t}\Big{(}\frac{y_{t}-y_{t,fit}}{\sigma_{y_{t}}}\Big{)}^{2}$
(8)

where $x\equiv$ RA, $y\equiv$ Dec, $(x,y)_{t}$ are the positions in the data,
$(x,y)_{t,fit}$ are the positions as derived from the linear motion fit, and
$\sigma_{(x,y),t}$ are the positional uncertainties at time $t$. The
distributions of $\chi^{2}_{x}$ and $\chi^{2}_{y}$ for the reference stars
detected in all epochs $N_{ep}$ are shown in Figure 9, with the expected
$\chi^{2}$ distribution overplotted on top. The distributions for the
positions in F814W and F606W are shown separately. The expected residual
distribution has $N_{ep,detect}-2$ degrees of freedom, as there are two free
parameters in the linear motion fit (initial position and proper motion). Note
that unlike Gaia, this linear model fit does not include parallax. Parallax is
only included when modeling the microlensing event.

Figure 9: Histogram of $\chi^{2}$ residual values to the linear fits with no
parallax for the reference stars of each target (Table 4.2.3). In each panel,
the left column shows the distributions for reference stars in F606W, while
the right column shows F814W. The number of reference stars is listed as $N$.
In each panel, the top row shows the $\chi^{2}$ distribution of residuals of a
linear fit to positions vs. time in $x\equiv$ RA and $y\equiv$ Dec. The bottom
row shows the $\chi^{2}$ distribution of residuals of a constant fit to
magnitude vs. time. The expected $\chi^{2}$ distributions are shown in black,
with the number of degrees of freedom listed as $N_{dof}$.  Table 5: Reference
star criteria

Note. — These are the criteria for the last pass.

####  4.2.4 OB110462 bias correction

There is a bright star (“the neighbor”) $\sim$10 pixels ($\sim$0.4 arcsec)
west of OB110462. The neighbor is $\sim$3 magnitudes brighter than OB110462 at
baseline (F814W = 16.7 mag, F606W = 19.0 mag). Because of its proximity and
high contrast, the neighbor’s PSF might “leak” onto OB110462 and alter its
astrometry and photometry. We perform injection and recovery tests to
ascertain the reliability of faint source extraction near a bright source, in
order to determine whether the astrometry and photometry of OB110462 as
determined in §4.2 and §4.3 is biased by the bright neighbor.

The methodology and results of the injection and recovery analysis are
detailed in Appendix B. In summary, the positional bias is negligible when the
target is highly magnified and of similar brightness to the neighbor. However,
in epochs where the target is no longer magnified, the bright star biases the
position of the target. In F814W, where the resolution is lower, the measured
position of the target is biased toward the neighbor by $\sim 0.4$ mas along
the target-neighbor separation vector. In F606W, where the resolution is
higher, the bias is less ($\sim 0.25$ mas) with the direction of bias more
randomly oriented. Similarly, the photometric bias is larger when the contrast
is large, with the bright neighbor causing the extracted photometry of
OB110462 to be brighter than the injected values. The effect is again more
severe in F814W than F606W because of the lower resolution.

Using the results of the injection and recovery analysis, we calculate a bias
correction to apply to OB110462 astrometry and photometry (Table B.1). The
values in the table are added to astrometry and photometry derived in §4.2 and
§4.3; the uncertainties are added in quadrature to the uncertainties in §4.2
and §4.3.

We only perform this analysis for OB110462, as it is the only faint target
near a bright companion. All the other targets are either bright with faint
companions, isolated, or both bright and isolated.

####  4.2.5 Astrometric color offset

As mentioned in §4.1 the data taken in F606W and F814W filters are treated as
independent measurements. For OB110037 and OB110462, the astrometric
measurements in F606W and F814W do not agree within the uncertainties.

For OB110037, although the 2011 and 2012 epochs show good agreement, the 2013
to 2017 epochs become increasingly discrepant as time goes on. We attribute
this difference to binarity (§7.4).

In contrast, for OB110462 the astrometry in the F606W and F814W are discrepant
in all datasets, but the difference appears to be a relatively constant offset
with time. This is true both before and after applying the bias correction in
§4.2.4. Because the nature of the color difference appears to be a constant
offset, we apply a constant shift to the OB110462 F606W astrometry

$\displaystyle\Delta_{x}$ $\displaystyle=\frac{\sum_{t}w_{x,t}(x_{{\rm
F814W},t}-x_{{\rm F606W},t})}{\sum_{t}w_{x,t}}$ (9) $\displaystyle\Delta_{y}$
$\displaystyle=\frac{\sum_{t}w_{y,t}(y_{{\rm F814W},t}-y_{{\rm
F606W},t})}{\sum_{t}w_{y,t}}$ (10)

where

$\displaystyle w_{x,t}$ $\displaystyle=(\sigma_{x,{\rm
F814W},t}^{2}+\sigma_{x,{\rm F606W},t}^{2})^{-1/2}$ (11) $\displaystyle
w_{y,t}$ $\displaystyle=(\sigma_{y,{\rm F814W},t}^{2}+\sigma_{y,{\rm
F606W},t}^{2})^{-1/2}.$ (12)

and $t$ indexes the observation times. Thus the modified astrometry for F606W
is

$\displaystyle x_{{\rm F606W},t}^{\prime}$ $\displaystyle=x_{{\rm
F606W},t}+\Delta_{x}$ (13) $\displaystyle y_{{\rm F606W},t}^{\prime}$
$\displaystyle=y_{{\rm F606W},t}+\Delta_{y}.$ (14)

The values of the offset are $\Delta_{x}=-0.57$ mas and $\Delta_{y}=0.39$ mas.
Note that these offsets are calculated using the bias-corrected
astrometry.999Note that even before the bias correction of §4.2.4 is applied,
this color offset is still present. In fact, it is slightly larger, with
$\Delta_{x}=-0.79$ mas and $\Delta_{y}=0.52$ mas. The bias correction is not
the source of the color dependent astrometric offset; rather, it helps to
slightly decrease the offset. See Appendix H for further details and
justification.

We also investigate stars nearby to determine whether any of them show similar
behavior. For the 70 stars within $3^{\prime\prime}$ of OB110462, we calculate
the average positional offset between F814W and F606W in RA and Dec using Eqs.
9 and 10. We then search for stars where the average positional uncertainty in
F814W and F606W (whichever is larger) is smaller than the average positional
offset to determine which differences are significant. There are 4 stars where
the average positional offset is greater than the average positional
uncertainty in RA, and an additional 4 stars the average offset is greater
than the average positional uncertainty in Dec. Hence, a total of 8 out of the
70 stars near the target also show these significant offsets. Thus, this
effect is seen for roughly 10% of the stars, and so is not very unusual.
Although we currently have no explanation for its significance, it appears
random, and thus include the astrometric offset when analyzing the data,
although it may be attributed to binarity (see §7.4 and §4.1 in the main
paper). See Appendix H for additional details.

####  4.2.6 Comparison to Gaia proper motions and parallaxes

Table 6: Gaia vs. HST Proper Motions

Note. — The source and lens proper motions here have been transformed into the
absolute Gaia proper motion frame, which is offset to the HST proper motion
frame as described in Appendix C. The uncertainties on $\mu_{L}$ and $\mu_{S}$
also reflect the uncertainty in the Gaia to HST proper motion transformation;
the standard error on the mean of that transformation was added in quadrature
to the uncertainties from the proper motion fits. For this reason, the
uncertainties for $\mu_{L}$ and $\mu_{S}$ in this table do not match those in
Tables 7.3 \- 7.5.

Figure 10:  Proper motions for stars in the field of the target. The proper
motions derived from F814W HST observations for stars within
$30^{\prime\prime}$ of the target are shown in gray with 1$\sigma$
uncertainties. Only stars with F814W $<23$ for MB09260, OB110310, OB110462 and
F814W $<22$ for MB10364, OB110037 are shown. For OB110462 there are two models
depending on the likelihood used (default weighted “DW” or equal weighted
“EW”, see §5.3 for details). As the alignment procedure places the stars in a
reference frame where the relative motion is zero, a constant offset must be
added to obtain proper motions in the original Gaia absolute reference frame.
This offset is calculated by matching the stellar positions in HST to those in
Gaia with astrometric_excess_noise_sig$<2$, and then calculating the 3$\sigma$
clipped average weighted by the uncertainty in their difference. The offset
value is given in the title for each field. For targets in Gaia where a
single-star proper motion is estimated (MB10364, OB11037, OB110310), they are
plotted as black squares. The lens and source proper motion as determined from
fitting the HST data with a microlensing model is shown in blue and red
$1-2-3\sigma$ contours, respectively. Note that the red contours are extremely
small.  Figure 11: Comparison of Gaia parallaxes with the fit lens and source
parallax.  Figure 12: Comparison of $b_{SFF}$ inferred from fitting
microlensing events, as compared to upper limits inferred from the high-
resolution HST images. The typical seeing disk _radius_ ranges are shown in
light gray for MOA (MB09260, MB10364) and OGLE (OB110037, OB110310, OB110462).
The inferred value of $b_{SFF}$ from fitting each microlensing event is shown
as a black line, with uncertainties in dark gray. The upper limits on blending
$b_{SFF,max}$ as a function of aperture radius, as estimated using the method
described in §7.7, are shown as the solid red line (_F814W_) and the dotted
blue line (_F606W_). For the three OGLE targets, at the relevant seeing disk
radii, the inferred value of $b_{SFF}$ from the fit are at or below
$b_{SFF,max}$ as inferred from the HST images. For the two MOA targets, the
inferred values of $b_{SFF}$ are higher than the maximum value as inferred
from the HST images. In the case of MB09260, this is likely due to the bright
star around $1.0^{\prime\prime}$ from the target which causes an abrupt change
in $b_{SFF}$ around that radius; proper convolution with a PSF would likely
resolve the difference. In the case of MB10364, the difference is minimal, and
likely just due to the imperfect nature of the comparison (for example, lack
of proper PSF convolution, and F814W not exactly matching the red filter in
MOA/OGLE).

For the three targets with astrometric solutions in Gaia EDR3 (MB10364,
OB110037, and OB110310; Table 3.5), the Gaia proper motions and parallaxes are
compared to the fit proper motions and parallaxes presented in §7 and Tables
7.3 \- 7.5. Direct comparisons are made in Table 4.2.6 and Figures 10 and 11.
Note that we fit a proper motion and parallax to the source and lens along
with an astrometric microlensing model. On the other hand, Gaia fits a proper
motion and parallax to the source and lens (the “target”) as they are
unresolved, and assumes the target is a single star with parallax. The effect
of astrometric lensing in Gaia on the proper motions is negligible (Appendix
D). Additional details about the various Gaia metrics discussed are in
Appendix E and Lindegren et al. (2021a).

To make proper motion comparisons between Gaia and HST, the proper motions
from the HST frame in which the fitting was performed need to be transformed
into the absolute Gaia frame, as the iterative astrometric alignment procedure
described in §4.2 produces a reference frame that is at rest with the average
proper motion of the aligned stars (Lu et al., 2016). See Appendix C on how
the proper motion offset between the HST and Gaia frames is calculated. The
source and lens proper motion in the Gaia frame, as well as the Gaia target
proper motions, are listed in Table 4.2.6 for MB10364, OB110037, and OB110310.
The vector point diagram for all stars from the HST observations transformed
to the Gaia frame, along with the $1-2-3\sigma$ contours from the source and
lens fit, are shown in Figure 10. For MB10364, OB110037, and OB110310 the
target proper motion from Gaia is also included.

##### MB10364

MB10364’s proper motion in Gaia is ($-7.43\pm 0.08$, $-6.80\pm 0.05$) mas/yr.
The fit from HST for the lens is ($-5.11_{-1.10}^{+1.62}$,
$-7.78_{-0.89}^{+0.58}$) mas/yr, inconsistent with Gaia in RA and Dec at $\sim
2\sigma$. The fit from HST for the source is ($-7.56_{-0.12}^{+0.12}$,
$-6.49_{-0.11}^{+0.11}$) mas/yr, inconsistent with Gaia in Dec at $\sim
2\sigma$. MB10364 has parallax $\pi=0.43\pm 0.08$ in Gaia. The source and lens
parallax from the MB10364 fits are $\pi_{S}=0.11^{+0.02}_{-0.02}$,
$\pi_{L}=0.24^{+0.08}_{-0.06}$, neither of which are consistent with the Gaia
value. The Gaia fit for MB10364 has a large renormalized unit weight error
(RUWE = 1.388), and a large astrometric excess noise ($\epsilon=0.406$ mas)
with high significance ($D=12.020$), indicating the single star model is not
providing a good fit. This mismatch is not due to astrometric microlensing
(Appendix D). The most likely explanation for the discrepancies are crowding–
there are several stars close to MB10364 that would cause confusion (Figures 1
and 12). As the source and lens are not resolvable, the fact that 15% of the
Image Parameters Determination (IPD) algorithm has identified a double peak is
likely due to confusion. This means the Gaia measurement is not reliable. See
Appendix E for further discussion.

##### OB110037

OB110037’s proper motion in Gaia is ($2.40\pm 0.13$, $-3.91\pm 0.09$) mas/yr.
The fit from HST for the lens is ($6.27_{-1.20}^{+1.27}$,
$-6.56_{-0.81}^{+0.77}$) mas/yr, inconsistent with Gaia in both RA and Dec at
$\sim 3\sigma$. The fit from HST for the source is ($2.19_{-0.24}^{+0.24}$,
$-3.87_{-0.20}^{+0.20}$) mas/yr, consistent with Gaia. The caveat to this is
that the fit to the HST F606W astrometry is poor (Figure 18). However, the
astrometric lensing model is mostly independent of the proper motion model,
hence the proper motion value can still be believed. OB110037 is very well
measured and behaved in Gaia (it is the only target with a 5-parameter
solution, see Table 3.5). Additionally from the image (Figures 1 and 12)
OB110037 is relatively bright and isolated. The fact that the fit source and
Gaia proper motions are consistent would indicate that the lens is dim in
comparison to the star in Gaia $G$ band. This makes sense as the astrometric
shift for OB110037 is undetectable by the time Gaia begins observing.

OB110037 is well modeled in Gaia by a single source with parallax. Although
the microlensing model is a poor fit to the astrometry due to a time-dependent
color-offset, the proper motions from our model are in good agreement with
Gaia. OB110037 has parallax $\pi=0.19\pm 0.13$ in Gaia. The source and lens
parallax from the OB110037 fits are $\pi_{S}=0.12^{+0.02}_{-0.02}$,
$\pi_{L}=0.58^{+0.14}_{-0.13}$. The source parallax is consistent with the
Gaia value. We note that the source may appear well-behaved in Gaia because
the astrometry is in a single filter. The multi-band HST astrometry may be
useful in identifying binary companions. See Appendix E for further
discussion.

##### OB110310

OB110310’s proper motion in Gaia is ($-2.08\pm 1.12$, $-6.75\pm 0.58$) mas/yr.
The fit from HST for the lens is ($-0.02_{-1.16}^{+1.93}$,
$-4.68_{-2.13}^{+2.39}$) mas/yr, consistent with Gaia due to large
uncertainties in both Gaia and the fit. The fit from HST for the source is
($-2.41_{-0.12}^{+0.12}$, $-7.26_{-0.08}^{+0.08}$) mas/yr, consistent with
Gaia. OB110310 has parallax $\pi=0.53\pm 1.16$ in Gaia (note the OB110310 zero
point correction is an extrapolation), consistent with non-detection of
parallax. The source and lens parallax from the OB110310 fit are
$\pi_{S}=0.10^{+0.02}_{-0.02}$, $\pi_{L}=0.22^{+0.12}_{-0.07}$, which are both
consistent within the very wide uncertainties of Gaia. The fact that the
source and Gaia proper motions are consistent would indicate that the lens is
dim in comparison to the star in Gaia $G$ band; blending in F814W and F606W
also suggest a dark lens (Table 7.5). The fact that in the IPD 55% of transits
have either truncation or multiple gates flagged in one or more windows
indicates likely contamination. OB110310 is not very bright and in a somewhat
crowded region (Figures 1 and 12). The astrometric noise is large
($\epsilon=0.894$ mas) but the value is insignificant ($D=0.332$). Together
these explain why the Gaia measurement does not produce very good constraints.
See Appendix E for further discussion.

###  4.3 HST Photometric Analysis

To obtain precise relative photometry, for each epoch a small constant
magnitude offset is applied to the stars. The offset is calculated by assuming
the reference stars have constant brightness in time, which we define as the
$3\sigma$ clipped mean. As with the positional uncertainties, for the
magnitude uncertainties we used the error on the mean, with an additional
additive error empirically determined during the astrometric alignment
process; details are in Appendix A.

Analogous to the positional transformation, we evaluate the magnitude
transformation by checking how well a constant magnitude describes the stars

$\chi^{2}_{m}=\sum_{i}\Big{(}\frac{m_{i}-m_{0}}{\sigma_{m_{i}}}\Big{)}^{2}$
(15)

where $m_{i}$ are the calibrated magnitudes, $m_{0}$ is the constant magnitude
fit, and $\sigma_{m_{i}}$ is the uncertainty on the calibrated magnitude. The
$\chi^{2}$ distributions for the magnitude residuals of the reference stars
detected in all $N_{ep}$ epochs are shown in the bottom panels of Figure 9.
The expected residual distribution has $N_{ep,detect}-1$ degrees of freedom,
as there is one free parameter in the constant magnitude fit.

Table 7: HST Calibrated Data for Each Target

Note. — Relative positions and magnitude of the target by epoch and filter.

##  5 Microlensing Modeling

The mass of the lens $M_{L}$ in a microlensing event is given by

$M_{L}=\frac{\theta_{E}}{\kappa\pi_{E}},$ (16)

where $\theta_{E}$ is the angular Einstein radius (Equation 17), $\pi_{E}$ is
the microlensing parallax (Equation 21), and $\kappa=4G/(1AU\cdot
c^{2})=8.14{\rm\,mas}/M_{\odot}$ is a constant. Densely sampled photometric
microlensing observations constrain $\pi_{E}$, while astrometric microlensing
observations constrain $\theta_{E}$.101010Theoretically, astrometric
microlensing observations should also be able to constrain $\pi_{E}$, but due
to the cadence of observations, this is currently unachievable.

To measure these quantities of interest, we simultaneously fit the ground-
based photometry and HST photometry and astrometry with a point-source point-
lens (PSPL) microlensing model with parallax. We do not consider models
involving either binary lenses or sources, nor higher-order effects beyond
parallax; these are beyond the scope of this work. Discussion of the need for
models more complex than PSPL with parallax are discussed in §7.

Throughout this section and the remainder of the paper, we report vector
quantities decomposed into their RA and Dec components, subscripted by
“$\alpha*$” and “$\delta$”. Because we work solely in the Equatorial
coordinate system, we also equivalently refer to RA as East and Dec as North,
subscripted by “E” or “N”, where RA increases to the East and Dec increases to
the North.

###  5.1 Microlensing Definitions

All microlensing quantities defined in the following section are in the
heliocentric reference frame.

By rearranging the terms in Equation 16, the Einstein radius, which sets the
angular scale of the microlensing event, can be written

$\theta_{E}=\sqrt{\kappa\pi_{rel}M_{L}}$ (17)

where $\pi_{rel}=\pi_{L}-\pi_{S}$ is the relative parallax of the lens and
source.

The Einstein crossing time $t_{E}$, the time it takes for the source to
traverse the angular radius of the lens and sets the timescale of the events,
is given by

$t_{E}=\frac{\theta_{E}}{\mu_{rel}}$ (18)

where $\mu_{rel}$ is the lens-source proper motion
$|\boldsymbol{\mu}_{S}-\boldsymbol{\mu}_{L}|$.

The source-lens separation on sky
$\boldsymbol{\theta}_{S}-\boldsymbol{\theta}_{L}$ normalized by the Einstein
radius is denoted $\boldsymbol{u}(t)$. The minimum separation is denoted
$\boldsymbol{u}(t_{0})=\boldsymbol{u}_{0}$. The impact parameter $u_{0}$ is
the scalar quantity associated with $\boldsymbol{u}_{0}$. We follow the
convention of Gould (2004) where if $u_{0,E}>0$, the source is to the East of
the lens, and $u_{0}>0$; if $u_{0,E}<0$, the source is to the West of the
lens, and $u_{0}<0$.111111Note that we define our coordinate system
differently than Gould (2004). Gould (2004) works in a system where the
position of the source relative to the lens is defined in a coordinate system
that is right-handed in the relative proper motion and minimum separation
vector. However, we work in a system where the coordinate system is consistent
on the sky; this means a coordinate system based on the relative proper motion
and minimum separation vector does not always preserve handedness. If
rectilinear motion is assumed, the lens-source separation is given by

$\boldsymbol{u}(t)=\boldsymbol{u}_{0}+\frac{t-t_{0}}{t_{E}}\boldsymbol{\hat{\mu}}_{rel}.$
(19)

However, an Earthly observer’s perspective of the lensing event is modulated
by the Earth’s motion around the Sun. For events with long duration
($t_{E}\gtrsim 3$ months), the Earth’s orbital motion violates this
rectilinear assumption and must be taken into account. This modifies
$\boldsymbol{u}(t)$ to

$\boldsymbol{u}(t)=\boldsymbol{u}_{0}+\frac{t-t_{0}}{t_{E}}\boldsymbol{\hat{\mu}}_{rel}-\pi_{E}\boldsymbol{P}(t)$
(20)

where

$\pi_{E}=\frac{\pi_{rel}}{\theta_{E}}$ (21)

is the microlensing parallax and $\boldsymbol{P}(t)$ is the parallax vector,
defined to be the Sun-Earth separation vector normalized by 1 AU. The
microlensing parallax vector $\boldsymbol{\pi}_{E}$ (not to be confused with
the parallax vector $\boldsymbol{P}(t)$) encodes the magnitude of the
microlensing parallax and the direction of the relative source-lens proper
motion:

$\boldsymbol{\pi}_{E}=\pi_{E}\boldsymbol{\hat{\mu}}_{rel}.$ (22)

The photometric brightening of the source is given by

$A(u)=\frac{u^{2}+2}{u\sqrt{u^{2}+4}}.$ (23)

where the total flux $F(t)$ in the telescope aperture is

$F(t)=A(t)F_{S}+F_{L}+F_{N}$ (24)

where $F_{S}$, $F_{L}$, and $F_{N}$ are the fluxes of the source, lens, and
neighboring un-lensed stars in the aperture of the telescope, respectively.
The source flux fraction is

$b_{SFF}=\frac{F_{S}}{F_{S}+F_{L}+F_{N}}$ (25)

and quantifies the fraction of light lensed in an observed microlensing event.
Note that $b_{SFF}$ depends on the observing wavelength and seeing/aperture.
The non-source flux $F_{L}+F_{N}$ is also called blend flux. Blend flux
decreases $b_{SFF}$ and dilutes the magnitude of both the photometric
amplification and astrometric shift.

A PSPL photometric microlensing event is characterized by five geometric
parameters: $t_{0}$, $u_{0}$, $t_{E}$, and $\boldsymbol{\pi}_{E}$. For each
telescope that observes this event, two additional parameters, $m_{base}$ and
$b_{SFF}$, are needed to describe each lightcurve; these depend on the seeing
and camera filter. As all photometric microlensing observable quantities are
normalized by $\theta_{E}$, the lens mass cannot be determined.

In contrast, astrometric microlensing provides a direct measurement of
$\theta_{E}$. The apparent position of the source, i.e. the centroid of the
lensed source images $\boldsymbol{\theta}_{S,c}$, is given by

$\boldsymbol{\theta}_{S,c}(u,\theta_{E})=\frac{(u^{2}+3)\boldsymbol{u}\theta_{E}}{u^{2}+2}.$
(26)

Assuming no blended light, the difference between the source’s apparent and
true positions, i.e. the astrometric shift, is given by

$\boldsymbol{\delta}_{c}(u,\theta_{E})=\frac{\boldsymbol{u}\theta_{E}}{u^{2}+2}.$
(27)

The astrometric shift $\delta_{c}$ is maximized at $u=\sqrt{2}$. This
corresponds to the value of the maximum astrometric shift

$\delta_{c,max}=\frac{\theta_{E}}{\sqrt{8}},$ (28)

which is directly proportional to the Einstein radius.

Note that Equation 27 is for an unblended event, i.e. $b_{SFF}=1$. If an event
has $b_{SFF}<1$, the non-source light would dilute the astrometric shift (see
Equations 11 and 12 in Lam et al. (2020)). For the five candidates analyzed
here, $b_{SFF}\sim 1$ in the HST filters, so assuming that the astrometry is
unblended in HST is valid.

###  5.2 Modeling framework

We perform parameter estimation using a Bayesian framework. Bayes’ theorem

$\pi(\boldsymbol{\Theta})\mathcal{L}(\boldsymbol{y}|\boldsymbol{\Theta})=\mathcal{Z}(\boldsymbol{y})\mathcal{P}(\boldsymbol{\Theta}|\boldsymbol{y})$
(29)

relates the prior $\pi$ and likelihood $\mathcal{L}$ to the evidence
$\mathcal{Z}$ and posterior $\mathcal{P}$. The goal of parameter estimation is
to calculate $\mathcal{P}$. The likelihood
$\mathcal{L}(\boldsymbol{y}|\boldsymbol{\Theta})$ is presented in Appendix F,
and the priors $\pi(\boldsymbol{\Theta})$ are discussed in Appendix G.

The data is fit using MultiNest (Feroz et al., 2009), an implementation of the
nested sampling algorithm (Skilling, 2004). Nested sampling produces an
estimate of the evidence
$\mathcal{Z}=\int\pi(\boldsymbol{\Theta})\mathcal{L}(\boldsymbol{\Theta})\,d\boldsymbol{\Theta}$,
and as a by-product, the posterior $\mathcal{P}$. In contrast to methods such
as MCMC, nested sampling is designed to better explore multimodal likelihood
spaces; however, care must still be taken to ensure that all local minima are
explored.

###  5.3 Likelihood weighting

There is a question of how best to combine the photometry and astrometry data
sets, as they represent two different types of measurements. In particular,
the question is how much weight each dataset should receive given that there
are several orders of magnitude more ground-based data points than HST data
points. Thus, the ground-based photometry has an outsize effect on the
likelihood. It can be argued that each data point should contribute equally to
the likelihood. We consider this to be “default weight” (hereafter abbreviated
as DW) likelihood, i.e.

$\log\mathcal{L}_{tot}=\log\mathcal{L}_{O,phot}+\log\mathcal{L}_{H,phot}+\log\mathcal{L}_{H,ast}$
(30)

where $\mathcal{L}_{O,phot}$, $\mathcal{L}_{H,phot}$, and
$\mathcal{L}_{H,ast}$ are the likelihoods of the OGLE or MOA photometry, HST
photometry, and HST astrometry respectively.

However, it could also be argued that each dataset is independent, and so
should each contribute equally to the likelihood. We call this the “equal
weight” (hereafter abbreviated as EW) likelihood, i.e.

$\log\mathcal{L}_{tot}=\frac{\log\mathcal{L}_{O,phot}}{n_{O,phot}}+\frac{\log\mathcal{L}_{H,phot}}{n_{H,phot}}+\frac{\log\mathcal{L}_{H,ast}}{n_{H,ast}}$
(31)

where $n_{O,phot}$, $n_{H,phot}$, and $n_{H,ast}$ are the number of data
points in in the OGLE photometry, HST photometry, and HST astrometry.

These different likelihoods are essentially giving the different datasets
different weights, in the case that they are inconsistent with each other. For
all targets, we fit using the default weight likelihood; for OB110462 we
additionally fit using the equal weight likelihood.

##  6 Constraining the nature of the lens

By modeling photometric and astrometric microlensing data as described in §5,
the lens’ mass can be measured. However, a mass measurement alone cannot
distinguish the difference between a $5M_{\odot}$ stellar lens and a
$5M_{\odot}$ BH. Additional information about the lens’ brightness is needed,
which can be estimated using the source flux fractions in the high-resolution
HST filters. By calculating the brightest star allowed by the inferred source
flux fractions, we can determine whether a luminous lens (i.e. stellar lens)
or dark lens (i.e. compact object lens) scenario is more likely.

We follow a procedure similar to Wyrzykowski et al. (2016) in order to
calculate the probability of a dark lens. The two main differences between the
analysis of Wyrzykowski et al. (2016) and ours is that they work with
photometric data only and must invoke a Galactic model to obtain lens masses
and distances, while in our case astrometric data allows us to fit lens mass
and distances. This greatly simplifies the analysis as we do not need to
calculate the Jacobian to transform the PDFs between different variables. In
addition, we use realistic stellar and Galactic models to determine the lens
luminosities, instead of relying on simple mass-luminosity scaling relations.

From fitting the microlensing data, we have posterior distributions for the
lens mass $M_{L}$, distance $d_{L}$, baseline magnitude $m_{\rm base}$, and
source flux fraction $b_{SFF}$. In the following analysis we reassign any fit
values with $b_{SFF}\geq 1$ to instead be equal to 0.999999. This is because
values of $b_{SFF}\geq 1$ would result in an infinite lens magnitude (Equation
32).121212One might worry this could bias the results. However, a visual check
comparing the posteriors (both individual and joint) for lens mass and
distance don’t show a change when excluding fits with $b_{SFF}>1$. In
particular, the lens mass and distances are nearly independent of the lens
magnitude. In addition, if the true value of $b_{SFF}$ is 1, we would expect
some scatter around that value, including unphysical $>1$ values.

We draw a random sample of 1000 stars from our posterior. For each star we
calculate 1) the brightest lens allowed by $b_{SFF}$ and $m_{base}$, and 2)
the brightest star allowed by $M_{L}$ and $d_{L}$.

To calculate 1), by assuming there are no contaminating neighbor stars,
$b_{SFF}=F_{S}/(F_{S}+F_{L})$, an upper limit can be placed on the brightness
of the lens:

$m_{L}=m_{\rm base}-2.5\log_{10}(1-b_{SFF}).$ (32)

We denote this $m_{L}(b_{SFF},m_{base})$.

To calculate 2), we use the simple stellar population synthesis code SPISEA
(Hosek et al., 2020) to generate a suite of simple stellar populations (SSPs)
to simulate the possible lens population. We use the MISTv1.2 solar
metallicity isochrones (Choi et al., 2016), get_merged_atmosphere atmosphere
model131313This is a combination of the ATLAS9 (Castelli & Kurucz, 2004),
PHOENIX v16 (Husser et al., 2013), BTSettl (Baraffe et al., 2015), abd
Koester10 (Koester, 2010) models. For further details see Appendix B of Hosek
et al. (2020)., Damineli et al. (2016) reddening law, and Kroupa (2001)
initial mass function (IMF) over the mass range $0.1M_{\odot}<M<120M_{\odot}$.
Each cluster is $10^{4}M_{\odot}$, in order to reduce stochastic effects in
the sampling of the IMF.

SSPs are generated at the distances spanned by each target’s lens distance
posteriors, sampled every 0.25 kpc. SSPs of ages 7.0 to 10.0 $\log_{10}$ years
in increments of 0.5 $\log_{10}$ years are simulated at each sampled distance
in order to cover the age range of disk and bulge stars. The stellar age
distribution as a function of distance is drawn from the Galaxia Milky Way
stellar simulation (Sharma et al., 2011), which implements a version of the
Besançon Galactic model of Robin et al. (2003). As our target sight lines are
toward the bulge, most of the stars are old, with $\gtrsim 85\%$ of stars
being at least 9.0 $\log_{10}$ years old, although younger stars tend to be at
closer distances as disk stars dominate. To calculate the simulated stars’
apparent magnitudes, we use the 3-D extinction map of Marshall et al. (2006),
accessed via the dustmaps software package (Green, 2018), to obtain an $A_{\rm
K_{s}}$ value for each distance.

For each sample from the posterior, we sample one of the SSPs at the
corresponding distance at a particular stellar age from the Galactic
distribution. All stars in the simulated SSP with masses within 10% of the
sampled posterior mass are identified, and the brightest apparent magnitude of
the star in this group is denoted $m_{L}(M_{L},D_{L},\bigstar)$. If no stars
are found within 10% of the sampled posterior mass, this indicates that stars
of that mass have all died and formed compact objects.

If $m_{L}(M_{L},D_{L},\bigstar)<m_{L}(b_{SFF},m_{base})$, a star with the
inferred stellar mass from stellar evolution models would be too bright to be
hidden in the blended light allowed by the fit. This means that the lens is
not a star. We take this to mean the lens is a compact object, and hence a
dark lens. Thus the lower limit on the probability of a dark lens in the
observed filter is the fraction of samples where
$m_{L}(M_{L},D_{L},\bigstar)<m_{L}(b_{SFF},m_{base})$ or no mass match is
found.

Objects with $m_{L}(M_{L},D_{L},\bigstar)>m_{L}(b_{SFF},m_{base})$ are samples
where a luminous stellar lens is consistent with the inferred amount of
blending. Note that $m_{L}(b_{SFF},m_{base})$ is an upper limit on the
brightness of the lens, as we assumed all the blended flux in Equation 32 is
due to the lens; it could be due to unresolved unrelated neighbor stars. Thus
the upper limit on the probability of a stellar lens in the observed filter is
thus the fraction of samples where
$m_{L}(M_{L},D_{L},\bigstar)>m_{L}(b_{SFF},m_{base})$.

For dark lenses (i.e. any objects where
$m_{L}(M_{L},D_{L},\bigstar)<m_{L}(b_{SFF},m_{base})$ or no star of the same
mass was found in the simulated SSP), we categorize them as brown dwarfs (BD),
white dwarfs (WD), neutron stars (NS), or BHs by simplistically sorting them
by their masses:

  * •

Brown dwarfs: $M<0.2M_{\odot}$

  * •

White dwarfs: $0.2M_{\odot}<M<1.2M_{\odot}$

  * •

Neutron stars: $1.2M_{\odot}<M<2.2M_{\odot}$

  * •

Black holes: $M>2.2M_{\odot}$.

In reality there is overlap between white dwarf and neutron star masses, and
the overlap between neutron star and BH masses (if they overlap at all) is
unknown. In addition, the maximum brown dwarf mass set by stellar physics is
$0.08M_{\odot}$, but we extend this up to $0.2M_{\odot}$ to have continuity
between the lowest mass WDs of around $0.2M_{\odot}$. Hence these values are
only approximate.

The above analysis is performed for both the HST F606W and F814W filters. The
reported probabilities are the joint constraint. A lens is dark if no masses
consistent with a stellar lens are found in either filter. A lens is also dark
if the maximum inferred lens flux is insufficient to hide a star (i.e.
$m_{L}(M_{L},D_{L},\bigstar)<m_{L}(b_{SFF},m_{base})$) in either filter. A
lens is luminous if the maximum inferred lens flux is sufficient to hide a
star (i.e. $m_{L}(M_{L},D_{L},\bigstar)>m_{L}(b_{SFF},m_{base})$) in both
filters. We do not perform this analysis for the OGLE or MOA photometry
parameters as the high resolution HST images show the seeing-limited apertures
have unrelated neighbor stars in the blend, and hence the limits will all be
weaker than using HST.

Figure 13:  Probabilities for dark (_gray_) and luminous lenses (_red_) , as a
function of lens mass for each target. Lenses that are dark because all stars
of that mass have already evolved and died (implying they are compact objects)
are subscripted with $A$ (_square hatch_), while lenses that are dark because
a star would be too bright to be hidden in the allowed lens flux are
subscripted with $B$ (_diamond hatch_). The probabilities for the luminous
lenses are upper limits, while the probabilities for the dark lenses are lower
limits, since the method described in §6 only places an upper limit on the
brightness of the lens. Note that there are two fits for OB110462, one with
equal weighting to the astrometry and photometry data (OB110462 EW) and one
with the default weighting of the astrometry and photometry data (OB110462
DW). See §5.3 for details.  Figure 14: OB110462 astrometry, using the equally
weighted likelihood. _Left column, top to bottom_ : RA vs. time with maximum
likelihood (MLE) unlensed source motion model subtracted; residuals to the MLE
model for RA vs. time fit. HST F814W astrometry data is shown in red; HST
F606W astrometry data is shown in blue. The MLE model is shown in black. Fifty
random draws from the posterior distribution are shown in light gray. _Middle
column, top to bottom_ : Same as left column, except Dec instead of RA. _Right
panel_ : astrometry as seen on-sky, in the barycentric frame. OB110462 shows a
strong $>1$ mas astrometric microlensing signal.  Figure 15: Same as Figure
14, but using the default weighted likelihood for OB110462.  Figure 16:
MB09260 astrometry. Same as Figure 14, but for MB09260. The astrometric signal
is small, and around the limit of the precision of the F814W measurements.
Figure 17: Same as Figure 14, but for MB10364. The astrometric signal is
small, at or below the limit of the precision of the F814W measurements.
Figure 18: Same as Figure 14, but for OB110037. The photometry and astrometry
seem to indicate this object is a binary.  Figure 19: Same as Figure 14, but
for OB110310. The astrometric signal is small, and around the limit of the
precision of the F814W measurements.

##  7 Results in Detail

Table 8: 99.73% credible intervals/upper limits

Note. — There are two entries for OB110462, one based on a fit with “default
weighting” (OB110462 DW) and one based on a fit with “equal weighting”
(OB110462 EW). See §5.3 for details.

Table 9: Lens type probabilities (%)

Note. — The Star probabilities are upper limits, while the brown dwarf (BD),
white dwarf (WD), neutron star (NS), and BH probabilities are lower limits.
This is because the luminous lens probabilities are upper limits, and the dark
lens probabilities are lower limits; see §6 for details. Note that there are
two entries for OB110462, one based on a fit with “default weighting”
(OB110462 DW) and one based on a fit with “equal weighting” (OB110462 EW). See
§5.3 for details.

A large ($>1$ mas) astrometric microlensing signal is detected in OB110462
(§7.1). A filter-dependent astrometric signal is also detected in OB110037
(§7.4); averaging the astrometry across the F814W and F606W filters shows a
small but significant 0.4 mas signal. The remaining three targets MB09260
(§7.2), MB10364 (§7.3), and OB110310 (§7.5) have astrometric signals that are
below HST’s detection threshold. Table 7 reports the lens masses, Einstein
radii, maximum astrometric shifts, and microlensing parallaxes either as
measured values or upper limits.

§7.1 \- §7.5, the posteriors of the joint photometry and astrometry
microlensing fits for the individual targets are presented. Tables 7.2 \- 7.5
and Tables 1 - 2 in Lam et al. (2022) list the median and $1\sigma$ (68%)
credible intervals for each parameter, in addition to the maximum a posteriori
(MAP) and maximum likelihood (MLE)141414For most of the targets, the MAP
solution is equal to the MLE solution, i.e. values of the parameters for the
mode of the posterior distribution are the same as those where the likelihood
function is maximized. solution.

Using the methodology described in §6, constraints are placed on the lens
types of each target. We report the relative probabilities of Star:BD:WD:NS:BH
in Table 7. Figure 13 shows the distribution of dark vs. luminous lens
probability as a function of mass for each target. Table 7 lists the upper
limit on the probability of a stellar lens and lower limits on the
probabilities of different dark lenses for each target. In §7.9, we compare
the yield of BHs from our search to that expected from simulations assuming
there are $\sim 10^{8}$ BHs in the Milky Way.

###  7.1 OB110462

The results for OB110462 are presented in the main paper, §4.1 of Lam et al.
(2022).

###  7.2 MB09260

The data and model for MB09260 are shown in Figures 3 (photometry) and 16
(astrometry), and the fit posteriors are summarized in Table 7.2. The inferred
Einstein crossing time $t_{E}$ is $142.64_{-2.87}^{+3.49}$ days, the
microlensing parallax $\pi_{E}$ is $0.09_{-0.01}^{+0.03}$, the Einstein radius
$\theta_{E}$ is $1.04_{-0.39}^{+0.42}$, and the lens mass is
$1.37_{-0.60}^{+0.74}M_{\odot}$.151515These are the values for the posteriors
before they are split into modes, which is why the values differ from what is
reported in Table 7.2. However, since these parameters are all globally
unimodal, their distributions nearly identical across modes within the
uncertainties, and neither mode is strongly favored, the values are not very
different and we report the global median and uncertainties.

The probability that MB09260 has a dark lens is at least 96%, and the
probability of a stellar lens is at most 4%. The relative probabilities of
Star:BD:WD:NS:BH are 4:0:38:44:14. Stellar lenses are only allowed below
$1M_{\odot}$. A white dwarf or neutron star is the most probable type of lens,
with black holes possible but less likely.

Table 10: MB09260 Fit Values, Modes

Note. — The columns list the median $\pm 1\sigma$ (68%) credible intervals,
maximum a posteriori (MAP) solution, and and maximum likelihood estimator
(MLE) solution for the microlensing parameters of MB09260. The posterior is
multimodal (primarily in $u_{0}$); it has been split and the parameters for
each mode reported separately. The fraction that each mode contributes to the
whole posterior ($\Sigma w_{i}$) and log evidence ($\log\mathcal{Z}$) are
listed for each mode at the end of the table.

###  7.3 MB10364

The data and model for MB10364 are shown in Figures 4 (photometry) and 17
(astrometry), and the fit posteriors are summarized in Table 7.3. The inferred
Einstein crossing time $t_{E}$ is $61.11^{+0.24}_{-0.24}$ days, the
microlensing parallax $\pi_{E}$ is $0.27^{+0.01}_{+0.01}$, the Einstein radius
$\theta_{E}$ is $0.46^{+0.31}_{-0.21}$, and the lens mass is
$0.21^{+0.14}_{-0.10}M_{\odot}$.

MB10364 is a low mass object, with the possibility of a neutron star or BH
lens ruled out. The relative probabilities of Star:BD:WD are 36:29:35.

Table 11: MB10364 Fit Values

Note. — Same as Table 7.2, but for MB10364. The solution is unimodal, hence
there is only one column and $\Sigma w_{i}$ and $\log\mathcal{Z}$ are not
reported.

###  7.4 OB110037

Table 12: OB110037 Fit Values

Note. — Same as Table 7.3, but for OB110037.

The data and model for OB110037 are shown in Figures 5 (photometry) and 18
(astrometry), the fit posteriors are summarized in Table 7.4. The inferred
Einstein crossing time $t_{E}$ is $92.78^{+2.63}_{-2.60}$ days, the
microlensing parallax $\pi_{E}$ is $\pi_{E}$ is $0.37_{-0.01}^{+0.01}$, the
Einstein radius $\theta_{E}$ is $1.24^{+0.36}_{-0.35}$, and the lens mass is
$0.41^{+0.12}_{-0.12}M_{\odot}$.

The probability that OB110037 has a dark lens is at least 26%, and the
probability of a stellar lens is at most 74%. The relative probabilities of
Star:BD:WD:NS:BH are 74:0:26:0:0. Stellar lenses are only allowed below
$0.6M_{\odot}$, and white dwarfs are the only type of compact objects allowed.

The lightcurve of OB110037 appears to have some type of perturbation at MJD
$\sim$ 55690\. This feature is also apparent in the MOA lightcurve, raising
our confidence that the lightcurve feature is real. This perturbation may be
attributed to a binary lens.

In addition, the astrometry fit, in particular for the F606W filter, is quite
poor (Figure 18). Although the first 5 observations from 2011-2012 seem to
agree between the two filters, a drastic difference that increases as time
goes on begins in 2013-2017. This may be attributed to a binary source.

It is curious that the photometry may be better explained by a binary lens
model, while the astrometry is likely better explained by a binary source
model. Re-analysis of this event with both types of binary models would be a
worthwhile and interesting pursuit, but is beyond the current scope of this
paper.

###  7.5 OB110310

The data and model for OB110310 are shown in Figures 6 (photometry) and 19
(astrometry), and the fit posteriors are summarized in Table 7.5. The inferred
Einstein crossing time $t_{E}$ is $82.64^{+2.18}_{-1.50}$ days, the
microlensing parallax $\pi_{E}$ is $0.13_{-0.04}^{+0.08}$, and Einstein radius
$\theta_{E}$ is $0.88^{+0.61}_{-0.42}$, and the lens mass is
$0.78^{+0.71}_{-0.39}M_{\odot}$.161616These are the values for the posteriors
before they are split into modes, which is why the values differ from what is
reported in Table 7.5. However, since these parameters are all globally
unimodal, their distributions nearly identical across modes within the
uncertainties, and neither mode is strongly favored, the values are not very
different and we report the global median and uncertainties.

The probability that OB110310 has a dark lens is at least 95%, and the
probability of a stellar lens is at most 5%. The relative probabilities of
Star:BD:WD:NS:BH are 5:3:65:22:5. Stellar lenses are only allowed below
$1M_{\odot}$. A white dwarf or neutron star is the most probable type of lens,
although brown dwarfs and black holes are still allowed at the low and high
mass ends, respectively.

Table 13: OB110310 Fit Values, Modes

Note. — Same as Table 7.2, but for OB110310.

###  7.6 Source properties inferred from CMDs

As there is very little blending for all the targets in the high resolution
F814W and F606W filters ($b_{SFF}\sim 1$), the difference between the target
and source on the CMD does not change much in color nor magnitude space
($\Delta$F814W $\lesssim 0.1$ mag and $\Delta$(F606W $-$ F814W) $\lesssim 0.1$
mag, Figure 2). Based on a CMD analysis, the source stars in MB09260, MB10364,
OB110037, and OB110310 are likely red giant stars in the bulge, as is typical
for microlensing events in this part of the sky.

The source of OB110462 in the CMD is around the main sequence turnoff (MSTO)
on the redder and more luminous side of the main sequence, suggesting it is
most likely a giant or sub-giant star. However, a main sequence source could
still be consistent.

The region of the CMD around the MSTO contains both foreground stars from the
disk as well as bulge stars. We compare the proper motions of OB110462’s
source to stars in the bright blue foreground as well as in the bulge red
giant branch to determine which population it most likely belongs to (Figure
20). The source is consistent with either population, although it falls within
the bulk of the bulge population and more on the edge of the disk population,
hence we consider it is most likely a bulge star. This is also consistent with
the source parallax $\pi_{S}=0.11\pm 0.02$ inferred from the fit which also
indicates the source is likely in the bulge.

MB09260 and OB110310 are in the the most highly reddened field, OB110462 and
OB110037 are in intermediately reddened fields, and MB10364 is in the least
reddened field. This serves as a reminder that within the bulge the amount of
extinction is highly variable, even over relatively small fields.

Figure 20: _Top left:_ CMD with the blue foreground and red giant branch (RGB)
stars marked. The source color and magnitude as inferred from the fit are
shown as the yellow star. _Top right:_ VPD of the foreground and RGB stars.
The source proper motion as inferred from the fit is shown as the yellow star.
_Bottom left:_ Histogram of the foreground and RGB RA proper motions. The
proper motion of the source is shown in the black line. _Bottom right:_ Same
as bottom left, but for Dec. The source proper motion is consistent with
either the foreground or RGB bulge population, although it falls closer to the
bulk of the RGB population, hence we consider it is most likely a bulge star.

###  7.7 Verifying fit results with $b_{SFF}$

The fitting results are validated by comparing the best-fit OGLE or MOA I-band
source flux fraction ($b_{SFF,O}$, $b_{SFF,M}$ in Tables 1 - 2 in the main
paper, and Tables 7.2 \- 7.5) to the high-resolution F814W HST images (Figure
12). We only compare the F814W images, as it is has a similar effective
wavelength to OGLE I-band. Assuming a seeing disk radius of $\sim
0.65^{\prime\prime}$ for OGLE and $\sim 1.25^{\prime\prime}$ for MOA, we add
up the flux from all stars detected by hst1pass within this radius around the
target. Next, by assuming all the flux from the microlensing target at
baseline as seen by HST is due to the source (i.e. no flux from the lens or
blended neighbor stars within HST’s diffraction limited aperture
$\theta_{HST}\sim 0.09^{\prime\prime}$), we can estimate an upper limit on the
source flux fraction for a given seeing disk radius
$\theta_{see}>\theta_{HST}$ using

$b_{SFF,HST-derived}\leq\frac{F_{event}}{\sum_{i}F_{i}}$ (33)

where $F_{i}$ are the fluxes of the $i$ stars within $\theta_{see}$ and
$F_{event}$ is the flux of the event at baseline in HST. This estimate on the
upper limit of the source flux fraction derived using HST, $b_{SFF,O,HST-
derived}$ or $b_{SFF,M,HST-derived}$, is compared to the values inferred from
the lightcurve fit, $b_{SFF,O}$ or $b_{SFF,M}$. This approach is sufficient to
cross check our results, considering other differences prevent a perfect
comparison (e.g. HST F814W is not a perfect match to the OGLE I filter).

From the photometry fits, MB10364 has $b_{SFF,M}\sim 0.93$, OB110037 has
$b_{SFF,O}\sim 0.90$, OB110310 has $b_{SFF,O}\sim 0.97$, and OB110462 has
$b_{SFF,O}\sim 0.05$, which are all below or no more than $\sim 0.05$ higher
than the upper limits inferred from the HST-derived values shown in Figure 12.

MB09260 has $b_{SFF,M}\sim 0.60$, which is higher than the upper limit of
$b_{SFF,HST-derived}\sim 0.4$ in F814W at $1.25^{\prime\prime}$ estimated from
the high resolution HST images. However, there is an abrupt transition from
$b_{SFF,HST-derived}\sim 1$ to $b_{SFF,HST-derived}\sim 0.4$ (indicating a
very bright star) at a radius of 1” from the source. The best-fit
$b_{SFF,M}\sim 0.60$ falls within these two values, indicating that properly
treating this situation by convolving the HST image with a MOA PSF would
result in a better match. Hence, MB09260 seems consistent with the high
resolution image.

The consistency between the source flux fraction inferred from the fits and
the high resolution imaging provides an extra degree of confidence in the
inferred fit values.

###  7.8 Prospects for resolving lens and source

By definition, microlensing means the source and lens are not resolvable.
However, after many years, the lens and source can separate far enough to be
resolved (e.g. Batista et al. (2015); Bennett et al. (2015)). In the case of
dark lenses, such as black holes, non-detections of the lens many years after
the event can be used to place constraints on its properties (e.g. Abdurrahman
et al. (2021)). For the five targets in this paper, we provide estimates to
determine if and when taking late-time follow-up data would enable such
analyses.

From the results of the fit, we can estimate the time necessary to resolve the
source and lens $t_{\rm res}$ via

$t_{\rm res}=\theta_{\rm res}/\mu_{\rm rel}$ (34)

where the relative proper motion $\mu_{\rm rel}$ comes from fitting the data,
and the minimum angular separation $\theta_{\rm res}$ can be estimated using
the Rayleigh criterion. For HST with a mirror diameter of 2.4m, $\theta_{\rm
res}$ corresponds to 63.53 mas in F606W and 85.35 mas in F814W.

Assuming there are no contaminating stars, $b_{SFF}=F_{S}/(F_{S}+F_{L})$. Then
the ratio of lens to source flux, or contrast, is
$F_{L}/F_{S}=(1-b_{SFF})/b_{SFF}$. Note that $b_{SFF}>1$ are allowed by some
of the fits which would lead to an unphysical negative contrast. $b_{SFF}>1$
is often referred to as “negative blending” because it means there is negative
non-source flux (see Equation 25). Negative blending can occur in ground-based
images if the background subtraction is imperfect (Park et al., 2004). For
HST, where the photometry is not derived from difference imaging but rather
PSF fitting, this is unlikely to be the case. However, when fitting a
microlensing model to data, negative blending can occur if $b_{SFF}\approx 1$,
simply due to normal photometric uncertainties. This is the most likely
explanation for HST; all fits with negative blending values have posterior
probabilities that encompass $b_{SFF}=1$. Thus, if a contrast value is
negative, we cap it at 0.171717This is technically not correct, analogous to
truncating negative parallaxes in Gaia; a proper treatment would involve a
full Bayesian analysis as explained in Luri et al. (2018). Such a treatment is
beyond the scope of this paper. We proceed with capping the contrast at 0, as
if $b_{SFF}\sim 1$ then the target is likely very dim or dark, and the
contrast unlikely to be detectable. In addition, the fit values for $\mu_{\rm
rel}$ and $b_{SFF}$ are nearly independent, and so these results for $t_{\rm
res}$ and $F_{L}/F_{S}$ can be thought of as independent.

The resolving time after photometric peak and maximum contrast in the F814W
and F606W HST filters for each target are listed in Table 7.8. Currently, only
OB110037 could potentially have its source-lens pair resolved in F606W. All
the other targets have lens-source resolving times at least 15 years post-
photometric peak or very low contrast. In addition, since they all have very
high source flux fractions $b_{SFF}$, their lens-source contrasts are very
low, which means even after enough time has passed for the lens and source to
separate, a luminous lens would be difficult to detect. Thus, high resolution
imaging of OB110037 in the near future could confirm the results of the fit if
the separating source and lens could be detected. For the other targets, the
absence of any lens detection would imply consistency with the results
presented here, but could not confirm them; however, any detection of a lens
would imply the fit results here are incorrect.

Table 14: Lens/source resolvability

Note. — There are two entries for OB110462, one based on a fit with “default
weighting” (OB110462 DW) and one based on a fit with “equal weighting”
(OB110462 EW). See §5.3 for details.

###  7.9 Number of detected BHs

Table 15: Fraction of expected BH detections vs. $t_{E}$ from PopSyCLE
simulation

Note. — For each target, the $t_{E}$ range is the median $\pm 3\sigma$. For
OB110462 there are two entries, one with equal weighting to the astrometry and
photometry data (OB110462 EW) and one with the default weighting of astrometry
and photometry data (OB110462 DW). See §5.3 for details.

Next, we compare our observed BH yield to the theoretical expectation
calculated using the PopSyCLE simulations. For a sample of simulated events
that would be observable by OGLE (see Table 4 of Lam et al. (2020)), we
calculate the fraction of those events due to BHs as a function of the
Einstein crossing time $t_{E}$. Assuming the OGLE observability criterion for
the MOA sample is not strictly correct; however, OGLE and MOA often observe an
overlapping set of events, so this approximation suffices.

We wish to calculate the probability of detecting $k$ BHs in our sample given
$n$ events, where the probability of detecting a BH in the $i$-th event is
$p_{i}$. This is described by a Poisson binomial distribution, which
characterizes a “success/no success” experiment with $n$ independent trials,
where the $i$-th trial has probability $p_{i}$ of success (Wang, 1993). The
probability of $k$ successes is given by

$P(K=k)=\sum_{A\in F_{k}}\prod_{i\in A}p_{i}\prod_{j\in A^{c}}(1-p_{j})$ (35)

where $F_{k}$ is the set of all subsets of $k$ integers that can be selected
from the set $\\{1,...,n\\}$, i.e.

$F_{k}=\\{A:A\subseteq\\{1,...,n\\},|A|=0,1,...,n\\}.$

$|A|$ is the number of elements in $A$, and $A^{c}$ is the complement of $A$.
In the limit where all $p_{i}$ are equal, the Poisson binomial distribution is
equivalent to the ordinary binomial distribution.

In our case, there are $n=5$ independent trials (i.e. microlensing targets).
We calculate the probability of success $p_{i}$ (i.e. BH detection) using
PopSyCLE. We define the success probability for the $i$-th target as the
fraction of BH lensing events in PopSyCLE over the range of $t_{E}$ inferred
from the fit

${\rm med}(t_{E})-3\sigma<t_{E}<{\rm med}(t_{E})+3\sigma.$

The probabilities of BH detection for each target are listed in Table
7.9.181818The fact that BH probability does not increase monotonically as a
function of $t_{E}$ appears to contradict the claims made in the main paper
that events with longer $t_{E}$ are more likely to have a BH lens. The issue
is that the uncertainty in BH probability rapidly increases at longer $t_{E}$,
as the number of events steeply drops off for $t_{E}\gtrsim 50$ days (Figure 4
in the main paper). For example, for both OB110462 solutions, the estimates
are based on $<10$ events. Because of the small number of events, the estimate
is highly sensitive to the particular $t_{E}$ range.

The results of evaluating Eq. 35 for $k=0,...,5$ BH detections assuming
success probabilities $p_{i}$ is presented in Figure 6 in the main paper.
PopSyCLE’s Galactic model contains $2\times 10^{8}$ BHs ranging from
$5-16M_{\odot}$ (Lam et al., 2020). There are no $2-5M_{\odot}$ mass gap NSs
or BHs in the simulation, and hence no exact OB110462 analogue. Thus, we
consider a mass gap detection as falling between 0 to 1 BH detections. The
probability of detecting 0 or 1 BHs are $\sim$25% and $\sim$40%, respectively.
This estimate is consistent with our single detection of a NS-BH mass gap
object.

##  8 Discussion in Detail

OB110462 is the first definitive detection of a compact object discovered with
astrometric microlensing. Depending on the fit, it is either a neutron star
(50% probability for the EW fit), a black hole (44% probability for the EW
fit, 100% probability for the DW fit), or a white dwarf (6% probability for
the EW fit). MB09260 and OB110310 are mostly likely white dwarfs or neutron
stars, although a NS-BH mass gap object cannot be ruled out. MB10364 and
OB110037 are definitively low mass objects; OB110037 is most likely a star or
white dwarf, while MB10364 is either a star, brown dwarf, or white dwarf.

Discussion of the BH yield as compared to theoretical expectations is found in
§5 of the main paper (Lam et al., 2022). Here, we discuss in more detail our
sample of BH candidates and the questions it raises about the Galactic BH
population (§8.1), additional observations of OB110462 (§8.2), potential
sources of systematics in OB110462 observations (§8.3), and the future of BH
microlensing searches (§8.4 and §8.5).

###  8.1 OB110462 in comparison to the BH population

####  8.1.1 Low-mass X-ray binaries

Several attempts have been made to determine the Milky Way BH mass function
using dynamical mass measurements of BHs in low-mass X-ray binaries (LMXBs).
Using a sample of 7 LMXBs, Bailyn et al. (1998) found the BH mass function to
be tightly centered around $7M_{\odot}$ with a dearth of systems between
$3-5M_{\odot}$; they argued this dearth was not due to observational selection
effects. Later work by Özel et al. (2010) and Farr et al. (2011) using samples
of 15-16 LMXBs found similar results, providing further support to the idea of
the NS-BH mass gap. However, Kreidberg et al. (2012) cautioned that systematic
errors in the analysis of LMXB systems could push their inferred masses high,
artificially creating the mass gap, and Jonker et al. (2021) identified
potential observational biases that prevent measurement of high-mass LXMB
systems. Additionally, LXMBs occupy a very small and specific part of BH
evolutionary parameter space, and the BHs found in those systems are likely
not representative of the Galactic BH population as a whole.

####  8.1.2 Filling the NS-BH mass gap

Measurements of BH masses in non-LMXB systems do not show evidence of a mass
gap. Gravitational wave searches have found mass gap objects both as the
merger remnant ($\sim 3M_{\odot}$ in GW170917 (Abbott et al., 2017) and $\sim
3.4M_{\odot}$ in GW190425 (Abbott et al., 2020a)), and in the merger
components ($\sim 2.6M_{\odot}$ in GW190814 (Abbott et al., 2020b)). Non-
interacting mass gap BHs of $\sim 3M_{\odot}$ with red giant companions have
also been detected in the Milky Way (Thompson et al., 2019). With a mass of
1.6 - 4.4 $M_{\odot}$, the lens of OB110462 is the first measured isolated
Galactic NS or mass gap BH. These detections of mass-gap objects will improve
our theoretical understanding of compact object formation channels.

####  8.1.3 A lack of higher-mass BH systems?

In order to gain a full understanding of the Galactic BH population, BHs must
be uncovered outside of closely interacting X-ray binary systems.

We first consider searches for isolated BHs using microlensing, as discussed
in this paper. From our sample of 5 events, we have a single detection of a
mass gap object; all other lenses are lower mass non-BH detections. In
addition, we have a single detection of a $>1$ mas astrometric shift; most of
the remaining detections are at the $\sim 0.5$ mas level, near the limit of
HST’s precision. As discussed in §5.1 of Lam et al. (2022), the low yield of
BHs in this sample is consistent with predictions by Galactic models. However,
this presents tentative evidence that Galactic BHs may be less massive than
the $\sim 10M_{\odot}$ expectation. If the BH mass function truly peaks at
8$M_{\odot}$, then selecting candidates via long duration microlensing events
should doubly bias us toward finding these high-mass lenses. First, the
lensing cross section $\sigma=\pi\theta_{E}^{2}$ is proportional to the mass
of the lens $M_{L}$, since $\theta_{E}\propto\sqrt{M_{L}}$. Thus, more massive
objects are more likely to be microlenses. Second, the Einstein crossing time
is proportional to the square root of the lens mass,
$t_{E}\propto\sqrt{M_{L}}$. Thus, long duration events are also more likely to
be due to massive microlenses.

Next, we consider searches for BHs in detached/non-interacting binary systems.
The mass function $f$ of a single lined spectroscopic binary is given by

$f=\frac{P_{orb}K^{3}(1-e^{2})^{3/2}}{2\pi
G}=\frac{M_{2}^{3}\sin^{3}i}{(M_{1}+M_{2})^{2}}$ (36)

where $P_{orb}$ is the orbital period, $K$ is the radial velocity (RV) semi-
amplitude, $e$ is the orbital eccentricity, $i$ is the orbital inclination,
$M_{1}$ is the mass of the visible component, and $M_{2}$ is the mass of the
unseen component. If measurements of $P_{orb}$, $K$, $e$, and $M_{1}$ can be
obtained, then a minimum mass on $M_{2}$ can be derived. If $M_{2}>5M_{\odot}$
without evidence of luminosity, the unseen system is inferred to be a BH.
Since systems with larger RV semi-amplitudes have larger mass functions, they
are most likely to host unseen BH companions. To date, searches for large RV
semi-amplitudes in spectroscopic catalogs have detected an object which fall
within the mass gap, but no $\sim 10M_{\odot}$ BHs, which suggest a paucity of
higher-mass systems (Thompson et al. (2019); El-Badry et al. (2022)).
Complementary searches using ellipsoidal variables (Rowan et al., 2021) also
suggest that higher-mass systems are rare.

Although both microlensing and RV searches should be biased toward finding
$10M_{\odot}$ objects more easily than $3M_{\odot}$ objects, only the latter
are being detected. It may be that the selection bias for massive objects is
cancelled out by the fact that the mass function of stars, WDs, NS, and BHs
sharply decreases from low to high mass. Additional work to quantify and
compare these two competing effects, combined with larger sample sizes, will
be needed to understand the Galactic BH mass function.

Finally, astrometric searches for detached binaries are also eagerly
anticipated with Gaia (e.g. Yamaguchi et al. (2018); Yalinewich et al. (2018);
Wiktorowicz et al. (2020)). It will be very fruitful to compare the results of
those searches to the X-ray transient, microlensing, and RV searches.

###  8.2 Auxiliary OB110462 observations

Additional observations of OB110462 may assist in ascertaining its NS or BH
nature. Here we describe planned observations of OB110462, as well as data
found in searches of archival catalogs.

The 2021 astrometric observations of OB110462 are crucial to the modeling as
they extend the temporal baseline of the original archival observations by
50%, from 6 years to 10 years. The remaining data from HST Cycle 29 program
GO-16760 to be taken in Fall 2022 (Lam & Lu, 2021) will further extend the
baseline by another year and improve the characterization of the astrometric
signal. The astrometric microlensing shift is a deflection with respect to the
unlensed position of source (Equation 27), and we can only have confidence in
our measurement of that shift if we also have confidence in our measurement of
the source’s unlensed proper motion. The astrometric shift due to lensing when
$u\gg 1$ can be approximated as $\delta_{c}\approx\theta_{E}t_{E}/(t-t_{0})$
(see Appendix D). For OB110462 ($t_{E}=280$ days and $\theta_{E}=4$ mas), the
microlensing astrometric shift did not dip below HST’s astrometric precision
of $\sim 0.3$ mas until 2021, 10 years after source-lens closest approach.
This calculation illustrates the importance of having a long temporal baseline
for OB110462 in order to properly measure the source’s unlensed proper motion
and characterize the astrometric microlensing signal.

Additional follow-up observations in the X-ray can place limits on accretion
from the ISM (Agol & Kamionkowski, 2002). For example, Maeda et al. (2005) and
Nucita et al. (2006) looked for X-rays at the location of BH microlensing
candidate MACHO-96-BLG-5 reported in Bennett et al. (2002), using ACIS on
Chandra and EPIC on XMM-Newton, respectively. Neither detected any X-rays.

We searched several X-ray catalogs that have observed at OB110462’s
coordinates to determine whether there are any coincident sources. OB110462
was not detected as an X-ray source in any of the following catalogs:

  * •

Chandra Source Catalog 2.0191919https://cxc.harvard.edu/csc/ (Evans et al.,
2019a, b). The limiting sensitivity provides an upper limit of $1.91\times
10^{-14}$ erg/s/cm2 at 0.5-7.0 keV.

  * •

XMM-Newton Science Archive202020http://nxsa.esac.esa.int/nxsa-web/#search
(Sarmiento et al., 2019). This provides an upper limit of $1.52\times
10^{-14}$ erg/s/cm2 and $<1.38\times 10^{-3}$ counts/s at 0.2 - 12.0 keV.

  * •

Swift XRT Point Source Catalogue212121https://www.swift.ac.uk/2SXPS/ (Evans et
al., 2020). This provides an upper limit of $2.4\times 10^{-3}$ counts/s at
0.3 - 10 keV.

OB110462’s coordinates are not in the eROSITA-DE Early Data Release
catalog222222https://erosita.mpe.mpg.de/edr/.

We also searched the Australia Telescope National Facility Pulsar
Catalogue232323https://www.atnf.csiro.au/research/pulsar/psrcat/ (Manchester
et al., 2005) Version 1.65 for any pulsars coincident with the target. There
are no coincident pulsars; the nearest pulsar is 0.55 deg away.

As mentioned in §3.5, OB110462 is not in the Gaia EDR3 catalog. OB114062’s
baseline magnitude is F606W $\approx 22$, while Gaia’s nominal magnitude limit
is $G\approx 20.7$, and will be even brighter in crowded regions (Fabricius et
al., 2021).

###  8.3 Potential systematics in OB110462 observations

We briefly discuss whether the tension between the photometry and astrometry
could be due to systematics in the data.

First we consider the OGLE photometry. Possible sources of systematic error
include differential refraction, proper motion, or low-level stellar
variabliliy are possibilities for such a long baseline. However, no trends due
to proper motion or stellar variability are seen in the lightcurve. In
addition, the F814W $-$ F606W color of the bright star within $\sim 0.3$
arcsec of OB110462 that are blended together in the OGLE images are within 0.1
of each other, and hence effects due to differential refraction would be
undetectable. We also inspect the lightcurves of several stars near OB110462
that have similar magnitude to the baseline magnitude of OB110462. No trends
can be seen that resemble the residuals in the EW fit. In addition, we
explored re-scaling the photometric uncertainties of the OGLE data; however,
even inflating all uncertainties by a factor of $\sim 2$ did not significantly
change the structure of the residuals, nor the inferred microlensing
parameters, including lens mass.

Potential sources of systematics in the HST astrometry are discussed in the
next section.

###  8.4 Improving experimental strategy and design

####  8.4.1 Multi-filter astrometry

In contrast to photometric observations, multi-filter astrometry is not
routinely obtained, as astrometric observations are expensive and facilities
with the requisite precision are rare. This is one of the first projects to
explore the impact of different filters on relative astrometry. The nature of
the difference between F606W and F814W filter astrometric observations of
OB110037 and OB110462 is an open question, whether it be astrophysical (e.g.
binaries sources with different color), systematic (e.g. uncorrected CTE), or
statistical. However, it demonstrates that multi-filter observations are worth
continuing to pursue in future astrometric microlensing studies. For example,
it could help break degeneracies between certain types of binaries lenses or
sources.

####  8.4.2 CTE correction

As mentioned in §4.1, the CTE correction in the _flc files is not perfect.
Future pursuits will explore other methods of correction, such as a re-
analysis of OB110462 that uses the newer and more accurate tabular correction
for CTE (Anderson, 2021). In addition, trying to fix CTE via a magnitude-
dependent astrometric alignment is another avenue that is being explored.

####  8.4.3 Observational strategy

For these precise astrometric measurements, taking good observations is
critical. The lens mass constraints for several of the targets are only upper
limits, as the astrometric shifts were so small as to be undetectable at the
precision of the measurements. A dominant source of astrometric uncertainty
with HST WFC3-UVIS observations is the undersampling of the PSF. It has been
shown that there is a floor in the astrometric precision that can be achieved,
even at high SNR, when only a few exposures are used (Hosek et al., 2015). As
the majority of the observations in each filter had 4 or less exposures, this
limited the achievable astrometric precision for several of the targets, in
particular OB110462. Increasing the number of exposures and implementing small
uncorrelated dithers to sample different pixel phases can reduce this floor as
$\sqrt{N_{dithers}}$.

In addition, the effects of CTE are worsening with time. Actively mitigating
CTE through careful planning of observations rather than trying to correct it
afterwards is even more important than before (e.g. as described in §7 of
Anderson et al. (2021)).

####  8.4.4 Event selection

Although all events presented in this work were selected to have $t_{E}>200$
days (Sahu, 2009), the inferred $t_{E}$ values for 4 of the 5 events did not
satisfy this criteria. As a result, the true $t_{E}$ range probed extended
down to $t_{E}=60$ days, and did not sample the $t_{E}$ range that maximizes
the expected yield of BHs ($t_{E}\gtrsim 100$ days). Only MB09260 and OB110462
had $t_{E}>100$ days, weakening the constraints on the BH fraction. We are
currently attempting to determine whether prediction of $t_{E}$ before the
photometric peak of the event could be improved. If possible, this would
enable improved selection of BH candidates for astrometric follow-up.

A secondary concern is the target field itself. A sufficient number of
reference stars is needed, hence the field in the immediate vicinity of the
target must be sufficiently crowded. However, the magnitude range of those
nearby stars must also be similar to the target. Because of the steepness of
the luminosity function, bright targets or targets with high magnification are
more difficult to analyze as they lack sufficient reference stars to perform
relative astrometric alignment. This is in tension with the need to have high
photometric precision in order to precisely measure the microlensing parallax.
Although bright and highly magnified stars should still be followed up if they
are long duration events, special care must be taken when designing
observations to ensure good astrometric alignment.

###  8.5 BH searches with the Roman Space Telescope

Although the initial idea and subsequent design requirements for the Roman
Space Telescope (hereafter _Roman_) microlensing survey are driven by
exoplanet searches (Penny et al., 2019), it also hails the next generation of
astrometric microlensing campaigns searching for BHs. Presently, each event
must be followed up individually, with only two facilities (HST and Keck)
capable of the precision in the near-infrared required to make such a
measurement. Such measurements are expensive (requiring an $\sim$orbit of HST
or $\sim$hour of Keck time per measurement), prohibiting dense astrometric
temporal sampling or a large sample of targets.

_Roman_ will change this with its ability to simultaneously obtain precise
photometry and astrometry over a wide field of view $100\times$ the area of
HST and astrometric precision almost an order of magnitude better than HST
(Spergel et al., 2015; WFIRST Astrometry Working Group et al., 2019). This
will also allow the masses of NS-BH mass gap objects to be precisely measured,
and allow a sample of 100-1000 BH candidates to be built up over the duration
of the survey (Lam et al., 2020). In addition, _Roman_ will probe a large
sample of shorter $t_{E}$ events which will place constraints on BH kicks.

To make _Roman_ as effective as possible for finding BHs, there are several
considerations that need to be addressed. Due to the placement of the
observatory’s solar panels, which dictates the available pointings of the
telescope, the Galactic Bulge can only be observed twice a year during a 72
day window centered on the vernal and autumnal equinoxes. In addition, for
_Roman_ ’s planned 5 year mission, only 6 of the 10 available Bulge seasons
are to be dedicated to the microlensing survey. These large temporal gaps are
generally not a concern for exoplanets searches, as the transient portion of
the lightcurve is nearly covered within the 72 day window. However, for long
duration events where the transient portion of the lightcurve is much longer
than the window and where a measurement of small microlensing parallax is
crucial, incomplete lightcurve coverage will mean the difference between a
confirmed BH mass measurement and only an upper limit. Observations filling in
these gap will be crucial. Collaboration with a ground-based telescope to
provide imaging during the gaps, or a smaller independent follow-up efforts
would be very important.

## Acknowledgements

We thank Dan Foreman-Mackey, Tharindu Jayasinghe, Tom Loredo, Greg Martinez,
and Jeff Andrews for helpful and interesting conversations. In addition, we
thank the referee for feedback that improved this paper. We also thank Kailash
Sahu, Howard Bond, Jay Anderson, Martin Dominik, Philip Yock, and Annalisa
Calamida for proposing and taking the archival HST observations used in this
work.

C.Y.L. and J.R.L. acknowledge support by the National Science Foundation under
Grant No. 1909641 and the National Aeronautics and Space Administration (NASA)
under contract No. NNG16PJ26C issued through the WFIRST (now Roman) Science
Investigation Teams Program. C.Y.L. also acknowledges support from NASA
FINESST grant No. 80NSSC21K2043. D.P.B. was supported by NASA grants
NASA-80NSSC18K0274 and 80GSFC17M0002. Ł.W. acknowledges support from the
Polish National Science Centre (NCN) grants Harmonia No. 2018/30/M/ST9/00311
and Daina No. 2017/27/L/ST9/03221 as well as the European Union’s Horizon 2020
research and innovation programme under grant agreement No 101004719 (OPTICON-
RadioNet Pilot, ORP) and MNiSW grant DIR/WK/2018/12.

Based on observations made with the NASA/ESA Hubble Space Telescope, obtained
from the data archive at the Space Telescope Science Institute (STSci)
operated by the Association of Universities for Research in Astronomy, Inc.
under NASA contract NAS 5-26555, and obtained from the the Hubble Legacy
Archive, a collaboration between STScI/NASA, the Space Telescope European
Coordinating Facility (ST-ECF/ESA) and the Canadian Astronomy Data Centre
(CADC/NRC/CSA).

This paper makes use of data obtained by the MOA collaboration with the 1.8
metre MOA-II telescope at the University of Canterbury Mount John Observatory,
Lake Tekapo, New Zealand. The MOA collaboration is supported by JSPS KAKENHI
(Grant Number JSPS24253004, JSPS26247023, JSPS23340064, JSPS15H00781,
JP16H06287, and JP17H02871) and the Royal Society of New Zealand Marsden Fund.
The MOA project has received funding from the Royal Society of New Zealand,
grant MAU1901 to IB.

This work presents results from the European Space Agency (ESA) space mission
Gaia. Gaia data are being processed by the Gaia Data Processing and Analysis
Consortium (DPAC). Funding for the DPAC is provided by national institutions,
in particular the institutions participating in the Gaia MultiLateral
Agreement (MLA).

This research has made use of data obtained from the Chandra Source Catalog,
provided by the Chandra X-ray Center (CXC) as part of the Chandra Data
Archive.

##  Appendix A Rescaling of uncertainties

For each epoch, hst1pass returns the RMS error of extracted source positions
and magnitudes over multiple frames $\sigma_{x}$, $\sigma_{y}$ and
$\sigma_{m}$, respectively. Following Hosek et al. (2015), instead of using
the RMS errors for our uncertainties, we use the error on the mean
$\sigma/\sqrt{N}$ where $N$ is the number of frames the star is detected in,
inflated with a empirical additive error. The empirical constant additive
error on the positions $\Delta_{xy}$ and magnitudes $\Delta_{m}$ is calculated
for each epoch and added in quadrature to the error on the mean. This produces
the final rescaled uncertainties used in the analysis

$\displaystyle\sigma_{x}^{\prime}$
$\displaystyle=\sqrt{\sigma_{x}^{2}/N+\Delta_{xy}^{2}}$ (A1)
$\displaystyle\sigma_{y}^{\prime}$
$\displaystyle=\sqrt{\sigma_{y}^{2}/N+\Delta_{xy}^{2}}$ (A2)
$\displaystyle\sigma_{m}^{\prime}$
$\displaystyle=\sqrt{\sigma_{m}^{2}/N+\Delta_{m}^{2}}.$ (A3)

To determine $\Delta_{xy}$ and $\Delta_{m}$, a sample of bright, unsaturated
stars are selected. The exact magnitude range constituting “bright” is roughly
saturation to 3-5 magnitudes fainter, with the exact range determined
empirically through the astrometric alignment process. In this sample, the
additive error is selected such that the $\chi^{2}$ distribution of the
reference stars position and magnitude fits is roughly consistent with the
expected distribution (Figure 9). The resulting values are listed in Table A.
Note that for MB09260 F606W, the same additive magnitude error was used across
all epochs.

As mentioned in §4.2.1, for the five microlensing targets, the RMS and
rescaled astrometric uncertainties are generally similar. However, adopting an
additive error makes the resulting astrometric uncertainties more uniform
across the field, particularly for bright stars. This is critical as these
uncertainties are used as weights in the astrometric alignment (§4.2.3 and
Figure 9). The reference stars’ uncertainties are used to determine how good
the reference frame transformation is, which is ultimately used to measure the
astrometry of the target.

Table 16: Additive errors for HST data

MB10364 2010-09-13 0.14 1.61 8.6 11.6 2010-10-26 0.14 0.49 4.5 8.6 2011-10-31
0.17 0.26 5.9 5.0 2012-09-25 0.18 0.26 8.5 6.1 2013-10-24 0.17 0.25 6.8 5.8
OB110037 2011-08-15 0.11 0.22 16.0 8.5 2011-09-26 0.11 0.18 4.1 14.0
2011-11-01 0.10 0.21 4.0 3.8 2012-05-07 0.10 0.26 12.9 12.2 2012-09-25 0.10
0.22 6.7 4.3 2013-10-21 0.10 0.23 3.9 4.9 2014-10-26 0.10 0.21 4.2 9.7
2017-03-13 0.10 0.21 4.9 5.1 2017-09-04 0.11 0.21 8.5 19.6 OB110310 2011-09-21
0.10 0.20 9.3 12.9 2011-10-31 0.12 0.30 4.2 4.9 2012-04-24 0.13 0.30 10.0 24.3
2012-09-24 0.14 0.29 10.2 10.9 2013-10-21 0.13 0.29 4.3 6.3 2017-03-14 0.13
0.31 4.1 18.8 2017-09-01 0.14 0.32 10.0 24.9 OB110462 2011-08-08 0.13 0.41
18.1 17.0 2011-10-31 0.19 0.37 3.9 5.4 2012-09-09 0.19 0.36 4.7 7.2 2012-09-25
0.19 0.36 6.8 7.2 2013-10-22 0.19 0.36 3.6 4.8 2014-10-26 0.19 0.39 3.8 10.7
2017-08-29 0.20 0.40 5.4 14.7 2021-10-01 0.17 0.36 7.5 12.9

##  Appendix B Injection and recovery tests

We use the ks2 software (Anderson et al. (2008); Bellini et al. (2018), see
also §3 of Sabbi et al. (2016) for a detailed description) to inject
artificial stars into the OB110462 HST images to determine how well we can
recover the magnitudes and positions of faint sources near bright sources.
Injection and recovery tests are performed for all OB110462 epochs in two
different manners:

  1. 1.

Sources are injected radially around the bright neighbor star at the radius of
OB110462 at a variety of azimuths excluding the azimuth of OB110462 itself, as
we cannot recover a source planted on top of a real star.

  2. 2.

Sources are injected near a star of similar brightness to the neighbor, in the
same radial and azimuthal configuration relative to the star as OB110462
relative to the bright neighbor.

The first test directly probes the region around the neighbor itself but
excludes the actual position of the OB110462, while the second test probes a
region around a star similar to the neighbor, at the same separation and angle
of OB110462 relative to the neighbor.

The star we dub “the neighbor” is $\sim$10 pixels ($\sim$0.4 arcsec) west of
OB110462. The star we dub “the neighbor-like star” in F814W is $\sim$75 pixels
($\sim$3 arcsec) northeast of OB110462; for F606W the neighbor-like star is
$\sim$40 pixels ($\sim$1.6 arcsec) north of OB110462 (Figure 21).

Figure 21: Injection and recovery for OB110462. _Black circle, center_ :
OB110462. _White circle_ : the bright neighbor. _Blue circle_ : neighbor-like
star for F606W. _Cyan circle_ : neighbor-like star for F814W. Positions where
stars are injected in reference to the neighbor-like star are shown as blue
and cyan stars (F606W and F814W, respectively).

In summary, both methods produce similar results (Figures 3.1 and 3.1). We use
the single-azimuth results from injection around the neighbor-like star in
order to capture any azimuthal dependence which would be lost by averaging
over multiple azimuths (Table B.1).

In F814W, the bias is negligible in the first epoch, when the target is
magnified and is of similar brightness to the neighbor. However, in later
epochs once the target is no longer magnified, the measured position of the
target is biased by $\sim 0.4$ mas toward the bright neighbor star along the
target-neighbor separation vector. The magnitude and direction of the
positional offset are comparable across the two approaches, except in 2014,
where the injection around the neighbor-like star leads to a larger difference
than when injecting around the neighbor. The uncertainties are larger in the
approach of planting around the neighbor, as it averages over more azimuths
and results in a larger uncertainty.

In F606W, the bias is smaller than F814W because the shorter wavelength
results in higher resolution. Like F814W, the bias is also negligible in the
first epoch, and around $\sim 0.25$ mas in later epochs. The uncertainties are
also larger when averaging across many azimuths. Unlike F814W, the bias is
mixed between radial and azimuthal components when injecting around the
neighbor-like star at a single azimuth.

###  B.1 Injection around the neighbor

For each epoch, we determine the separation $r$ of OB110462 and the bright
neighbor. We inject three rings of stars surrounding the bright neighbor, of
radii $r-0.2$, $r$, and $r+0.2$ pixels. Each ring consists of 24 evenly spaced
stars, resulting in one star every 15 degrees. Because we cannot recover stars
injected on top of the target itself, we do not attempt to recover injected
stars that fall within 4 pixels of the target itself; this excludes three of
the positions. We thus inject a total of $(24-3)\times 3=63$ stars per epoch.

Table 17: Bias correction derived from injection and recovery

F814W 2011-08-08 -0.049 $\pm$ 0.000 -0.014 $\pm$ 0.000 -0.001 $\pm$ 0.000
2011-10-31 -0.203 $\pm$ 0.001 -0.070 $\pm$ 0.001 -0.004 $\pm$ 0.000 2012-09-09
-0.278 $\pm$ 0.001 -0.150 $\pm$ 0.001 -0.016 $\pm$ 0.000 2012-09-25 -0.292
$\pm$ 0.001 0.025 $\pm$ 0.002 -0.007 $\pm$ 0.000 2013-10-22 -0.490 $\pm$ 0.001
-0.081 $\pm$ 0.001 -0.011 $\pm$ 0.000 2014-10-26 -0.665 $\pm$ 0.001 -0.138
$\pm$ 0.001 -0.013 $\pm$ 0.000 2017-08-29 -0.258 $\pm$ 0.001 -0.129 $\pm$
0.001 -0.013 $\pm$ 0.000 2021-10-01 -0.315 $\pm$ 0.000 -0.001 $\pm$ 0.001
-0.007 $\pm$ 0.000

Note. — Bias correction derived from injection/recovery around a star of
comparable brightness at the same separation, azimuth, and magnitude
difference as the target to its bright neighbor.

###  B.2 Injection around the neighbor-like star

The neighbor-like star we inject around is different for F814W and F606W,
because the surrounding stars do not have the same colors as the neighbor and
target. The neighbor-like stars were chosen to have similar magnitude and
saturation level to the neighbor. In F814W, the neighbor tended to be
saturated; the F814W neighbor-like star is also saturated. On the other hand,
in F606W the neighbor was not saturated, and the F606W neighbor-like star is
also not saturated. The F814W neighbor-like star is brighter than the neighbor
in F814W by $\sim 0.6$ mag, and the F606W neighbor-like star is fainter than
the neighbor in F606W by $\sim 0.5$ mag.

In each epoch, we inject three arcs of radii $r-0.2$, $r$, and $r+0.2$ pixels
centered on the neighbor-like star, where $r$ is the target-neighbor
separation ($\sim 10$ pixels). Each arc consists of 15 stars at the azimuth of
the target relative to the neighbor $\pm$ 0.2 pixels/(target - neighbor
separation in pixels), which corresponds to a subtended angle of approximately
2.2 degrees. This corresponds to a region of $\sim$0.04 pix2 where $3\times
3\times 15=135$ stars are injected.

At each position we inject stars of magnitude $m_{I}-0.1$, $m_{I}$, and
$m_{I}+0.1$, where $m_{I}$ is the magnitude that results in the same contrast
with the neighbor-like star as OB110462 to the neighbor. That is, if OB110462
has magnitude $m_{T}$ and the neighbor $m_{N}$, and the injected star is
$m_{I}$ around the neighbor-like star $m_{C}$, then $m_{T}-m_{N}=m_{I}-m_{C}$.

###  B.3 Recovery of injected sources

After planting fake stars into the image, we determine how well we can recover
the positions and magnitudes. To match the properties of our original dataset,
we consider stars to be recovered if they are detected in at least $N$ frames,
where $N$ is the number of frames that were used to calculate the position of
the source. Iterative 3-sigma clipping is performed to exclude outliers due to
confusion, e.g. from the diffraction spike mask. We then use the
transformation parameters derived for this epoch (as described in §4.2.1) to
convert from (x,y) pixel positions to (RA, Dec) coordinates.

We define a polar coordinate system with the origin located at the neighbor
star for the analysis in §B.1 and located at the neighbor-like star for the
analysis in §B.2. The azimuthal direction is measured counterclockwise from
the origin-OB110462 separation vector. The average offset in the radial and
azimuthal directions are $r$ and $\theta$ (Figures 3.1 and 3.1).

The color of the neighbor and OB110462 are very similar (F606W - F814W = 2.25
and 2.15, respectively). At baseline, the neighbor is about 3.1 magnitudes
brighter than OB110462. During magnification in the first epoch (2011-08-08),
the neighbor is only about 0.4 magnitudes brighter. Since the resolution is
higher at shorter wavelengths, it is not unexpected that the positional bias
is less in F606W than F814W, since the separation is the same in both filters.
It is also not surprising that the bias is smallest in the first epoch when
the magnitude difference between OB110462 and the neighbor is small, and
larger in the remaining epochs when the magnitude difference is large.

Figure 22: Positional bias in F814W (_top panel_) and F606W (_bottom panel_)
as calculated from injection and recovery tests, injecting around a neighbor-
like star (_left column_) or around the neighbor (_right column_). The
positional bias (position recovered minus the true position input) is shown as
a function of radial $r$ and azimuthal $\theta$ components (_top row_), as
well as total positional bias (_bottom row_).

Figure 23: Magnitude bias in F814W (_top panel_) and F606W (_bottom panel_) as
calculated from injection and recovery tests, injecting around a neighbor-like
star (_left column_) or around the neighbor (_right column_). The magnitude
bias is defined as the magnitude recovered minus the true magnitude input.

##  Appendix C Absolute proper motion reference frame

The HST astrometry in §4.2 is derived in a procedure where the average motion
of the stars is at rest with respect to the reference frame. To interpret the
lens’ proper motion or transverse velocity, we must place the relative
astrometry into an absolute reference frame. To do this, we calculate the
proper motion offset between stars in the relative HST frame and the absolute
Gaia frame. We match all stars in the Gaia EDR3 catalog within 1 arcmin of the
target that have astrometric_excess_noise_sig $<2$ and a proper motion
measurement (i.e. 5-parameter and 6-parameter solutions) to the bright stars
in our HST proper motion catalog (F814W $<$ 22 for MB10364 and OB110037; F814W
$<$ 23 for MB09260, OB110310, and OB110462). The 1-iteration $3\sigma$-clipped
uncertainty-weighted average difference in the proper motion between the
cross-matched stars is calculated, then applied to the relative HST astrometry
to place it into the absolute Gaia proper motion frame. The values to convert
between the HST and Gaia frames for each target’s field are listed in Table
3.1.

The vector point diagram of proper motion differences between cross-matched
sources in the Gaia and HST F814W catalogs, after applying the proper motion
offset to place the HST catalog into the Gaia reference frame, is shown in
Figure 3.1. In general, the proper motions of bright stars in Gaia are
inconsistent with those derived using HST. For fainter stars, the
uncertainties are much larger, so there is more consistency between Gaia and
HST; however, there is substantial scatter between the measurements of the two
catalogs. The discrepancies between the HST and Gaia proper motions could
indicate that the uncertainties are underestimated in one or both catalogs, or
that there are higher-order distortions between the two reference frames that
cannot be captured by a constant offset. However, the most likely explanation
is that the Gaia proper motions are not accurate, as it is clear that the Gaia
EDR3 astrometry is unreliable in crowded regions like the Galactic Bulge (see
Appendix E and Rybizki et al. (2022)).

Figure 24: Vector point diagram of proper motion differences between cross-
matched sources in the Gaia and HST F814W catalogs within $60^{\prime\prime}$
of the targets, after applying the constant proper motion offset derived in
Appendix C to place the HST catalog into the Gaia reference frame. The color
of the point denotes the cross-matched source’s magnitude in Gaia G-band.
MB10364, OB110037, and OB110310 are sources in Gaia; they are shown as the
black-outlined squares and error bars.

Table 18: HST to Gaia Proper Motion Reference Frame Offset

Note. — $\Delta\mu_{\alpha*}=\mu_{{\rm Gaia},\alpha*}-\mu_{{\rm
HST\;F814W},\alpha*}$ and $\Delta\mu_{\delta}=\mu_{{\rm
Gaia},\delta}-\mu_{{\rm HST\;F814W},\delta}$.

##  Appendix D Astrometric lensing in Gaia

As discussed in §4.2.6, there are several differences between the proper
motion modeling of microlensing targets in this paper as compared to Gaia.
First, we simultaneously model the proper motions of the lens and source with
parallax within a microlensing model. In contrast, Gaia models the proper
motion of a single star with parallax. As the source and lens are not resolved
by Gaia during the observations, if the lens is luminous, the proper motion of
the source and lens are being treated as a single object. Depending on the
relative source-lens motion and colors, this can produce different inferred
proper motions. Second, the temporal baseline of the data used to calculate
the proper motions differ. The HST time baseline spans 2010-2013 for MB10364
and 2011-2017 for OB110037 and OB110310, in contrast to 2014-2017 in Gaia. As
the astrometric lensing signal is time dependent, temporal baseline
differences can also lead to different inferred proper motions.

With regard to the second point, we estimate the amount by which astrometric
microlensing could affect the proper motion measurement. The lens-source
separation $\boldsymbol{u}(t)$ (Equation 20) sets the time evolution of the
astrometric shift $\delta_{c}(t)$ (Equation 27). Three terms contribute to
$\boldsymbol{u}(t)$: the impact parameter $\boldsymbol{u}_{0}$, the source-
lens proper motion $(t-t_{0})\boldsymbol{\hat{\mu}}_{rel}/t_{E}$, and the
microlensing parallax $\pi_{E}\boldsymbol{P}(t)$. As the Earth’s orbit around
the Sun in nearly circular, $|\boldsymbol{P}(t)|\approx 1$ to within 3%. For
typical microlensing events like MB10364, OB110037, and OB110310,
$|\boldsymbol{u}_{0}|<1$ and $\pi_{E}<1$. Thus, long after source-lens closest