# Primordial Black Hole Dark Matter Simulations Using PopSyCLE

Kerianne Pruett Lawrence Livermore National Laboratory, Livermore, CA 94550
William Dawson Lawrence Livermore National Laboratory, Livermore, CA 94550
Michael S. Medford Department of Astronomy, University of California,
Berkeley, CA 94720  Lawrence Berkeley National Laboratory, Berkeley, CA 94720
Jessica R. Lu Department of Astronomy, University of California, Berkeley, CA
94720  Casey Lam Department of Astronomy, University of California, Berkeley,
CA 94720  Scott Perkins Lawrence Livermore National Laboratory, Livermore, CA
94550  Peter McGill Lawrence Livermore National Laboratory, Livermore, CA
94550  Nathan Golovich Lawrence Livermore National Laboratory, Livermore, CA
94550  George Chapline  Lawrence Livermore National Laboratory, Livermore, CA
94550

###### Abstract

Primordial black holes (PBHs), theorized to have originated in the early
universe, are speculated to be a viable form of dark matter. If they exist,
they should be detectable through photometric and astrometric signals
resulting from gravitational microlensing of stars in the Milky Way.
Population Synthesis for Compact-object Lensing Events, or PopSyCLE, is a
simulation code that enables users to simulate microlensing surveys, and is
the first of its kind to include both photometric and astrometric microlensing
effects, which are important for potential PBH detection and characterization.
To estimate the number of observable PBH microlensing events we modify
PopSyCLE to include a dark matter halo consisting of PBHs. We detail our PBH
population model, and demonstrate our PopSyCLE \+ PBH results through
simulations of the OGLE-IV and Roman microlensing surveys. We provide a proof-
of-concept analysis for adding PBHs into PopSyCLE, and thus include many
simplifying assumptions, such as $f_{\text{DM}}$, the fraction of dark matter
composed of PBHs, and $\bar{m}_{\text{PBH}}$, mean PBH mass. Assuming
$\bar{m}_{\text{PBH}}=30$ $M_{\odot}$, we find $\sim$ 3.65$f_{\text{DM}}$
times as many PBH microlensing events than stellar evolved black hole events,
a PBH average peak Einstein crossing time of $\sim$ 91.4 days, estimate on
order of $10^{2}f_{\text{DM}}$ PBH events within the 8 year OGLE-IV results,
and estimate Roman to detect on the order of $10^{3}f_{\text{DM}}$ PBH
microlensing events throughout its planned microlensing survey.

gravitational lensing: micro – black hole physics – dark matter – Galaxy:halo

tablenum SIXtablenum

##  1 Introduction

Approximately 85% of the matter in our universe is composed of dark matter,
which is not well-understood. Dark matter cannot be observed via
electromagnetic mechanisms (hence “dark”) and is only observable through
gravitational interactions, making it challenging to detect. Primordial black
holes (PBHs), hypothesized in the 1960s (Zel’dovich & Novikov, 1967), are a
theoretical type of black hole believed to be a viable dark matter candidate.
Unlike black holes that are formed through stellar evolution, PBHs are thought
to have formed in the early universe (within 1 second of the Big Bang, e.g.,
Carr & Silk, 2018), via gravitational collapse in over-dense regions. PBHs are
theorized to have formed during the radiation dominated era, before Big Bang
nucleosynthesis, thus they are considered non-baryonic and are speculated to
behave like “cold” dark matter. Our current understanding of the universe
places a $5\%$ constraint on the baryonic matter fraction (Collaboration et
al., 2014), meaning that stellar evolved black holes can never account for
more than $5\%$ of the energy density in our universe. Thus, while we know
that stellar evolved black holes exist and may contribute to the dark matter
fraction, we know they cannot explain the majority of the missing dark matter
in our universe.

Because PBHs formed at roughly the cosmological horizon (Chapline, 1975a),
their initial masses should correspond to the horizon mass:

$M\sim\frac{c^{3}t}{G},$

where $c$ is the speed of light, $t$ is the time since the big bang, and $G$
is the gravitational constant (Carr, 1975; Carr & Kuhnel, 2021). This means
that PBHs should have masses between $\sim$$10^{-8}$ $kg$ (if formed at Planck
time $t=10^{-43}$ $s$) and $\sim$$10^{5}$ $M_{\odot}$ (if formed at $t=1$ $s$)
(Carr, 1975). PBHs are considered “free-floating” in the late universe, and
should retain roughly their initial masses (Chapline, 1975a). The initial
conjecture of PBHs drove Hawking (1971) to study quantum gravitational effects
of black holes, and led to the discovery that PBHs formed before $t=10^{-23}$
$s$ ($M\sim 10^{12}$ $kg$) should have evaporated by present day (Hawking,
1974).

In 2016, the Laser Interferometer Gravitational-Wave Observatory (LIGO)
detected gravitational waves from two merging $\sim$30 $M_{\odot}$ black holes
(Abbott et al., 2016), which were potentially primordial in origin (Bird et
al., 2016). The idea of PBHs as a viable form of dark matter had nearly died
out over the last few decades due to insufficient new evidence or microlensing
results (e.g., Alcock et al., 1993; Aubourg et al., 1993; Croon et al., 2020;
Udalski et al., 1994; Wyrzykowski et al., 2011, etc.), however, our inability
to detect dark matter (of any variety) throughout that time, combined with
recent scientific discoveries, has reinvigorated scientists to explore PBHs as
a potential dark matter candidate. These recent explorations have placed very
tight constraints on $f_{\text{DM}}$ covering many decades in mass (see
summary in Bird et al., 2022).

If PBHs exist, they should be present within our Milky Way galaxy (Chapline,
1975b) and observable via gravitational microlensing signals from local stars
(e.g., Chapline & Frampton, 2016). Microlensing is the effect when a
foreground object (the lens) passes in front of a luminous background object
(the source), relative to the observer, in such a way that causes
magnification of the source flux. While microlensing can cause magnification
on the same order as strong gravitational lensing (i.e. multiple images), it
happens on a smaller angular scale, such that lensed images are too small to
be resolved.

Microlensing imparts a characteristic photometric signal on the measured light
curves of lensed stars, which can be used to detect and characterize
foreground objects irrespective of their luminosity, or whether they are
gravitationally bound to a companion. Unlike other probes, such as
gravitational waves or X-ray binaries, microlensing can detect free-floating
primordial black holes, in addition to binary black hole systems, and enables
us to infer properties of the lens without observing the lens itself. Because
primordial black holes interact gravitationally and not electromagnetically,
microlensing is the most direct means of probing PBHs.

The angular Einstein radius of a microlensing event, $\theta_{E}$, sets the
angular scale of the microlensing event, and is given by:

$\theta_{E}=\sqrt{\frac{4GM}{c^{2}}\Bigg{(}\frac{1}{d_{L}}-\frac{1}{d_{S}}\Bigg{)}},$
(1)

where ${d_{S}}$ is the observer-source distance, and ${d_{L}}$ is the
observer-lens distance.

Microlensing parallax, $\pi_{E}$, causes asymmetries in the photometric light
curve and is caused by differences in the relative location of the observer
with respect to the source-lens axis, as can happen when Earth orbits the Sun
during long duration events:

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

where $\pi_{rel}$ is the relative parallax:

$\pi_{rel}=1\mathrm{AU}\bigg{(}\frac{1}{d_{L}}-\frac{1}{d_{S}}\bigg{)}.$

The Einstein crossing time can then be determined using the magnitude of the
relative source-lens proper motion, $\mu_{rel}$:

$t_{E}=\frac{\theta_{E}}{\mu_{rel}},$

which describes the characteristic timescale of the event (i.e., the time it
takes for the source to transit the Einstein radius of the lens).

In addition to photometric measurements, one can obtain astrometric
measurements for microlensing events using the location and movement of stars.
In the absence of lensing and parallax a source star would appear to make a
straight line on the sky over most observable timescales, however, in the
presence of lensing, the source appears to shift non-linearly throughout the
duration of the microlensing event. Thus, by measuring the change in location
of the source image centroid throughout the microlensing event, one can
determine the astrometric shift, $\vec{\delta_{c}}$:

$\vec{\delta_{c}}=\frac{\theta_{E}}{u^{2}+2}\vec{u}.$ (3)

Then, using the angular on-sky positions of the source and lens ($\varphi_{S}$
and $\varphi_{L}$, respectively), one can determine the source-lens separation
vector:

$\vec{u}=\frac{\vec{\varphi_{S}}-\vec{\varphi_{L}}}{\theta_{E}}.$

By comparing the observed distribution of microlensing events to simulated
distributions in $\pi_{E}$, $t_{E}$, $\bar{\delta}_{c}$, and $\mu_{rel}$
space, both with and without PBHs, we can ultimately place constraints on PBH
dark matter abundance in the Milky Way. We can then estimate expected
statistical signatures for each population by simulating the stellar and dark
matter components of the Milky Way, and can analyze resulting microlensing
properties after being filtered by a simulated observational survey.

In this work, we describe how PopSyCLE is used to generate a simulated
microlensing survey in the Milky Way (§2), the steps used to inject PBHs into
PopSyCLE (§3), and analysis on the resulting microlensing events (§4). As a
proof-of-concept, we then consider our simulation results in the context of
current and future surveys (§5).

##  2 PopSyCLE Pipeline

Population Synthesis Code for Compact Object Microlensing Events, or
PopSyCLE111https://github.com/jluastro/PopSyCLE, is a Python package that
enables the user to simulate synthetic microlensing surveys in the Milky Way.
PopSyCLE is a resolved microlensing simulation code, and the first of its kind
to include both photometric and astrometric microlensing signals, and
numerically derived compact object models (Lam et al., 2020, henceforth
Lam2020).

PopSyCLE relies on first using Galaxia222http://galaxia.sourceforge.net
(Sharma et al., 2011) to generate a synthetic model of the Milky Way populated
with stars, and SPISEA333https://github.com/astropy/SPISEA  (Hosek Jr. et al.,
2020) to generate and inject compact objects into the Milky Way model.
PopSyCLE then uses the resulting outputs to simulate a microlensing survey
with user-specified parameters, enabling the user to control survey length,
observing frequency, and various observational parameters. We briefly
summarize these tools in the remainder of this section (additional detail in
Lam2020).

###  2.1 Galaxia

Galaxia is first used to populate a stellar model of the Milky Way, using the
Besançon analytic model (Robin et al., 2004). We adopt changes to the bulge
kinematics (following the “v3” modifications from Lam2020 Appendix A) and
reddening laws (Lam2020 Appendix B) in the standard Galaxia package in an
attempt to eliminate discrepancies with observations, and to more accurately
represent the event rate predictions for the scope of this study.

For bulge kinematics, we modify the Milky Way parameters in Galaxia to use a
pattern speed of $\Omega=40$ $km$ $s^{-1}$ $kpc^{-1}$, bulge velocity
dispersion of $\sigma_{R}=\sigma_{\phi}=100$ $km$ $s^{-1}$, bar angle of
$\alpha=$$$, and bar length of $x_{0}=0.7$ $kpc$. The default Galaxia
extinction maps are not accurate in the high extinction regions of the inner
bulge, so we modify our output in accordance to the reddening law from
Damineli et al. (2016).

###  2.2 SPISEA

The Stellar Population Interface for Stellar Evolution and Atmospheres,
SPISEA, is a Python package used to generate star clusters (single-age and
single-metallicity populations) using various input parameters, including age,
mass, metallicity, extinction, atmospheric models, and initial mass functions
(Hosek Jr. et al., 2020). SPISEA allows for compact object generation by
supporting user control of the initial-final mass relation (IFMR) (Rose et
al., 2022). We follow the SPISEA modifications outlined in Lam2020 §2.2, which
uses the Raithel18 (Raithel et al., 2018) IFMR from recent simulations, along
with the zero-age main sequence (ZAMS) mass of each star, in order to build a
population of stellar evolved black holes (BH), neutron stars (NS), and white
dwarfs (WD).

Table 1: PBH Simulation Parameters  Parameter |  Description |  Value
---|---|---
$\bar{m}_{\text{PBH}}$ |  Mean PBH mass |  30 $M_{\odot}$
$f_{\text{DM}}$ |  PBH Dark matter fraction |  1.0
$v_{\text{esc}}$ |  Milky Way escape velocity |  550 $km$ $s^{-1}$
$\rho_{0}$ |  Characteristic density of Milky Way dark matter halo |  0.0093 $M_{\odot}$ $pc^{-3}$
$r_{s}$ |  Scale radius of Milky Way dark matter halo |  18.6 $kpc$
$r_{gc}$ |  Distance between the Sun and Galactic center |  8.3 $kpc$
$r_{\text{max}}$ |  Max distance from Earth to populate PBHs ($2\times r_{gc}$) |  16.6 $kpc$
$\gamma$ |  Milky Way Halo inner slope |  1, 0.5, or 0.25

