# Observing the Galactic Underworld: Predicting photometry and astrometry from
compact remnant microlensing events

David Sweeney1,2,3, Peter Tuthill1, Alberto Krone-Martins2, Antoine Mérand3,
Richard Scalzo4, Marc-Antoine Martinod1
1Sydney Institute for Astronomy (SIfA), The University of Sydney, Physics
Road, Sydney 2050, Australia
2Donald Bren School of Information and Computer Sciences, University of
California, Irvine, CA 92697, USA
3European Southern Observatory, Karl-Schwarzschild-Straße 2, 85748 Garching,
Germany
4CSIRO Data61, Clayton, VIC 3168, Australia E-mail:
<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Isolated black holes (BHs) and neutron stars (NSs) are largely undetectable
across the electromagnetic spectrum. For this reason, our only real prospect
of observing these isolated compact remnants is via microlensing; a feat
recently performed for the first time. However, characterisation of the
microlensing events caused by BHs and NSs is still in its infancy. In this
work, we perform N-body simulations to explore the frequency and physical
characteristics of microlensing events across the entire sky. Our simulations
find that every year we can expect $88_{-6}^{+6}$ BH, $6.8_{-1.6}^{+1.7}$ NS
and $20^{+30}_{-20}$ stellar microlensing events which cause an astrometric
shift larger than 2 mas. Similarly, we can expect $21_{-3}^{+3}$ BH,
$18_{-3}^{+3}$ NS and $7500_{-500}^{+500}$ stellar microlensing events which
cause a bump magnitude larger than 1 mag. Leveraging a more comprehensive
dynamical model than prior work, we predict the fraction of microlensing
events caused by BHs as a function of Einstein time to be smaller than
previously thought. Comparison of our microlensing simulations to events in
Gaia finds good agreement. Finally, we predict that in the combination of Gaia
and GaiaNIR data there will be $14700_{-900}^{+600}$ BH and
$1600_{-200}^{+300}$ NS events creating a centroid shift larger than 1 mas and
$330_{-120}^{+100}$ BH and $310_{-100}^{+110}$ NS events causing bump
magnitudes $<-1$. Of these, $<10$ BH and $5_{-5}^{+10}$ NS events should be
detectable using current analysis techniques. These results inform future
astrometric mission design, such as GaiaNIR, as they indicate that, compared
to stellar events, there are fewer observable BH events than previously
thought.

###### keywords:

gravitational lensing: micro – astrometry – stars: black holes – stars:
neutron, methods: statistical – software: simulations

††pubyear: 2023††pagerange: Observing the Galactic Underworld: Predicting
photometry and astrometry from compact remnant microlensing events–F

##  1 Introduction

Gravitational microlensing occurs when there is a close alignment between a
background source, a foreground object and an observer. The mass of the
foreground object warps nearby spacetime and thus the light of the background
source is “bent” as it travels through this warped region. The end result is
that when we observe the background source, the signal has been altered by the
presence of the foreground object, i.e. the background source has been
_lensed_ by the foreground object. This can manifest in several ways: the
apparent brightening of the background source, a shift in the apparent
position of the background source or in multiple observable images of the
background source.

As two celestial bodies participate in microlensing, we can use events to gain
information on either or both. Since the strength of the lensing depends on
the closeness of the alignment between the lens and the background source and
the mass of the lens, microlensing events can be used to ascertain not only
the presence of a lens, but also its mass. This has been done with white
dwarfs (Sahu et al., 2017) and was recently done for an isolated BH of about 8
M⊙ (Lam et al., 2022; Sahu et al., 2022; Mróz et al., 2022; Lam & Lu, 2023).
Microlensing events caused by isolated BHs and NSs are of particular interest
to modern astronomy because these objects are virtually invisible to
traditional methods of observation due to their lack of detectable
electromagnetic signatures.

The population of compact remnants, BHs and NSs, has not been observationally
characterised; a problem that is particularly acute for isolated remnants for
which few are confirmed in the literature. Isolated pulsars are readily
observable in the radio for an interval after their birth, but with spin-down
they fade from view on a (cosmically) brief timescale. Once these pulsars
darken, NSs become nearly as invisible as BHs. With the singular exception
mentioned above, all confirmed stellar mass BHs have been identified through
their interactions in binaries: as a component in massive X-Ray binaries,
radial velocity measurements (El-Badry et al., 2023; Nagarajan et al., 2023)
or through inspiral signals recorded in gravitational wave detectors (Abbott
et al., 2016). Outside of binary interactions, which obviously do not occur
for isolated remnants, the only current technology capable of detecting these
remnants is the use of microlensing measurements. There have been a number of
large-scale surveys studying microlensing events: MACHO (Bennett et al.,
1993), EROS (Aubourg et al., 1993), OGLE (Udalski et al., 1992) and KMTNet
(Kim et al., 2016), but while these surveys are sensitive to NS and BH
microlensing events, they cannot distinguish them from stellar events because
light-curve-only observations cannot break model ambiguity in parameters
describing the lens mass, distance and relative proper motion of the event.

The follow-up required for full characterisation of these events is demanding
— the only successfully confirmed BH detection required the use of the Hubble
Space Telescope (HST) over 6 years (Lam et al., 2022; Sahu et al., 2022). Even
with HST data, several reanalyses were required before agreement was reached
on the mass of the lensing BH (Mróz et al., 2022; Lam & Lu, 2023). Utilising
interferometric follow up is a promising avenue for obtaining astrometric
information on these events without the need for years of HST observations,
however this method is yet to post a successful detection. Either way, it is
crucial that we understand the characteristics of these events so that we may
best target them with a coherent, integrated and multi-facility observing
strategy (see Gould, 2023).

To do this we must model BH and NS microlensing events and compare them to
stellar events. Previous work has used various stellar population models,
which improve with every iteration, to construct a population of compact
remnants and calculate microlensing events (Han & Gould, 2003; Wood & Mao,
2005; Osłowski et al., 2008; Dai et al., 2015; Lam et al., 2020). The latest
of these, Lam et al. (2020) and subsequent Rose et al. (2022), advanced the
field by simulating realistic extinction for high mass stars.

In this work we take advantage of an improved population synthesis model for
compact remnants and use the Gaia catalogue as our source for background stars
so that our selection function is more realistic than previous modelling. Our
improved lens population synthesis model comes from the work by Sweeney et al.
(2022) which, in contrast to PopSyCLE (Lam et al., 2020) and the work by
Olejak et al. (2020), evolves BHs and NSs post-supernovae. This evolution
substantially alters the BH and NS distribution in the Galaxy due to the large
kicks ($\sim 450$ km/s for NSs, $\sim 80$ km/s for BHs) received by NSs and
BHs during their birth in supernovae. We use these populations to calculate
microlensing events across the entire sky and compare our results to the
previous characterisation by Lam et al. (2020) and the microlensing events
observed in Gaia (Wyrzykowski et al., 2023). We also make predictions about
microlensing events detectable in the combination of Gaia and GaiaNIR (Hobbs
et al., 2016) data.

The first confirmed detection of a stellar mass BH via microlensing (Lam et
al., 2022; Sahu et al., 2022; Mróz et al., 2022; Lam & Lu, 2023) along with
the handful of ambiguous detections of compact objects (Mróz & Wyrzykowski,
2021) demonstrate that we are on the cusp of uncovering the large scale
distribution of compact remnants. In order to efficiently build upon this
discovery, it is crucial that observing strategies be informed by the accurate
parameterisation of the population of compact remnant lensing events. The
impending opening of a new observational window on this previously hidden
galactic population motivates re-characterisation of microlensing events
involving compact remnants, now leveraging the improved distributions of
Sweeney et al. (2022). In turn, distributions of events may be confronted with
new data from Gaia and inform forthcoming global astrometry mission proposals
such as GaiaNIR.

##  2 Data

Figure 1: Comparison between the probability density map of the population of
NSs + BHs and a synthetic visible Milky Way created by GALAXIA. Because of
their greater number, NSs dominate the BH + NS distribution. The probability
density is estimated by a gaussian kernel density estimator such that the
integral over the whole $x$ – $z$ plane is 1 and logarithmically spaced
contour lines are placed on top of the plot. This plot has been recreated from
Sweeney et al. (2022).

There are three sets of data used in this work. The first is the data used to
describe the distribution of BGSs which is taken from the Gaia catalogue (Gaia
Collaboration et al., 2016, 2023b). The sources which populate the catalogue
have already been affected by extinction due to dust, observational biases
toward luminous objects and crowded field effects. By utilising the Gaia
catalogue as BGSs our simulations inherit this selection function, ensuring
that the microlensing events simulated are from plausibly monitored sources.
Null values in the Gaia catalogue are discussed in Appendix A. Using the Gaia
catalogue has the added advantage that it allows comparisons to be easily made
between the predictions given in this work and the observed rate of
microlensing in Gaia.

