[a]Hao Chen

# The pole counting rule and XYZ states

Hong-Rong Qi

###### Abstract

The pole counting rule is a powerful and model-independent method to
distinguish a confining state from a hadronic molecule. It has been applied to
the explorations of $X(6900)$, $X_{1}(2900)$ as well as $Z_{c}(3900)$,
$X(3872)$, $X(4660)$, etc.. For $X(6900)$, both a confining state and a
molecular state are not excluded, because lacking of enough data. For
$X_{1}(2900)$, the analysis shows that it should be a $\bar{D}_{1}K$ molecule,
with $J^{P}=1^{-}$ and an iso-singlet interpretation is much more favorable.
Finally, it is noted that almost all XYZ particles with exotic quantum numbers
can be interpreted as hadronic molecules. The $X(3872)$ is, however, more like
a charmonium.

##  1 Introduction

In the past 20 years, many exotic hadrons (also known as “XYZ" states) have
been observed. These states cannot be explained by ordinary quark model for
mesons and baryons which have valence quarks of $\bar{q}q$ and $qqq$,
respectively. So, a natural idea could be that these states have multiquark
constituents, i.e., they may have valence quarks of $\bar{q}\bar{q}qq$,
$\bar{q}qqqq$. An explanation for these multiquark states is that they are
compact or elementary states, whose building blocks are quarks and anti-quarks
and they interact with each other by exchanging gluons. However, there is
another explanation that these XYZ states are molecules of ordinary hadrons.
From this point of view, the building block of these states are ordinary
mesons or baryons and they interact with each other by exchanging color
neutral forces. On the one hand, multiquark states are not forbidden by our
present understanding of QCD, the underlying theory of strong interaction and
this indicates that the compact state explanation is possible. On the other
hand, these states are always observed near the thresholds of two-hadron
channels which indicates that the hadronic molecule explanation is also
possible. So, it is important to distinguish hadronic molecules from compact
states and this has become one of the most active topics in hadron physics.

Many methods are proposed to do this task, for example, QCD sum rule, quasi-BS
equation, diquark model, etc [1, 2]. In this paper, the pole counting rule
(PCR), a model-independent method will be introduced and some applications on
XYZ states such as $X(3872),~{}Z_{c}(3900),~{}X(4660)$ as well as
$X(6900),~{}X_{1}(2900)$ will be reviewed. For $X(3872)$, the result from PCR
suggests that it is more like a charmonium, i.e., not a hadronic molecule [4,
3]. The same situation also happens on $X(4660)$ [5]. As for $Zc(3900)$ and
$X_{1}(2900)$, they are regarded as hadronic molecules of $\bar{D}^{*}D$ [6,
7] and $\bar{D}_{1}K$ [8], respectively. Furthermore, due to lack of
experimental data on $X(6900)$, the fit results which are regarded as compact
states as well as the results which are regarded as hadronic molecules
according to PCR have the same fit goodness. So, one cannot confirm the nature
of $X(6900)$ from current data [9].

This talk will be organized as follows, firstly, the statements of PCR will be
briefly reviewed. Then, the applications of PCR on
$X(3872),~{}Z_{c}(3900),~{}X(4660)$ as well as $X(6900),~{}X_{1}(2900)$ will
be introduced. Finally, some conclusions will be made.

##  2 The pole counting rule

In this section, the statements of the pole counting rule will be introduced
briefly. This method was proposed by D. Morgan in 1992 [10]. Consider an
$S$-wave scattering amplitude:

$T=(M-ik)^{-1},$ (1)

where $k$ is the momentum of scattering particles, and the pole position is
determined by $M-ik=0$. If it is a potential scattering, then $M$ can be
written in the form of effective range expansion:

$M=-\frac{1}{a_{0}}+\frac{r_{0}}{2}k^{2}+O\left(k^{4}\right),$ (2)

where $a_{0}$ is the scattering length and $r_{0}$ is the effective range.
Substitute (2) into the equation $M-ik=0$, a quadratic equation with
respective to pole position is obtained. Two roots of this equation satisfy:

$\left|k_{1}\right|+\left|k_{2}\right|\geq\frac{2}{\left|r_{0}\right|},$ (3)

where effective range of strong interaction usually takes the value of about
(200 MeV)-1, so the r.h.s. of the inequality (3) is a quite large quantity. On
the one hand, in order to satisfy this inequality, there is only one pole can
appear around the threshold and another pole should located far enough from
the threshold to make the l.h.s. of (3) big enough. On the other hand, the
potential scattering scenario means that the pole is dynamically generated,
i.e., it represents that the structure should be regarded as a hadronic
molecule.

Moreover, if the pole hardly couples to the hadronic channel in which the
structure is observed, then the pole will behave more like a CDD pole in the
N/D method. In other words, the pole is more "elementary" and cannot be
generated from hadron interactions. It should be regarded as a compact state.
In this situation, $M$ can be written as,

$M=\frac{k^{2}-k_{0}^{2}}{g^{2}},$ (4)

two roots satisfy,

$\left|k_{1}\right|+\left|k_{2}\right|\geq g^{2}.$ (5)

If the coupling constant $g$ is small, two poles can be found around the
threshold in this case.

So, the pole counting rule provides a model-independent method to distinguish
hadronic molecules from compact states. In the $S$-wave scattering, if there
is only one pole near the threshold, then the scattering state should be
regarded as a hadronic molecule. If there are two poles around the threshold,
the scattering state should be regarded as a compact state. Usually, the poles
will be searched for in the complex energy plane, but the criteria of
molecular states and compact states are the same. It is noted that if there
are more than one couple channels, the poles will have imaginary parts since
the resonance can decay into lower channels and the pole positions will move
to the complex energy plane. Furthermore, the PCR is still a qualitative
method, and the efforts for a quantitative version can be found in [11].

##  3 Applications

In this section, the applications of PCR on XYZ states such as
$X(3872),~{}Z_{c}(3900),~{}X(4660)$ as well as $X(6900),~{}X_{1}(2900)$ will
be introduced.

###  3.1 $X(3872)$

$X(3872)$ was firstly observed in 2003 in the $J/\psi\pi^{+}\pi^{-}$ final
state [12] and then it was observed in $\bar{D}^{0}D^{0}\pi^{0}$ final state
and near the threshold of $\bar{D}^{0*}D^{0}$ [13]. So, it is natural to think
$X(3872)$ to be a molecule of $\bar{D}^{*}D$ or a compact state, e.g., a
charmonium or a tetraquark state.

In order to study the nature of $X(3872)$, Ref. [4] considered the coupling
between $X(3872)$ and $\bar{D}^{*0}D^{0},~{}D^{*-}D^{+}$, as well as
$J/\psi\pi^{+}\pi^{-},~{}J/\psi\pi^{+}\pi^{-}\pi^{0}$. Using Flatté-like
parametrization, the data sets corresponding to the final states of
$\bar{D}^{*0}D^{0}$, $\bar{D}^{0}D^{0}\pi^{0}$, $J/\psi\pi^{+}\pi^{-}$ were
fitted. During the fitting, the parameter $\Gamma_{c}$ which represents the
decay width of some modes far away from the threshold such as
$\psi^{\prime}\gamma$ can provide important information on the nature of
$X(3872)$. If $X(3872)$ is a charmonium, then its most important hidden-charm
decay mode is the inclusive light hadronic decay. In this case,
$\Gamma_{c}\sim\mathcal{O}$(1) MeV [14]. But if $X(3872)$ is a
$\bar{D}^{*0}D^{0}$ molecule, the most important hidden-charm decay mode is
decaying to $\chi_{cJ}\pi\pi$. In this case, $\Gamma_{c}$ is expected to be
smaller than $\Gamma_{J/\psi\pi\pi}\sim$0.05 MeV [15].

Besides $\Gamma_{c}$, the parameter $\mathcal{B}$ representing the branch
ratio of $B\to X(3872)K$ is also important. From the fit results, it is found
that if $\mathcal{B}$ takes a reasonable value, there always exist two poles
near the threshold and $\Gamma_{c}$ will be at $\mathcal{O}(1)$ MeV which
supports the charmonium statement. However, the fit result corresponding to
pure molecular state where there is only one pole near the threshold can also
be found. But in this case, the value of $\mathcal{B}$ is too large to be
realistic and $\Gamma_{c}$ is of $\mathcal{O}(1)$ MeV which is inconsistent
with the analyses in Ref. [15]. So, with the help of PCR, it is suggested in
Ref. [4] that $X(3872)$ is more like a charmonium rather than a pure hadronic
molecule, but the coupling to $\bar{D}^{*0}D^{0}$ is also important, since the
pole of $X(3872)$ is located around the threshold of $\bar{D}^{*0}D^{0}$.

Furthermore, in Ref. [3] the dynamically generated $X(3872)$ which is produced
by $\bar{D}^{*}D$ bubble chain and the explicitly introduced $X(3872)$ as well
as the competition between these two schemes, were under consideration. It is
found that a pure explicitly introduced $X(3872)$ (called Fit I in Ref. [3])
can fit the data better than a pure dynamically generate $X(3872)$ (called Fit
II in Ref. [3]) as shown in Fig. 1. As the results, Fit I gives more than one
pole around the threshold of $\bar{D}^{*}D$, which indicates that $X(3872)$ is
more like a charmonium in accordance with the PCR.

Figure 1: Fitting results of two methods in Ref. [3].

###  3.2 $Z_{c}(3900)$

$Z_{c}(3900)$ is another well known exotic state besides $X(3872)$. It was
firstly observed in the $J/\psi\pi$ final state by BESIII [16], with the mass
$(3899.0\pm 3.6)$ MeV and the width $(46\pm 10)$ MeV determined by Breit-
Wigner fit. Soon after, $Z_{c}(3900)$ was also observed in the $\bar{D}^{*}D$
final state [17] where the peak is very close to the threshold. So, a natural
assumption is regarding it as the hadronic molecule of $\bar{D}^{*}D$, but the
scenario that $Z_{c}(3900)$ is a compact state, e.g. a tetraquark state cannot
be excluded.

Ref. [6] constructed the amplitudes of processes $e^{+}e^{-}\to
J/\psi\pi\pi,~{}\bar{D}D^{*}\pi,~{}h_{c}\pi\pi$ and fitted the invariant mass
spectrum of $\pi\pi,~{}J/\psi\pi,~{}\bar{D}^{*}D$ simultaneously. Two
mechanisms are involved, one is dynamically generated by $\bar{D}^{*}D$
contact interaction (corresponding to Fit I) and another is exchanging
explicit $Z_{c}$ in $s$ channel (corresponding to Fit II). These two
mechanisms are quite different, but they can lead to the same result as shown
in Ref. [6]. Only one pole can be found around the threshold of $\bar{D}^{*}D$
in both Fit I and Fit II. In Fit I, the pole is located at
$(3.880-0.004i)~{}$GeV and in Fit II, the pole is located at
$(3.979-0.001i)~{}$GeV. According to the PCR, it can be concluded that
$Z_{c}(3900)$ is a hadronic molecule of $\bar{D}^{*}D$.

###  3.3 $X(4660)$

The PCR was also applied to study the properties of $X(4660)$. $X(4660)$ was
firstly observed by Belle [18] in 2007. Later, it was observed in the final
state $\Lambda_{c}\bar{\Lambda}_{c}$ [19] and the peak is near the threshold
of $\Lambda_{c}\bar{\Lambda}_{c}$. However, in 2017, BESIII released new data
on $X(4660)$ and it was shown that the new data conflicted with previous data.

As shown in Ref. [5], this conflict between the old and the new data can be
explained by a virtual pole near the $\Lambda_{c}\bar{\Lambda}_{c}$ threshold
rather than the threshold enhancement of electromagnetic force. The
$\Lambda_{c}\bar{\Lambda}_{c}$ contact interactions, as well as a $X(4660)$
Breit-Wigner resonance are combined to explain the odd line-shape. The cross
sections of $e^{+}e^{-}\to\Lambda_{c}\bar{\Lambda}_{c}$ and
$e^{+}e^{-}\to\psi^{\prime}\pi\pi$ are fitted. A K-matrix formula which mixes
the contact interaction and $s$ channel $X(4660)$ exchange is used to describe
$\Lambda_{c}\bar{\Lambda}_{c}$ final state interactions :

$A_{K}=\frac{1}{1-i\rho K}=\frac{1}{1-i\rho(T_{X}+T_{\text{ci}})},$ (6)

where $\rho$ is the two-body phase space factor of
$\Lambda_{c}\bar{\Lambda}_{c}$ and $T_{X}$, $T_{\text{ci}}$ are the tree
diagrams of $s$ channel $\Lambda_{c}\bar{\Lambda}_{c}\to
X(4660)\to\Lambda_{c}\bar{\Lambda}_{c}$ and $\Lambda_{c}\bar{\Lambda}_{c}$
contact vertex as final states interactions. From the pole analyses, there is
a stable virtual pole with mass of $M_{V}=(4.566\pm 0.003)$ GeV. Furthermore,
$X(4660)$ is a compact state, if exists, since there are two nearby poles.

###  3.4 $X(6900)$ and $X(4260)$

In 2020, LHCb observed a structure around 6.9 GeV, named as $X(6900)$, in the
di-$J/\psi$ invariant mass spectrum [21]. This state is the first observed
state which has four (anti-) charm quarks as valence quarks. Since it is
difficult to produce $\bar{c}c$ from the vaccum, so $X(6900)$ is hardly a
excited charmonium. For this reason, most studies suggested that $X(6900)$ is
a tetraquark state [22, 23] or a hadronic molecule [24, 25]. In Ref. [9],
Flatté-like parametrization was involved to describe the structure around 6.9
GeV as well as a higher structure around 7.2 GeV, called $X(7200)$. Using this
parametrization, the di-$J\psi$ invariant mass spectrum was fitted and the
nature of these structures were analysed according to the PCR.

Three types of fits are carried out in Ref. [9] called S-S coupling, P-P
coupling, S-P (or P-S) coupling, where the first and the second character
represent the angular momentum quantum number of channels coupled to $X(6900)$
and $X(7200)$, respectively. From the results, it is found that both the
solutions of compact states and the solutions of molecular states can be found
and these two types of solutions have almost the same fit goodnesses (an
example is shown in Fig. 2). As an example, the pole positions of S-S
couplings are listed in Tab. 1. Here, the solution I stands for compact state
and solution II stands for hadronic molecule. Three cases for the $X(6900)$
coupling are under consideration: Case I: $J/\psi J/\psi$ and
$J/\psi\psi(3770)$, Case II: $J/\psi J/\psi$ and $J/\psi\psi_{2}(3823)$, Case
III:$J/\psi J/\psi$ and $J/\psi\psi_{3}(3842)$. For $X(7200)$, the $J/\psi
J/\psi$ and the near-threshold $J/\psi$ $\psi(4160)$ channels are used in the
couple channel analysis.

Figure 2: Fit results of case I in S-S coupling in Ref. [9]: (a): compact state, (b): molecular state.  | Case | State | Sheet II | Sheet III
---|---|---|---|---
Sol. I | I | $X(6900)$ | $6885.4-68.0i$ | $6874.4-80.0i$
$X(7200)$ | $7202.2-16.6i$ | $7187.1-18.0i$
II | $X(6900)$ | $6947.6-172.0i$ | $6810.4-274.0i$
$X(7200)$ | $7220.8-31.0i$ | $7220.8-31.0i$
III | $X(6900)$ | $6845.2-117.0i$ | $6789.2-138.0i$
$X(7200)$ | $7221.9-28.0i$ | $7221.9-28.0i$
Sol. II | I | $X(6900)$ | $6937.9-97.0i$ | $6527.3-323.0i$
$X(7200)$ | $7210.7-27.5i$ | $7037.3-47.5i$
II | $X(6900)$ | $6933.9-111.0i$ | $6443.8-275.0i$
$X(7200)$ | $7218.9-24.0i$ | $7067.9-41.5i$
III | $X(6900)$ | $6933.3-113.0i$ | $6452.3-275.0i$
$X(7200)$ | $7221.9-23.0i$ | $7073.7-41.0i$
Table 1: Summary of pole positions (written in MeV) for $S-S$ couplings. In
practice, it is adopted that if a pole is more than three times of line-shape
width away from the threshold, it can be seen as a faraway pole. So, for every
case, there is only one nearby pole in solution II. The Riemann sheets are
defined by changing signs of the phase space factors for the first and the
second couple channels where $(+,+),~{}(-,+),~{}(-,-),~{}(+,-)$ stand for the
sheet I, II, III, IV, respectively.

According to the pole analysis in Ref. [9], the nature of $X(6900)$ as well as
$X(7200)$ cannot be determined using current data. So, more data are required
to study their properties better.

It is worthwhile to mention that in Refs. [26, 27], similar methods are
applied to study $X(4260)$, another exotic state. The couplings to two nearby
channels, $\omega\chi_{c0},~{}D_{s}^{*}\bar{D}_{s}^{*}$ are carefully
considered. The results shows that with a reasonable production widths of
$e^{+}e^{-}\to X(4260)$, the couplings of $X(4260)$ to these channels are not
strong enough to form a hadronic molecule, i.e., $X(4260)$ is more like a
compact state.

###  3.5 $X_{1}(2900)$

In 2020, LHCb observed a spin-1 resonance around 2.9 GeV in the $D^{-}K^{+}$
invariant mass spectrum, labeled as $X_{1}(2900)$ [28]. This resonance is near
the threshold of $\bar{D}_{1}K$ and $\bar{D}^{*}K^{*}$, so it is more likely
to be a hadronic molecule of these two hardron pairs. However, $X_{1}(2900)$
has $J^{P}=1^{-}$ given by the experiment, and an $S$-wave $\bar{D}^{*}K^{*}$
molecule and $\bar{D}_{1}K$ molecule has $J^{P}=1^{+},~{}1^{-}$, respectively.
So, if $X_{1}(2900)$ is a hadronic molecule, it is more likely to be a
molecule of $\bar{D}_{1}K$.

In Ref. [8], two schemes of the production of $X_{1}(2900)$ are taken into
consideration. One is to produce $X_{1}(2900)$ dynamically and another is to
introduce $X_{1}(2900)$ explicitly. In these two schemes, channel $\bar{D}K$
and $\bar{D}_{1}K$ are involved. The dynamically generated $X_{1}(2900)$ is
achieved by the couple channel K-matrix formula,

$T=\mathcal{K}\cdot[1-i\rho(s)\mathcal{K}]^{-1},~{}\mathcal{K}=\left(\begin{array}[]{cc}\mathcal{T}_{\bar{D}K\to\bar{D}K}&\mathcal{T}_{\bar{D}K\to\bar{D}_{1}K}+\mathcal{P}_{12}\\\
\mathcal{T}_{\bar{D}_{1}K\to\bar{D}K}+\mathcal{P}_{12}&\mathcal{T}_{\bar{D}_{1}K\to\bar{D}_{1}K}\end{array}\right)^{IJ},$
(7)

where $\rho=\mathrm{diag}\\{\rho_{\bar{D}K},\rho_{\bar{D}_{1}K}\\}$, is the
diagnoal matrix of phase space factors, $\mathcal{T}$s are the tree level
diagrams with special isospin and angular momentum, $\mathcal{P}_{12}$ is a
first order polynomial of c.m. energy to simulate the effects from higher
order interactions such as vector meson exchanging diagrams in physical
region. The explicitly introduced $X_{1}(2900)$ is achieved by the Flatté-like
parametrization.

The isospin symmetry is considered in the method of couple channel K-matrix,
however, the values of parameter $\zeta_{1}$ which is involved in the tree-
level amplitude $\mathcal{T}_{\bar{D}K\to\bar{D}_{1}K}$ are quite different in
iso-singlet and iso-triplet solutions. This parameter comes from the
interaction lagrangian [29],

$\mathcal{L}=-i\zeta_{1}\left\langle\bar{H}_{a}^{(\bar{Q})}\left(\mathcal{V}_{\mu}-\rho_{\mu}\right)_{ab}T_{b}^{(\bar{Q})\mu}\right\rangle+\mathrm{h.c.},$
(8)

where $\bar{H}^{(\bar{Q})}$ is the field operator with $J^{P}=0^{-}$ such as
$\bar{D}$ and $T^{(\bar{Q})}$ is the field operator with $J^{P}=1^{+}$ such as
$\bar{D}_{1}$. $\rho_{\mu}$ is the vector meson octet. So, in principle, the
parameter $\zeta_{1}$ should be estimated from the decay $D_{1}\to
D+\rho(\omega)$, but this process is forbidden kineically. In Ref. [30], it is
estimated from the process $K_{1}\to K\rho$ and takes the value of about $\pm
0.16$ . Since the corresponding values of $\zeta_{1}$ given by the iso-triplet
solution and iso-singlet solution are $5.85\pm 2.00$ and $-0.9\pm 0.48$, so it
indicates that $X_{1}(2900)$ is an iso-singlet.

Figure 3: Fit results in Ref. [8] for (a): dynamically generated $X_{1}(2900)$
with $IJ^{P}=01^{-}$, (b): explicitly introduced $X_{1}(2900)$. These two fits
have almost the same goodness and only one pole is found near the threshold of
$\bar{D}_{1}K$ in both cases.

In both two schemes, there is only one pole found near the threshold of
$\bar{D}_{1}K$. For dynamically generated $X_{1}(2900)$, the pole is located
at $(2.928-0.034i)$ GeV and for explicitly introduced $X_{1}(2900)$, the pole
is located at $(2.910-0.039i)$ GeV. According to the PCR, this suggests that
$X_{1}(2900)$ is a hadronic molecule of $\bar{D}_{1}K$ and at this time, two
different schemes lead to the same conclusion, just like the study of
$Z_{c}(3900)$ in Ref. [7].

##  4 Summary and conclusion

In this talk, it is showen that PCR provides a model-independent method to
distinguish hadronic molecules from compact states. But it still has a
disadvantage that only $S$-wave amplitudes can be used as inputs.
Nevertheless, almost all the known molecular states are $S$-wave couplings, so
this method is well practical.

In this talk, the applications of PCR in studying some XYZ states such as
$X(3872),~{}Z_{c}(3900),$ $X(4660),~{}X(6900),~{}X_{1}(2900)$ are reviewed.
With the help of PCR, $X(3872)$ is regarded as charmoniums, i.e., compact
states and $Z_{c}(3900)$, $X_{1}(2900)$ are suggested to be hadronic molecules
of $\bar{D}^{*}D$ and $\bar{D}_{1}K$, respectively. For $X(6900)$, its nature
cannot be well determined using current data, so more experimental information
are required. Furthermore, it is mentionable that in some works, molecule
solutions for $X(3872)$ can be found, e.g., in Ref. [31] where it is argued
that $\bar{D}^{*}D$ molecule. However, the analysis for the decay width to
lower channels $\Gamma_{c}$ and the $B$ decay branch ratio $\mathcal{B}(B\to
X(3872)K)$ as shown in Refs. [3, 4] suggest that $X(3872)$ is more like a
charmonium despite the strong coupling to $\bar{D}^{*}D$.

From the applications listed in this paper, it is suggested that all the XYZ
states which have exotic quantum number are more likely to be hadronic
molecules and the state like $X(3872)$, is compact states since they have
$\bar{c}c$ quantum numbers. For more discussion about methods in studying
exotic hadron states, one is referred to the review [32].

## Acknowledgements

The authors would like to thank Han-Qing Zheng and Qin-Fang Cao for helpful
discussions. This work is supported in part by National Nature Science
Foundations of China (NSFC) under contract numbers 11975028 and 10925522.

## References

  * [1] A. Hosaka, T. Iijima, K. Miyabayashi, Y. Sakai and S. Yasui, PTEP 2016, 062C01 (2016)
  * [2] N. Brambilla, S. Eidelman, C. Hanhart, A. Nefediev, C. P. Shen, C. E. Thomas, A. Vairo and C. Z. Yuan, Phys. Rept. 873, 1 (2020)
  * [3] C. Meng, J. J. Sanz-Cillero, M. Shi, D. L. Yao and H. Q. Zheng, Phys. Rev. D 92, 034020 (2015)
  * [4] O. Zhang, C. Meng and H. Q. Zheng, Phys. Lett. B 680, 453 (2009)
  * [5] Q. F. Cao, H. R. Qi, Y. F. Wang and H. Q. Zheng, Phys. Rev. D 100, 054040 (2019)
  * [6] Q. R. Gong, J. L. Pang, Y. F. Wang and H. Q. Zheng, Eur. Phys. J. C 78, 276 (2018)
  * [7] Q. R. Gong, Z. H. Guo, C. Meng, G. Y. Tang, Y. F. Wang and H. Q. Zheng, Phys. Rev. D 94, 114019 (2016)
  * [8] H. Chen, H. R. Qi and H. Q. Zheng, Eur. Phys. J. C 81, 812 (2021)
  * [9] Q. F. Cao, H. Chen, H. R. Qi and H. Q. Zheng, Chin. Phys. C 45, 103102 (2021)
  * [10] D. Morgan, Nucl. Phys. A 543, 632 (1992)
  * [11] Z. Xiao and Z. Y. Zhou, J. Math. Phys. 58 (2017) 072102
  * [12] S. K. Choi et al. [Belle], Phys. Rev. Lett. 91, 262001 (2003)
  * [13] B. Aubert et al. [BaBar], Phys. Rev. D 77, 011102 (2008)
  * [14] C. Meng and K. T. Chao, Phys. Rev. D 75, 114002 (2007)
  * [15] S. Fleming and T. Mehen, Phys. Rev. D 78, 094019 (2008)
  * [16] M. Ablikim et al. [BESIII], Phys. Rev. Lett. 110, 252001 (2013)
  * [17] M. Ablikim et al. [BESIII], Phys. Rev. Lett. 112, 022001 (2014)
  * [18] X. L. Wang et al. [Belle], Phys. Rev. Lett. 99, 142002 (2007)
  * [19] G. Pakhlova et al. [Belle], Phys. Rev. Lett. 101, 172001 (2008)
  * [20] M. Ablikim et al. [BESIII], Phys. Rev. Lett. 120, 132001 (2018)
  * [21] R. Aaij et al. [LHCb], Sci. Bull. 65, 1983 (2020)
  * [22] H. X. Chen, W. Chen, X. Liu and S. L. Zhu, Sci. Bull. 65, 1994 (2020)
  * [23] C. Deng, H. Chen and J. Ping, Phys. Rev. D 103, 014001 (2021)
  * [24] C. Gong, M. C. Du, Q. Zhao, X. H. Zhong and B. Zhou, Phys. Lett. B 824, 136794 (2022)
  * [25] Z. R. Liang, X. Y. Wu and D. L. Yao, Phys. Rev. D 104, 034034 (2021)
  * [26] L. Y. Dai, M. Shi, G. Y. Tang and H. Q. Zheng, Phys. Rev. D 92, 014020 (2015)
  * [27] Q. F. Cao, H. R. Qi, G. Y. Tang, Y. F. Xue and H. Q. Zheng, Eur. Phys. J. C 81, 83 (2021)
  * [28] R. Aaij et al. [LHCb], Phys. Rev. D 102, 112003 (2020)
  * [29] G. J. Ding, Phys. Rev. D 79, 014001 (2009)
  * [30] X. K. Dong, Y. H. Lin and B. S. Zou, Phys. Rev. D 101, 076003 (2020)
  * [31] F. K. Guo, C. Hanhart, U. G. Meißner, Q. Wang, Q. Zhao and B. S. Zou, Rev. Mod. Phys. 90, 015004 (2018)
  * [32] D. L. Yao, L. Y. Dai, H. Q. Zheng and Z. Y. Zhou, Rept. Prog. Phys. 84, 076201 (2021)


approach, i.e. $(t-t_{0})/t_{E}\gg 1$, the contribution by the impact
parameter and microlensing parallax terms to $\boldsymbol{u}(t)$ are
subdominant to that of the source-lens relative proper motion, and the source-
lens separation can be approximated

$u(t)\approx\frac{t-t_{0}}{t_{E}}\gg 1.$ (D1)

The astrometric shift for large source-lens separations $u\gg 1$ can be
approximated as

$\delta_{c}(t)\approx\frac{\theta_{E}}{u(t)}.$ (D2)

Putting these two equations together, we can approximate how much a proper
motion measurement based on observations made at times $t_{2}$ and $t_{1}$
will change due to astrometric microlensing as

$\displaystyle\Delta_{PM}$
$\displaystyle\approx\frac{\delta_{c}(t_{2})-\delta_{c}(t_{1})}{t_{2}-t_{1}}$
(D3)
$\displaystyle\approx\frac{\theta_{E}}{t_{2}-t_{1}}\Bigg{(}\frac{t_{E}}{t_{2}-t_{0}}-\frac{t_{E}}{t_{1}-t_{0}}\Bigg{)},$
(D4)

where $(t_{2}-t_{0})/t_{E}$, $(t_{1}-t_{0})/t_{E}\gg 1$. For MB10364,
OB110037, and OB110310, for $t_{1}$ and $t_{2}$ corresponding to the start and
end of Gaia EDR3 observations, $\Delta_{PM}\approx 0.003$, 0.02, and 0.01
mas/yr, respectively. In other words, the effect of lensing on proper motions
should be negligible in Gaia. Thus, if the lens is dark, the proper motions of
the source as measured by the HST data and the proper motion of the target
reported by Gaia should be the same within the uncertainties.

##  Appendix E Gaia diagnostics

Table 19: Distribution of Gaia parameters towards fields of interest

Parameter | MB09260 | MB10364 | OB110037 | OB110310
---|---|---|---|---
$N_{\rm stars}$ | 611 | 1157 | 1056 | 884
$f_{D>2}$ | 0.44 | 0.51 | 0.59 | 0.55
$f_{\epsilon=0}$ | 0.81 | 0.77 | 0.85 | 0.86
$f_{{\rm multi}=0}$ | 0.79 | 0.41 | 0.53 | 0.70
med(multi>0) | 4.50 | 14.00 | 9.00 | 5.00
multi>0[95%]) | 37.55 | 54.00 | 42.30 | 43.00
$f_{{\rm odd}=0}$ | 0.76 | 0.61 | 0.64 | 0.64
med(odd>0) | 25.50 | 30.00 | 26.00 | 27.50
odd>0[95%]) | 83.55 | 78.20 | 80.00 | 78.70
$\langle\log_{10}(\rm{amp})\rangle$ | -1.24 | -0.89 | -1.03 | -1.13
$\sigma_{\log_{10}(\rm{amp})}$ | 0.45 | 0.39 | 0.58 | 0.44
$\langle\log_{10}(\rm{ruwe})\rangle$ | 0.04 | 0.12 | 0.08 | 0.06
$\sigma_{\log_{10}(\rm{ruwe})}$ | 0.10 | 0.16 | 0.11 | 0.11
$\langle\log_{10}(C)\rangle$ | 0.23 | 0.22 | 0.25 | 0.26
$\sigma_{\log_{10}(C)}$ | 0.12 | 0.16 | 0.15 | 0.15