###  2.3 PopSyCLE

Galaxia assigns objects to four populations: thin disk ($150$ $Myr$ to $10$
$Gyr$), bulge ($10$ $Gyr$), thick disk ($11$ $Gyr$), and stellar halo ($14$
$Gyr$). PopSyCLE then splits the thin disk population into 29 age bins and
changes the age of the stellar halo population to $13.8$ $Gyr$ (Lam2020 §3).

Once SPISEA has generated the compact objects444We use PopSyCLE at commit
#018af97 which considers only singular stars and compact objects., they get
injected back into the stellar population555PopSyCLE is being updated to
support binary populations, updated IFMRs, and a NS mass spectrum, however,
these updates are not reflected in the results of this work.. Following
Lam2020, we assign initial “kick” velocities (in addition to stellar
velocities) to NSs and BHs, assigning values of 350 $km$ $s^{-1}$ and 100 $km$
$s^{-1}$, respectively. Compact object (BHs, WDs, NSs) position and velocity
follow that of the stellar population in Galaxia, and are considered “dark”
sources (with the exception of luminous white dwarfs) (Lam2020).

Microlensing events are assumed to be “point source, point lens” (meaning both
the source and lens are treated as point sources), and event calculation
follows that described in Lam2020 §4. Initial event candidates are determined
using a cut on the separation between each object and its nearest on-sky
neighbor, $\Delta\theta\leq\theta_{blend}$, in which $\theta_{blend}$ is the
survey specific seeing-disk radius, and defines the angular projected distance
from the lens in which objects are considered blended.

Event detectability is determined by calculating the Einstein radius,
$\theta_{E}$, following Equation 1. The list of candidates is then narrowed
down to contain only those events which satisfy the following equation:

$\frac{\Delta\theta}{\theta_{E}}<u_{0,max},$

where $u_{0,max}$ is the impact parameter threshold for the survey. Events are
then cut again to ensure that they satisfy the user-specified selection
criteria, related to the survey duration and cadence. Lastly, events with dark
objects as sources are removed, as they are non-observable events. Every event
candidate remaining is considered an event in the simulation and the user is
left with catalogs containing photometric and astrometric microlensing
parameters for each event.

##  3 Primordial Black Hole Injection

Before this work, PopSyCLE supported the compact object populations of NSs,
WDs, and BHs, but had no implementation for PBHs. We update PopSyCLE to
include the addition of an optional PBH
population666https://github.com/jluastro/PopSyCLE/tree/pbh using the
parameters and assumptions given in Table 1.

Including a PBH population in PopSyCLE requires injecting PBHs into our
projected survey line-of-sight (LOS) by first determining how many PBHs to
simulate for the given field (detailed in §3.2). PBHs have spatial and
velocity distributions following that of the dark matter halo, rather than the
distributions for stellar remnant compact objects which were born in the
galactic disk. For this reason, we then have to assign positions (detailed in
§3.3) and velocities (detailed in §3.4) to each PBH based on the halo density
and speed profiles. Three-dimensional positions and velocities are assigned by
sampling the radial components (with respect to the Galactic center) from the
relevant distributions, and randomly sampling the angular components over the
possible ranges. Once a PBH population for a given field has been generated,
we inject the objects into the synthetic Milky Way generated from Galaxia and
SPISEA, thus enabling PopSyCLE to simulate a microlensing survey with PBHs
present.

###  3.1 PBH Distribution

We assume that Milky Way dark matter follows the halo density profile in
McMillan (2016):

$\rho\left(r\right)=\frac{\rho_{0}}{\left(\frac{r}{r_{s}}\right)^{\gamma}\left(1+\frac{r}{r_{s}}\right)^{\left(3-\gamma\right)}},$
(4)

where $r$ is the distance from the galactic center (in the galactocentric
coordinate frame), $\gamma$ is the halo inner density slope, $r_{s}$ is the
scale radius of the halo, and $\rho_{0}$ is the halo characteristic density.
An inner slope of $\gamma=1$ corresponds to the Navarro-Frenk-White (NFW)
profile (Navarro et al., 1996), and is used for all simulation results in this
work (though PopSyCLE + PBH supports $\gamma$ values of 1, 0.5, and 0.25). The
upper panel of Figure 1 demonstrates how the halo density changes as a
function of radius and $\gamma$. We choose $r_{gc}=8.3$ $kpc$ to be the
distance between the Sun and Galactic center (which is typical of the
literature, e.g., Gillessen et al., 2009; Chatzopoulos et al., 2015),
$r_{s}=18.6$ $kpc$ to be the halo scale radius, $\rho_{0}=0.0093$ $M_{\odot}$
$pc^{-3}$ to be the halo characteristic density (which are the median $r_{s}$
and $\rho_{0}$ values determined by McMillan, 2016, given their mass model),
and $r_{\text{max}}\equiv 2\;r_{gc}$ to be the maximum galactocentric radius
(towards the galactic center) to populate objects out to (as we do not expect
to detect a significant number of microlensing events beyond that).

###  3.2 PBH Number

To determine the number of PBHs to include in a given simulation, we first
construct a LOS extending from $r=0$ (the Sun) to the center of the field
being simulated ($l$ and $b$, at a heliocentric distance of
$r=r_{\text{max}}$). We then assume that PBHs will populate a heliocentric LOS
cylinder (shown in Figure 2, described in §3.3) whose center axis is the LOS,
and radius, $R$, corresponds to the projected survey radius at the distance
$r_{\text{max}}$ (see Figure 2). Next, we split the LOS into 100,000 equally
spaced points for numerically integrating the LOS dark matter density
(accounting for galactocentric to LOS coordinate transformations), enabling us
to approximate the total mass within the LOS cylinder, $M_{\text{LOS}}$777This
approximation becomes less accurate as $R$ increases, however, we expect this
error to be sub-dominant to other approximations, as we simulate at most 0.34
$deg^{2}$..

Using the fraction of dark matter composed of PBHs, $f_{\text{DM}}$, and mean
PBH mass, $\bar{m}_{\text{PBH}}$, we can then determine the number of PBHs to
inject in a given field, $N_{\text{PBH}}$:

$N_{\text{PBH}}=f_{\text{DM}}\left(\frac{M_{\text{LOS}}}{\bar{m}_{\text{PBH}}}\right)$
(5)

Figure 1: Upper: Halo density as a function of radius, for varying halo inner
density slopes. Lower: A recreation of the Eddington inversion curves from
Figure 11 of Lacroix et al. (2018), demonstrating how the mean speed profile
changes as a function of inner slope. The standard halo model is shown by the
solid line. In both panels the stars show that the injected PBHs trace the
intended profiles.

For the scope of this work we assume $f_{\text{DM}}=1$ and assign PBHs a
monochromatic mass of $\bar{m}_{\text{PBH}}=30$ $M_{\odot}$, as we find this
to be roughly consistent with inferred masses of LIGO detected black holes
(e.g., Abbott et al., 2016, 2017a, 2017b; Lu et al., 2019). We acknowledge
that these assumptions contradict existing $f_{\text{DM}}$ constraints (e.g.,
Blaineau et al., 2022; Alcock et al., 2001; Tisserand et al., 2007;
Wyrzykowski et al., 2011), and expectations of an extended mass distribution
(e.g., Chapline, 1975a, b; Chapline & Frampton, 2016; Chapline & Barbieri,
2018; Carr, 1975; Carr & Silk, 2018; Carr & Kuhnel, 2021), however, our
primary focus is on demonstrating the addition of a PBH population into the
PopSyCLE simulation code. For this reason, we believe these choices are
acceptable, as they enable readers to easily transform the results to other
values of $f_{\text{DM}}$ and $\bar{m}_{\text{PBH}}$, while reducing Poisson
noise.

Figure 2: PBH position assignment for a single PBH in a simulation _Left_ : A
cylinder is constructed with radius $R$ (the survey radius projected at
$r_{max}$), and centered on the LOS (dotted line extending from the Sun to the
field center). PBHs are assigned a distance, $r$, between [0, $r_{max}$], by
uniformly sampling the cumulative distribution function of Equation 4 using
standard inverse transform sampling methods. _Middle_ : A cross-sectional view
of the LOS cylinder for the location assignment of a single PBH. Three-
dimensional positions within the LOS cylinder are assigned by randomly
sampling a radial distance $\hat{R}$ between [0, $R$] (proportional to the
cross-sectional area, $\pi\hat{R}^{2}$), and uniformly sampling an angle
between [0, 2$\pi$]. This is repeated until each PBH has a 3D position within
the cylinder. _Right_ : PBHs that fall outside of the survey light cone are
masked out.

###  3.3 PBH Positions

After determining the number of PBHs that should be simulated within the LOS
for a given field, $N_{\text{PBH}}$, we assign each PBH a distance along the
LOS, $r$, using inverse transform sampling (e.g., Press et al., 1988) to
assign each PBH a galactocentric radius based on a numerical approximation of
the cumulative distribution function of the normalized dark matter density
function. Each PBH galactocentric radius is then converted to a heliocentric
distance along the LOS, $r$, which falls directly on the LOS (see left panel
of Figure 2 for a single PBH example).

Each PBH then needs two additional dimensional quantities for a final 3D
position within the LOS cylinder. We begin by projecting each PBH onto a
circular cross-section of the LOS cylinder at sampled distance $r$, then,
assuming an isotropic distribution in projection, we can sample radial and
angular component values for each PBH (see middle panel of Figure 2 for a
single PBH example). For each PBH, a radial distance (from the center axis of
the LOS cylinder to the projected field radius, $R$) is sampled proportional
to the area of the cylindrical cross-section ($\pi\hat{R}^{2}$), while an
angular component is uniformly sampled between [0, 2$\pi$].

The result of this sampling process888This cumulative integration technique
neglects correlation between PBH mass and PBH location in the dark matter halo
(as might be introduced through dynamical friction forces on an extended PBH
mass spectrum). is a cylinder along the LOS populated with the number of PBHs
determined by Equation 5, following that of the halo profile in Equation 4. We
performed convergence checks on the discritized space and find that as you
pass over the galactic center, distances within $\approx 0.157$ $kpc$ of the
center cause the density to converge, so for our simulations we do not go
below this distance radially from the galactic center. Before velocity
assignment, we mask out PBHs that cannot be observed by cutting out PBHs which
fall outside of the survey light cone (i.e., solid angle). An example of this
cut is shown in the right panel of Figure 2.

###  3.4 PBH Velocities

Each PBH is assigned a mean speed $\langle\text{v}\rangle$ given its
galactocentric radius, using the Eddington formalism (Eddington, 1916) for the
halo velocity as derived in Lacroix et al. (2018). The Eddington inversion
model (which relates total gravitational potential and mass density to the
dark matter phase-space distribution) deviates from the standard halo model
(SHM; in which velocities follow a truncated Gaussian), especially at smaller
Galactic radii (see “SHM” curve on lower panel of Figure 1). However, Lacroix
et al. (2018) think the Eddington model more representative of PBH velocities,
thus, we use the NFW parameterization of this relationship (illustrated by the
$\gamma=1$ panel on Figure 1).

Next, we assume a Maxwellian velocity distribution, and for each PBH we use
its assigned $\langle\text{v}\rangle$ as the distribution mean, to get the
Maxwell distribution parameter:

$a=\frac{1}{2}\langle\text{v}\rangle\sqrt{\frac{\pi}{2}},$

then calculate the probability density function:

Figure 3: Velocity magnitude distributions for selected PBH radii, at three
different $\gamma$ values.

$f(v)=\sqrt{\frac{2}{\pi}}\frac{v^{2}e^{-\frac{v^{2}}{2a^{2}}}}{a^{3}},$ (6)

in which $v$ is the root-mean square (RMS) velocity. We uniformly sample a
single value from the cumulative distribution function (again, using inverse
transform sampling):

$F(v)=\int
f(v)=\text{erf}\left(\frac{v}{a\sqrt{2}}\right)-\sqrt{\frac{2}{\pi}}\frac{ve^{\frac{-v^{2}}{2a^{2}}}}{a},$
(7)

where _erf_ is the error function. Each PBH is assigned the corresponding
sampled RMS velocity.

We let $v$ take values between 0 and the escape velocity $v_{esc}$ = 550 $km$
$s^{-1}$, which is consistent with current estimates, such as those from the
Radial Velocity Experiment (RAVE) survey (Smith et al., 2007; Piffl et al.,
2014). For PBHs with $\langle\text{v}\rangle<250$ $km$ $s^{-1}$, this
$v_{esc}$ cut filters out $<$ 1% of the PBH distribution, while PBHs with mean
speeds closer to the maximum, $\langle\text{v}\rangle$ $\sim{300}$ $km$
$s^{-1}$, equates to a loss of $\sim$3.5%. Because the number of PBHs in the
LOS cylinder was previously selected to match the total mass of the cylinder,
the number of PBH lensing events will be biased slightly low. Figure 3
illustrates $f(v)$ and $F(v)$ distribution examples (in the main figures and
subfigures, respectively) for a few selected galactocentric PBH radii and
$\gamma$ values.