The second set of data is the set of BH and NS lenses which comes from the
paper by Sweeney et al. (2022). This distribution is reached by starting with
a synthetic stellar distribution of the Milky Way generated by GALAXIA (Sharma
et al., 2011). GALAXIA formerly only modelled stellar populations: objects
still burning nuclear fuel, but was modified to output stars that evolve past
their nuclear burning lifespan. Critically, GALAXIA synthesises a realistic
population which evolves through cosmic time, resulting in a galaxy which has
a record of when each remnant formed, allowing the distribution to be modified
to account for remnant-specific evolution (e.g. natal kicks). From this
distribution of exhausted stars, NS and BH progenitors were identified based
on the initial mass of each progenitor star: 8–$25$ M⊙ were evolved into NSs
and $>25$ M⊙ became BHs (Heger et al., 2003; Smartt, 2009). Remnants were then
given a natal kick drawn from the bimodal Maxwellian distribution published by
Igoshev (2020) which was derived from the observed velocity of pulsars. This
natal kick was added to the kinematics the remnant inherited from its main
sequence progenitor. They key difference from previous work is that each
remnant was then evolved through the Galactic potential (MWPotential2014 from
galpy; Bovy, 2015) for a duration equal to the time since the remnant formed.
This process resulted in a remnant population which has a dramatically
different spatial distribution compared to the visible Galaxy, as shown in
Figure 1. It should be noted that the data files from Sweeney et al. (2022)
undersample the true galactic population by a factor of 1000 and we correct
for this during our analysis.

The BH and NS population makes a number of simplifying assumptions, summarised
as follows:

  * •

The distribution generated by GALAXIA is correct. GALAXIA was used to generate
data based on the Besançon Galaxy (Robin et al., 2003) model with some
modifications. This model matches observational constraints well but has no
history of Galactic mergers. GALAXIA also doesn’t generate multiple star
systems. This is not thought to be a major simplification if considering only
isolated BHs and NSs, as we have done here.

  * •

Compact remnants are split into BHs and NSs based on the initial mass of their
progenitor. Progenitors with initial masses of 8–25M☉ become NSs and
progenitors with initial masses >25M☉ become BHs. This simplifying assumption
is frequently made (e.g. Fryer, 1999; Belczynski & Taam, 2008; Fryer et al.,
2012) and has some observation and theoretical support (Smartt, 2015; Adams et
al., 2017). Metalicity likely also plays a role, but averaged over the entire
population we expect this effect to be small.

  * •

NSs receive a natal kick following the bimodal Maxwellian derived in Igoshev
(2020) from pulsar observations. While the exact formulation of NS natal kicks
is not settled, the distribution is thought to be similar to Igoshev (2020)
(Lyne & Lorimer, 1994; Hobbs et al., 2005; Katsuda et al., 2018). BHs receive
a natal kick with the same momentum, which, due to their larger mass, launches
them at lower velocity. BHs with progenitors with initial masses >40M☉ were
said to undergo direct collapse and were not kicked. Little is known about BH
natal kicks, so this is not the only formulation possible.

  * •

BHs all have a mass of 7.8M☉ and NSs all have a mass of 1.35M☉. While there is
not significant variance in NS masses (Özel et al., 2012; Postnov & Yungelson,
2014; Pejcha & Thompson, 2015; Sukhbold et al., 2016), there is thought to be
a large range of masses for stellar mass BHs. The range and distribution of BH
masses is not known and would (in this formulation) impact the natal kicks
that BHs receive. 7.8M☉ was chosen because it is a widely accepted median
value (Özel et al., 2010; Spera et al., 2015; Pejcha & Thompson, 2015;
Sukhbold et al., 2016).

  * •

That the MWPotential2014 Galactic potential is correct. In practice this is
likely not the case as it does not contain a bar which affects microlensing
observations toward the bulge (Han & Gould, 1995), however aside from this
deficiency, the MWPotential2014 potential is a good fit to measured Galactic
properties (Bovy, 2015). The stellar population does not have to be evolved
through the Galactic potential and so still exhibits the bar generated by
GALAXIA.

For more details see the original paper (Sweeney et al., 2022).

The third set of data describes the population of stellar lenses which is also
generated by GALAXIA. The standard release version of GALAXIA was used to
generate a realistic population of stars in our Galaxy. The stellar data
undersamples the true population by a factor of $10^{6}$, which is also
corrected for during our analysis. To convert the absolute magnitudes provided
by GALAXIA to apparent magnitudes we use the conversion to CTIO V and R
filters described in the GALAXIA documentation and then use the Johnson-
Cousins relationships specified in the Gaia documentation to convert to
GaiaDR3 G magnitudes (Busso et al., 2022):

$\displaystyle V_{\text{CTIO}}$
$\displaystyle=V_{\texttt{GALAXIA}}+5\ln\left[\frac{100d}{\text{kpc}}\right]+3.24E^{Schlegel}_{B-V}$
(1) $\displaystyle R_{\text{CTIO}}$
$\displaystyle=R_{\texttt{GALAXIA}}+5\ln\left[\frac{100d}{\text{kpc}}\right]+2.634E^{Schlegel}_{B-V}$
(2) $\displaystyle G_{\text{Gaia}}$
$\displaystyle=V_{\text{CTIO}}-0.03088-0.04653(V_{\text{CTIO}}-R_{\text{CTIO}})$
$\displaystyle\;\;\;\;\;-0.8794(V_{\text{CTIO}}-R_{\text{CTIO}})^{2}+0.1733(V_{\text{CTIO}}-R_{\text{CTIO}})^{3}$
(3)

Where $d$ is the star’s distance from Earth and $E^{Schlegel}_{B-V}$ is the 3D
extinction modelled by GALAXIA based on the Schlegel maps (Schlegel et al.,
1998).

Note that we do not include a population of white dwarfs as they were not
modelled in Sweeney et al. (2022). Modelling microlensing events caused by
white dwarfs was done by Lam et al. (2020) who found they yield signatures
easily distinguishable from BH events, especially when considering
microlensing events on the $t_{E}$–$\pi_{E}$ plane.

##  3 Methods

###  3.1 An introduction to microlensing

Microlensing, and gravitational lensing more broadly, is a prediction from
General Relativity (Einstein, 1936). It describes an effect which occurs when
there is a close alignment between a massive foreground lens and a background
source (BGS). In the case of a point mass, this effect creates two images of
the BGS, a bright major image near to the “true" position of the BGS and a dim
minor image on the opposite side of the lens. In the case of an extremely
close alignment between the lens and BGS an Einstein ring can appear around
the lens instead of a pair of images. While very rare, this Einstein ring
gives rise to some natural quantities which describe a lensing event. Firstly
there is the Einstein angle which is the angular radius of the Einstein ring
(Paczynski, 1986). The Einstein angle is typically used as a unit of the scale
for microlensing events and is defined as:

$\displaystyle\theta_{E}=\sqrt{\frac{4GM_{L}}{c^{2}}\frac{D_{S}-D_{L}}{D_{S}\cdot
D_{L}}}=\sqrt{\frac{4GM_{L}}{c^{2}\cdot 1\text{
pc}}\frac{\bar{\omega}_{L}-\bar{\omega}_{S}}{1\text{ arcsec}}}$ (4)

Where $G$ is the gravitational constant, $M_{L}$ is the mass of the lens, $c$
is the speed of light, $D_{L}$ and $D_{S}$ are the distance to the lens and
BGS and $\bar{\omega}_{L}$ and $\bar{\omega}_{S}$ are the parallax of the lens
and BGS, respectively.

The Einstein angle then defines a number of other quantities, such as the
separation between the lens and BGS as measure in Einstein angles:

$\displaystyle\boldsymbol{u}=\frac{\boldsymbol{\theta}}{\theta_{E}}$ (5)

Where $\boldsymbol{\theta}$ is the two dimensional vector of the unlensed
lens–BGS angular separation. As a note, we follow the convention that bold
symbols (e.g. $\boldsymbol{u}$) represent vectors and their unbold variants
represent their magnitudes (e.g. $u$). Also note that both $\boldsymbol{\mu}$
and $\boldsymbol{\theta}$ are time dependant since the lens–BGS separation
changes with time; sometimes this is made explicit in notation, e.g. by
writing $\boldsymbol{u}(t)$. By convention, quantities sub-scripted with 0
(e.g. $u_{0}$) are the value this quantity takes when the lens and BGS are
closest together and the lensing event is strongest.

The timescale of microlensing events is usually related in Einstein times, the
time taken for the lens to travel across the Einstein angle:

$\displaystyle t_{E}$ $\displaystyle=\frac{\theta_{E}}{\mu_{\text{rel}}}$ (6)

Where $\mu_{\text{rel}}$ is the absolute value of the relative proper motion
between the lens and BGS.

The creation of two images during a microlensing event leads to two observable
quantities of interest: an astrometric shift and a photometric magnification.
The astrometric shift of the major and minor images (i.e. the difference in
position between the major/minor image and the true position of the BGS) is
given by:

$\displaystyle\delta\boldsymbol{\theta_{\pm}}=\frac{\pm\sqrt{u^{2}+4}-u}{2}\cdot\frac{\boldsymbol{u}}{u}\cdot\theta_{E}$
(7)