Note. — $N_{\rm stars}$: Number of Gaia sources within 1 deg2 of the target.

$f_{D>2}$: Fraction of sources with astrometric_excess_noise_significance$>2$.

$f_{\epsilon=0}$: Fraction of sources with astrometric_excess_noise$>2$.

$f_{{\rm multi}=0}$: Fraction of stars with ipd_frac_multi_peak = 0.

med(multi>0): Median of nonzero ipd_frac_multi_peak values.

multi>0[95%]): 95th quantile of nonzero ipd_frac_multi_peak values.

$f_{{\rm odd}=0}$: Fraction of stars with ipd_frac_odd_win= 0.

med(odd>0): Median of nonzero ipd_frac_odd_winvalues.

odd>0[95%]): 95th quantile of nonzero ipd_frac_odd_win values.

$\langle\log_{10}(\rm{amp})\rangle$: Mean of
$\log_{10}$(ipd_gof_harmonic_amplitude).

$\sigma_{\log_{10}(\rm{amp})}$: Standard deviation of
$\log_{10}$(ipd_gof_harmonic_amplitude).

$\langle\log_{10}(\rm{ruwe})\rangle$: Mean of $\log_{10}$(ruwe).

$\sigma_{\log_{10}(\rm{ruwe})}$: Standard deviation of $\log_{10}$(ruwe).