Following similar projection methods as described in §3.3, we take our
velocity magnitude $v$, and randomly sample two angular components for
latitudinal and longitudinal velocity components. The longitudinal component
is sampled uniformly between [0, $2\pi$], while the latitudinal component is
sampled uniformly on the unit sphere (i.e., uniform in $\cos({\rm latitude})$
between [$-1$, $1$]). It is worth noting that we assume a Maxwellian
distribution for PBH velocities as a computational shortcut, but use mean
speeds following the Eddington inversion results from Lacroix et al. (2018),
which assume velocities that deviate from the truncated Maxwellian of the SHM.
However, we do not expect this simplifying assumption to be a dominant source
of uncertainty in our analysis.

###  3.5 Other PBH Parameters

In addition to six-dimensional kinematics, PopSyCLE requires additional
parameters for microlensing survey simulation. Because PBHs form in the early
universe and not through astrophysical formation channels, they do not have a
ZAMS mass. However, to ensure smooth integration into PopSyCLE, the PBH mass
at the time of the simulation is assigned as the ZAMS mass. We add a
population type for PBHs dubbed “dark matter halo” (in addition to those
supplied by Galaxia) and assign an age of 13.8 Gyr. In addition to the
PopSyCLE remnant IDs of 0 (star), 101 (WD), 102 (NS), and 103 (BH), a new
remnant ID of 104 is assigned to PBHs. Magnitudes and photometric information
does not need to be considered, as PBHs are “dark”.

##  4 Microlensing Simulations

Each of our mock surveys closely resemble observational criteria for the
survey being simulated. The PopSyCLE survey parameters used are listed in
Table 2 and include: duration (length of the survey), cadence (time between
each observation), $\theta_{blend}$ (seeing-disk radius) and
$u_{0,\text{max}}$ (maximum impact parameter) as described in §2.3, area
(surface area of the survey field-of-view projected on the sky), limiting
magnitude of the survey, optical filter for the survey, $\Delta_{mag}$ (“bump
magnitude”; the difference between baseline magnitude and peak magnitude, used
to remove low-amplitude events), and $t_{e}$ (the Einstein crossing time
described in §1). For this work we assume an NFW dark matter profile, and use
$\gamma=1$ when generating the PBH population.

PopSyCLE cadence and actual observational cadence do not exactly correlate, as
PopSyCLE spreads out observations over the user-specified duration, and with
the user-specified cadence, not necessarily accounting for daytime or seasonal
gaps. Our OGLE cadences are derived from the reported observations per field
each night, and Roman cadences are determined by spreading the expected
microlensing survey over the duration of telescope operation, meaning we
neglect these daytime and seasonal gaps as well (discussed further in §4.1 and
§4.2). This choice of cadence may lead to underestimating the number of events
with timescales shorter than the specified cadence, but we find this affect
negligible.

Table 2: PopSyCLE Survey Parameters

| Mock OGLE-IV | Mock
---|---|---
Parameter | EWS Cut | Mróz19 Cut | Roman
Duration | 2920 days | 2920 days | 1825 days
Cadence∗ |  $\sim$2.5 days† |  $\sim$2.5 days† |  $\sim$3 days
$\theta_{blend}$ | 0.65” | 0.65” | 0.09”
$u_{0,\text{max}}$ |  $\leq$ 2 |  $\leq$ 1 |  $\leq$ 2
Area | 0.34 $\text{deg}^{2}$ | 0.34 $\text{deg}^{2}$ | 0.16 $\text{deg}^{2}$
Min. Mag.‡ |  $\leq$ 21 |  $\leq$ 21 |  $\leq$ 26
Filter | I | I | H
$\Delta_{mag}$ |  $\geq$ 0.1 | – |  $\geq$ 0.1
$t_{E}$ | – | $0.5\leq t_{E}\leq 300$ | –

Note. — ∗PopSyCLE defined cadence, which does not account for daytime gaps or
survey gaps.

†Fields BLG-500 and BLG-611 use a PopSyCLE survey cadence of 10 days.

‡Mock EWS and Mock Roman use minimum baseline magnitude. Mock Mróz19 uses
minimum source magnitude.

For each survey, we specify a set of fields to run PopSyCLE on. Fields and
locations (with the area from Table 2 drawn to scale) are detailed in Figure
4. Each field in Figure 4 is run one time at minimum, where light pink denotes
each simulated OGLE-IV field that is run only once, dark pink denotes each
OGLE-IV field that is run over three random seeds each, and cyan denotes each
simulated Roman field that is run over three random seeds each. The dark pink
and cyan fields are 9 fields taken from Lam2020, in which each field is
centered at a Galactic longitude and Galactic latitude ($l$ and $b$,
respectively, noted on the right panel of Figure 4) that either corresponds to
the center of a bulge field from the Optical Gravitational Lensing Experiment
(OGLE) survey (Udalski et al., 1992) or the Microlensing Observations in
Astrophysics (MOA) survey (Muraki et al., 1999). Only a 0.34 $\text{deg}^{2}$
or 0.16 $\text{deg}^{2}$ (for OGLE-IV and Roman, respectively) portion of each
survey footprints is simulated, then extrapolation is done (when necessary)
with the assumption that stellar density is constant over each field. Fields
with a label in bold on Figure 4 are “high-cadence” OGLE-IV fields (as defined
in Mróz et al. (2019)), while everything else is considered a low cadence
field.

Figure 4: _Left:_ Simulated survey areas (drawn to scale) plotted over the
OGLE-IV survey footprint. The background color scale represents the number of
stars in a given field, and the number on each field is the OGLE bulge (BLG)
field number. Each light pink field was simulated with one seed, while the
dark pink fields were simulated with three seeds each. The fields that have
their numbers in bold are the high cadence OGLE-IV fields (as defined in Mróz
et al. (2019)), and everything else is considered a low cadence field.
_Right:_ The name and Galactic coordinates for each field that was simulated
three times (each field is centered on the given coordinates, and has a
simulated area specified in Table 2).
∗BLG-500 is both a 0.34 $\text{deg}^{2}$ Mock OGLE-IV field and 0.16
$\text{deg}^{2}$ Mock Roman field.
†OGLE-2014-BLG-0613a is centered on a microlensing event reported by OGLE EWS.

###  4.1 Mock OGLE-IV Simulations

The OGLE survey (Udalski et al., 1992) is designed to detect microlensing
events in the Milky Way’s bulge and Magellanic Clouds by frequently surveying
the sky to look for microlensing signatures. OGLE phase 1 (OGLE-I) began in
1992, and recently finished the fourth phase (OGLE-IV), making observations
with the 1.3 $m$ Warsaw telescope at the Las Campanas Observatory in Chile
over an 8 year period (Udalski et al., 2015).

Our “Mock OGLE-IV” surveys use the corresponding values from Table 2 to
simulate the portions of the OGLE microlensing survey shown in light pink and
dark pink on Figure 4. We adopt these field parameters from the “Mock
EWS999EWS = Early Warning System. The EWS reduces OGLE data in real time to
alert on potential microlensing events. We base parameter estimates on
reported numbers from 2016-2018.” and “Mock Mróz19” surveys from Lam2020. The
Mróz19 cut returns OGLE-IV events that might be detectable, and that are also
of sufficient signal-to-noise such that their lensing parameters could be
reasonably estimated. The EWS cut includes additional events that could
conceivably be detected with OGLE, but that one may not be able to get
reliable lensing parameter estimates for.

Using Table 6 from Mróz et al. (2019) we calculate the average cadence for
each OGLE field we simulate, then, averaging over all fields, we decide on a
simulated survey cadence of 1 observation every $\sim$2.5 days (with the
exception of BLG-500 and BLG-611 which use a cadence of 10 days due to
computational constraints). Cadence affects the sensitivity of shorter scale
events (in which $t_{E}$ is less than the time between observations) because
the event will fall entirely between observations. While our simulated
cadences are longer than some of the actual OGLE-IV field cadences, Lam2020
show that a cadence of $\leq$ 10 days is sufficient to provide an accurate
prediction of the $t_{E}$ distribution peak.

All quantities reported from our PopSyCLE simulations are considered as
corrected, as PopSyCLE does not account for events that fall between gaps,
noisy observations, etc.

To calculate event parameters with and without PBHs, we first determine the
number of detectable stars in each field, $N_{s}$, by cutting the full object
catalog from Galaxia on the limiting magnitude of the survey, $I=21$. To
estimate the number of microlensing events (not including PBHs) in each field
for each OGLE-IV cut, including relevant scaling and correcting for the actual
survey efficiencies, we use the following equation:

$N_{\text{ev}}=N^{sim}_{\text{ev}}\left(\frac{1.4}{0.34}\right)\left(\frac{N^{M19}_{star}}{N^{M19}_{s}}\right)\langle\epsilon\rangle$
(8)

where $N^{sim}_{\text{ev}}$ is the number of simulated microlensing events
(excluding PBHs) output from PopSyCLE (after the cuts from Table 2 are
applied), (1.4/0.34) scales our simulated number of microlensing events up to
the full field area, $N^{M19}_{star}$ is the number of stars being monitored
in that particular OGLE-IV field (values taken from Table 6 of Mróz et al.
(2019)), $N^{M19}_{s}$ is the number of sources in a field that are brighter
than $I=21$ (with blends accounted for; values from Table 7 of Mróz et al.
(2019)), and $\langle\epsilon\rangle$ is the mean efficiency averaged over all
of our simulated microlensing events in that field (i.e., all events in
$N^{sim}_{\text{ev}}$). The quantity ($N^{M19}_{stars}$/$N^{M19}_{s}$) serves
as a luminosity completeness correction factor.

Then, to estimate the number of PBH microlensing events in each field, with
relevant survey and efficiency corrections included, for each OGLE-IV cut, we
use the following equation:

$\begin{multlined}N_{\text{ev, PBH}}=\\\ N^{sim}_{\text{ev,
PBH}}\left(\frac{N^{obs}_{\text{ev}}}{N^{sim}_{\text{ev}}\langle\epsilon\rangle+N^{sim}_{\text{ev,
PBH}}\langle\epsilon^{\prime}\rangle}\right)\langle\epsilon^{\prime}\rangle\end{multlined}N_{\text{ev,
PBH}}=\\\ N^{sim}_{\text{ev,
PBH}}\left(\frac{N^{obs}_{\text{ev}}}{N^{sim}_{\text{ev}}\langle\epsilon\rangle+N^{sim}_{\text{ev,
PBH}}\langle\epsilon^{\prime}\rangle}\right)\langle\epsilon^{\prime}\rangle$
(9)

where $N^{sim}_{\text{ev, PBH}}$ is the number of simulated PBH microlensing
events output from PopSyCLE (after the cuts from Table 2 are applied),
$N^{obs}_{\text{ev}}$ is the number of detected events listed in Table 7 of
Mróz et al. (2019) or the number of events reported between the 2011 and 2018
EWS seasons (depending on which cut we are determining $N_{\text{ev, PBH}}$
for), $N^{sim}_{\text{ev}}$ is the number of simulated microlensing events
(excluding PBHs) output from PopSyCLE (after the cuts from Table 2 are
applied), $\langle\epsilon\rangle$ is the mean efficiency averaged over all of
our simulated microlensing events (excluding PBHs) in that field (i.e., all
events in $N^{sim}_{\text{ev}}$), and $\langle\epsilon^{\prime}\rangle$ is the
mean efficiency averaged over each of our simulated PBH microlensing events in
that field (i.e., all events in $N^{sim}_{\text{ev, PBH}}$).
$N^{sim}_{\text{PBH}}$ and ($N^{sim}_{\text{ev}}+N^{sim}_{\text{ev, PBH}}$)
both have a dependence on field area and corrections due to luminosity
completeness, meaning that both the area and luminosity correction ratios
cancel out and do not need to be considered for these calculations.

Efficiencies are calculated using each event $t_{E}$, with the published
$t_{E}$ vs. efficiency binned values from the OGLE-IV
survey101010https://www.astrouw.edu.pl/ogle/ogle4/
microlensing_maps/eff/. For simulated Mróz19 fields that do not have PBH
events, $\langle\epsilon^{\prime}\rangle$ is taken to be the average
efficiency value across all other fields that do have PBHs,
$\langle\epsilon^{\prime}\rangle=0.268$. For simulated Mróz19 fields that do
not have any events (i.e., the number of NS, BH, WD, Star, and PBH events are
all zero), $\langle\epsilon\rangle$ is taken to be the average efficiency
value across all other fields that have non-zero events,
$\langle\epsilon\rangle=0.172$.