Where the astrometric shift of the major/minor image is denoted with
$\delta\boldsymbol{\theta_{+}}$/$\delta\boldsymbol{\theta_{-}}$, respectively.

In practice this shift is hard to measure as the major and minor images
typically are not resolved. Instead we observe the astrometric shift of the
centroid due to the combination of the two images:

$\displaystyle\delta\boldsymbol{\theta_{c}}=\frac{\boldsymbol{u}}{u^{2}+2}\cdot\theta_{E}$
(8)

This astrometric shift is further complicated by blending between the source
and a luminous lens (Dominik & Sahu, 2000). Defining the flux ratio between
the lens and source, $g=F_{L}/F_{S}$, we can write the blended centroid shift
as:

$\displaystyle\delta\boldsymbol{\theta_{c,LL}}=\frac{\theta_{E}}{1+g}\cdot\frac{1+g(u^{2}-u\sqrt{u^{2}+4}+3)}{u^{2}+2+gu\sqrt{u^{2}+4}}\cdot\boldsymbol{u}$
(9)

Note that in the case of a non-luminous lens $g=0$ and this reduces to
Equation 8.

The photometric magnification of the BGS is derived by summing the
magnifications of the major and minor image:

$\displaystyle\mu=\frac{u^{2}+2}{u\sqrt{u^{2}+4}}$ (10)

In the case of a luminous lens, the source-lens blending reduces the
magnification observed in a photometric light curve, as the luminous lens
increases the observed baseline magnitude. To account for this we define a
bump magnitude which is the difference between the peak magnitude and the
baseline magnitude:

$\displaystyle\Delta m=-2.5\log_{10}\left(\frac{\mu_{0}+g}{1+g}\right)$ (11)

For a detailed derivation of these quantities see the review by Mao (2008).

One final quantity is the microlens parallax:

$\displaystyle\pi_{E}=\frac{\bar{\omega}_{L}-\bar{\omega}_{S}}{\theta_{E}}$
(12)

For long microlensing events ($t_{E}>100$ days) the microlens parallax can
usually be measured from photometry, as the Earth’s orbit around the sun
imprints an asymmetry into the light curve.

###  3.2 Calculating microlensing events

To calculate the number of microlensing events caused by compact remnants
(ignoring remnants further than 20 kpc from the Galactic centre) we identified
the stars in Gaia which were close enough to each remnant to potentially cause
a major image shift larger than 1 $\mu$as. Similar to the approach taken by
Klüter et al. (2022) for white dwarfs, this was done by calculating the arc of
sky traversed by each remnant during the $n$ years of “observations”. A circle
was then drawn which encapsulates this arc with the following buffer term
added to the radius:

$\displaystyle\text{buffer}=\frac{\theta_{E}^{2}}{1\,\mu\text{as}}+\bar{\omega}_{L}+1.2\,\text{mas}+10.5n\,\text{mas}$
(13)

Where the first term on the right-hand side is an approximation of Equation 8
for the maximum angular separation between lens and BGS to result in a 1
$\mu$as astrometric shift. $\theta_{E}$ is calculated with a BGS parallax of
$-1.2$ mas, chosen because 95% of queried Gaia sources had a larger parallax
than this value. The second and third terms are added to encapsulate the
parallax motion of the lens and most BGSs. The final term is added to include
high proper motion BGSs which might move near to the path. The value of this
term is 10.5$n$ mas as 95% of queried Gaia sources had a proper motion smaller
than 10.5 mas/yr. An illustration of the resulting circle is shown in Figure
2. It should be noted that because of the chosen parallax and proper motion
thresholds our results will slightly underestimate the number of events which
occur as we will miss consideration of some BGSs. Testing revealed this
underestimation is <1% for events with a lensing magnification of >1.0001
(i.e. $u\lesssim 20$).

Figure 2: Diagram depicting the motion of a lensing object and the region in
which Gaia stars are selected. For each star in this region the microlensing
effect caused by the lens will be calculated.

The Gaia catalogue was then queried for stars inside this circle. The returned
sources were then filtered to remove objects with a parallax value which was
negative by more than 5 times the parallax uncertainty (i.e.
$\bar{\omega}_{S}<0$ and $|\text{parallax over error}|>5$) and those with a
parallax larger than the lens. Negative parallaxes are unphysical, but
removing them entirely has been shown to bias the data (Luri et al., 2018), so
we use this filter as a middle ground.

The position, neglecting parallax effects, of the remnant and each star was
then calculated over the span of $n$ years. The minimum separation between the
two objects was calculated via the haversine formula to within a day. A cut on
the size of the microlensing event is then applied:

$\displaystyle\left|\frac{\theta_{E}^{2}}{(\theta-(\bar{\omega}_{L}-\bar{\omega}_{S}))}\right|\geq
0.8\,\mu\text{as},$ (14)
$\displaystyle\theta-(\bar{\omega}_{L}-\bar{\omega}_{S})<0$ (15)

If neither condition is true then the BGS was discarded, otherwise the minimum
separation was recalculated with parallax taken into account. The calculation
is broken up like this as the minimum separation between two objects without
parallax effects is a convex function, and so requires much less computational
power to minimise.

The parallax for each object was increased by 1% to emulate the effect for the
Gaia telescope at L2. Taking into account these parallax effects, the minimum
separation is then calculated by calculating the minimum separation every day
for 100 years either side of the previously found closest approach. If the new
minimum separation value is found within 1% of either extreme of this checked
period the period is doubled (respecting the observation bounds) and this
process repeats. Finally, the minimum separation is recomputed, checking every
14.4 minutes (0.01 days) within a day either side of the closest approach.
This process results in the minimum separation being calculated to within
$\sim 15$ minutes.

The minimum separation is then used to calculate the astrometric shifts,
photometric magnification and bump magnitude using Equations 7–11. If the
astrometric shift of the major image is smaller than 1 $\mu$as the event is
discarded, otherwise details of the remnant, BGS, minimum separation, Einstein
angle, Einstein time, major image shift, centroid shift, blended centroid
shift, lensing magnification, bump magnitude and the time at which the event
occurred are recorded. Since the centroid shift (Equation 8) is maximised at
$u=\pm\sqrt{2}$, if the minimum separation has $u_{0}<\sqrt{2}$ the centroid
shift is calculated with $u=\sqrt{2}$ as the lens will cause this shift at
some point in time. The $u$ which maximises the blended centroid shift depends
on $g$. This $u$ is calculated for each event and if $u_{0}$ is smaller than
this value, the maximum blended centroid shift is used. In the case of a
remnant that does not cause any major image shifts larger than 1 $\mu$as, its
details are recorded with null values in the remaining columns. The 1 $\mu$as
major image astrometric cut was chosen to be well below the threshold of
observable events.

This results in a list of microlensing events with major image astrometric
shifts larger than 1 $\mu$as caused by the compact remnant population derived
by Sweeney et al. (2022). In practice this procedure is parallelised (this is
easily done as the effects of each remnant are independent) and the results
are combined at the end.

###  3.3 Different timespans

In this work three different simulations, each with a distinct timespan, were
performed:

_Simulation 1_ spanned 10 000 years for which BHs and NSs constituted the
lensing population. As noted in Section 2, the compact remnant population is
undersampled by a factor of 1000, so this simulation provides an effective 10
years of compact remnant lensing in the Milky Way.

_Simulation 2_ spanned 100 000 years, modelling the scenario where stars form
the lensing population. Similar to the BH and NS data, the stellar data is
undersampled, now by a factor of $10^{6}$; therefore this simulation is
equivalent to 1/10th of a year of real microlensing observation.

_Simulation 3_ modelled a 32 year timespan, once again using BHs and NSs as
lenses. This timespan consisted of 12 years of “observation", where the
microlensing events were calculated, a 15-year gap and then 5 years of further
observation. This format was chosen to mimic the lifespan of Gaia and then a
hiatus of observation before a potential 5-year future GaiaNIR mission. To
deliver robust statistics despite the undersampled population, 101 different
32-year timespans were modelled. This was achieved by using the publicly
available code from Sweeney et al. (2022) to evolve the compact remnant
distribution from its nominal initial (“present day”) configuration by 20 000
years into the future, saving the distribution every 200 years. Each timestep
forms an effectively independent lens population for the purposes of the
process described in Section 3.2 as 99.9% of lensing events have Einstein
times less than 37 years.

The first and final year of any observation period exhibits an increased
number of microlensing event “detections” due to events which peak outside of
the observing period. This is an important effect for Simulation 3 as these
events have been detected in Gaia data (Wyrzykowski et al., 2023), but is
undesirable for Simulations 1 and 2 which are used to calculate yearly event
rates. To remove this issue the first and last year was removed for
Simulations 1 and 2; the slightly reduced observational periods taken into
account in all further analysis. For convenience and brevity we continue to
refer to these simulation timespans as 10 000 or 100 000 years of data.