$\langle\log_{10}(C)\rangle$ : Mean of
$\log_{10}(\texttt{phot\\_bp\\_rp\\_excess\\_factor})$.

$\sigma_{\log_{10}(C)}$: Standard deviation of
$\log_{10}(\texttt{phot\\_bp\\_rp\\_excess\\_factor})$.

We would ideally like to incorporate information from Gaia EDR3 into our
target analysis (e.g. using the reported parallaxes and proper motions to
inform the prior on the source parallax in the fit, or to compare to the
posteriors as a cross-check to validate the results). However, it is known
that there are many as-of-yet-unresolved systematics in the Gaia EDR3
astrometry, especially toward the Bulge, and extra care must be taken to
verify if a proper motion or parallax for a particular Gaia source is
reliable. Gaia EDR3 is much better than Gaia DR2 in terms of photometry, but
the astrometry still has issues that need to be worked out (Fabricius et al.,
2021). We investigate several different metrics for the Gaia solutions for our
targets to determine the reliability of the reported parallax and proper
motions. A brief summary of the meaning of relevant Gaia statistics is
presented in Table 3.5; we refer the reader to the Gaia Early Data Release 3
documentation (van Leeuwen et al., 2021) for details.

###  E.1 Solution type

Sources in Gaia have varying amounts of information available, described by
the astrometric_params_solved parameter. A value of 3 signifies a 2-parameter
solution (position), a value of 31 signifies a 5-parameter solution (position,
parallax, proper motion), and a value of 95 signifies a 6-parameter solution
(position, parallax, proper motion, astrometrically estimated effective
wavenumber). 5-parameter solutions are generally the most accurate, followed
by 6-parameter, then 2-parameter solutions. 6-parameter solutions are worse
than 5-parameter solutions because it means an assumption had to be made about
the source color. This usually indicates an issue with determining the
properties of the source, which reduces the accuracy of the solution. This can
happen in very crowded regions, like the Bulge. MB09260 has a 2-parameter
solution, MB10364 and OB110310 have 6-parameter solutions, and only OB110037
has a 5-parameter solution. Based on the high-resolution HST images (Figure 1)
this is not surprising, as OB110037 is the brightest object in its vicinity,
while the other targets are near comparably bright stars.