EWS efficiencies are unknown, therefore, we assume
$\langle\epsilon\rangle=\langle\epsilon^{\prime}\rangle$. However, when
calculating the number of events (excluding PBHs) for our EWS simulations, we
repeat the process of finding $\langle\epsilon\rangle$ using the reported
OGLE-IV efficiency values, and fields that do not have any simulated events
get assigned the average efficiency value across all other fields with non-
zero events: 0.169.

The event rate, $\Gamma$ (which excludes PBH events), is calculated by
dividing the number of events without PBHs, $N_{\text{ev}}$, by the number of
stars that fall above the limiting magnitude for the given survey, $N_{s}$,
and then divide by the given survey duration. We report our simulated OGLE-IV
results in Table 3, which compares values from our analysis, with relevant
quantities (i.e. $N_{s}$ and $N_{\text{ev}}$) scaled up from 0.34
$\text{deg}^{2}$ to the full field area, 1.4 $deg^{2}$, for equivalent
comparison to their observational counterparts. We list $N_{s}$, $\Gamma$,
$N_{\text{ev}}$, and $N_{\text{ev, PBH}}$ for our mock Mróz19 and Mock EWS
cuts, and list $N_{s}$, $\Gamma$, and $N_{\text{ev}}$ observational results
reported by Mróz et al. (2019) or the EWS website. $N_{s}$ and $\Gamma$ are
completeness corrected values, while $N_{\text{ev}}$ and $N_{\text{ev, PBH}}$
are not, but are appropriately scaled by the various factors listed in
Equations 8 and 9 to make non-completeness corrected values to compare to
their OGLE-IV observational counterparts. Three values are listed for fields
that are run over three seeds each, and Poisson noise is assumed for all
simulated values (this is because the exact noise profile is hard to
characterize, and our sample of seeds follows Poisson uncertainty).

Table 3: OGLE-IV Event Statistic Comparison