Uncertainties were calculated by computing a 95% confidence interval via
bootstrapping. To perform this bootstrapping, Simulations 1 and 2 were broken
into yearly observations and the various event rates (see Section 4) were
computed. These quantities were then sampled and combined to provide $10^{5}$
bootstrap samples. The 2.5th and 97.5th percentiles was then calculated for
each statistic, becoming the confidence interval. In the case of Simulation 3,
it already exists as 101 32-year timespans so these were sampled from instead
of the yearly divisions made for Simulations 1 and 2.

##  4 Results

###  4.1 Expected microlensing events in a year

Figure 3: Expected yearly distribution of blended astrometric shifts, bump
magnitudes and photometric magnifications caused by microlensing events where
the lenses are stars, BHs and NSs. The shaded regions indicate $3\sigma$
uncertainty or the undersampling factor if the histogram bin is empty. Note
that in the magnification plot the $x$ axis is the magnification minus one
(enabling a log scale can be used). Also note that the granularity of the
simulation affects the plots at the right-hand side — the stellar data only
spans 1/10th of a year and so (after correction with the multiplicative factor
of 10) the line falls as it reaches $\sim$10 events per year. A more fine-
grained simulation would show the stellar relationship smoothly continuing
below the 101 threshold. Similarly, the BH and NS lines fall at $\sim$0.1
events per year as the data accounts for 10 years of microlensing events.

The expected annual distributions of astrometric shifts, bump magnitudes and
photometric magnifications caused by microlensing events is given in Figure 3.
These plots were generated by taking the data set reached by Simulation 1 and
Simulation 2 and transforming them into a yearly distribution based on the
undersampling factor of the original distribution — 103 for NSs and BHs, 106
for stars. These results indicate that every year we can expect $88_{-6}^{+6}$
BH, $6.8_{-1.6}^{+1.7}$ NS and $20^{+30}_{-20}$ stellar microlensing events
which cause a blended astrometric shift larger than 2 mas ($231_{-9}^{+10}$,
$55_{-5}^{+5}$ and $690_{-160}^{+170}$ larger than 1 mas, respectively).
Similarly, we expect $21_{-3}^{+3}$ BH, $18_{-3}^{+3}$ NS and
$7500_{-500}^{+500}$ stellar events causing bump magnitudes $>1$ mag as well
as $5.5_{-1.4}^{+1.5}$, $4.3_{-1.2}^{+1.3}$ and $1600_{-200}^{+300}$ events
with bump magnitudes larger than 2.5. In terms of unblended photometric
magnification $26_{-3}^{+3}$ BH, $23_{-3}^{+3}$ NS and $12800_{-700}^{+700}$
stellar events cause magnifications of $>2$ (i.e. a doubling in brightness) as
well as $5.5_{-1.4}^{+1.5}$, $4.3_{-1.2}^{+1.3}$ and $2300_{-300}^{+300}$
events with magnifications larger than 10.

Figure 4: Expected yearly distribution of dimensionless separations (top
panel) and angular separations (bottom panel) caused by microlensing events
where the lenses are stars, BHs and NSs. The shaded regions indicate $3\sigma$
uncertainty or the undersampling factor if the histogram bin is empty. Note
that the stellar data only accounts for 1/10th of a year and so is multiplied
by 10, which is why the line falls as it reaches $\sim$10 events per year. The
stellar line likely continues below the 101 threshold and is just not picked
up by our simulation. This figure includes events down to a major image shift
of $1\mu$as, many of which will not be observable.

These results can be generalised to other all-sky surveys with similar
selection functions to Gaia, by scaling the number of stars surveyed in
comparison to Gaia, which encompasses approximately 2 billion stars.

###  4.2 Characterisation of microlensing events

Figure 5: Expected yearly distribution of microlensing events as a function of
Einstein time and type of lens. The largest bin/point in this figure was set
to be $6\,000$ days so the larger stellar undersampling doesn’t mislead the
viewer. BH, NS and some stellar events continue beyond this point. This figure
includes events down to a bump magnitude of $-10^{-10}$ mag, many of which
will not be observable.

The characterisation of microlensing events, particularly by Einstein time, is
important for identifying which events to follow-up to efficiently probe the
population of BHs and NSs. In particular, identifying events likely to be
caused by BHs is of great interest and furthermore their more extreme
properties renders them easier to detect. Figure 5 shows the distribution of
microlensing events as a function of Einstein time split by the type of lens.

These figures may be directly contrasted with Figure 7 from Lam et al. (2020)
which, using a different methodology, modelled microlensing events toward the
Galactic Bulge. Our results show that BHs make up a significantly smaller
fraction of long Einstein time events than the 50% of $t_{E}>120$ days that
was predicted in Lam et al. (2020). The discrepancy between the two studies
likely arises due to differences in the dynamical model: in the application of
the natal kick and subsequent evolution of the remnant. In addition, the two
methodologies take divergent approaches to arrive at their populations of
BGSs: an augmented GALAXIA distribution for the case of Lam et al. (2020) and
Gaia sources in the case of this work.

In modelling lensing BHs, the Lam et al. (2020) work applies natal kicks of
100 km/s with no subsequent evolution through the Galactic gravitational
potential over cosmic time. By contrast, the Sweeney et al. (2022) work
applies natal kicks which follow a bimodal Maxwellian distribution which, in
the case of BHs, has expected value of $\sim 77$ km/s. The BHs are then
evolved through the Galactic gravitational potential for the timespan from
birth to the present day. This acts to reduce the average speed of the BHs
while expanding their distribution to larger radii. The end result is that the
BHs from Lam et al. (2020) have a much larger average speed and are more
concentrated towards the galactic centre than Sweeney et al. (2022): features
which decrease the Einstein time and increases the event rates due to BHs.

This comparatively lower fraction of BHs for long Einstein time events renders
the strategy of identifying likely BH events using photometric lightcurves
alone much more difficult. Analysis of photometric surveys which leveraged the
Lam et al. (2020) finding that 50% of events with $t_{E}\geq 120$ days are BHs
may therefore significantly overestimate the number of events ascribed to BH
lensing (e.g. Golovich et al. (2022)’s reanalysis of OGLE-III and -IV events).
The detection of an isolated BH via microlensing (Lam et al., 2022; Sahu et
al., 2022) with an Einstein time of $\sim 250$ days constitutes a reasonably
unlikely outcome under our modelling assuming a small number of attempts.
However, without knowing the number of events monitored or the selection
criteria used it is impossible to draw any conclusions on how this detection
should impact the modelling. The characteristics of the event itself were also
quite unusual: the magnification was $\sim 400$ which is much larger than the
typical event, and the lens distance was found to be $\sim 1.6$ kpc, which is
unusually close (as shown in Figure 20).

In the case of NSs, Lam et al. (2020) apply a natal kick of 350 km/s as
opposed to the bimodal Maxwellian distribution with expectation of $\sim 450$
km/s in Sweeney et al. (2022). The evolution through the Galactic potential in
Sweeney et al. (2022) leads to a major reduction in speed and an increase in
the scale height of the NSs, as shown in Figure 1. In light of these
differences, it’s quite surprising that the distributions of microlensing
events caused by NSs are relatively similar between the two studies. The
similarity may be partially explained by the opposing nature of the two
effects: the reduced speed increases the Einstein time of events, but the
inflated distance to NSs decreases it.

The NSs’ migration to larger scale heights and slower speeds should
dramatically reduce the number of events compared to that of Lam et al.
(2020). This study finds $\sim 3$ times fewer NS events in ratio to BH events
compared to $\sim$an order of magnitude fewer from the Lam et al. (2020) work.
This cannot be explained by differing numbers of BHs and NSs in the two
underlying populations. The Lam et al. (2020) population contains $2\times$
more NSs than BHs. In comparison, the Sweeney et al. (2022) start with
$4.8\times$ more NSs, but the evolution of the remnants through the Galactic
potential (which is not done in Lam et al., 2020) acts to kick 40% of NSs out
of the Galaxy entirely. The remaining population then contains $2.3\times$
more NSs than BHs, if remnants $>20$ kpc of the Galactic centre are excluded.
This agreement in population make-up combined with the increased speed of the
Lam et al. (2020) NS population should result in many more NS events being
recorded in their counts. However, their low rate of NS events is likely due
to a separate model factor entirely: their simulations lack sensitivity to
short events which are much more likely due to NSs than BHs. This can be
traced to the Lam et al. (2020) simulated “observational cadence” which is
only every 10 days: a span that can easily miss many of these short events,
and most particularly those additional events generated by the large
velocities of their NS population.

Figure 6: Expected fractional distribution of microlensing events in a
scenario where BHs have a mass of $20$ M⊙. The largest bin/point in this
figure was set to be $6\,000$ days so the larger stellar undersampling doesn’t
mislead the viewer. BH, NS and some stellar events continue beyond this point.
This figure includes events down to a bump magnitude of $-10^{-10}$ mag, many
of which will not be observable.

Unsurprisingly, the fractional contributions seen in Figure 5 depend on the
masses of the lenses involved. For example, if BHs have a mass of $20$ M⊙
(instead of $7.8$ M⊙) but still receive the same natal kick velocity then the
fractional contribution as a function of Einstein time is shown in Figure 6.
This plot shows that the increased BH mass causes microlensing events to
become dominated by BHs at long Einstein times, but even with this extreme
mass BHs are only a few percent of events with $t_{E}\approx 250$ days. Event
rates are given for this example in Appendix E.