###  E.2 Image parameter determination parameters

The Gaia EDR3 solution assumes each source is a single star. Gaia’s image
parameter determination parameters (IPD) can be used to determine whether this
assumption is valid. Different types of binaries can be identified through
large of IPD values: ipd_gof_harmonic_amplitude (indicates partially resolved
binaries, asymmetric images), ipd_frac_multi_peak (indicates resolved, close
binaries), ipd_frac_odd_win (indicates another bright source, or observation
window conflicts for wide pairs) (Lindegren et al., 2021b). Given that the
targets all have lightcurves well described by point lenses (although see §7.4
regarding indications of binarity in OB110037), we do not expect any of these
parameters to be unusual for our targets; unusual values likely indicates
problems with the astrometric solution itself due to crowding or
contamination.

The ipd_gof_harmonic_amplitude values for MB09260, MB10364, OB110037, and
OB11030 (0.089, 0.064, 0.036, 0.042, respectively) are somewhat large compared
to the median values reported in Tables 4-6 of Lindegren et al. (2021b).
However, this is to be expected toward the Bulge. Compared to other sources
within 1 deg2 (Table 19) of these targets, the values are not
unusual.242424Note the values reported in Table 19 are the logarithms of
ipd_gof_harmonic_amplitude. We choose to report the logarithms as they are
approximately normally distributed.

MB10364 has ipd_frac_multi_peak = 15. The high-resolution HST images do not
indicate MB10364 is a resolved binary. The last HST image is from 2014, and
Gaia EDR3 observations span 2014-2017. A possibility is that the lens and
source separation after 2014 was large enough to resolve with Gaia EDR3.
However, this is ruled out by considering the lens-source separation criteria
(Table 7.8) and that Gaia’s resolution is worse than HST’s. The most likely
explanation is confusion in Gaia due to stars near MB10364 since the crowding
within 1 arcsec of MB10364 is high. The other targets have ipd_frac_multi_peak
= 0.

MB09260 has ipd_frac_odd_win = 18 and OB110310 has ipd_frac_odd_win = 55.
Similar to the case of MB10364 in the previous paragraph, this is likely due
to bright stars within 1 arcsec of these targets that led to confusion in
Gaia. The other targets have ipd_frac_odd_win = 0.

###  E.3 Astrometric goodness-of-fit statistics

Gaia also presents several relevant goodness-of-fit (GOF) statistics from the
astrometric fit: the renormalized unit weight error (RUWE; ideally 1), the
extra noise required per observation to explain the residual in the
astrometric fit of the source (astrometric_excess_noise; ideally 0), and the
significance of this source noise (astrometric_excess_noise_sig; insignificant
if $<2$).

MB09260 does not have a RUWE because it is not calculated for 2-parameter
solutions. OB110037 and OB110310 both have RUWE $\sim 1$ (0.971 and 0.981,
respectively). MB10364 has a large RUWE (1.388), although not unusual compared
to the sources nearby (Table 19).252525Note the values reported in Table 19
are the logarithms of RUWE. We choose to report the logarithms as they are
approximately normally distributed.

MB09260 and MB10364 both have significant astrometric excess noise (1.241 mas
and 0.406 mas, respectively). OB110037 does not have astrometric excess noise
and the noise for OB110310 is not significant.

Table 19 lists the distributions of $D=$ astrometric excess noise significance
for the stars in Gaia within 1 deg2 of our targets. Gaia documentation notes
that $D=0$ for roughly half the sources and $D>2$ for a few percent of sources
with well-behaved astrometric solutions (van Leeuwen et al., 2021). Near
MB09260, 84% sources have $D>0$ and 44% have $D>2$. Near MB10364, 79% sources
have $D>0$ and 51% with have $D>2$. Near OB110037, 86% sources have $D>0$ and
59% have $D>2$. Near OB110310, 87% sources have $D>0$ and 54% have $D>2$. This
indicates the astrometric solution is not well-behaved. This is not
surprising, as there still exist many systematics in the astrometry,
especially toward the Bulge. Figure 25 shows astrometric excess noise
significance plotted against astrometric excess noise for systems where
astrometric excess noise is nonzero.

###  E.4 Color excess

Fabricius et al. (2021) use $C$ = phot_bp_rp_excess_factor as a proxy for
contamination due to crowding. Stars with large excess ($C>5$) tend to have
underestimated proper motion uncertainties by a factor of $\sim 1.7$. All of
our targets have $C<5$; less than 1% of stars within 1 deg2 of the targets
have $C>5$. This metric does not seem to capture crowding toward the Bulge.