OGLE |  $N_{s}$ [$10^{6}$] |  $\Gamma$ [$10^{-6}$ $\text{star}^{-1}$ $\text{year}^{-1}$] | ${N_{\text{ev}}}$ | ${N_{\text{ev, PBH}}}$
---|---|---|---|---
Field | Mróz19 | M/E | Mróz19 | M | EWS | E | Mróz19 | M | EWS | E | M | E
BLG-500 | 6.78 | 4.71 | 23.9 | 23.5 | 110.9 | 134.42 | 164 | 111.52 | 752.0 | 646.23 | 7.01 | 20.83
|  |  |  | 23.59 |  | 134.77 |  | 112.04 |  | 648.33 | 3.56 | 21.89
|  |  |  | 26.16 |  | 130.46 |  | 123.97 |  | 626.27 | 2.43 | 14.84
BLG-501 | 13.31 | 12.24 | 24.1 | 33.0 | 89.9 | 123.6 | 317 | 253.3 | 1197.0 | 947.3 | 13.0 | 40.9
BLG-502 | 10.02 | 4.67 | 11.0 | 10.3 | 31.8 | 39.7 | 171 | 43.9 | 319.0 | 160.1 | 17.8 | 17.6
BLG-503 | 10.7 | 5.25 | 5.5 | 7.8 | 14.2 | 24.3 | 91 | 33.6 | 152.0 | 97.9 | 11.7 | 13.4
BLG-504 | 11.86 | 4.34 | 16.9 | 18.4 | 77.9 | 98.6 | 225 | 70.6 | 924.0 | 383.9 | 9.4 | 27.0
BLG-505 | 19.02 | 9.88 | 22.2 | 24.2 | 71.2 | 92.5 | 441 | 143.2 | 1354.0 | 550.1 | 9.5 | 37.1
BLG-506 | 12.85 | 5.36 | 16.5 | 20.7 | 57.0 | 88.3 | 247 | 81.7 | 732.0 | 357.5 | 9.5 | 20.1
BLG-507 | 11.22 | 5.27 | 12.3 | 15.4 | 33.3 | 59.3 | 216 | 67.1 | 374.0 | 259.3 | 17.0 | 14.2
BLG-508 | 9.53 | 5.0 | 6.8 | 10.4 | 22.9 | 33.4 | 119 | 45.4 | 218.0 | 145.1 | 9.4 | 14.4
BLG-509 | 9.82 | 4.72 | 4.6 | 6.1 | 13.3 | 20.6 | 79 | 23.8 | 131.0 | 79.0 | 11.0 | 7.8
BLG-510 | 7.67 | 3.47 | 3.6 | 3.0 | 6.6 | 10.2 | 38 | 7.8 | 51.0 | 27.9 | 6.4 | 4.1
BLG-511 | 13.45 | 5.11 | 13.5 | 16.8 | 56.1 | 64.0 | 204 | 64.0 | 755.0 | 244.2 | 4.5 | 23.1
BLG-512 | 17.48 | 8.34 | 14.0 | 18.0 | 53.0 | 63.9 | 280 | 108.4 | 926.0 | 384.9 | 9.9 | 25.2
BLG-513 | 13.71 | 7.87 | 8.5 | 14.6 | 26.8 | 54.0 | 213 | 78.1 | 368.0 | 287.5 | 8.7 | 14.1
BLG-514 | 10.27 | 5.44 | 6.2 | 10.2 | 18.9 | 41.8 | 108 | 45.9 | 194.0 | 179.9 | 9.8 | 8.0
BLG-515 | 8.1 | 4.32 | 5.7 | 8.1 | 19.6 | 24.8 | 80 | 38.3 | 159.0 | 115.1 | 1.7 | 10.7
BLG-516 | 7.61 | 3.76 | 4.6 | 7.0 | 7.6 | 15.0 | 48 | 17.8 | 58.0 | 41.0 | 0.0 | 1.5
BLG-517 | 5.96 | 3.18 | 2.6 | 4.5 | 0.3 | 11.3 | 12 | 9.6 | 2.0 | 23.1 | 2.1 | 0.2
BLG-518 | 12.34 | 4.72 | 7.3 | 11.1 | 25.2 | 36.1 | 160 | 38.9 | 311.0 | 126.7 | 4.4 | 12.6
BLG-519 | 12.93 | 6.3 | 9.3 | 10.5 | 29.3 | 36.1 | 223 | 51.4 | 379.0 | 177.3 | 11.3 | 21.1
BLG-520 | 11.63 | 5.94 | 5.6 | 9.1 | 20.3 | 32.1 | 124 | 54.5 | 236.0 | 187.7 | 13.8 | 9.2
BLG-521 | 9.18 | 4.95 | 4.5 | 9.2 | 12.2 | 26.4 | 69 | 40.9 | 112.0 | 127.0 | 0.0 | 7.8
BLG-522 | 8.62 | 4.24 | 7.4 | 7.8 | 19.4 | 21.7 | 108 | 40.5 | 167.0 | 107.8 | 10.4 | 12.9
BLG-523 | 10.4 | 4.84 | 6.0 | 7.0 | 16.5 | 21.1 | 99 | 26.6 | 172.0 | 76.3 | 10.0 | 9.8
BLG-524 | 9.66 | 4.49 | 6.0 | 6.5 | 7.9 | 20.2 | 75 | 17.0 | 76.0 | 50.7 | 9.4 | 4.9
BLG-525 | 7.84 | 4.04 | 3.2 | 4.3 | 6.5 | 15.4 | 36 | 16.4 | 51.0 | 55.1 | 7.9 | 5.6
BLG-526 | 5.25 | 3.34 | 2.3 | 5.1 | 4.4 | 13.5 | 20 | 34.0 | 23.0 | 88.8 | 3.2 | 1.7
BLG-527 | 6.36 | 2.85 | 5.5 | 2.89 | 5.3 | 8.30 | 48 | 10.63 | 34.0 | 33.64 | 8.44 | 1.42
|  |  |  | 2.89 |  | 9.58 |  | 10.63 |  | 38.76 | 16.70 | 5.40
|  |  |  | 3.06 |  | 9.55 |  | 11.29 |  | 38.76 | 8.03 | 1.24
BLG-528 | 7.39 | 2.99 | 5.5 | 4.1 | 5.7 | 11.2 | 50 | 12.8 | 42.0 | 33.1 | 0.0 | 0.6
BLG-529 | 7.79 | 3.2 | 5.6 | 4.7 | 3.6 | 11.9 | 39 | 7.0 | 28.0 | 15.3 | 1.9 | 1.8
BLG-530 | 7.37 | 3.24 | 8.5 | 2.4 | 4.9 | 9.9 | 54 | 3.5 | 36.0 | 14.2 | 10.6 | 2.7
BLG-531 | 6.85 | 3.05 | 5.1 | 3.0 | 4.8 | 9.5 | 34 | 6.4 | 33.0 | 20.9 | 6.9 | 5.0
BLG-532 | 6.02 | 2.82 | 2.7 | 2.6 | 2.0 | 6.9 | 23 | 8.2 | 12.0 | 24.7 | 4.8 | 1.4
BLG-533 | 1.02 | 0.21 | 7.1 | 0.0 | 0.0 | 0.0 | 1 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0
BLG-534 | 9.06 | 4.51 | 17.2 | 19.8 | 72.3 | 101.1 | 176 | 93.2 | 655.0 | 481.0 | 3.6 | 19.3
BLG-535 | 7.57 | 3.58 | 15.2 | 14.4 | 43.7 | 67.1 | 189 | 46.0 | 331.0 | 217.5 | 9.4 | 9.0
BLG-536 | 5.21 | 2.47 | 3.8 | 1.9 | 0.8 | 4.8 | 15 | 0.7 | 4.0 | 1.7 | 2.1 | 0.6
BLG-539 | 4.52 | 2.37 | 4.0 | 2.8 | 0.4 | 6.1 | 13 | 8.3 | 2.0 | 18.4 | 3.9 | 0.4
BLG-543 | 4.76 | 1.62 | 5.9 | 2.9 | 7.1 | 10.8 | 36 | 5.0 | 34.0 | 18.9 | 9.3 | 2.8
BLG-544 | 5.66 | 2.87 | 6.1 | 5.2 | 6.0 | 15.0 | 40 | 11.7 | 34.0 | 33.2 | 0.0 | 1.9
BLG-545 | 6.47 | 3.05 | 6.1 | 6.4 | 14.1 | 17.6 | 57 | 14.8 | 91.0 | 40.1 | 7.7 | 8.0
BLG-546 | 4.25 | 2.47 | 3.4 | 3.8 | 0.5 | 7.7 | 12 | 7.5 | 2.0 | 14.6 | 0.0 | 0.2
BLG-547 | 3.67 | 2.07 | 1.3 | 1.99 | 0.3 | 5.46 | 3 | 2.57 | 1.0 | 7.80 | 0.88 | 0.19
|  |  |  | 2.23 |  | 6.21 |  | 2.89 |  | 8.86 | 1.07 | 0.14
|  |  |  | 1.74 |  | 4.47 |  | 2.25 |  | 6.38 | 0.00 | 0.18
BLG-566 | 3.03 | 1.28 | 0.9 | 2.4 | 0.0 | 4.0 | 2 | 5.0 | 0.0 | 7.8 | 0.3 | 0.0
BLG-573 | 7.16 | 3.83 | 3.7 | 8.1 | 4.9 | 19.5 | 33 | 14.5 | 35.0 | 35.3 | 4.9 | 3.5
BLG-580 | 8.47 | 3.19 | 8.0 | 9.2 | 29.2 | 31.9 | 116 | 22.6 | 247.0 | 76.7 | 7.2 | 9.6
BLG-588 | 5.39 | 3.19 | 2.7 | 6.5 | 3.9 | 12.4 | 23 | 13.6 | 21.0 | 27.9 | 3.6 | 2.4
BLG-597 | 4.28 | 2.57 | 4.0 | 2.2 | 3.7 | 7.4 | 27 | 5.3 | 16.0 | 17.5 | 2.9 | 1.9
BLG-599 | 7.8 | 3.69 | 3.9 | 4.5 | 3.8 | 11.3 | 41 | 14.6 | 30.0 | 34.5 | 5.0 | 2.4
BLG-600 | 5.29 | 2.5 | 2.9 | 3.5 | 2.3 | 6.0 | 18 | 9.0 | 12.0 | 14.5 | 1.4 | 0.8
BLG-603 | 10.02 | 5.24 | 10.6 | 11.1 | 29.8 | 36.6 | 167 | 50.2 | 299.0 | 168.8 | 5.8 | 19.5
BLG-604 | 9.8 | 4.42 | 7.0 | 6.4 | 12.4 | 20.3 | 100 | 19.7 | 122.0 | 65.6 | 7.0 | 9.1
BLG-605 | 6.13 | 2.6 | 3.9 | 1.8 | 4.2 | 8.1 | 29 | 3.9 | 26.0 | 15.8 | 12.8 | 3.8
BLG-606 | 4.48 | 2.32 | 3.8 | 1.8 | 0.2 | 3.8 | 11 | 2.8 | 1.0 | 6.2 | 1.6 | 0.1
BLG-609 | 6.77 | 3.05 | 8.4 | 7.9 | 21.6 | 36.6 | 81 | 17.9 | 146.0 | 92.5 | 5.2 | 6.4
BLG-610 | 5.78 | 2.86 | 6.2 | 5.9 | 9.3 | 19.6 | 38 | 12.3 | 54.0 | 42.3 | 0.0 | 1.9
BLG-611 | 6.93 | 5.14 | 16.2 | 18.79 | 68.0 | 75.07 | 158 | 154.64 | 471.0 | 624.8 | 8.34 | 23.25
|  |  |  | 18.70 |  | 81.09 |  | 153.82 |  | 674.72 | 5.69 | 17.88
|  |  |  | 16.88 |  | 78.52 |  | 139.01 |  | 653.92 | 7.27 | 20.63
BLG-612 | 7.34 | 4.14 | 10.1 | 14.3 | 21.0 | 56.5 | 87 | 30.2 | 154.0 | 117.9 | 4.1 | 5.2
BLG-613 | 7.99 | 3.96 | 12.4 | 15.1 | 24.8 | 68.0 | 119 | 40.4 | 198.0 | 181.3 | 11.4 | 11.1
BLG-614 | 6.32 | 3.33 | 5.5 | 10.4 | 13.3 | 31.7 | 43 | 23.5 | 84.0 | 72.2 | 3.0 | 5.4
BLG-615 | 8.09 | 3.7 | 6.8 | 13.5 | 19.0 | 40.8 | 74 | 31.0 | 154.0 | 94.1 | 4.8 | 5.6
BLG-616 | 5.08 | 2.75 | 3.4 | 3.7 | 0.2 | 11.4 | 12 | 5.5 | 1.0 | 14.5 | 1.5 | 0.1
BLG-617 | 6.32 | 3.08 | 5.2 | 6.20 | 10.6 | 16.75 | 43 | 18.79 | 67.0 | 50.63 | 5.45 | 6.64
|  |  |  | 3.68 |  | 13.90 |  | 11.17 |  | 42.02 | 2.47 | 7.20
|  |  |  | 5.19 |  | 14.57 |  | 15.75 |  | 44.05 | 3.43 | 3.64
BLG-619 | 4.47 | 2.57 | 3.3 | 6.8 | 0.0 | 19.6 | 3 | 1.9 | 0.0 | 5.2 | 0.7 | 0.0
BLG-621 | 5.15 | 2.24 | 6.8 | 5.1 | 8.3 | 11.0 | 44 | 11.5 | 43.0 | 22.9 | 8.4 | 4.1
BLG-622 | 4.76 | 2.43 | 3.6 | 5.3 | 6.7 | 15.7 | 25 | 13.0 | 32.0 | 35.4 | 3.6 | 3.1
BLG-624 | 6.02 | 2.54 | 4.3 | 3.0 | 6.3 | 8.7 | 31 | 6.7 | 38.0 | 21.0 | 2.9 | 4.0
BLG-625 | 6.47 | 2.6 | 5.5 | 5.1 | 7.9 | 15.4 | 46 | 11.6 | 51.0 | 30.7 | 10.6 | 5.3
BLG-626 | 8.32 | 4.3 | 8.3 | 9.2 | 15.9 | 32.9 | 92 | 31.1 | 132.0 | 112.8 | 7.3 | 4.2
BLG-629 | 4.56 | 2.09 | 3.4 | 2.71 | 0.4 | 7.65 | 11 | 4.71 | 2.0 | 12.75 | 3.9 | 0.37
|  |  |  | 2.96 |  | 7.15 |  | 5.14 |  | 11.93 | 0.00 | 0.19
|  |  |  | 2.96 |  | 7.15 |  | 5.14 |  | 11.93 | 1.85 | 0.24
BLG-630 | 5.11 | 2.21 | 4.4 | 5.4 | 7.4 | 11.2 | 30 | 11.9 | 38.0 | 25.5 | 8.4 | 5.4
BLG-631 | 7.09 | 2.99 | 7.1 | 6.9 | 8.7 | 20.3 | 61 | 16.2 | 62.0 | 47.6 | 8.9 | 3.5
BLG-632 | 9.59 | 5.43 | 9.0 | 10.2 | 22.8 | 33.1 | 128 | 45.8 | 219.0 | 155.7 | 5.6 | 10.2
BLG-633 | 6.3 | 2.78 | 13.1 | 12.2 | 23.3 | 55.5 | 106 | 27.4 | 147.0 | 119.9 | 4.8 | 6.1
BLG-636 | 4.26 | 2.26 | 5.1 | 5.5 | 0.5 | 11.9 | 14 | 8.1 | 2.0 | 16.5 | 0.9 | 0.1
BLG-637 | 5.84 | 2.56 | 5.2 | 7.4 | 8.6 | 19.5 | 40 | 14.9 | 50.0 | 41.4 | 1.6 | 2.0
BLG-638 | 7.96 | 3.66 | 6.9 | 10.1 | 11.1 | 32.2 | 64 | 22.4 | 88.0 | 73.9 | 5.1 | 3.7
BLG-639 | 5.08 | 1.89 | 10.4 | 8.2 | 19.5 | 34.8 | 70 | 13.4 | 99.0 | 54.2 | 3.2 | 3.7
BLG-641 | 5.33 | 3.05 | 4.9 | 9.6 | 4.5 | 25.8 | 31 | 22.7 | 24.0 | 60.0 | 2.4 | 1.2
BLG-642 | 3.22 | 1.37 | 8.7 | 4.1 | 8.7 | 15.8 | 36 | 5.0 | 28.0 | 19.6 | 0.0 | 1.3
BLG-643 | 2.11 | 1.53 | 5.7 | 4.7 | 8.5 | 16.1 | 17 | 8.4 | 18.0 | 29.8 | 0.0 | 0.7
BLG-644 | 1.02 | 0.34 | 0.0 | 0.0 | 0.0 | 0.0 | 0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0
BLG-645 | 3.33 | 0.66 | 13.6 | 5.4 | 32.1 | 13.9 | 63 | 3.4 | 107.0 | 9.4 | 0.0 | 5.6
BLG-646 | 5.09 | 0.62 | 8.3 | 8.3 | 19.8 | 19.9 | 57 | 4.6 | 101.0 | 9.5 | 0.0 | 0.0
BLG-647 | 0.68 | 0.21 | 0.0 | 0.0 | 0.0 | 2.4 | 0 | 0.0 | 0.0 | 0.1 | 0.0 | 0.0
BLG-648 | 2.86 | 1.74 | 18.3 | 8.0 | 42.3 | 45.5 | 67 | 12.1 | 121.0 | 67.7 | 6.9 | 4.5
BLG-649 | 0.93 | 0.18 | 0.0 | 0.0 | 0.0 | 0.0 | 0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0
BLG-650 | 0.91 | 0.18 | 0.0 | 0.0 | 0.0 | 0.0 | 0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0
BLG-651 | 1.56 | 0.77 | 0.0 | 8.0 | 0.0 | 20.8 | 0 | 1.1 | 0.0 | 2.9 | 0.0 | 0.0
BLG-652 | 6.55 | 3.62 | 14.2 | 16.2 | 22.3 | 84.5 | 100 | 33.9 | 146.0 | 177.9 | 7.3 | 5.0
BLG-653 | 5.37 | 4.56 | 20.1 | 20.1 | 39.9 | 109.4 | 127 | 64.1 | 214.0 | 362.3 | 4.7 | 5.0
BLG-654 | 4.92 | 2.35 | 16.7 | 15.4 | 48.8 | 84.9 | 110 | 27.8 | 240.0 | 159.4 | 5.6 | 9.0
BLG-655 | 2.13 | 0.25 | 9.5 | 0.0 | 0.0 | 2.1 | 3 | 0.0 | 0.0 | 0.1 | 0.0 | 0.0
BLG-657 | 4.43 | 2.38 | 2.6 | 1.9 | 0.7 | 7.1 | 6 | 3.4 | 3.0 | 11.0 | 0.9 | 0.5
BLG-659 | 1.05 | 0.73 | 0.0 | 8.5 | 0.0 | 21.1 | 0 | 0.9 | 0.0 | 1.6 | 0.0 | 0.0
BLG-660 | 4.63 | 2.07 | 23.7 | 12.2 | 35.4 | 50.6 | 94 | 15.7 | 164.0 | 72.4 | 10.6 | 10.6
BLG-661 | 4.22 | 2.78 | 15.2 | 7.6 | 25.6 | 44.6 | 80 | 14.7 | 108.0 | 88.9 | 0.0 | 5.5
BLG-662 | 6.16 | 1.79 | 12.2 | 12.1 | 20.0 | 56.0 | 87 | 12.8 | 123.0 | 66.2 | 10.4 | 7.7
BLG-663 | 2.15 | 1.35 | 3.7 | 6.9 | 2.3 | 23.3 | 4 | 6.6 | 5.0 | 22.4 | 0.0 | 0.0
BLG-664 | 1.85 | 0.64 | 5.7 | 1.6 | 4.3 | 8.8 | 6 | 1.1 | 8.0 | 3.1 | 2.1 | 0.7
BLG-665 | 2.72 | 1.34 | 1.2 | 5.8 | 0.4 | 16.9 | 1 | 1.2 | 1.0 | 3.3 | 0.1 | 0.0
BLG-666 | 2.04 | 1.9 | 9.8 | 7.6 | 0.0 | 34.6 | 3 | 2.1 | 0.0 | 12.3 | 0.3 | 0.0
BLG-667 | 2.41 | 1.72 | 11.3 | 11.4 | 15.8 | 43.8 | 27 | 17.8 | 38.0 | 61.7 | 2.4 | 1.3
BLG-668 | 1.93 | 0.48 | 8.8 | 4.3 | 5.2 | 10.8 | 15 | 1.4 | 10.0 | 3.7 | 0.0 | 0.0
BLG-670 | 1.65 | 1.5 | 10.7 | 4.5 | 10.9 | 14.8 | 19 | 3.8 | 18.0 | 21.5 | 0.0 | 1.2
BLG-672 | 2.45 | 2.88 | 5.6 | 5.36 | 8.6 | 17.14 | 17 | 12.03 | 21.0 | 34.56 | 3.05 | 1.43
|  |  |  | 4.11 |  | 20.2 |  | 9.22 |  | 40.68 | 3.77 | 1.06
|  |  |  | 5.36 |  | 18.6 |  | 12.03 |  | 37.44 | 3.05 | 2.49
BLG-675 | 5.64 | 5.5 | 26.5 | 22.2 | 53.2 | 118.5 | 160 | 91.9 | 300.0 | 491.8 | 9.3 | 9.8
BLG-676 | 4.58 | 4.27 | 3.9 | 21.7 | 0.0 | 136.9 | 1 | 5.9 | 0.0 | 38.3 | 0.1 | 0.0
BLG-677 | 0.96 | 0.84 | 0.0 | 5.5 | 0.0 | 15.9 | 0 | 0.3 | 0.0 | 1.1 | 0.0 | 0.0
BLG-680 | 3.57 | 2.0 | 5.0 | 3.3 | 0.8 | 11.9 | 8 | 1.8 | 3.0 | 9.5 | 1.2 | 0.3
BLG-683 | 4.28 | 3.51 | 26.2 | 16.0 | 1.6 | 74.2 | 21 | 7.6 | 7.0 | 37.7 | 1.2 | 0.3
BLG-705 | 1.1 | 0.65 | 0.0 | 0.0 | 0.0 | 0.0 | 0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0
BLG-706 | 0.94 | 0.54 | 3.4 | 0.0 | 0.0 | 0.0 | 1 | 0.0 | 0.0 | 0.0 | 1.0 | 0.0
BLG-707 | 0.81 | 0.5 | 0.0 | 0.0 | 0.0 | 1.0 | 0 | 0.0 | 0.0 | 0.5 | 0.0 | 0.0
BLG-708 | 0.84 | 0.42 | 0.0 | 1.2 | 0.0 | 1.2 | 0 | 0.5 | 0.0 | 0.5 | 0.0 | 0.0
BLG-709 | 0.61 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0
BLG-710 | 0.59 | 0.0 | 2.6 | 0.0 | 0.0 | 0.0 | 1 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0
BLG-711 | 0.48 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0
BLG-714 | 5.84 | 2.84 | 11.0 | 11.9 | 7.4 | 57.2 | 54 | 14.8 | 43.0 | 74.8 | 5.3 | 1.6
BLG-715 | 3.42 | 1.61 | 10.2 | 6.1 | 14.9 | 30.1 | 36 | 8.0 | 51.0 | 31.1 | 8.0 | 3.5
BLG-717 | 2.06 | 0.77 | 3.9 | 2.7 | 1.0 | 8.6 | 6 | 1.6 | 2.0 | 3.5 | 1.7 | 0.1