Figure 7: Fraction of events caused by BHs as a function of Einstein time for
a bump magnitude and an astrometric cut. The orange line has the same cut
applied as the orange line in Figure 5 while the red line instead applies an
astrometric cut.

The fractional contributions strongly depend on the selection criteria for
events. Figure 7 shows how the BH contribution varies for a photometric or
astrometric cut. The threshold of this cut does not significantly alter the
shape, as shown in Appendix B. Whether the events are located inside or
outside the bulge only makes a small difference to the BH fraction, as shown
in Appendix C. For microlensing events identified via photometric surveys, the
flatter curve of BH contribution is expected. The fraction of BH events is
much larger if a selection is made based on astrometric criteria, however
surveying the sky for small deviations caused by microlensing events is not
currently performed.

Figure 8: Scatter plot showing the splitting of stellar, BH and NS events
across the Einstein time–Einstein parallax plane. Note that there is a factor
of 100 difference in the undersampling between the stellar events and the BH
and NS events which is not accounted for in this figure. This figure only
includes events with a bump magnitude $<-2$ mag.

Instead of solely relying on the Einstein time as a predictor for lens type,
it is much better to predict based on both the Einstein time and the Einstein
parallax. We plot the scatter of lensing events with a bump magnitude $<-2$ in
Figure 8. If natal kicks are larger than those applied in this work then the
BH and NS events would migrate to smaller Einstein times, compressing the
$x$-axis of the figure. As discussed by Lam et al. (2020) and Gould (2023),
the $t_{E}$–$\pi_{E}$ basis splits lensing events by the mass of the lens,
with equal mass contours being approximately diagonal. Given the difficulty in
identifying BH microlensing events by Einstein time alone, shown in Figure 5,
these results suggest it is crucial to also measure the microlens parallax. It
should be noted that while BHs and NSs are reasonably distinct in Figure 8,
they are significantly less separate than in the scatter published by Lam et
al. (2020).

Figure 9: Expected yearly distribution of microlensing events as a function of
Einstein angle and type of lens. The largest bin in this figure was set to be
$15$ mas so the larger stellar undersampling doesn’t mislead the viewer, BH
events continue beyond this point. This figure includes events down to a bump
magnitude of $-10^{-10}$ mag, many of which will not be observable.

The distribution of Einstein angles for microlensing events is shown in Figure
9. We can see from this figure that microlensing events with Einstein angles
larger than $\sim$8 mas are mostly caused by BHs.

A curious result is that the lens-BGS position angle (i.e. the angle, east of
north, between the lens and BGS) a year in advance of the microlensing event
is not uniformly distributed, as shown in Figure 10. We suggest the non-
uniformity in this distribution is due to the net proper motion observed in
the sky as a result of our motion through the Galaxy.

In Appendix F we provide plots of the lens masses ($M_{L}$) and lens distances
($D_{L}$) broken down by type of lens. These are included so that they may be
used as priors for observation analysis.

Figure 10: Plot showing the lens-BGS position angle a year in advance of the
microlensing event, split by the type of lens.This figure includes events down
to a bump magnitude of $-10^{-10}$ mag, many of which will not be observable.

###  4.3 Comparison to Gaia

Figure 11: Comparison between the distribution of Einstein times (top panel)
and dimensionless separations (bottom panel) found in Gaia data (Wyrzykowski
et al., 2023) and those reached in this work. The distributions of Einstein
times and dimensionless separation show good agreement between this work and
the microlensing events found in the Gaia data. Microlensing events due to
stars dominate the distributions from this work. The distribution from this
work only includes events with a bump magnitude $<-0.4$ mag in order to
emulate a photometric selection criteria. The Gaia data plotted is the “Level
1” solutions derived in Wyrzykowski et al. (2023). Some outliers from the Gaia
data are beyond the bounds of this plot. The curves were generated using a
Gaussian kernel density estimator, so each curve will integrate to 1.

We can also compare the predictions of this work to the microlensing events
observed in Gaia photometry. Wyrzykowski et al. (2023) analyse the Gaia
photometry data from 2014 to 2017 for evidence of microlensing events. They
find evidence of 363 microlensing events, of which 163 are identified based on
Gaia data alone and the other 200 are confirmed to be present in the Gaia data
after a reexamination due to the events being reported in the OGLE-IV survey
(Mróz et al., 2019, 2020) or the All-Sky Automated Survey for Supernovae
(ASAS-SN; Shappee et al., 2014). This inclusion biases the sample due to the
increased detection efficiency towards the bulge where OGLE-IV survey fields
lie. Gaia struggles in crowded regions, such as the bulge, and so this
positive bias acts to offset the reduced native detection efficiency due to
Gaia’s sampling.

To compare the number of events in Gaia data with our simulation we impose
cuts to emulate Wyrzykowski et al. (2023)’s selection criteria. We filter both
sets of data to only look at the bulge (Galactic longitude $|l|<10$ degrees,
where OGLE assists the most with candidate selection), only consider $<19$ G
magnitude BGSs and events with Einstein times 20 days $<t_{E}<$ 250 days. With
this criteria, we find that there are $100$ microlensing events per year with
bump magnitudes $<-3$ mag, of which $<0.1$ are due to BHs and $0.1$ due to
NSs. In comparison, Wyrzykowski et al. (2023)’s “Level 1” analysis of Gaia
data finds there to be $\sim 50$ microlensing events per year meeting the
above criteria, in reasonable agreement with our simulations. This suggests
that it is quite unlikely that many BHs can be found amongst the events
identified by Wyrzykowski et al. (2023). Note that the cuts we apply don’t
take into account binaries, variable stars or blending with neighbouring
stars, all of which would make microlensing detections much more difficult.

We can also compare the distribution of Einstein times and dimensionless
separations between the bright events (bump magnitude $<-0.4$) in this work
and the distributions from events in Gaia. This comparison is shown in Figure
11. It shows that the distribution of Einstein times and dimensionless
separations in this work and Gaia data are in good agreement. Our bump
magnitude threshold was chosen so that the “plateau” of the dimensionless
separations from this work approximately matches peak width from the Gaia
data. Decreasing or increasing this threshold simply shortens or lengthens the
plateau as events with greater separations are removed or included. We also
plot the distribution of Einstein times for the bulge/non-bulge population in
Appendix C.

Figure 12: Comparison between the distribution of microlensing events over
galactic longitude (top panel) and galactic latitude (bottom panel) found in
Gaia data (Wyrzykowski et al., 2023) and those reached in this work. Overall,
the distribution of event locations shows good agreement between the
distributions. Microlensing events due to stars dominate the distributions
from this work. The distribution from this work only includes events with a
bump magnitude $<-0.4$ mag in order to emulate a photometric selection
criterion. The curves were generated using a Gaussian kernel density
estimator, so each curve will integrate to 1.

The locations of the events detected in Gaia data can also be compared to the
locations of events from this work. Figure 12 shows this comparison between
distributions. It shows that the locations of events largely agree, but that
this work predicts more events toward the bulge than are seen in the Gaia
data. This could be due to Gaia’s inefficiencies in making measurements in
crowded stellar fields, such as the bulge, and suggests losses are not
entirely offset by previously mentioned gains in candidate selection from
photometric surveys such as OGLE-IV. Leaving aside difficulties in crowded
stellar fields, Gaia’s detections are also impacted by its scanning law (Gaia
Collaboration et al., 2016), with the bulge being one of the areas least
frequently revisited. The scanning law also explains the positive-negative
longitude asymmetry for the Gaia data. During the period considered by
Wyrzykowski et al. (2023), the region between $-15$ and $-50$ longitude was
typically scanned $2$–$3\times$ more frequently than its positive counterpart.
As discussed in Wyrzykowski et al. (2023), a reduction in the number of data
points and sampling frequency reduces the number of detections that can be
made on the basis of Gaia photometry, explaining the larger number of events
identified by Wyrzykowski et al. (2023) at negative longitudes.

The lower panel of Figure 12 shows that our simulations yield two distinct
peaks whereas the Gaia data has one major and one minor peak. The peak at a
Galactic latitude of $\sim-2$ is a close match between both sets of data. Both
data sets then dip at 0 degrees, likely due to Gaia failing to observe stars
in the dusty Galactic plane. At $\sim 2$ degrees our data exhibits a peak
mirroring that at $\sim-2$, which makes sense — intuition would have the
microlensing event distribution be approximately symmetrical around the
Galactic plane. The Gaia data exhibits a much smaller peak in this region.
Wyrzykowski et al. (2023)’s analysis is heavily reliant on candidates provided
from OGLE IV to detect events in near the Galactic plane. OGLE IV samples the
$\sim 2$ latitude region much less frequently than the $\sim-2$ region, so
Wyrzykowski et al. (2023)’s asymmetric peaks are likely a selection effect
arising from their reliance on OGLE.