###  E.5 Summary

Compared to the high resolution images from HST, the various metrics reported
by Gaia make sense.

  * •

OB110037 is the brightest source in its vicinity, and hence has a good
astrometric solution. Its astrometric solution (parallax and proper motion)
are likely to be reliable. The HST source proper motion agrees with the Gaia
proper motion measurement.

  * •

MB10364 is bright but is in a crowded field. OB110310 is somewhat isolated but
is faint. Their astrometric solutions are likely to be unreliable.

  * •

MB09260 does not have enough visibility periods to generate an astrometric
solution (which requires at least 9 visibility periods, while MB09260 only has
8).

In conclusion, although Gaia is a dedicated astrometric mission, it is not
optimized for the crowded and extincted Bulge, and the astrometric parameters
are likely to be untrustworthy there (Fabricius et al., 2021; Rybizki et al.,
2022). Although Gaia EDR3 is also an improvement over DR2, those improvements
are in the photometry and not the astrometry. Placing too much weight on the
Bulge astrometry in analyses (especially the uncertainties) should be done
with caution. To make use of Gaia data for these targets, we will need to wait
for future data releases with improved astrometry in crowded fields as well as
per-epoch astrometry.

Figure 25: Astrometric excess noise vs. astrometric excess noise significance
as a function of $G$ magnitude, for sources with non-zero astrometric excess
noise in the 1 deg2 fields surrounding the microlensing targets found in Gaia.
The targets are marked as stars; OB110037 is not marked as it has astrometric
excess noise 0. The dotted horizontal line denotes
astrometric_excess_noise_sig $=2$; values $<2$ are not significant.

##  Appendix F Gaussian Process

The photometric microlensing survey data contains uncharacterized systematics
in the lightcurves, which may be due to unaccounted correlated noise from
astrophysical processes or data acquisition and reduction. Correlated noise
can bias the results of parameter estimation. Golovich et al. (2022) fit a set
of publicly released OGLE-III and OGLE-IV microlensing lightcurves using
gaussian processes (GPs) to account for correlated noise; we follow that
approach here. We use the celerite package (Foreman-Mackey et al., 2017) to
implement GPs in our microlensing model.

For a thorough reference about GPs and their application to machine learning
and inference, the reader may consult sources such as Rasmussen & Williams
(2006). In short, a GP is composed of two parts: the noise (the stochastic
component) and the mean function (the deterministic component). The properties
of the stochastic component are governed by the covariance matrix, also known
as the kernel, of the GP. The notation

$y\sim\mathcal{GP}(\mu_{\boldsymbol{\theta}}(t),K_{\boldsymbol{\alpha}}(t,t^{\prime}))$
(F1)

indicates a GP with mean function $\mu$ with parameters $\boldsymbol{\theta}$
and kernel $K$ with parameters $\boldsymbol{\alpha}$.

When trying to infer some set of parameters $\boldsymbol{\theta}$ from data

$\boldsymbol{y}=(y_{1}(t_{1})\quad...\quad y_{N}(t_{N}))^{T}$ (F2)

the Gaussian log-likelihood is

$\mathrm{log}\mathcal{L}(\boldsymbol{y}|\boldsymbol{\theta},\boldsymbol{\alpha})=-\frac{1}{2}\boldsymbol{r}_{\boldsymbol{\theta}}^{T}K_{\boldsymbol{\alpha}}^{-1}\boldsymbol{r}_{\boldsymbol{\theta}}-\frac{1}{2}\log|K_{\boldsymbol{\alpha}}|-\frac{N}{2}\log(2\pi),$
(F3)

where $N$ is the number of data points, and $\boldsymbol{r}_{\theta}$ is the
residual vector

$\boldsymbol{r}_{\boldsymbol{\theta}}=(y_{1}-\mu_{\boldsymbol{\theta}}(t_{1})\quad...\quad
y_{N}-\mu_{\boldsymbol{\theta}}(t_{N}))^{T}.$ (F4)

In the residual vector, $\mu_{\boldsymbol{\theta}}$ is the mean model, which
in our case is the microlensing model. The kernel or covariance matrix
$K_{\boldsymbol{\alpha}}$ describes the correlated errors. If the errors are
not correlated, $K_{\boldsymbol{\alpha}}$ is diagonal and the likelihood
reduces to the familiar form

$\displaystyle\mathrm{log}\mathcal{L}(\boldsymbol{y}|\boldsymbol{\theta})$
$\displaystyle=-\frac{1}{2}\sum_{i}\Bigg{(}\mathrm{log}[2\pi\sigma_{i}^{2}]+\frac{(y_{i}-\mu_{\boldsymbol{\theta}}(t_{i}))^{2}}{\sigma_{i}^{2}}\Bigg{)}$
(F5)
$\displaystyle=-\frac{1}{2}\sum_{i}\mathrm{log}[2\pi\sigma_{i}^{2}]-\frac{1}{2}\chi^{2},$

where $\sigma_{i}$ is the uncertainty on data point $y_{i}$ and $\chi^{2}$ is
the “goodness of fit”.

Following Golovich et al. (2022), a damped driven simple harmonic oscillator
(DDSHO) kernel $\kappa_{DDSHO}$ added to a Matérn-3/2 kernel $\kappa_{M3/2}$
is used to model the correlated noise in the photometric microlensing survey
lightcurves. Both these kernels are stationary, as they are a function of the
differences of the times only:

$\tau_{ij}=|t_{i}-t_{j}|.$ (F6)

The kernel is given by

$K_{\boldsymbol{\alpha}}(\tau_{ij})=\kappa_{DDSHO}(\tau_{ij})+\kappa_{M3/2}(\tau_{ij})+\delta_{ij}\sigma_{i}^{2}.$
(F7)

Qualitatively, the DDSHO kernel models smooth variations, while the Matérn-3/2
captures more irregular variations.

The DDSHO kernel is given by

$\kappa_{DDSHO}(\tau_{ij})=S_{0}\omega_{0}e^{-\omega_{0}\tau_{ij}/\sqrt{2}}\cos{\Bigg{(}\frac{\omega_{0}\tau_{ij}}{\sqrt{2}}-\frac{\pi}{4}\Bigg{)}}$
(F8)

where $S_{0}$ controls the amplitude of the deviation from the mean model and
$\omega_{0}$ controls the variation frequency. This kernel has been used in
asteroseismic modeling (Li et al. (2019) and references therein).

The Matérn-3/2 kernel is given by

$\kappa_{M3/2}(\tau_{ij})=\sigma^{2}e^{-\sqrt{3}\tau_{ij}/\rho}\Bigg{(}1+\frac{\sqrt{3}\tau_{ij}}{\rho}\Bigg{)}$
(F9)

where $\sigma$ determines the amount of deviation allowed from the mean model,
and $\rho$ is the characteristic coherence scale. The Matérn-3/2 kernel has
been used to model correlated noise in the lightcurves of transiting and
eclipsing exoplanets (Gibson et al., 2013; Evans et al., 2015), and in
particular is appropriate for modeling non-smooth behaviors (Gilbertson et
al., 2020). For numerical reasons (see §4 of Foreman-Mackey et al., 2017), the
Matérn-3/2 kernel is approximated

$\kappa_{M3/2}(\tau_{ij})=\sigma^{2}\Bigg{[}\Big{(}1+\frac{1}{\epsilon}\Big{)}e^{-\sqrt{3}(1-\epsilon)\tau_{ij}/\rho}\Big{(}1-\frac{1}{\epsilon}\Big{)}e^{-\sqrt{3}(1+\epsilon)\tau_{ij}/\rho}\Bigg{]}.$
(F10)

In the limit $\epsilon\rightarrow 0$, this is exactly the Matérn-3/2 kernel.
We implement the approximation with $\epsilon=0.01$.

The ground-based OGLE and MOA photometry are fit using the Gaussian likelihood
with a full covariance matrix (Equation F3). The HST photometry and astrometry
are fit using the Gaussian likelihood assuming a diagonal covariance matrix
(Equation F5).

Note for MB10364, instead of fitting the MOA lightcurve using a GP, we instead
fit an additive error on the ground-based photometry. During nested sampling,
the GP showed some numerical instability.

##  Appendix G Priors

The distributions for the priors $\pi$ are described in this sections.
$\mathcal{N}(\mu,\sigma)$ denotes a normal distribution with mean $\mu$ and
standard deviation $\sigma$.
$\mathcal{N}_{T}(\mu,\sigma,l_{\sigma},u_{\sigma})$ denotes a normal
distribution with a low end truncation at $\mu+\sigma l_{\sigma}$ and a high
end truncation at $\mu+\sigma u_{\sigma}$. $\mathcal{U}(a,b)$ denotes a
uniform distribution from $a$ to $b$. $\Gamma^{-1}(\alpha,\beta)$ is the
inverse gamma distribution

$\Gamma^{-1}(x;\alpha,\beta)=\frac{\beta^{\alpha}}{\Gamma(\alpha)}x^{-\alpha-1}\exp[-\beta/x].$
(G1)

The prior distributions for each target are summarized in Table G.

Table 20: Priors

Note. — For definitions of the different variables, see Appendix G. There are
two fits for OB110462 depending on the likelihood used, “equal weighting”
(OB110462 EW) or “default weighting” (OB110462 DW). See §5.3 for details.

###  G.1 Photometry priors

The five microlensing parameters in a PSPL with parallax fit are $t_{0}$,
$u_{0}$, $t_{E}$, $\pi_{E,E}$, and $\pi_{E,N}$.

The prior on $t_{0}$ is a normal distribution centered on the time of peak
magnification in the geocentric frame, with a spread of 75 days