Note. — Microlensing estimates for the following, all scaled to the full OGLE-
IV field area of 1.4 $\text{deg}^{2}$: (1) “Mróz19”: observational results
from Mróz et al. (2019), (2) “M”: our “Mróz19 Cut” simulation results, (3)
“EWS”: Early Warning System events triggered between the 2011 and 2018 OGLE-IV
seasons, and (4) “E”: our “EWS Cut” simulation results. $N_{s}$ is the
completeness corrected number of stars brighter than the magnitude limit of
the survey ($I=21$), $\Gamma$ is the completeness corrected event rate
(excluding PBH events), $N_{\text{ev}}$ is the appropriately scaled number of
events excluding PBHs (using Equation 8), and $N_{\text{ev, PBH}}$ is the
number of scaled simulated PBH microlensing events (using Equation 9). Columns
with three values were run over three seeds each, while only one seed is run
for each remaining OGLE-IV estimates. Poisson uncertainty is assumed for all
simulated values.

Figure 5: OGLE-IV PBH microlensing event ($N_{\text{ev, PBH}}$) estimation as
a function of dark matter fraction ($f_{\text{DM}}$). _Upper Left:_ PBH event
estimation per OGLE-IV field, for our Mock Mróz19 simulations. High cadence
and low cadence OGLE-IV BLG fields are shown with bold colored lines, and thin
gray lines, respectively. _Upper Right:_ Cumulative PBH event estimate over
all fields, with the $f_{\text{DM}}$=1 estimate reported. _Lower Left:_ same
as upper left panel but for our EWS cut. _Lower Right:_ same as upper right
panel but for our EWS cut.

Table 3 assumes that 100% of dark matter consists of PBHs. Figure 5 then steps
further to illustrate the breakdown of the number of predicted PBH events per
field, as a function of $f_{\text{DM}}$. The high-cadence OGLE-IV fields (from
Mróz et al. (2019), labeled in bold on the left panel of Figure 4) are plotted
in color on Figure 5, while each low-cadence field is plotted as a gray line.
The right panels of Figure 5 show the cumulative number of PBH events over
each field, giving the number of PBH events estimated per survey. In all
panels, the number of PBH events is determined using Equation 9.

We combine results across all simulated OGLE-IV fields with the EWS cut and
plot the resulting $\pi_{E}$, $t_{E}$ distribution for all events (see left
panel of Figure 6). We then compare the $t_{E}$ distributions between PBHs and
all populations except PBHs, and then compare the fractional contributions
between compact object populations. The top panel of Figure 7 shows our
simulated EWS $t_{E}$ results. The upper row shows the $t_{E}$ distributions
for all populations combined, all populations except PBHs, and PBHs alone,
while the bottom row shows the relative contribution of each object type to
the total number of microlensing events, as a function of $t_{E}$. The columns
of each panel correspond to $f_{\text{DM}}=0.01$, 0.1, and 1, respectively.
Both Figure 6 and Figure 7 include values from one seed run per field, and are
generated before applying any of the scaling factors given in Equations 8 or 9
(i.e., these are the total subset of events output from PopSyCLE with only the
cuts listed in Table 2).

Figure 6: Simulated microlensing parallax ($\pi_{E}$) and Einstein crossing
time ($t_{E}$) distributions. _Left_ : Mock OGLE-IV simulation distribution
after applying the Mock EWS observational cuts from Table 2. One seed run from
each field is represented on this panel. _Right_ : Same as left panel, this
time for Mock Roman. The arrows on the bottom of the left panel demonstrate
the direction that points on this plot would shift, as a function of
increasing the quantities specified. Only a small portion of the Roman survey
footprint is simulated, and therefore only our simulated regions are
represented within the right panel. This panel reflects events from all three
seed runs per field. Neither panel considers any scaling or correction
factors, and are considered the number of PopSyCLE detections given only the
cuts from Table 2. This figure should be used for distribution demonstration
purposes only.

###  4.2 Mock Roman Simulations

The Nancy Grace Roman Space Telescope (abbreviated as “Roman” and previously
named WFIRST; Spergel et al., 2015), is a space-based telescope set to launch
in 2027. During Roman’s $\sim$5 year operation, Roman will conduct a
microlensing survey towards the galactic bulge consisting of a seven-field
mosaic, over six 72 day seasons (Penny et al., 2019; Johnson et al., 2020).
During each microlensing survey season, images are estimated to be taken every
$\sim$15 minutes from one of the seven fields in the “Cycle 7” configuration
(Penny et al., 2019; Johnson et al., 2020). The survey footprint of this
arrangement (See Figure 7 of Penny et al., 2019) spans the region between
$l\approx(-.5,1.5)$ and $b\approx(-.5,-2)$).

Figure 7: $t_{E}$ distributions and fractional contributions across various
compact object populations. _Upper_ : Mock OGLE-IV simulations (EWS cut)
results, from one seed per field. _Lower_ : Mock Roman survey results,
including all three seed runs for each of the three simulated Roman fields.
The left, middle, and right columns of each panel correspond to
$f_{\text{DM}}=1$, $f_{\text{DM}}=0.1$, and $f_{\text{DM}}=0.01$,
respectively. The upper row of each panel has three $t_{E}$ distributions: all
compact object populations, all populations except PBHs, and only the PBH
population. Median PBH $t_{E}$ values are shown by the dotted vertical lines.
The lower row of each panel shows the fractional $t_{E}$ contribution from
each population, as a function of Einstein crossing time. This figure should
be used for distribution demonstration purposes only.

Our Mock Roman survey uses the corresponding values from Table 2 to simulate
the Mock Roman labeled fields in Figure 4, which cover a subset of the
projected Roman microlensing survey footprint. The values in this table match
the “Mock WFIRST” values from Lam2020, which come from a combination of Penny
et al. (2019) and OGLE EWS values, and are described by Lam2020 as “reasonable
criteria expected for the mission”. Each field is simulated three times, each
with a different random seed, with all calculated quantities reported in this
work using the average of all simulated results from that field. We assume
constant stellar density over each simulated field whenever interpolating
results.

The Roman microlensing survey will be divided up into seasons throughout Roman
operation (exact dates are not yet settled on), where expected Roman operation
time is $\sim$5 years. Due to the microlensing survey seasons not being fully
defined, we choose to spread the total microlensing survey over the estimated
operation time, and to scale observing frequencies accordingly. Given the six
planned seasons, each 72 days in length, with $\sim$96 observations each day,
we estimate 41,472 observations total. Then, spreading the number of
observations equally over a 5 year estimated telescope operation time works
out to $\sim$23 observations per day. Because PopSyCLE cadence does not
account for various time gaps, we evenly spread the observations over the
telescope duration, ignoring daytime gaps and survey gaps. Because we do not
expect any PBH events with $t_{E}<10$ days, and we want to improve our
computational feasibility, we choose to slow the cadence to 1 observation
every 3 days. This enables us to stay within our computational boundaries,
while still accurately representing the $t_{E}$ distribution peak for all
other compact objects.

Slowing the cadence this way can lead to underestimating the number of
potentially detectable events, however, our choice of cadence should not
affect the BH and PBH populations. Additionally, the actual Roman survey
layout (i.e. seasons throughout telescope operation) will lead to event
overestimation, as some events will fall entirely within observing gaps
between seasons, or will be missed by only detecting one side of the light
curve, etc. The actual Roman survey may be more sensitive to short duration
events, but fail to detect a fraction of events with timescales that fall
between seasons (which can be $>$ 200 days if seasons are evenly spaced out
over operation), which are detected in our PopSyCLE simulations. We estimate
these gross cadence and duration approximations to be sufficient for an order
of magnitude estimate of the number of PBH microlensing events expected with
Roman.

The right panel of Figure 6 illustrates the $\pi_{E}$, $t_{E}$ distribution
for all events from all object populations, and the lower panel of Figure 7
shows the $t_{E}$ distributions and fractional $t_{E}$ contribution across
populations for our Mock Roman simulations. Note that this plot reflects
events from all three seed runs, from only a portion of the full Roman survey
footprint, and has no scaling applied beyond the Table 2 cuts, as they are
strictly for distribution demonstration.

Due to its resolving power and diffraction limited seeing, Roman will be
capable of obtaining simultaneous photometric and astrometric measurements
from observed microlensing events, with an estimated astrometric precision
(for isolated BHs) of $\sim$ 0.05 $mas$ (Sanderson et al., 2019). Blending
between the source and lens needs to be considered for luminous lenses, but
this work primarily examines dark lenses, thus we do not account for source-
lens blending. We also choose to ignore blending from neighboring stars as we
expect the scatter in including this blending to be subdominate to other
simplifying assumptions111111Considering blending from neighboring stars will
cause scatter in $\delta_{c,max}$, and generally lead to smaller astrometric
signals.. Therefore, we define the maximum astrometric shift as
$\delta_{c,max}=\delta_{c}$ (where $\delta_{c}$ is given in Equation 3),
calculated at the point in which the maximum astrometric shift occurs,
$u=\pm{\sqrt{2}}$. This yields the maximum distance that the image will
deviate from the expected path of the star with no lensing present. The left
panel of Figure 8 illustrates $\delta_{c,max}$ as a function of $\pi_{E}$, and
includes a dashed line for the 3$\sigma$ astrometric precision required for
Roman ($\sim$0.15 $mas$) to detect isolated black holes.

Figure 8: Simulated Roman survey distributions. _Left:_ Maximum astrometric
shift ($\delta_{c,max}$) vs. microlensing parallax ($\pi_{E}$). _Right:_
Maximum astrometric shift ($\delta_{c,max}$) vs. H-band source magnitude. The
horizontal dashed line represents the 3$\sigma$ astrometric precision required
for Roman ($\sim$0.15 mas) to detect isolated black holes (Sanderson et al.,
2019), the vertical dotted line represents the approximate limiting magnitude
for Roman (H$\sim$26), and the two diagonal lines represent the 1 exposure
(dash-dotted line) and 100 exposure (solid line) estimated precision curves.
The solid point denotes the fixed 1 exposure estimate, where
$\delta_{c,max}$=1.0 mas for a H=21.6 source (Sanderson et al., 2019). Roman
should be able to detect events that fall above the solid black line, and that
fall below the limiting magnitude of the survey. Note that this plot includes
events from simulations on only a small portion of the full planned survey
footprint, all three seed runs per simulated field are included in this
figure, and events represented here include only the cuts from Table 2.

Additionally, due to the assumed delta-function mass distribution for PBHs and
NSs, they appear highly correlated in this space. However, extending the mass
distributions will extend their footprint in this space, likely causing
overlap with the BH population. Resulting $\delta_{c,max}$ vs. $\pi_{E}$
distributions and $\delta_{c,max}$ vs. H-band source magnitude distributions
for our Mock Roman simulation are shown on the left and right panels of Figure
8, respectively.

We estimate the number of PBH microlensing events that may be observed by
Roman in three ways:

  1. 1.

Least conservative cut: estimate the number of PBH events in each of our
simulated 0.16 $\text{deg}^{2}$ fields by first averaging over the three seed
runs per field, then scale the results up from $3\times 0.16$ $\text{deg}^{2}$
to the full Roman survey area of 1.97 $\text{deg}^{2}$

  2. 2.

Realistic cut: count up the number of events that fall both above the solid
100 exposure precision line on Figure 8 and to the left of the dotted $H=26$
limiting magnitude limit line, divide results by three (as each of the three
seeds is represented on that plot), then scale up from $3\times 0.16$
$\text{deg}^{2}$ to the full Roman survey area of 1.97 $\text{deg}^{2}$

  3. 3.

Most conservative cut: multiply the result from method 1 by the survey duty
cycle:

$\text{duty cycle}=\frac{6\text{ seasons}\times 72\text{ days}}{5\text{
years}\times 365.24\frac{\text{ days}}{\text{ year}}}=0.24$ (10)

For determining our realistic event estimate cut, we first relate the H-band
apparent source magnitude, $H$, to the number of detected photons, $F$, using
Equation 11, and then set the signal-to-noise, SNR, approximately equal to
$\sqrt{F}$ using Equation 12:

$H\approx-2.5log_{10}(F),$ (11)
$\text{SNR}\approx\sqrt{10^{\frac{\text{H}}{-2.5}}}.$ (12)

Assuming $\delta_{c}\propto\frac{1}{SNR}$ we solve for the proportionality
constant $A$, using the estimated Roman precision of $\delta=1.0$ $mas$ at
$H=21.6$ (for a single exposure, Sanderson et al., 2019). Then, we calculate
$\delta_{c,max}$ as a function of H-band magnitude with
$\delta_{c,max}=\frac{A}{SNR}$ (over the full range of H values), to get the 1
exposure curve on the right panel of Figure 8. For the 100 exposure limit
curve we assume that the SNR scales as a factor of 10, and repeat using
$\delta_{c,max}=\frac{A}{\sqrt{100}SNR}$.