###  4.4 Predictions for future data releases

Photometric microlensing events should continue to be found in Gaia data at
the same rate as they were detected in Wyrzykowski et al. (2023). This rate
changes throughout the study, but using the penultimate year of observation
(to avoid edge effects) we predict 126 microlensing events per year. Thus, in
Gaia DR4 (which will contain 5 years of data) we expect 640 microlensing
events and in a potential Gaia DR5 (containing 10.5 years of data) we expect
1300 microlensing events.

We can estimate the uncertainty of raw positions that will appear in DR4/5
using the recent Gaia focused product release which contains the observed raw
position for quasars (Gaia Collaboration et al., 2023a). By calculating the
centroid of each source we reach a standard deviation $\sim 50$ mas in each
individual measurement. A typical Gaia source is observed $\sim 20$ times per
year for a total of 210 measurements over the full 10.5 years of observation.
Assuming the uncertainty scales with the square root of the number of
observations, the uncertainty of a typical object analysed in isolation will
be $\sim 3$ mas — much too large to detect astrometric microlensing events
which are typically 5 mas at most. This result highlights the necessity of
including a microlensing model into the astrometric solution of the Gaia data
analysis pipeline. With a microlensing model included, astrometric
microlensing events should be detectable below the 3 mas level, as
demonstrated by the tight astrometric solutions the existing Gaia pipeline
provides for sources.

A similar analysis applied to the Roman Galactic Bulge Time-Domain Survey
which is currently expected to achieve 1 mas uncertainty on individual
astrometric measurements with $\sim 7000$ measurements made over 6 observing
seasons results in astrometric microlensing events being detectable as low as
20 $\mu$as. If achieved, this would result in colossal numbers of microlensing
events being detected.

###  4.5 Predictions for GaiaNIR

Figure 13: Expected distribution for GaiaNIR of astrometric shifts, bump
magnitudes and lensing magnifications caused by microlensing events where the
lenses are BHs and NSs. The shaded regions indicate $3\sigma$ uncertainty.
Note that in the magnification plot the $x$ axis is (magnification - 1),
permitting a log scale. The relationships will continue below the 101
threshold, but are below the sensitivity of our simulation.

In addition to comparing our results to existing Gaia data, we can make
predictions for microlensing events observed by the combination of Gaia and
GaiaNIR data. As described in Section 3.3, Simulation 3 contains 101
“observation runs” based on an observing window of 12 years (to mimic the Gaia
observation window), then a break of 15 years followed by another observing
window of 5 years intended to emulate GaiaNIR observations (assuming GaiaNIR
begins observing in 2041). By performing 101 such simulations we undersample
the true distribution by a factor of $\sim$10, since the distribution from
Sweeney et al. (2022) undersamples by a factor of 1 000.

The predicted distributions of centroid shifts, bump magnitudes and
photometric magnifications are shown in Figure 13. These results show that the
combination of Gaia and GaiaNIR data is likely to contain $4200_{-400}^{+400}$
BH and $190_{-80}^{+90}$ NS microlensing events with an astrometric shift
larger than 2 mas ($14700_{-900}^{+600}$ BH and $1600_{-200}^{+300}$ NS events
larger than 1 mas). Similarly, there are likely to be $330_{-120}^{+100}$ BH
and $310_{-100}^{+110}$ NS events with bump magnitudes $<-1$ mag and
$90_{-60}^{+60}$ and $60_{-40}^{+50}$ BH and NS events larger with bump
magnitudes $<-2.5$ mag. If we use half the number of events with bump
magnitude $\geq 3$ mag (lensing magnification $\geq 15$), G magnitude $<19$
and $t_{E}>20$ days as a proxy for Wyrzykowski et al. (2023)’s selection
criteria then the same analysis is likely to detect $<10$ events from BHs and
$5_{-5}^{+10}$ from NSs in the combination of Gaia and GaiaNIR data. Note that
these numbers are optimistic as the approximation of the Wyrzykowski et al.
(2023) selection function is for their optimal detection efficiency, not the
entire sky.

By contrast, our yearly results (Section 4.1) would estimate only a fraction
of the astrometric event counts. This apparent surplus of microlensing events
in our combination of Gaia and GaiaNIR data is due to our simulation
“detecting” many events which do not peak during the observation windows. In
the analysis by Wyrzykowski et al. (2023) they were able to detect events
which peaked outside their observation period so we expected many of these
events with longer Einstein times will be recoverable. The observation windows
also lead to more 1 mas deflection BH events being detected in comparison to
$1$ mas deflection NSs events, $9.2\times$ more as opposed to $4.2\times$ more
for events peaking in a typical year. This increase in the fraction of BHs
being detected is because of their generally longer Einstein times, so there
are more “heads” and “tails” of BH events occurring at the start/end of any
given observation period.

Figure 14: Figure showing the fraction of the BH or NS population which
experiences at least one microlensing event larger than an astrometric or
photometric threshold in the combination of Gaia and GaiaNIR data.

The fraction of the BH or NS population which experiences at least one
microlensing event is a function of both the observing period and the
detection criteria for a microlensing event. Figure 14 shows this lensing
fraction for BHs and NSs as a function of centroid shift and bump magnitude
for the combination of Gaia and GaiaNIR data. From this figure we can see that
in order to probe a reasonable fraction of the isolated BH or NS population we
must be able to identify lensing events which cause small astrometric shifts.

##  5 Conclusions

This paper has explored the rates at which microlensing events occur across
the entire sky. This was modelled on both a yearly basis and for the ongoing
Gaia mission combined with a future GaiaNIR mission. Microlensing event rates
were calculated by performing an N-body simulations using the BH and NS
population published by Sweeney et al. (2022) and a stellar population
generated by GALAXIA as the lensing populations. The population of BGSs was
drawn from the Gaia catalogue so the selection function of stars monitored for
microlensing events would be highly realistic.

The key methodological advancements in this work are (1) the population of NSs
and BHs modelled with realistic natal kicks and subsequently evolved through
the Galactic potential and (2) the use of Gaia stars as BGSs to ensure a
realistic selection function. The results of these simulations were compared
to microlensing events detected in Gaia data (Wyrzykowski et al., 2023) and
were found to be consistent with these observations. Our main findings are as
follows:

  * •

On a yearly basis we can expect $88_{-6}^{+6}$ BH, $6.8_{-1.6}^{+1.7}$ NS and
$20^{+30}_{-20}$ stellar microlensing events which cause a blended astrometric
shift larger than 2 mas. For events causing an astrometric shift larger than 1
mas, we expect $231_{-9}^{+10}$ BH events, $55_{-5}^{+5}$ NS events and
$690_{-160}^{+170}$ stellar events. Similarly, we can expect $21_{-3}^{+3}$
BH, $18_{-3}^{+3}$ NS and $7500_{-500}^{+500}$ stellar microlensing events
with bump magnitudes larger than 1 mag and $5.5_{-1.4}^{+1.5}$,
$4.3_{-1.2}^{+1.3}$ and $1600_{-200}^{+300}$ events with bump magnitudes
larger than 2.5.

  * •

Our results update the BH and NS event statistics from Lam et al. (2020) who
found that events with $t_{E}>120$ days have a 50% probability of being a BH.
We find that when the BHs are evolved post-natal-kick to the present day the
fraction of BH events at $t_{E}\approx 120$ days is much lower, due to both an
increase in the Einstein times of BH events as well as a decrease in the
number of BH events. Here we find that 28% of $t_{E}>1\,000$ day events are
caused by BHs if events are selected by a bump magnitude cut. This fraction is
larger if events are selected by an astrometric cut.

  * •

We also consider the scenario where BHs have a mass of 20 M⊙, but receiving
the same natal kick. This leads to BHs causing $\sim$50% more photometric
microlensing events, still mostly at very long Einstein times.

  * •

In agreement with Lam et al. (2020) and Gould (2023), we find that
measurements of both $\pi_{E}$ and $t_{E}$ are required to identify BH events.

  * •

We provide distributions of lens masses and lens distances broken down by
species that may be valuable as informative priors for observational data
analysis.

  * •

The microlensing events found in Gaia data (Wyrzykowski et al., 2023) are
consistent with our data. Considering mid-length ($20<t_{E}<250$) bulge events
($|l|<10$) for bright ($G<19$ mag) Gaia stars, emulating the filters in
Wyrzykowski et al. (2023)’s analysis, we predict 100 yearly events with bump
magnitudes $<-3$ mag compared to Wyrzykowski et al. (2023)’s $\sim 50$ yearly
detected events. We find that the distribution of Einstein times and their
longitude also show very good agreement. Our data is symmetrical around a 0
Galactic latitude, whereas Wyrzykowski et al. (2023)’s analysis is not. We
conclude that this is due to the selection bias in the OGLE-IV survey upon
which their analysis depends on for candidate selection.

  * •