$\pi(t_{0})\sim\mathcal{N}(t_{peak,\oplus},75\textrm{ days}).$ (G2)

Note that the time at peak magnification in the heliocentric frame $t_{0}$ is
not necessarily the same as in the geocentric frame $t_{peak,\oplus}$, hence
the large amount of spread in the prior.

The prior on $u_{0}$ is a Gaussian with mean 0 and standard deviation 0.5

$\pi(u_{0})\sim\mathcal{N}(0,0.5)$ (G3)

which takes into account that events with smaller $|u_{0}|$ are more likely to
be detected, and that events with $|u_{0}|>1.5$ are not robustly detectable
with current ground-based surveys.262626$|u_{0}|=1.5$ corresponds to a
brightening of no more than around 0.1 mag. When selecting microlensing
events, those with a brightening less than 0.1 mag are generally excluded in
survey samples to prevent contamination from low-amplitude variables (e.g.
Mróz et al. (2019)).

The prior on $t_{E}$ is a Gaussian centered at 200 days with a large spread of
$\sigma=100$ days. The distribution is truncated at $-1.8\sigma$ and $3\sigma$
(20 and 500 days, respectively):

$\pi(t_{E})\sim\mathcal{N}_{T}(200,50,-1.8,10)\textrm{ days}.$ (G4)

The priors on the microlensing parallax are estimated from the population of
bulge microlensing events from the PopSyCLE simulation:

$\displaystyle\pi(\pi_{E,E})$ $\displaystyle\sim\mathcal{N}(-0.02,0.12)$ (G5)
$\displaystyle\pi(\pi_{E,N})$ $\displaystyle\sim\mathcal{N}(-0.03,0.13).$ (G6)

For each dataset filter, $b_{SFF}$ and $m_{base}$ are also fit. For the
ground-based photometry, we use a prior

$b_{SFF,ground}\sim\mathcal{U}(0,1.1)$ (G7)

where the negative blend flux implied by $b_{SFF}>1$ allows for some extra
noise such as imperfect background subtraction. Similarly for the HST
astrometry, we use a uniform prior on $b_{SFF}$

$b_{SFF,HST}\sim\mathcal{U}(0,1.05).$ (G8)

$m_{base}$ is a normal distribution

$m_{base}\sim\mathcal{N}(\overline{m}_{base},\sigma_{\overline{m}_{base}})$
(G9)

where $\overline{m}_{base}$ is the average magnitude during the un-magnified
seasons, weighted by the measurement uncertainties, and
$\sigma_{\overline{m}_{base}}$ is 0.1 for OGLE, 0.2 for MOA, and 0.05 for HST.

###  G.2 Gaussian Process hyperparameter priors

The ground-based photometry includes correlated noise we fit. We follow a very
similar parametrization to Golovich et al. (2022) for the GP priors. The main
difference is that fit in magnitude space instead of flux space, and so our
priors are also in magnitudes instead of fluxes.

For $\sigma$, we use the prior

$\log(\sigma/\textrm{mag})\sim\mathcal{N}(0,5)$ (G10)

which allows a wide range of lightcurve amplitude variability.

For $\rho$, we use the prior

$\rho\sim\Gamma^{-1}(a,b)$ (G11)

where $a$ and $b$ are the constants that satisfy the relation

$\displaystyle 0.01$ $\displaystyle=\int_{0}^{\textrm{med}(\Delta
t)}\Gamma^{-1}(x;a,b)\;dx$ (G12) $\displaystyle 0.01$
$\displaystyle=1-\int_{0}^{\Delta T}\Gamma^{-1}(x;a,b)\;dx$ (G13)

where $\textrm{med}(\Delta t)$ is the median duration between observations and
$\Delta T$ is the duration of full dataset. This helps suppress values at
extremely short or long timescales that might lead to ill-behaved
models.272727See betanalpha.github.io/assets/case_studies/gaussian_processes

For $S_{0}$ and $\omega_{0}$ we use the priors

$\displaystyle\log S_{0}\omega_{0}^{4}$
$\displaystyle\sim\mathcal{N}(\textrm{med}(\sigma_{m}^{2}),5)$ (G14)
$\displaystyle\log\omega_{0}$ $\displaystyle\sim\mathcal{N}(0,5).$ (G15)

###  G.3 Astrometry priors

The prior on the Einstein radius $\theta_{E}$ is a lognormal distribution
estimated from PopSyCLE for events with $t_{E}>50$ or $t_{E}>120$ days as

$\displaystyle\pi(\log_{10}(\theta_{E}/\textrm{mas}))$
$\displaystyle\sim\mathcal{N}(-0.2,0.3)$ (G16)
$\displaystyle\pi(\log_{10}(\theta_{E}/\textrm{mas}))$
$\displaystyle\sim\mathcal{N}(0,0.5),$ (G17)

respectively. We use the prior from $t_{E}>120$ days for OB110462 and
$t_{E}>50$ days for the other 4 targets.

The prior on the source parallax $\pi_{S}$ is estimated from the population of
bulge microlensing events from the PopSyCLE simulation

$\pi(\pi_{S})\sim\mathcal{N}_{T}(0.1126,0.0213,-2.9390,90.0)\textrm{ mas}$
(G18)

which corresponds to source distances ranging from 0.5 to 20 kpc.

The prior on the source proper motion $\mu_{S,E}$ and $\mu_{S,N}$ are uniform
distributions

$\displaystyle\pi(\mu_{S,E})$
$\displaystyle\sim\mathcal{U}(\overline{\mu}_{S,E}-f\sigma_{\overline{\mu}_{S,E}},\overline{\mu}_{S,E}+f\sigma_{\overline{\mu}_{S,E}})$
(G19) $\displaystyle\pi(\mu_{S,N})$
$\displaystyle\sim\mathcal{U}(\overline{\mu}_{S,N}-f\sigma_{\overline{\mu}_{S,N}},\overline{\mu}_{S,N}+f\sigma_{\overline{\mu}_{S,N}})$
(G20)

where $\overline{\mu}_{S,E}$, $\overline{\mu}_{S,N}$ are the proper motions
inferred from assuming straight-line motion (no parallax) from the F814W data,
$\sigma_{\overline{\mu}_{S,E}}$, $\sigma_{\overline{\mu}_{S,N}}$ are the
uncertainties to that fit, and $f$ is an inflation factor. To allow a wide
range of proper motions we use $f=100$.

The prior on the source position at $t_{0}$, $x_{0_{S},E}$ and $x_{0_{S},N}$,
is

$\displaystyle\pi(x_{0_{S},E})$
$\displaystyle\sim\mathcal{U}(\textrm{min}(x_{E})-f\sigma_{x_{E}},\textrm{max}(x_{E})+f\sigma_{x_{E}})$
(G21) $\displaystyle\pi(x_{0_{S},N})$
$\displaystyle\sim\mathcal{U}(\textrm{min}(x_{N})-f\sigma_{x_{N}},\textrm{max}(x_{N})+f\sigma_{x_{N}})$
(G22)

where $x_{E}$, $x_{N}$ are the positions in the F814W data of the target,
$\sigma_{x_{E}}$, $\sigma_{x_{N}}$ is the standard deviation, and $f$ is an
inflation factor. We use $f=5$.

##  Appendix H Astrometric color analysis

For some stars, the astrometric measurements taken in the F814W and F606W
filter are discrepant at the level of the reported uncertainties. We explore
this discrepancy specifically for OB110462, but this issue is also seen in
OB110037 and reference stars for all targets.

First, we must quantify the discrepancy. We consider several ways to measure
the total offset between the astrometry of the two filters across all epochs.
$\overline{dx}$ is the average of the offsets in RA across all epochs for a
particular star, and can be thought of measuring the amount of translation
between F814W and F606W. $\overline{|dx|}$ is the average of the magnitude of
the offsets in RA across all epochs for a particular star, and can be thought
of measuring the absolute amount of translation between F814W and F606W.
$|\overline{dx}|$ is the absolute value of the average of the magnitude of the
offsets in RA across all epochs for a particular star, and can be thought of
as measuring the total amount of deviation between F814W and F606W. Note that
$|\overline{dx}|$ is distinct from $\overline{|dx|}$. The definitions for
$\overline{dy}$, $\overline{|dy|}$, and $|\overline{dy}|$ are analogous,
except they are the offsets in Dec. We also consider these quantities in units
of sigma, where the differences in each epoch $dx_{i}$ and $dy_{i}$ are
normalized by the positional uncertainties $\sigma_{x,i}$ and $\sigma_{y,i}$.

In Figure 26, we show the distributions of these quantities as a function of
magnitude for stars within $30^{\prime\prime}$ of OB110462. While not falling
in the bulk of the distribution, OB110462 is not an extreme outlier.
Considering how large the variation in positional differences is, especially
for fainter stars, OB110462 seems well within the other positions. For this
reason, we assume the positional differences are a systematic we can correct
empirically. We apply a constant positional offset to the F606W OB110462
observations to match the positions of the F814W observations as described in
§4.2.5. However, further investigation to determine whether the source of the
filter dependent astrometry of OB110462 and other stars may actually be
astrophysical is worth pursuing, as are more observational programs to study
filter dependence on astrometry.

Figure 26:  Difference in astrometry between observations taken in F814W vs.
F606W for stars in the OB110462 images, aggregated across all epochs.
Definitions of $\overline{dx}$, $\overline{|dx|}$, and $|\overline{dx}|$ are
given in text of Appendix H. The squares mark OB110462. The dashed lines
indicate the magnitude range of the reference stars in the OB110462 images.

##  Appendix I Directly confronting the photometry and astrometry tension