The right panel of 8 shows maximum astrometric shift, $\delta_{c,max}$, vs.
H-band source magnitude, for each of our individual microlensing events. Also
marked on this figure are Roman’s estimated limiting magnitude (H $\approx$
26), the calculated 1 exposure and 100 exposure estimated precision curves,
and a fixed solid point estimate for the 1 exposure precision, at
$\delta_{c,max}$ = 1.0 $mas$ for a H = 21.6 source (Sanderson et al., 2019).
Roman should be able to detect events that lie above the solid black line, and
that fall below the limiting magnitude of the survey (which we estimate here
as the dashed line at H = 26). Figure 8 includes the full distribution of
events over all three seeds, over all three simulated Roman fields, and assume
not scaling ratios or cuts other than those from Table 2.

##  5 Results

Adding a PBH model to PopSyCLE serves as a capable means for simulating PBH
dark matter, enabling us to estimate the number and statistical distribution
of detectable PBH microlensing events for various microlensing surveys, and to
constrain the abundance of PBH dark matter in the Milky Way. Throughout the
remainder of this section we discuss the results from our Mock OGLE-IV and
Mock Roman simulations, comparing estimated photometric and astrometric
measurement distributions between PBHs and other compact object populations,
using similar analysis techniques as used by Lam2020.

Lam2020 find that analyzing microlensing events in $\pi_{E}$, $t_{E}$ space
proves to be a compelling way of disentangling black holes from other types of
lenses, given only photometric information. While $\pi_{E}$, $t_{E}$ space can
enable ensemble statistical constraints, it does have real-time alert
limitations. Lam2020 notes that this analysis cannot be criteria for
astrometric follow-up selection, as $\pi_{E}$ is difficult to measure
precisely at or before the time of the photometric peak of the light curve
occurs. Fortunately, the next wave of microlensing surveys (e.g., Roman) are
expected to obtain simultaneous photometric and astrometric observations, and
Lam2020 further finds that combining photometric analysis with astrometric
measurements has the ability to further increase separability between
populations, making PopSyCLE a powerful tool for potentially differentiating
PBHs from other compact objects in these parameter spaces.

###  5.1 Mock OGLE-IV

Given the resulting PBH event distributions from Figure 6, we find that 30
$M_{\odot}$ PBH events tend towards higher $t_{E}$ and lower $\pi_{E}$, but
that the population significantly overlaps the stellar evolved BH
distribution. This overlap is likely due to our choice of PBH mass, as we
would only expect PBHs to shift in this space by a factor of
$1/\sqrt{\bar{m}_{\text{BH}}}$ in $\pi_{E}$, and by a factor of
$\sqrt{\bar{m}_{\text{BH}}}$ in $t_{E}$. Additionally, the median $\mu_{rel}$
is slightly higher for PBHs than BHs, with $\mu_{rel,\text{PBH}}$ $\sim$ 12.49
and $\mu_{rel,\text{BH}}$ $\sim$ 7.14 (averaged between the Mróz19 and EWS
cuts), causing an additional shift to lower $t_{E}$, undoing some of the
offset shift to higher $t_{E}$ due to the higher PBH masses.

Comparing the $t_{E}$ distributions with and without PBHs in Figure 7 reveals
that PBH microlensing events peak at longer duration $t_{E}$ than other
compact objects, with a peak of $\sim$ 91.7 days for both our Mróz19 and EWS
cuts. Additionally, analyzing the fractional contribution of PBH to BH events
in our EWS cut simulations (using the full distribution of events after only
the cuts from Table 2), we find $\sim$ 4.3$f_{\text{DM}}$ times more PBH
events than BH events. Note that our results use a higher cadence than the
Lam2020 “v3” results, catching more events with long duration. However, when
running our simulations without PBHs using the exact same simulation
parameters as Lam2020, we find identical results to those reported in Lam2020
Appendix A.1 and by Rose et al. (2022).

We find a small shift in Einstein crossing time peak values between our
simulations and Mróz et al. (2019) (consistant with Rose et al., 2022), likely
due to field-of-view extrapolation and the absence of binary lenses and binary
sources, which should not have a significant impact on our predicted long-
duration PBH events. Similar to Rose et al. (2022) and Lam2020, our star
counts and event rates can be off from Mróz et al. (2019) and EWS estimates by
a factor of $\approx$ 2 (see Table 3). The star count differences are not due
to stellar binarity or confusion (Abrams et al. in prep., Rose et al., 2022),
and the event rate discrepancies may be due to the galactic model chosen for
our simulations and extinction being variable closer to the galactic bulge
(Lam2020). It is worth noting that our EWS predicted PBH event rates will be
artificially low, as we have not accounted for spurious and anomalous events,
beyond our general scaling factor between our simulated results without PBHs
vs. the observational EWS results.

Given Table 3 results, and Figure 5, we find that 681 PBH events may be within
the 8 year OGLE-IV dataset, in which each PBH may conceivably be detected, but
reliable lensing parameter estimates may not be attainable for (i.e., EWS
cuts). We additionally estimate 543 PBH events may be detectable and also have
sufficient SNR such that lensing estimates can be reasonably estimated for
(i.e., Mróz19 cuts).

###  5.2 Mock Roman

Figure 6 illustrates the distribution results that we expect to be detectable
in the upcoming Roman microlensing survey. We find a significant overlap
between the PBH and BH populations which we expect to be for the same reason
as described in §5.1, which is mainly due to our choice of
$\bar{m}_{\text{PBH}}$ and the resulting $\mu_{rel}$ distributions, with
median values of $\mu_{rel,\text{PBH}}$ $\sim$ 11.4 and $\mu_{rel,\text{BH}}$
$\sim$ 7.1. We find $\sim$ 3.0 $f_{\text{DM}}$ times more PBH microlensing
events than we do BH microlensing events, and find a peak $t_{E}$ of $\sim$
91.1 days. The estimated number of PBH events to BH events is lower than that
of our simulated OGLE survey, but that is expected, as PBH events tend towards
longer timescale events, and OGLE is a longer survey.

Given our simplified PBH model and assumed Roman microlensing survey
structure, we estimate the number of potentially detectable PBH microlensing
events expected with Roman by using the seed-averaged and survey corrected
values per field over our 0.48 $\text{deg}^{2}$ simulated Roman survey area,
and scale up to the full Roman survey area of 1.97 $\text{deg}^{2}$
change(accounting for the small portion of overlap between two of the
simulated fields, where the overlap fraction = 0.19%). Assuming constant
stellar density over the Roman microlensing survey footprint,
$\bar{m}_{\text{PBH}}=30$ $M_{\odot}$, and that there are no survey gaps, we
estimate $\sim$ $10^{3}$$f_{\text{DM}}$ PBH microlensing events to be
detectable during the 5 year telescope operation.

For our least conservative estimate (scaling method 1, that scales our
predicted event rate up to the full survey area of 1.97 $\text{deg}^{2}$), we
estimate 5,422 events to be detectable with Roman. Our realistic cut estimate
(scaling method 2, which includes only the events that would fall above the
solid 100 exposure precision line, and to the left of the limiting magnitude
dashed line on Figure 8, and scaled to the full survey size) results in 4,835
PBH events detectable by Roman. Lastly, our most conservative estimate
(scaling 3, which multiplies the result of estimate 1 by the Roman duty cycle
given in equation 10, and then scales that result up to the full survey area
of 1.97 $\text{deg}^{2}$) results in 1,306 PBH events that might be detectable
by the future Roman survey.

Ignoring potentially important survey selection effects, the actual number of
expected PBH lensing events will scale roughly proportional to $f_{\text{DM}}$
and approximately constant as a function of mass. Because we slow the cadence
from $\sim$ 23 observations a day to 1 observation every 3 days (following
calculations described in §4.2), we underestimate the signal-to-noise of
microlensing events by a factor of $\sim$ 8, and ultimately underestimate the
number of microlensing events that may be detectable, however, we then
overestimate the predicted number due to our survey gap assumptions. For these
reasons, we expect our event estimates to hold to $\sim
10^{3}$$f_{\text{DM}}$.

Figure 9: $\pi_{E}$, $t_{E}$, $\mu_{rel}$ distribution of compact objects for
our simulated Roman survey (using all three seeds from each of the three
simulated fields, including only the cuts from Table 2). This figure
demonstrates that while there is significant overlap of the BH and PBH
populations in $\pi_{E}$, $t_{E}$ space (Figure 6), there is potential to
disentangle the populations with the addition of a $\mu_{rel}$ component.

Analyzing event distributions in $\delta_{c,max}$, $\pi_{E}$ space (left panel
of Figure 8) finds separability between the PBH and BH populations, more so
than in $\pi_{E}$, $t_{E}$ space. Note that because we use a delta-function
PBH mass distribution, the astrometric shift appears as a perfect linear
function in this space, however, exploring broader mass spectra will cause
scatter about $\delta_{c,max}$, likely causing overlap with the BH population.
We then explore the compact object event distributions in three-dimensional
$\pi_{E}$, $t_{E}$, $\mu_{rel}$ space, and find complete separability between
the PBH and BH populations (Figure 9). Again, this complete separability is
likely due to our choice of a constant mass, however, more analysis is needed
to confirm how much overlap would appear in this space given an extended PBH
mass distribution.

In general, decreasing $\bar{m}_{\text{PBH}}$ will result in PBH microlensing
events overlapping in $\pi_{E}$, $t_{E}$ space with compact object
microlensing events with similar masses, although different spatial and
velocity distributions will not make this exact. Additionally, our simulations
assume no blending, however, due to the sensitivity of Roman, binary sources
and flux from neighboring background stars will dilute the astrometric signal,
causing Figure 8 and Figure 9 to change, introducing scatter in the
$\delta_{c,max}$ and $\mu_{rel}$ dimensions. Penny et al. (2019) finds that
the mean blend fraction of detected microlensing events is $\sim$ 0.2 meaning
that the astrometric lensing signal could be diluted by a factor of $\sim$ 5\.
Regardless, our results demonstrates that accurate $\delta_{c,max}$ and
$\mu_{rel}$ estimates will play an important role in separating PBH and BH
events, making precise astrometric measurements with Roman a priority in the
next decade of microlensing.

##  6 Discussion and Conclusion

Dark matter makes up $\sim$85% of the matter content in the universe, and
though not well-understood, may be predominately comprised of PBHs. Recent
explorations have placed very tight constraints on $f_{\text{DM}}$ covering
many decades in mass (see summary in Bird et al., 2022). Due to cosmological
constraints, stellar evolved black holes can only account for a small portion
of dark matter (Collaboration et al., 2014), making detection and
characterization of PBHs an important key in understanding the structure and
formation of our universe.

Throughout this work we have detailed and demonstrated the utility of adding a
PBH population model into the microlensing simulation code PopSyCLE, enabling
us to simulate microlensing surveys in the Milky Way, serving as a powerful
tool for analyzing the estimated astrometric and photometric measurements of
potential PBH microlensing events. Because this work is first and foremost a
proof-of-concept for the addition of PBHs to PopSyCLE simulations, we make
many simplifying assumptions about PBHs (see section §3) and microlensing
survey parameters (see section §4). We acknowledge our assumption of
$f_{\text{DM}}=1$ and monochromatic PBH mass of $\bar{m}_{\text{PBH}}=30$
$M_{\odot}$ do not necessarily align with recent constraints on
$f_{\text{DM}}$ (e.g., Blaineau et al., 2022; Alcock et al., 2001; Tisserand
et al., 2007; Wyrzykowski et al., 2011) or the expectation of an extended mass
distribution, however they enable us to demonstrate our PopSyCLE \+ PBH
simulations, and allows the reader to scale our results to other values of
$f_{\text{DM}}$ and $\bar{m}_{\text{PBH}}$ as they see fit. Reported results
scale roughly proportional to $f_{\text{DM}}$ and approximately constant as a
function of mass, modulo potentially important survey selection effects.

Given our simplified PBH model and survey assumptions, we find the following:

  * •

Assuming a monochromatic PBH mass of $\bar{m}_{\text{PBH}}=30$ $M_{\odot}$,
PBH microlensing events tend towards longer duration ($t_{E}$), smaller
parallax ($\pi_{E}$), larger maximum astrometric shift ($\delta_{c,max}$), and
larger source-lens proper motion ($\mu_{rel}$) than BHs, but in most cases
PBHs overlap with the stellar BH population.

  * •