We predict 640 microlensing events will be detected in the photometry of Gaia
DR4 and 1300 microlensing events in Gaia DR5. We use a recent Gaia focused
product release (Gaia Collaboration et al., 2023a) to estimate the number of
astrometric events which will be recoverable from the raw Gaia positions to be
released in DR4/5 and conclude that the uncertainties will be too large to
detect any events unless a microlensing model is included into the Gaia data
analysis pipeline. We briefly consider astrometric events which will be
detectable in the Roman Galactic Bulge Time Domain Survey and conclude that if
they achieve the expected 1 mas uncertainty on individual measurements then
microlensing events as small as 20 $\mu$as may be detectable.

  * •

We predict $4200_{-400}^{+400}$ BH and $190_{-80}^{+90}$ NS microlensing
events will be present in the combination of Gaia and GaiaNIR data with
centroid shifts larger than 2 mas ($14700_{-900}^{+600}$ BH and
$1600_{-200}^{+300}$ NS events larger than 1 mas). Similarly, we predict
$330_{-120}^{+100}$ BH and $310_{-100}^{+110}$ NS events with bump magnitudes
$<-1$ mag ($90_{-60}^{+60}$ BH and $60_{-40}^{+50}$ NS events $<-2.5$ mag).
Using half the number of events with bump magnitudes $<-3$ mag, G magnitude
$<19$ and $t_{E}>20$ days as a proxy for Wyrzykowski et al. (2023)’s selection
function, we predict that $<10$ BH and $5_{-5}^{+10}$ NS events will be
detectable in the combination of Gaia and GaiaNIR data if the same analysis is
performed.

The recent detection of an isolated BH (Lam et al., 2022; Sahu et al., 2022)
has demonstrated that it is possible to peer into the hidden world of dark
compact remnants. Observing these electromagnetically invisible objects
provides a probe for massive binary evolution and supernovae through Galactic
history. We characterise these signatures as precisely as modern theory allows
and make predictions against which the observational community can test future
measurements. In particular, the BH and NS events are strongly dependant on
the natal kicks received during supernovae and the masses of the compact
remnants. Sweeney et al. (2022) established that 40% of NSs are ejected from
the Galaxy due to these natal kicks, which significantly reduces the number of
microlensing events expected from NSs. Sufficient statistics from future
studies would allow constraints to be placed on these natal kicks and masses,
both of which are uncertain — especially in the case of BHs.

## Acknowledgements