Figure 27: _Left:_ Comparison of the source-lens relative proper motion
inferred from the photometry and astrometry using the default weighted
likelihood (_DW_), the photometry and astrometry using the equal weighted
likelihood (_EW_), from the HST astrometry alone (_HST ast_), and from the
OGLE photometry alone (_OGLE phot_). The median is shown as the thick line;
the shaded region denotes the $1\sigma$ uncertainties. _Right_ : 1 and
$2\sigma$ contours of $\pi_{E}$ and direction of the source-lens relative
proper motion $\phi$, for the 4 fits presented in the left panel. the EW fit
is consistent with the fit to the HST astrometry, although the astrometry
alone has very little constraining power. The DW and OGLE phot fits are
consistent with each other, as expected.

To try and elucidate the tension between the photometry and astrometry, we fit
the OGLE photometry alone, and separately fit the HST astrometry alone. The
results are shown in Figure 27. The HST astrometry alone does not have very
much constraining power– the uncertainties on parameters such as $t_{0}$,
$t_{E}$, and $u_{0}$ are so wide that the fit itself is not useful. However,
the results are consistent with those of the EW fit. The OGLE photometry has
much more constraining power, but alone cannot constrain the lens mass. The
results are consistent with those of the DW fit.

## References

  * Abbott et al. (2017) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017, ApJ, 851, L16, doi: 10.3847/2041-8213/aa9a35
  * Abbott et al. (2020a) —. 2020a, ApJ, 892, L3, doi: 10.3847/2041-8213/ab75f5
  * Abbott et al. (2020b) Abbott, R., Abbott, T. D., Abraham, S., et al. 2020b, ApJ, 896, L44, doi: 10.3847/2041-8213/ab960f
  * Abdurrahman et al. (2021) Abdurrahman, F. N., Stephens, H. F., & Lu, J. R. 2021, ApJ, 912, 146, doi: 10.3847/1538-4357/abee83
  * Agol & Kamionkowski (2002) Agol, E., & Kamionkowski, M. 2002, MNRAS, 334, 553, doi: 10.1046/j.1365-8711.2002.05523.x
  * Anderson (2016) Anderson, J. 2016, Empirical Models for the WFC3/IR PSF, Space Telescope WFC Instrument Science Report
  * Anderson (2021) —. 2021, Table-Based CTE Corrections for flt-Format WFC3/UVIS, Space Telescope WFC Instrument Science Report
  * Anderson et al. (2021) Anderson, J., Baggett, S., & Kuhn, B. 2021, Updating the WFC3/UVIS CTE model and Mitigation Strategies, Space Telescope WFC Instrument Science Report
  * Anderson & King (2006) Anderson, J., & King, I. R. 2006, PSFs, Photometry, and Astronomy for the ACS/WFC, Instrument Science Report ACS 2006-01
  * Anderson et al. (2008) Anderson, J., Sarajedini, A., Bedin, L. R., et al. 2008, AJ, 135, 2055, doi: 10.1088/0004-6256/135/6/2055
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Bailyn et al. (1998) Bailyn, C. D., Jain, R. K., Coppi, P., & Orosz, J. A. 1998, ApJ, 499, 367, doi: 10.1086/305614
  * Baraffe et al. (2015) Baraffe, I., Homeier, D., Allard, F., & Chabrier, G. 2015, A&A, 577, A42, doi: 10.1051/0004-6361/201425481
  * Batista et al. (2015) Batista, V., Beaulieu, J. P., Bennett, D. P., et al. 2015, ApJ, 808, 170, doi: 10.1088/0004-637X/808/2/170
  * Bellini et al. (2011) Bellini, A., Anderson, J., & Bedin, L. R. 2011, PASP, 123, 622, doi: 10.1086/659878
  * Bellini et al. (2018) Bellini, A., Libralato, M., Bedin, L. R., et al. 2018, ApJ, 853, 86, doi: 10.3847/1538-4357/aaa3ec
  * Bennett et al. (2002) Bennett, D. P., Becker, A. C., Quinn, J. L., et al. 2002, ApJ, 579, 639, doi: 10.1086/342225
  * Bennett et al. (2015) Bennett, D. P., Bhattacharya, A., Anderson, J., et al. 2015, ApJ, 808, 169, doi: 10.1088/0004-637X/808/2/169
  * Bond et al. (2001) Bond, I. A., Abe, F., Dodd, R. J., et al. 2001, MNRAS, 327, 868, doi: 10.1046/j.1365-8711.2001.04776.x
  * Bond et al. (2017) Bond, I. A., Bennett, D. P., Sumi, T., et al. 2017, MNRAS, 469, 2434, doi: 10.1093/mnras/stx1049
  * Buchner et al. (2014) Buchner, J., Georgakakis, A., Nandra, K., et al. 2014, A&A, 564, A125, doi: 10.1051/0004-6361/201322971
  * Castelli & Kurucz (2004) Castelli, F., & Kurucz, R. L. 2004, ArXiv Astrophysics e-prints
  * Choi et al. (2016) Choi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102, doi: 10.3847/0004-637X/823/2/102
  * Damineli et al. (2016) Damineli, A., Almeida, L. A., Blum, R. D., et al. 2016, MNRAS, 463, 2653, doi: 10.1093/mnras/stw2122
  * Eddington (1919) Eddington, A. S. 1919, The Observatory, 42, 119
  * El-Badry et al. (2022) El-Badry, K., Seeburger, R., Jayasinghe, T., et al. 2022, MNRAS, 512, 5620, doi: 10.1093/mnras/stac815
  * Evans et al. (2019a) Evans, I. N., Allen, C., Anderson, C. S., et al. 2019a, in American Astronomical Society Meeting Abstracts, Vol. 233, American Astronomical Society Meeting Abstracts #233, 379.01
  * Evans et al. (2019b) Evans, I. N., Allen, C., Anderson, C. S., et al. 2019b, in AAS/High Energy Astrophysics Division, Vol. 17, AAS/High Energy Astrophysics Division, 114.01
  * Evans et al. (2020) Evans, P. A., Page, K. L., Osborne, J. P., et al. 2020, ApJS, 247, 54, doi: 10.3847/1538-4365/ab7db9
  * Evans et al. (2015) Evans, T. M., Aigrain, S., Gibson, N., et al. 2015, MNRAS, 451, 680, doi: 10.1093/mnras/stv910
  * Fabricius et al. (2021) Fabricius, C., Luri, X., Arenou, F., et al. 2021, A&A, 649, A5, doi: 10.1051/0004-6361/202039834
  * Farr et al. (2011) Farr, W. M., Sravan, N., Cantrell, A., et al. 2011, ApJ, 741, 103, doi: 10.1088/0004-637X/741/2/103
  * Feroz et al. (2009) Feroz, F., Hobson, M. P., & Bridges, M. 2009, MNRAS, 398, 1601, doi: 10.1111/j.1365-2966.2009.14548.x
  * Foreman-Mackey et al. (2017) Foreman-Mackey, D., Agol, E., Ambikasaran, S., & Angus, R. 2017, AJ, 154, 220, doi: 10.3847/1538-3881/aa9332
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1, doi: 10.1051/0004-6361/201629272
  * Gaia Collaboration et al. (2021) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 649, A1, doi: 10.1051/0004-6361/202039657
  * Gennaro (2018) Gennaro, M. e. a. 2018, WFC3 Data Handbook, Version 4.0, Baltimore: STScI
  * Gibson et al. (2013) Gibson, N. P., Aigrain, S., Barstow, J. K., et al. 2013, MNRAS, 428, 3680, doi: 10.1093/mnras/sts307
  * Gilbertson et al. (2020) Gilbertson, C., Ford, E. B., Jones, D. E., & Stenning, D. C. 2020, ApJ, 905, 155, doi: 10.3847/1538-4357/abc627
  * Golovich et al. (2022) Golovich, N., Dawson, W., Bartolić, F., et al. 2022, ApJS, 260, 2, doi: 10.3847/1538-4365/ac5969
  * Gould (2004) Gould, A. 2004, ApJ, 606, 319, doi: 10.1086/382782
  * Green (2018) Green, G. 2018, The Journal of Open Source Software, 3, 695, doi: 10.21105/joss.00695
  * Groth (1986) Groth, E. J. 1986, AJ, 91, 1244, doi: 10.1086/114099
  * Hearnshaw et al. (2006) Hearnshaw, J. B., Abe, F., Bond, I. A., et al. 2006, in The 9th Asian-Pacific Regional IAU Meeting, ed. W. Sutantyo, P. W. Premadi, P. Mahasena, T. Hidayat, & S. Mineshige, 272.  https://arxiv.org/abs/astro-ph/0509420
  * Hog et al. (1995) Hog, E., Novikov, I. D., & Polnarev, A. G. 1995, A&A, 294, 287
  * Hosek et al. (2020) Hosek, Matthew W., J., Lu, J. R., Lam, C. Y., et al. 2020, AJ, 160, 143, doi: 10.3847/1538-3881/aba533
  * Hosek et al. (2015) Hosek, Jr., M. W., Lu, J. R., Anderson, J., et al. 2015, ApJ, 813, 27, doi: 10.1088/0004-637X/813/1/27
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Husser et al. (2013) Husser, T.-O., Wende-von Berg, S., Dreizler, S., et al. 2013, A&A, 553, A6, doi: 10.1051/0004-6361/201219058
  * Jonker et al. (2021) Jonker, P. G., Kaur, K., Stone, N., & Torres, M. A. P. 2021, ApJ, 921, 131, doi: 10.3847/1538-4357/ac2839
  * Kains et al. (2017) Kains, N., Calamida, A., Sahu, K. C., et al. 2017, ApJ, 843, 145, doi: 10.3847/1538-4357/aa78eb
  * Koester (2010) Koester, D. 2010, Mem. Soc. Astron. Italiana, 81, 921
  * Kreidberg et al. (2012) Kreidberg, L., Bailyn, C. D., Farr, W. M., & Kalogera, V. 2012, ApJ, 757, 36, doi: 10.1088/0004-637X/757/1/36
  * Kroupa (2001) Kroupa, P. 2001, MNRAS, 322, 231, doi: 10.1046/j.1365-8711.2001.04022.x
  * Kuhn & Anderson (2021) Kuhn, B., & Anderson, J. 2021, WFC3/UVIS: New FLC External CTE Monitoring 2009-2020, Space Telescope WFC Instrument Science Report
  * Lam & Lu (2021) Lam, C., & Lu, J. R. 2021, First detection of an isolated stellar mass black hole with astrometric microlensing, HST Proposal
  * Lam et al. (2020) Lam, C. Y., Lu, J. R., Hosek, Matthew W., J., Dawson, W. A., & Golovich, N. R. 2020, ApJ, 889, 31, doi: 10.3847/1538-4357/ab5fd3
  * Lam et al. (2022) Lam, C. Y., Lu, J. R., Udalski, A., et al. 2022, ApJL, arXiv:2202.01903.  https://arxiv.org/abs/2202.01903
  * Li et al. (2019) Li, S. S., Zang, W., Udalski, A., et al. 2019, MNRAS, 488, 3308, doi: 10.1093/mnras/stz1873
  * Lindegren et al. (2021a) Lindegren, L., Bastian, U., Biermann, M., et al. 2021a, A&A, 649, A4, doi: 10.1051/0004-6361/202039653
  * Lindegren et al. (2021b) Lindegren, L., Klioner, S. A., Hernández, J., et al. 2021b, A&A, 649, A2, doi: 10.1051/0004-6361/202039709
  * Lu et al. (2016) Lu, J. R., Sinukoff, E., Ofek, E. O., Udalski, A., & Kozlowski, S. 2016, ApJ, 830, doi: 10.3847/0004-637X/830/1/41
  * Luri et al. (2018) Luri, X., Brown, A. G. A., Sarro, L. M., et al. 2018, A&A, 616, A9, doi: 10.1051/0004-6361/201832964
  * Maeda et al. (2005) Maeda, Y., Kubota, A., Kobayashi, Y., et al. 2005, ApJ, 631, L65, doi: 10.1086/491698
  * Manchester et al. (2005) Manchester, R. N., Hobbs, G. B., Teoh, A., & Hobbs, M. 2005, AJ, 129, 1993, doi: 10.1086/428488
  * Mao et al. (2002) Mao, S., Smith, M. C., Woźniak, P., et al. 2002, MNRAS, 329, 349, doi: 10.1046/j.1365-8711.2002.04986.x
  * Marshall et al. (2006) Marshall, D. J., Robin, A. C., Reylé, C., Schultheis, M., & Picaud, S. 2006, A&A, 453, 635, doi: 10.1051/0004-6361:20053842
  * McGill et al. (2020) McGill, P., Everall, A., Boubert, D., & Smith, L. C. 2020, MNRAS, 498, L6, doi: 10.1093/mnrasl/slaa118
  * Miyamoto & Yoshii (1995) Miyamoto, M., & Yoshii, Y. 1995, AJ, 110, 1427, doi: 10.1086/117616
  * Mróz et al. (2019) Mróz, P., Udalski, A., Skowron, J., et al. 2019, ApJS, 244, 29, doi: 10.3847/1538-4365/ab426b
  * Nucita et al. (2006) Nucita, A. A., De Paolis, F., Ingrosso, G., et al. 2006, ApJ, 651, 1092, doi: 10.1086/507784
  * Özel et al. (2010) Özel, F., Psaltis, D., Narayan, R., & McClintock, J. E. 2010, ApJ, 725, 1918, doi: 10.1088/0004-637X/725/2/1918
  * Paczynski (1986) Paczynski, B. 1986, ApJ, 304, 1, doi: 10.1086/164140
  * Park et al. (2004) Park, B. G., DePoy, D. L., Gaudi, B. S., et al. 2004, ApJ, 609, 166, doi: 10.1086/420926
  * Penny et al. (2019) Penny, M. T., Gaudi, B. S., Kerins, E., et al. 2019, ApJS, 241, 3, doi: 10.3847/1538-4365/aafb69
  * Poindexter et al. (2005) Poindexter, S., Afonso, C., Bennett, D. P., et al. 2005, ApJ, 633, 914, doi: 10.1086/468182
  * Rasmussen & Williams (2006) Rasmussen, C. E., & Williams, C. K. I. 2006, Gaussian Processes for Machine Learning
  * Robin et al. (2003) Robin, A. C., Reylé, C., Derrière, S., & Picaud, S. 2003, A&A, 409, 523, doi: 10.1051/0004-6361:20031117
  * Rowan et al. (2021) Rowan, D. M., Stanek, K. Z., Jayasinghe, T., et al. 2021, MNRAS, 507, 104, doi: 10.1093/mnras/stab2126
  * Rybicki et al. (2018) Rybicki, K. A., Wyrzykowski, Ł., Klencki, J., et al. 2018, MNRAS, 476, 2013, doi: 10.1093/mnras/sty356
  * Rybizki et al. (2022) Rybizki, J., Green, G. M., Rix, H.-W., et al. 2022, MNRAS, 510, 2597, doi: 10.1093/mnras/stab3588
  * Sabbi et al. (2016) Sabbi, E., Lennon, D. J., Anderson, J., et al. 2016, ApJS, 222, 11, doi: 10.3847/0067-0049/222/1/11
  * Sahu (2009) Sahu, K. 2009, Detecting Isolated Black Holes through Astrometric Microlensing, HST Proposal
  * Sahu (2012) —. 2012, Detecting Isolated Black Holes through Astrometric Microlensing, HST Proposal
  * Sahu et al. (2017) Sahu, K. C., Anderson, J., Casertano, S., et al. 2017, Science, 356, 1046, doi: 10.1126/science.aal2879
  * Sarmiento et al. (2019) Sarmiento, M. H., Arévalo, M., Arviset, C., et al. 2019, in Astronomical Society of the Pacific Conference Series, Vol. 521, Astronomical Data Analysis Software and Systems XXVI, ed. M. Molinaro, K. Shortridge, & F. Pasian, 104
  * Sharma et al. (2011) Sharma, S., Bland-Hawthorn, J., Johnston, K. V., & Binney, J. 2011, ApJ, 730, 3, doi: 10.1088/0004-637X/730/1/3
  * Skilling (2004) Skilling, J. 2004, in American Institute of Physics Conference Series, Vol. 735, Bayesian Inference and Maximum Entropy Methods in Science and Engineering: 24th International Workshop on Bayesian Inference and Maximum Entropy Methods in Science and Engineering, ed. R. Fischer, R. Preuss, & U. V. Toussaint, 395–405, doi: 10.1063/1.1835238
  * Skowron et al. (2016) Skowron, J., Udalski, A., Kozłowski, S., et al. 2016, Acta Astron., 66, 1.  https://arxiv.org/abs/1604.01966
  * Speagle (2020) Speagle, J. S. 2020, MNRAS, 493, 3132, doi: 10.1093/mnras/staa278
  * Spergel et al. (2015) Spergel, D., Gehrels, N., Baltay, C., et al. 2015, arXiv e-prints, arXiv:1503.03757.  https://arxiv.org/abs/1503.03757
  * Sumi (2008) Sumi, T. 2008, in Manchester Microlensing Conference, ed. E. Kerins, S. Mao, N. Rattenbury, & L. Wyrzykowski, 25
  * Thompson et al. (2019) Thompson, T. A., Kochanek, C. S., Stanek, K. Z., et al. 2019, Science, 366, 637, doi: 10.1126/science.aau4005
  * Udalski et al. (2015) Udalski, A., Szymański, M. K., & Szymański, G. 2015, Acta Astron., 65, 1\.  https://arxiv.org/abs/1504.05966
  * van der Walt et al. (2011) van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science & Engineering, 13, 22, doi: 10.1109/MCSE.2011.37
  * van Leeuwen et al. (2021) van Leeuwen, F., de Bruijne, J., Babusiaux, C., et al. 2021, Gaia EDR3 documentation, Gaia EDR3 documentation
  * Virtanen et al. (2019) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2019, arXiv e-prints, arXiv:1907.10121.  https://arxiv.org/abs/1907.10121
  * Walker (1995) Walker, M. A. 1995, ApJ, 453, 37, doi: 10.1086/176367
  * Wang (1993) Wang, Y. H. 1993, Statistica Sinica, 3, 295.  http://www.jstor.org/stable/24304959
  * WFIRST Astrometry Working Group et al. (2019) WFIRST Astrometry Working Group, Sanderson, R. E., Bellini, A., et al. 2019, Journal of Astronomical Telescopes, Instruments, and Systems, 5, 044005, doi: 10.1117/1.JATIS.5.4.044005
  * Wiktorowicz et al. (2020) Wiktorowicz, G., Lu, Y., Wyrzykowski, Ł., et al. 2020, ApJ, 905, 134, doi: 10.3847/1538-4357/abc699
  * Wozniak (2000) Wozniak, P. R. 2000, Acta Astron., 50, 421.  https://arxiv.org/abs/astro-ph/0012143
  * Wyrzykowski et al. (2016) Wyrzykowski, Ł., Kostrzewa-Rutkowska, Z., Skowron, J., et al. 2016, MNRAS, 458, 3012, doi: 10.1093/mnras/stw426
  * Yalinewich et al. (2018) Yalinewich, A., Beniamini, P., Hotokezaka, K., & Zhu, W. 2018, MNRAS, 481, 930, doi: 10.1093/mnras/sty2327
  * Yamaguchi et al. (2018) Yamaguchi, M. S., Kawanaka, N., Bulik, T., & Piran, T. 2018, ApJ, 861, 21, doi: 10.3847/1538-4357/aac5ec
  * Zurlo et al. (2018) Zurlo, A., Gratton, R., Mesa, D., et al. 2018, MNRAS, 480, 236, doi: 10.1093/mnras/sty1805