We predict Roman will have the ability to detect $\sim$ 1,306$f_{\text{DM}}$
(most conservative estimate, determined by area scaling ratios between our
simulations and the Roman survey, along with the Roman duty cycle of 0.24) to
$\sim$ 5,422$f_{\text{DM}}$ (least conservative estimate, determined by area
scaling ratios alone) PBH events, with a more realistic estimate of $\sim$
4,835$f_{\text{DM}}$ events (when considering events that would fall above the
100 exposure precision threshold for H-band source magnitude, and fall below
the limiting magnitude of the survey, then scaling to the full survey area).
These reported values assume that there is constant stellar density across
Roman fields, allowing us to scale up to the full Roman survey area of 1.97
$\text{deg}^{2}$, and that $\bar{m}_{\text{PBH}}=30$ $M_{\odot}$ lies within a
mass range detectable by Roman with moderate efficiency. We expect this to
hold to the order of $10^{3}f_{\text{DM}}$ after taking into account the
multiple overestimations and underestimations described throughout this work,
which are due to PBH population and survey assumptions.

  * •

We find that $\sim$ 355 PBH events may be within the 8 year OGLE-IV dataset,
in which each PBH may conceivably be detected, but reliable lensing parameter
estimates may not be attainable for. We additionally estimate $\sim$ 306 PBH
events may be detectable and also have sufficient SNR such that lensing
estimates can be reasonably estimated for each PBH. We expect this to hold on
the order of $10^{2}f_{\text{DM}}$.

  * •

With the assumption of $\bar{m}_{\text{PBH}}=30$, we find $\sim$
4.3$f_{\text{DM}}$ and $\sim$ 3.0$f_{\text{DM}}$ times more PBH microlensing
events than BH microlensing events for our Mock OGLE-IV (EWS cut) and Roman
simulations, respectively. The fraction of events is higher for the OGLE
survey as it has a longer survey duration, and PBH microlensing events tend
toward longer timescales. This averages out to 3.65$f_{\text{DM}}$ more PBHs
than BHs.

  * •

Using simulated PBH distributions with $f_{\text{DM}}=1$ and
$\bar{m}_{\text{PBH}}=30$, we find median $t_{E}$ values for our OGLE-IV EWS
cut and Roman simulations of $\sim$ 91.7 days and $\sim$ 91.1 days,
respectively, where peak $t_{E}$ scales approximately proportional to
$\sqrt{\bar{m}_{\text{PBH}}}$. This results in an average PBH estimated peak
Einstein crossing time of $\sim$ 91.4 days.

Roman will have the ability to conduct simultaneous photometric and
astrometric observations, with an astrometric sensitivity that enables
significant astrometric lensing measurement for the majority of events where
$M_{\text{PBH}}\geq 1$ $M_{\odot}$ (see Figure 8). Therefore, we expect Roman
to be an important instrument in the next decade of microlensing science (and
for potential PBH detection), in that it will enable better disambiguation
between the PBH and BH event populations in the $\pi_{E}$, $t_{E}$,
$\delta_{c,max}$ and $\mu_{rel}$ parameter spaces. Due to the overlap between
PBHs and BHs in these various spaces, and our average result of $\sim$
3.65$f_{\text{DM}}$ times more 30 $M_{\odot}$ PBH microlensing events than BH
events, surveying the Magellanic clouds may provide improved PBH dark matter
fraction constraints, as the lower stellar density should result in a higher
relative number of PBHs to BHs, enabling us to better constrain photometric
and astrometric PBH parameters. To be able to simulate this, a Magellanic
Cloud model would need to be added to PopSyCLE.

With all simplifying assumptions aside, we have demonstrated and validated our
addition of a PBH model into PopSyCLE that can be used to analyze astrometric
and photometric microlensing measurement predictions for PBHs, serving as a
powerful tool in the next decade of microlensing science and astronomical
surveys.

##  7 Acknowledgements

This work was performed under the auspices of the U.S. Department of Energy by
Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344 and
was supported by the LLNL-LDRD Program under Project Numbers 17-ERD-120 and
22-ERD-037, and by LLNL HEP Program Development. M.S.M acknowledges support
from the University of California Office of the President for the UC
Laboratory Fees Research Program In-Residence Graduate Fellowship (Grant ID:
LGF-19-600357). C.Y.L. and J.R.L. acknowledge support by the National Science
Foundation under grant No. 1909641, the National Aeronautics and Space
Administration (NASA) under contract No. NNG16PJ26C issued through the WFIRST
(now Roman) Science Investigation Teams Program, and the Heising-Simons
Foundation under Grant No. 2022-3542. C.Y.L. also acknowledges support from
NASA FINESST grant No. 80NSSC21K2043. LLNL-JRNL-834877.

## References

  * Abbott et al. (2016) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016, Physical Review Letters, 116, doi: 10.1103/physrevlett.116.061102
  * Abbott et al. (2017a) —. 2017a, Physical Review Letters, 118, doi: 10.1103/physrevlett.118.221101
  * Abbott et al. (2017b) —. 2017b, Physical Review Letters, 119, doi: 10.1103/physrevlett.119.141101
  * Alcock et al. (1993) Alcock, C., Akerlof, C. W., Allsman, R. A., et al. 1993, Nature, 365, 621, doi: 10.1038/365621a0
  * Alcock et al. (2001) Alcock, C., Allsman, R. A., Alves, D. R., et al. 2001, ApJ, 550, L169, doi: 10.1086/319636
  * Aubourg et al. (1993) Aubourg, E., Bareyre, P., Bréhin, S., et al. 1993, Nature, 365, 623, doi: 10.1038/365623a0
  * Bird et al. (2016) Bird, S., Cholis, I., Muñoz, J. B., et al. 2016, Physical Review Letters, 116, 201301, doi: 10.1103/PhysRevLett.116.201301
  * Bird et al. (2022) Bird, S., Albert, A., Dawson, W., et al. 2022, Snowmass2021 Cosmic Frontier White Paper:Primordial Black Hole Dark Matter, arXiv, doi: 10.48550/ARXIV.2203.08967
  * Blaineau et al. (2022) Blaineau, T., Moniez, M., Afonso, C., et al. 2022, arXiv, doi: 10.48550/ARXIV.2202.13819
  * Carr & Kuhnel (2021) Carr, B., & Kuhnel, F. 2021, Primordial Black Holes as Dark Matter Candidates.  https://arxiv.org/abs/2110.02821
  * Carr & Silk (2018) Carr, B., & Silk, J. 2018, Monthly Notices of the Royal Astronomical Society, 478, 3756, doi: 10.1093/mnras/sty1204
  * Carr (1975) Carr, B. J. 1975, ApJ, 201, 1, doi: 10.1086/153853
  * Chapline & Barbieri (2018) Chapline, G., & Barbieri, J. 2018, LHEP, 1, 17, doi: 10.31526/LHEP.1.2018.04
  * Chapline (1975a) Chapline, G. F. 1975a, Phys. Rev. D, 12, 2949, doi: 10.1103/PhysRevD.12.2949
  * Chapline (1975b) —. 1975b, Nature, 253, 251, doi: 10.1038/253251a0
  * Chapline & Frampton (2016) Chapline, G. F., & Frampton, P. H. 2016, Journal of Cosmology and Astroparticle Physics, 2016, 042–042, doi: 10.1088/1475-7516/2016/11/042
  * Chatzopoulos et al. (2015) Chatzopoulos, S., Fritz, T. K., Gerhard, O., et al. 2015, MNRAS, 447, 948, doi: 10.1093/mnras/stu2452
  * Collaboration et al. (2014) Collaboration, P., Ade, P. A. R., Aghanim, N., et al. 2014, Astronomy & Astrophysics, 571, A16, doi: 10.1051/0004-6361/201321591
  * Croon et al. (2020) Croon, D., McKeen, D., Raj, N., & Wang, Z. 2020, Physical Review D, 102, doi: 10.1103/physrevd.102.083021
  * Damineli et al. (2016) Damineli, A., Almeida, L. A., Blum, R. D., et al. 2016, MNRAS, 463, 2653, doi: 10.1093/mnras/stw2122
  * Eddington (1916) Eddington, A. S. 1916, Monthly Notices of the Royal Astronomical Society, 76, 572, doi: 10.1093/mnras/76.7.572
  * Gillessen et al. (2009) Gillessen, S., Eisenhauer, F., Trippe, S., et al. 2009, ApJ, 692, 1075, doi: 10.1088/0004-637x/692/2/1075
  * Hawking (1971) Hawking, S. 1971, MNRAS, 152, 75, doi: 10.1093/mnras/152.1.75
  * Hawking (1974) Hawking, S. W. 1974, Nature, 248, 30, doi: 10.1038/248030a0
  * Hosek Jr. et al. (2020) Hosek Jr., M. W., Lu, J. R., Lam, C. Y., et al. 2020, The Astronomical Journal, 160, 143, doi: 10.3847/1538-3881/aba533
  * Johnson et al. (2020) Johnson, S. A., Penny, M., Gaudi, B. S., et al. 2020, AJ, 160, 123, doi: 10.3847/1538-3881/aba75b
  * Lacroix et al. (2018) Lacroix, T., Stref, M., & Lavalle, J. 2018, Journal of Cosmology and Astroparticle Physics, 2018, 040–040, doi: 10.1088/1475-7516/2018/09/040
  * Lam et al. (2020) Lam, C. Y., Lu, J. R., Hosek, M. W., Dawson, W. A., & Golovich, N. R. 2020, The Astrophysical Journal, 889, 31, doi: 10.3847/1538-4357/ab5fd3
  * Lu et al. (2019) Lu, J. R., Lam, C. Y., Medford, M., Dawson, W., & Golovich, N. 2019, Research Notes of the AAS, 3, 58, doi: 10.3847/2515-5172/ab1421
  * McMillan (2016) McMillan, P. J. 2016, Monthly Notices of the Royal Astronomical Society, 465, 76–94, doi: 10.1093/mnras/stw2759
  * Mróz et al. (2019) Mróz, P., Udalski, A., Skowron, J., et al. 2019, The Astrophysical Journal Supplement Series, 244, 29, doi: 10.3847/1538-4365/ab426b
  * Muraki et al. (1999) Muraki, Y., Sumi, T., Abe, F., et al. 1999, Progress of Theoretical Physics Supplement, 133, 233, doi: 10.1143/PTPS.133.233
  * Navarro et al. (1996) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563, doi: 10.1086/177173
  * Penny et al. (2019) Penny, M. T., Gaudi, B. S., Kerins, E., et al. 2019, The Astrophysical Journal Supplement Series, 241, 3, doi: 10.3847/1538-4365/aafb69
  * Piffl et al. (2014) Piffl, T., Scannapieco, C., Binney, J., et al. 2014, Astronomy & Astrophysics, 562, A91, doi: 10.1051/0004-6361/201322531
  * Press et al. (1988) Press, W. H., Bellini, A., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1988, Numerical Recipes: The Art of Scientific Computing
  * Raithel et al. (2018) Raithel, C. A., Sukhbold, T., & Özel, F. 2018, The Astrophysical Journal, 856, 35, doi: 10.3847/1538-4357/aab09b
  * Robin et al. (2004) Robin, A. C., Reylé, C., Derriére, S., & Picaud, S. 2004, Astronomy & Astrophysics, 416, 157–157, doi: 10.1051/0004-6361:20040968
  * Rose et al. (2022) Rose, S., Lam, C. Y., Lu, J. R., et al. 2022, The Astrophysical Journal, 941, 116, doi: 10.3847/1538-4357/aca09d
  * Sanderson et al. (2019) Sanderson, R. E., Bellini, A., Casertano, S., et al. 2019, Astrometry with the Wide-Field InfraRed Space Telescope.  https://arxiv.org/abs/1712.05420
  * Sharma et al. (2011) Sharma, S., Bland-Hawthorn, J., Johnston, K. V., & Binney, J. 2011, The Astrophysical Journal, 730, 3, doi: 10.1088/0004-637x/730/1/3
  * Smith et al. (2007) Smith, M. C., Ruchti, G. R., Helmi, A., et al. 2007, Monthly Notices of the Royal Astronomical Society, 379, 755–772, doi: 10.1111/j.1365-2966.2007.11964.x
  * 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
  * Tisserand et al. (2007) Tisserand, P., Guillou, L. L., Afonso, C., et al. 2007, Astronomy & Astrophysics, 469, 387, doi: 10.1051/0004-6361:20066017
  * Udalski et al. (1992) Udalski, A., Szymanski, M., Kaluzny, J., Kubiak, M., & Mateo, M. 1992, Acta Astron., 42, 253
  * Udalski et al. (2015) Udalski, A., Szymański, M. K., & Szymański, G. 2015, Acta Astron., 65, 1\.  https://arxiv.org/abs/1504.05966
  * Udalski et al. (1994) Udalski, A., Szymanski, M., Stanek, K. Z., et al. 1994, The Astrophysical Journal
  * Wyrzykowski et al. (2011) Wyrzykowski, L., Skowron, J., Kozłowski, S., et al. 2011, Monthly Notices of the Royal Astronomical Society, 416, 2949, doi: 10.1111/j.1365-2966.2011.19243.x
  * Zel’dovich & Novikov (1967) Zel’dovich, Y. B., & Novikov, I. D. 1967, Soviet Ast., 10, 602