We would like to thank Subo Dong, Tim Bedding, James Crowley, Jessica Lu,
Natasha Abhrams, Casey Lam and Lukasz Wyrzykowski for their helpful
conversations. This collaboration work was in part made possible by ESO’s
early-career scientific visitor programme. AKM acknowledges partial support
from the Portuguese Fundação para a Ciência e a Tecnologia project EXPL/FIS-
AST/1368/2021. This work has made use of data from the European Space Agency
(ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia
Data Processing and Analysis Consortium (DPAC,
https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has
been provided by national institutions, in particular the institutions
participating in the Gaia Multilateral Agreement. In addition to those
mentioned in the main text, we would also like to acknowledge other python
packages used as part of this work: NumPy (Harris et al., 2020), Matplotlib
(Hunter, 2007), seaborn (Waskom, 2021), pandas (Wes McKinney, 2010; pandas
development team, 2020), Astropy (Astropy Collaboration et al., 2013, 2018,
2022), astroquery (Ginsburg et al., 2019) and ray (Moritz et al., 2017).

## Data Availability

The data and code underlying this article are available on Zenodo (Sweeney et
al., 2024), as well as at the following GitHub url: https://github.com/David-
Sweeney/Incidence.

## References

  * Abbott et al. (2016) Abbott B. P., et al., 2016, Phys. Rev. Lett., 116, 061102
  * Adams et al. (2017) Adams S. M., Kochanek C. S., Gerke J. R., Stanek K. Z., Dai X., 2017, MNRAS, 468, 4968
  * Astropy Collaboration et al. (2013) Astropy Collaboration et al., 2013, A&A, 558, A33
  * Astropy Collaboration et al. (2018) Astropy Collaboration et al., 2018, AJ, 156, 123
  * Astropy Collaboration et al. (2022) Astropy Collaboration et al., 2022, ApJ, 935, 167
  * Aubourg et al. (1993) Aubourg E., et al., 1993, The Messenger, 72, 20
  * Belczynski & Taam (2008) Belczynski K., Taam R. E., 2008, The Astrophysical Journal, 685, 400
  * Bennett et al. (1993) Bennett D. P., et al., 1993, in Akerlof C. W., Srednicki M. A., eds, Vol. 688, Texas/PASCOS 1992: Relativistic Astrophysics and Particle Cosmology. p. 612 (arXiv:astro-ph/9304014), doi:10.1111/j.1749-6632.1993.tb43945.x
  * Bovy (2015) Bovy J., 2015, The Astrophysical Journal Supplement Series, 216, 29
  * Busso et al. (2022) Busso G., et al., 2022, Gaia DR3 documentation Chapter 5: Photometric data, Gaia DR3 documentation, European Space Agency; Gaia Data Processing and Analysis Consortium. Online at <A href=“https://gea.esac.esa.int/archive/documentation/GDR3/index.html”>https://gea.esac.esa.int/archive/documentation/GDR3/index.html</A>, id. 5
  * Dai et al. (2015) Dai S., Smith M. C., Lin M. X., Yue Y. L., Hobbs G., Xu R. X., 2015, ApJ, 802, 120
  * Dominik & Sahu (2000) Dominik M., Sahu K. C., 2000, ApJ, 534, 213
  * Einstein (1936) Einstein A., 1936, Science, 84, 506
  * El-Badry et al. (2023) El-Badry K., et al., 2023, MNRAS, 521, 4323
  * Fryer (1999) Fryer C. L., 1999, The Astrophysical Journal, 522, 413
  * Fryer et al. (2012) Fryer C. L., Belczynski K., Wiktorowicz G., Dominik M., Kalogera V., Holz D. E., 2012, The Astrophysical Journal, 749, 91
  * Gaia Collaboration et al. (2016) Gaia Collaboration et al., 2016, A&A, 595, A1
  * Gaia Collaboration et al. (2023a) Gaia Collaboration et al., 2023a, arXiv e-prints, p. arXiv:2310.06295
  * Gaia Collaboration et al. (2023b) Gaia Collaboration et al., 2023b, A&A, 674, A1
  * Ginsburg et al. (2019) Ginsburg A., et al., 2019, AJ, 157, 98
  * Golovich et al. (2022) Golovich N., et al., 2022, ApJS, 260, 2
  * Gould (2023) Gould A., 2023, arXiv e-prints, p. arXiv:2310.19164
  * Han & Gould (1995) Han C., Gould A., 1995, ApJ, 447, 53
  * Han & Gould (2003) Han C., Gould A., 2003, ApJ, 592, 172
  * Harris et al. (2020) Harris C. R., et al., 2020, Nature, 585, 357
  * Heger et al. (2003) Heger A., Fryer C. L., Woosley S. E., Langer N., Hartmann D. H., 2003, ApJ, 591, 288
  * Hobbs et al. (2005) Hobbs G., Lorimer D. R., Lyne A. G., Kramer M., 2005, MNRAS, 360, 974
  * Hobbs et al. (2016) Hobbs D., et al., 2016, arXiv e-prints, p. arXiv:1609.07325
  * Hunter (2007) Hunter J. D., 2007, Computing in Science & Engineering, 9, 90
  * Igoshev (2020) Igoshev A. P., 2020, Monthly Notices of the Royal Astronomical Society, 494, 3663
  * Katsuda et al. (2018) Katsuda S., et al., 2018, ApJ, 856, 18
  * Kim et al. (2016) Kim S.-L., et al., 2016, Journal of Korean Astronomical Society, 49, 37
  * Klüter et al. (2022) Klüter J., Bastian U., Demleitner M., Wambsganss J., 2022, AJ, 163, 176
  * Lam & Lu (2023) Lam C. Y., Lu J. R., 2023, ApJ, 955, 116
  * Lam et al. (2020) Lam C. Y., Lu J. R., Hosek Matthew W. J., Dawson W. A., Golovich N. R., 2020, ApJ, 889, 31
  * Lam et al. (2022) Lam C. Y., et al., 2022, The Astrophysical Journal Letters, 933, L23
  * Luri et al. (2018) Luri X., et al., 2018, A&A, 616, A9
  * Lyne & Lorimer (1994) Lyne A. G., Lorimer D. R., 1994, Nature, 369, 127
  * Mao (2008) Mao S., 2008, arXiv e-prints, p. arXiv:0811.0441
  * Moritz et al. (2017) Moritz P., et al., 2017, arXiv e-prints, p. arXiv:1712.05889
  * Mróz & Wyrzykowski (2021) Mróz P., Wyrzykowski Ł., 2021, Acta Astron., 71, 89
  * Mróz et al. (2019) Mróz P., et al., 2019, ApJS, 244, 29
  * Mróz et al. (2020) Mróz P., et al., 2020, ApJS, 249, 16
  * Mróz et al. (2022) Mróz P., Udalski A., Gould A., 2022, ApJ, 937, L24
  * Nagarajan et al. (2023) Nagarajan P., et al., 2023, arXiv e-prints, p. arXiv:2312.05313
  * Olejak et al. (2020) Olejak A., Belczynski K., Bulik T., Sobolewska M., 2020, A&A, 638, A94
  * Osłowski et al. (2008) Osłowski S., Moderski R., Bulik T., Belczynski K., 2008, A&A, 478, 429
  * Özel et al. (2010) Özel F., Psaltis D., Narayan R., McClintock J. E., 2010, ApJ, 725, 1918
  * Özel et al. (2012) Özel F., Psaltis D., Narayan R., Santos Villarreal A., 2012, ApJ, 757, 55
  * Paczynski (1986) Paczynski B., 1986, ApJ, 301, 503
  * Pejcha & Thompson (2015) Pejcha O., Thompson T. A., 2015, ApJ, 801, 90
  * Postnov & Yungelson (2014) Postnov K. A., Yungelson L. R., 2014, Living Reviews in Relativity, 17, 3
  * Robin et al. (2003) Robin A. C., Reylé C., Derrière S., Picaud S., 2003, A&A, 409, 523
  * Rose et al. (2022) Rose S., Lam C. Y., Lu J. R., Medford M., Hosek M. W., Abrams N. S., Ramey E., Vasylyev S. S., 2022, ApJ, 941, 116
  * Sahu et al. (2017) Sahu K. C., et al., 2017, Science, 356, 1046
  * Sahu et al. (2022) Sahu K. C., et al., 2022, ApJ, 933, 83
  * Schlegel et al. (1998) Schlegel D. J., Finkbeiner D. P., Davis M., 1998, ApJ, 500, 525
  * Shappee et al. (2014) Shappee B. J., et al., 2014, ApJ, 788, 48
  * Sharma et al. (2011) Sharma S., Bland-Hawthorn J., Johnston K. V., Binney J., 2011, ApJ, 730, 3
  * Smartt (2009) Smartt S. J., 2009, ARA&A, 47, 63
  * Smartt (2015) Smartt S. J., 2015, Publications of the Astronomical Society of Australia, 32, e016
  * Spera et al. (2015) Spera M., Mapelli M., Bressan A., 2015, MNRAS, 451, 4086
  * Sukhbold et al. (2016) Sukhbold T., Ertl T., Woosley S. E., Brown J. M., Janka H. T., 2016, ApJ, 821, 38
  * Sweeney et al. (2022) Sweeney D., Tuthill P., Sharma S., Hirai R., 2022, MNRAS, 516, 4971
  * Sweeney et al. (2024) Sweeney D., Tuthill P., Krone-Martins A., Mérand A., Scalzo R., Martinod M.-A., 2024, Observing the Galactic Underworld: Code supplement ot predicting photomertric and astrometric microlensing events from compact remnants, doi:10.5281/zenodo.10840694, https://doi.org/10.5281/zenodo.10840694
  * Udalski et al. (1992) Udalski A., Szymanski M., Kaluzny J., Kubiak M., Mateo M., 1992, Acta Astron., 42, 253
  * Waskom (2021) Waskom M. L., 2021, Journal of Open Source Software, 6, 3021
  * Wes McKinney (2010) Wes McKinney 2010, in Stéfan van der Walt Jarrod Millman eds, Proceedings of the 9th Python in Science Conference. pp 56 – 61, doi:10.25080/Majora-92bf1922-00a
  * Wood & Mao (2005) Wood A., Mao S., 2005, MNRAS, 362, 945
  * Wyrzykowski et al. (2023) Wyrzykowski Ł., et al., 2023, A&A, 674, A23
  * pandas development team (2020) pandas development team T., 2020, pandas-dev/pandas: Pandas, doi:10.5281/zenodo.3509134, https://doi.org/10.5281/zenodo.3509134

##  Appendix A Gaia null values

The Gaia catalogue has null values scattered throughout, as you would expect
from a large survey. If our simulation queries the Gaia catalogue and receives
a null value we replace it with a value of 0 if the null value is a parallax
or radial velocity. If the null value is a proper motion in right ascension
then the value is set to $-1.885932688110818$; a proper motion in declination
is set to $-3.9745764127030273$. These proper motion values were chosen as
they were the median values of some early trial runs. We do not replace null
values in the G magnitude column. Where they exist the null values propagate
through to the blended centroid shift and bump magnitude columns (making these
values null).

##  Appendix B BH event fraction

Plotted in Figure 15 and Figure 16 are the fraction of events which are caused
by BHs as a function of time for various selection criteria applied to the
population. We can see that the threshold of the cut does not affect the shape
of the contribution curve, if the noise at high Einstein times is ignored
(where the number of stellar events is extremely low due to the undersampling
of the stellar population). Instead, the shape of the curve is mainly
determined by the type of cut applied: photometric or astrometric. For
microlensing events collected from photometric surveys the flatter curve of BH
contribution is expected.

Figure 15: Fraction of events caused by BHs as a function of Einstein time for
an astrometric cut and various bump magnitude cuts. These lines are the
counterpart of the orange line in Figure 5 depending on the selected events.
Figure 16: Fraction of events caused by BHs as a function of Einstein time for
various astrometric cuts. These lines are the counterpart of the orange line
in Figure 5 depending on the selected events.

##  Appendix C Einstein time variation in and outside the bulge

Defining the bulge as $|l|<30$ deg and $|b|<10$ deg, we can plot the fraction
of event which have BH lenses as a function of Einstein time for the bulge,
non-bulge and total population, as shown in Figure 17. This figure shows that
there is only a slight difference in the BH fraction inside/outside the bulge.
For completeness we also plot the results from (Lam et al., 2020), ignoring
their events due to white dwarfs (which they calculate while we don’t).

Figure 17: Fraction of yearly microlensing events caused by BHs as a function
of Einstein time, split by population and compared to the results from Lam et
al. (2020). The largest bin/point in this figure was set to be $6\,000$ days
so the larger stellar undersampling doesn’t mislead the viewer. BH, NS and
some stellar events continue beyond this point. This figure includes events
down to a bump magnitude of $-10^{-10}$ mag, many of which will not be
observable.

We also plot the distribution of Einstein times for events inside and outside
the bulge, as shown in Figure 18. We overplot this on Figure 11 which shows
the distribution for our entire population and the distribution from
Wyrzykowski et al. (2023).

Figure 18: Comparison between the distribution of Einstein times found in Gaia
data (Wyrzykowski et al., 2023) and those reached in this work, split by
whether they were/were not in the bulge. Microlensing events due to stars
dominate the distributions from this work. The distribution from this work
only includes events with a bump magnitude $<-0.4$ mag in order to emulate a
photometric selection criteria. The Gaia data plotted is the “Level 1”
solutions derived in Wyrzykowski et al. (2023). Some outliers from the Gaia
data are beyond the bounds of this plot. The curves were generated using a
Gaussian kernel density estimator, so each curve will integrate to 1.

##  Appendix D Distribution of microlens parallax

Plotted in Figure 19 is the yearly distribution of microlens parallax. We can
see from the figure that stellar events dominate the distribution at all
values of $|\pi_{E}|$, but that BH and NS events are most present at small
values of $|\pi_{E}|$. Note the uncertainty in the Gaia measurement for
parallax often results in a negative parallax value. For this reason we plot
the absolute value of $\pi_{E}$.

Figure 19: Expected yearly distribution of microlens parallax, shown for
events down to a bump magnitude of $-2$ mag. The largest bin in this figure
was set to be 2 mas, some events do occur beyond this point.

##  Appendix E Event rates for 20 M⊙ black holes

If the population of BHs have a mass of 20 M⊙ then they cause
$670_{-19}^{+19}$ and $268_{-11}^{+11}$ events with centroid shifts larger
than $1$ mas and $2$ mas, respectively. They cause $32_{-4}^{+4}$ and
$8.1_{-1.7}^{+1.8}$ events with bump magnitudes $<-1$ and $<-2.5$,
respectively. For comparison, with a mass of 7.8 M⊙ the respective numbers are
$231_{-9}^{+10}$, $88_{-6}^{+6}$, $21_{-3}^{+3}$ and $5.5_{-1.4}^{+1.5}$ (in
the same order). Thus, the $2.6\times$ increase in BH mass ($1.6\times$
increase in $\theta_{E}$) leads to a $\sim 3\times$ increase in astrometric
signals and a $\sim 1.5\times$ increase in photometric signals.

##  Appendix F Distribution of lens distances and masses

Figure 20: Distribution of lens distances for microlensing events caused by
different populations. The distribution does not change significantly for
different astrometric or photometric cuts. This figure includes events down to
a major image shift $>1$ $\mu$as, many of which will not be observable. Figure
21: Distribution of lens masses for stellar lenses. There is a mild dependence
on the photometric cut applied, so this figure includes events with bump
magnitude $<-0.01$ mag. The distribution is not significantly different for
events with bump magnitude $<-0.1$ mag. The mass distribution for BHs and NSs
is not plotted as they were all taken to have a mass of 7.8 M⊙ and 1.35 M⊙,
respectively.