# Analysis of the Full Spitzer Microlensing Sample. I. Dark Remnant Candidates and Gaia Predictions

Krzysztof A. Rybicki, Yossi Shvartzvald, Jennifer C. Yee, Sebastiano Calchi Novati, Eran O. Ofek, Ian A. Bond, Charles Beichman, Geoff Bryden, Sean Carey, Calen Henderson, et al.

► To cite this version:

Krzysztof A. Rybicki, Yossi Shvartzvald, Jennifer C. Yee, Sebastiano Calchi Novati, Eran O. Ofek, et al.. Analysis of the Full Spitzer Microlensing Sample. I. Dark Remnant Candidates and Gaia Predictions. The Astrophysical Journal, 2024, 975, 10.3847/1538-4357/ad7bb1. insu-04851492

HAL Id: insu-04851492

https://insu.hal.science/insu-04851492v1

Submitted on 21 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Analysis of the Full Spitzer Microlensing Sample. I. Dark Remnant Candidates and Gaia Predictions

Krzysztof A. Rybicki¹, Yossi Shvartzvald¹, Jennifer C. Yee², Sebastiano Calchi Novati³, Eran O. Ofek¹, Ian A. Bond⁴,
and

Charles Beichman³, Geoff Bryden³, Sean Carey³, Calen Henderson³, Wei Zhu⁵, Michael M. Fausnaugh⁶, Benjamin Wibking⁷
(The Spitzer Team),

Andrzej Udalski⁸, Radek Poleski⁸, Przemek Mróz⁸, Michał K. Szymański⁸, Igor Soszyński⁸, Paweł Pietrukowicz⁸,
Szymon Kozłowski⁸, Jan Skowron⁸, Krzysztof Ulaczyk⁸,⁹, Patryk Iwanek⁸, Marcin Wrona⁸

(The OGLE Collaboration),

Yoon-Hyun Ryu¹,¹⁰, Michael D. Albrow¹,¹¹, Sun-Ju Chung¹,¹⁰, Andrew Gould¹,¹²,¹³, Cheongho Han¹,¹⁴, Kyu-Ha Hwang¹,¹⁰,¹⁰
Youn Kil Jung¹,¹⁰,¹⁵, In-Gu Shin², Hongjing Yang¹,¹⁶, Weicheng Zang², Sang-Mok Cha¹,¹⁰,¹⁷, Dong-Jin Kim¹,¹⁰
Hyoun-Woo Kim¹, Seung-Lee Kim¹,¹⁰, Chung-Uk Lee¹,¹⁰, Dong-Joo Lee¹,¹⁰, Yongseok Lee¹,¹⁰,¹⁷, Byeong-Gon Park¹,¹⁰
Richard W. Pogge¹,¹⁸,¹⁹

(The KMTNet Collaboration),

and

Fumio Abe²⁰, Richard Barry²¹, David P. Bennett²¹,²², Aparna Bhattacharya²¹,²², Akihiko Fukui²³,²⁴, Ryusei Hamada²⁵,
Shunya Hamada²⁵, Naoto Hamasaki²⁵, Yuki Hirao²⁶, Stela Ishitani Silva²¹,²⁷, Yoshitaka Itow²⁰, Rintaro Kirikawa²⁵,
Naoki Koshimoto²⁵, Yutaka Matsubara²⁰, Shota Miyazaki²⁸, Yasushi Muraki²⁰, Tutumi Nagai²⁵, Kansuke Nunota²⁵,
Greg Olmschenk²¹, Clement Ranc²⁹, Nicholas J. Rattenbury³⁰, Yuki K. Satoh³¹, Takahiro Sumi²⁵, Daisuke Suzuki²⁵,
Paul J. Tristram³², Aikaterini Vandorou²¹,²², Hibiki Yama²⁵
(MOA collaboration),

Łukasz Wyrzykowski⁸, Kornel Howil⁸, and Katarzyna Kruszyńska³³

¹Department of Particle Physics and Astrophysics, Weizmann Institute of Science, Rehovot 76100, Israel

²Center for Astrophysics | Harvard & Smithsonian, 60 Garden St., MS-15, Cambridge, MA 02138, USA
³IPAC, Mail Code 100-22, Caltech, 1200 E. California Blvd., Pasadena, CA 91125, USA

⁴Institute of Natural and Mathematical Sciences, Massey University, Auckland 0745, New Zealand

⁵Department of Astronomy, Tsinghua University, Beijing 10084, People's Republic of China

⁶Department of Physics & Astronomy, Texas Tech University, Lubbock TX, 79410-1051, USA

⁷Michigan State University, 426 Auditorium Road, East Lansing, MI 48824, USA

⁸Astronomical Observatory, University of Warsaw, Al. Ujazdowskie 4, 00-478, Warszawa, Poland

⁹Department of Physics, University of Warwick, Coventry CV4 7 AL, UK

¹⁰Korea Astronomy and Space Science Institute, Daejeon 34055, Republic of Korea

¹¹University of Canterbury, School of Physical and Chemical Sciences, Private Bag 4800, Christchurch 8020, New Zealand

¹²Max-Planck-Institute for Astronomy, Königstuhl 17, 69117 Heidelberg, Germany

¹³Department of Astronomy, Ohio State University, 140 W. 18th Ave., Columbus, OH 43210, USA

¹⁴Department of Physics, Chungbuk National University, Cheongju 28644, Republic of Korea

¹⁵National University of Science and Technology (UST), Daejeon 34113, Republic of Korea

¹⁶Department of Astronomy, Tsinghua University, Beijing 100084, People's Republic of China

¹⁷School of Space Research, Kyung Hee University, Yongin, Kyeonggi 17104, Republic of Korea

¹⁸Department of Astronomy, Ohio State University, 140 West 18th Ave., Columbus, OH 43210, USA

¹⁹Center for Cosmology and AstroParticle Physics, Ohio State University, 191 West Woodruff Ave., Columbus, OH 43210, USA

²⁰Institute for Space-Earth Environmental Research, Nagoya University, Nagoya 464-8601, Japan

²¹Code 667, NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA

²²Department of Astronomy, University of Maryland, College Park, MD 20742, USA

²³Department of Earth and Planetary Science, Graduate School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan

²⁴Instituto de Astrofísica de Canarias, Via Lactea s/n, E-38205 La Laguna, Tenerife, Spain

²⁵Department of Earth and Space Science, Graduate School of Science, Osaka University, Toyonaka, Osaka 560-0043, Japan

²⁶Institute of Astronomy, Graduate School of Science, The University of Tokyo, 2-21-1 Osawa, Mitaka, Tokyo 181-0015, Japan
²⁷Department of Physics, The Catholic University of America, Washington, DC 20064, USA

²⁸Institute of Space and Astronautical Science, Japan Aerospace Exploration Agency, 3-1-1 Yoshinodai, Chuo, Sagamihara, Kanagawa 252-5210, Japan

²⁹Sorbonne Université, CNRS, UMR 7095, Institut d'Astrophysique de Paris, 98 bis bd Arago, 75014 Paris, France

³⁰Department of Physics, University of Auckland, Private Bag 92019, Auckland, New Zealand

³¹College of Science and Engineering, Kanto Gakuin University, 1-50-1 Mutsuurahigashi, Kanazawa-ku, Yokohama, Kanagawa 236-8501, Japan

³²University of Canterbury Mt. John Observatory, P.O. Box 56, Lake Tekapo 8770, New Zealand

³³Las Cumbres Observatory, 6740 Cortona Drive, Suite 102, Goleta, CA 93117, USA

Received 2024 May 7; revised 2024 August 15; accepted 2024 September 4; published 2024 November 5

<figure><img src="image_3.png" /></figure>
Abstract

In the pursuit of understanding the population of stellar remnants within the Milky Way, we analyze the sample of ~950 microlensing events observed by the Spitzer Space Telescope between 2014 and 2019. In this study we focus on a subsample of nine microlensing events, selected based on their long timescales, small microlensing parallaxes, and joint observations by the Gaia mission, to increase the probability that the chosen lenses are massive and the mass is measurable. Among the selected events we identify lensing black holes and neutron star candidates, with potential confirmation through forthcoming release of the Gaia time-series astrometry in 2026. Utilizing Bayesian analysis and Galactic models, along with the Gaia Data Release 3 proper-motion data, four good candidates for dark remnants were identified: OGLE-2016-BLG-0293, OGLE-2018-BLG-0483, OGLE-2018-BLG-0662, and OGLE-2015-BLG-0149, with lens masses of $3.0_{-1.3}^{+1.8} M_\odot$, $4.7_{-2.1}^{+3.2} M_\odot$, $3.15_{-0.64}^{+0.66} M_\odot$ and $1.40_{-0.55}^{+0.75} M_\odot$, respectively. Notably, the first two candidates are expected to exhibit astrometric microlensing signals detectable by Gaia, offering the prospect of validating the lens masses. The methodologies developed in this work will be applied to the full Spitzer microlensing sample, populating and analyzing the timescale ($t_E$) versus parallax ($\pi_E$) diagram to derive constraints on the population of lenses in general and massive remnants in particular.

*Unified Astronomy Thesaurus concepts: Astrometric microlensing effect (2140); Gravitational microlensing (672); Satellite microlensing parallax (2148); Microlensing parallax (2144); Black holes (162); Neutron stars (1108); Stellar remnants (1627); Gaia (2360); Stellar masses (1614)*

1. Introduction

Detecting and characterizing stellar remnants are vital for our understanding of the evolution of stars and populations in the Milky Way. Black holes, in particular, have been of great interest, due to their role in the growth and formation of galaxies. Neutron stars offer unique insights into stellar evolution, extreme matter, and fundamental physics. In addition, both can play a role in the distribution of dark matter in galaxies, which is still one of the biggest unresolved mysteries in astrophysics.

Detecting stellar remnants poses a challenging problem owing to their typically dim and elusive nature. Historically, efforts to identify these remnants have relied heavily on indirect methods that require the remnant to be in a binary system. For instance, the discovery of binary systems with compact objects, such as X-ray binaries, has provided crucial evidence of their existence. Notably, the detection of gravitational waves, pioneered by the Laser Interferometer Gravitational-Wave Observatory (LIGO) and the Virgo collaboration (e.g., B. P. Abbott et al. 2016), marked a monumental advancement in the direct detection of stellar remnants, particularly black holes and neutron stars.

Still, single, nonaccreting, stellar-mass black holes, as well as aged, isolated neutron stars, are practically inaccessible to date. Consequently, our knowledge about these celestial objects remains limited, although they hold valuable information related to stellar formation and evolution. The only practical way to observe them is through microlensing, i.e., by detecting their gravitational influence on the light from another object. Measuring the mass of the lensing object and constraining the flux that it emits allows us to assess whether it could be a stellar remnant. The only known isolated, stellar-mass black hole was identified using this technique (C. Y. Lam et al. 2022; P. Mróz et al. 2022; K. C. Sahu et al. 2022; C. Y. Lam & J. R. Lu 2023).

Microlensing is inherently limited in terms of mass measurement, due to its reliance on obtaining the Einstein radius $\theta_E$ from

$$M_L = \frac{\theta_E}{\kappa\pi_E}, \quad \pi_E = \frac{\pi_{rel}}{\theta_E}, \quad \pi_{rel} = \frac{1}{D_L} - \frac{1}{D_s} \quad (1)$$

where $\kappa = 8.144 \text{ mas}/M_\odot$, $\pi_E$ is the magnitude of the microlensing parallax vector (see, e.g., A. Gould 2004), $D_L$ is the distance to the lens, and $D_s$ is the distance to the source.

Direct measurements of $\theta_E$ are challenging to obtain in practice. In some cases it can be done through observations of caustic crossings, resolving the images using interferometry (F. Delplancke et al. 2001; N. Dalal & B. F. Lane 2003; N. J. Rattenbury & S. Mao 2006; A. Cassan & C. Ranc 2016) or astrometric microlensing (e.g., M. Dominik & K. C. Sahu 2000). All these approaches are somewhat suboptimal for studying large populations of objects. The first requires very special circumstances to occur: a caustic crossing of a binary lens or, less likely, the central passage of a single lens in front of the source. The second, resolving the images with high-precision interferometers, is a promising avenue. Two such detections have been reported to date: the first one (S. Dong et al. 2019; W. Zang et al. 2020) for the Kojima event (A. A. Nucita et al. 2018; A. Fukui et al. 2019) and the second one for Gaia19bld (E. Bachelet et al. 2022; A. Cassan et al. 2022; K. A. Rybicki et al. 2022). Although this technique is limited to brighter targets, we should expect an increasing number of interferometric observations of microlensing events thanks to developments in the field (Gravity+ Collaboration et al. 2022; A. Gould 2023). The third approach to direct $\theta_E$ measurement is through astrometric microlensing, which requires a sub-milliarsecond astrometric precision. This limitation will also slowly be overcome, especially thanks to the development of advanced space satellites like Gaia (e.g., K. A. Rybicki et al. 2018; J. Klüter et al. 2020) and Roman (e.g., S. Sajadian & K. C. Sahu 2023; J. Fardeen et al. 2024). It could also be possible to investigate seeing-limited data sets like OGLE or KMTNet (N. Segev et al. 2024, in preparation), where the number of measurements might help to overcome limited astrometric accuracy. However, up to now, success in observing astrometric microlensing has been limited. There are ongoing attempts to detect it using adaptive optics (e.g., J. R. Lu et al. 2016), but only limits were obtained from these kinds of studies. The signal has been detected only in a handful of events using the Hubble Space Telescope and under very special circumstances (K. C. Sahu et al. 2017; A. Zurlo et al. 2018; K. C. Sahu et al. 2022; P. McGill et al. 2023).<fnref n="34" /> As of today, it is still a technique that can only be applied to specific cases.

<footnote n="34">³⁴ Gaia preliminary astrometric time series also confirmed the light centroid deviations in the Gaia16aye event; see https://www.cosmos.esa.int/web/gaia/iow_20210924.</footnote>
An alternative method for estimating the mass of the lens in microlensing events involves using the timescale of the event and microlensing parallax, while making assumptions about the distribution of proper motions, because $\theta_E = \mu_{rel}t_E$. Then, one would calculate the mass as

$$M_L = \frac{\mu_{rel}t_E}{\kappa\pi_E} = 1.35M_{\odot} \left[ \frac{\mu_{rel}}{4 \text{ mas yr}^{-1}} \right] \left[ \frac{t_E}{100 \text{ days}} \right] \left[ \frac{0.1}{\pi_E} \right]. \quad (2)$$

However, this approach is reliant on assumptions about the proper-motion distribution and Galactic model, which can introduce uncertainties and biases in the final mass estimates. It is important to note that with this approach one does not directly measure the mass of the lens and only recovers statistical information about this parameter, depending on the assumed Galactic model. As mentioned before, it is also necessary to constrain the flux coming from the lens to be able to claim that it is a dark stellar remnant.

Several studies have already employed this mass estimation approach, where statistical information about $\mu_{rel}$ is applied. For example, Ł. Wyrzykowski et al. (2016) applied this method to the OGLE-III sample of 59 events exhibiting a parallax signal. Later, additional analysis was conducted for the same sample of events, but it also implemented the source proper motion from the Gaia Data Release 3 (GDR3) catalog, to tighten the prior on the relative proper-motion distribution (Ł. Wyrzykowski & I. Mandel 2020). In the following studies, P. Mróz & Ł. Wyrzykowski (2021) refined and improved the technique, while P. Mroz et al. (2021) applied it again to the OGLE-III data set.

The crucial element of this approach hinges on the detection or, at the very least, constraint of the microlensing parallax signal within the selected sample's events. Unfortunately, the ground-based-only photometric measurements usually do not provide strong constraints on the microlensing parallax, as it requires the Earth's accelerated motion around the Sun to be significant.

In addition, the more pronounced (and thus easier to measure) the microlensing parallax signatures, the higher the value of the $\pi_E$ parameter. This bias poses an even greater challenge in the examination of massive stellar remnants, as events with smaller microlensing parallaxes tend to favor more massive lenses (see Equation (1)). One way to avoid such bias and also measure smaller $\pi_E$ values would be to rely on simultaneous space satellite observations to identify the microlensing parallax signal, instead of Earth's orbital motion.

To utilize this idea, we reviewed the sample of ~950 events that have been a part of the Spitzer microlensing campaign, which was conducted in the 2014–2019 seasons. The Spitzer campaign was directed specifically toward the goal of extrasolar planet characterization. Nonetheless, the procedure of target selection and the observing strategy for the campaign are well-defined (J. C. Yee et al. 2015), meaning that it is a controlled sample, which allows drawing conclusions about the general stellar remnant population. In this study, we do not explore the whole population of stellar remnants based on the Spitzer sample, but rather select particular candidates with longer timescales and smaller values of the $\pi_E$ parameter, which are likely to be caused by a more massive lens. Then, after assessing the amount of light that is emitted by the blend and the lens, one can judge whether the lens is a good candidate for a dark stellar remnant or not. Furthermore, we select and analyze only those events that were

observed by the Gaia mission, as they might be verifiable in the near future, thanks to high-precision astrometric measurements from Gaia that could be used to measure or constrain $\theta_E$. It is worth mentioning here that the first candidate for stellar remnant lens from Spitzer has been found by S. Dong et al. (2007), in the first ever detection of the space parallax signal in a microlensing event. The best-fit solution suggested that the lens is likely a binary black hole system with a total mass of $10 M_{\odot}$.

This paper is organized as follows: In Section 2.4 we present an initial review and analysis of the whole Spitzer microlensing sample, which is performed to populate an initial $t_E-\pi_E$ diagram, necessary to select candidates for stellar remnant lenses. In Section 3 we perform an in-depth analysis of the nine events that were selected, including the derivation of lens mass and distance distributions, which utilizes our light-curve analysis and priors on lens-source proper motions based on the Milky Way models. In Section 4 we simulate realistic Gaia astrometry for a range of possible $\theta_E$ values, including the most probable ones resulting from the analysis presented in Section 3. We summarize and give conclusions in Section 5.

## 2. Populating the $t_E-\pi_E$ Plane

To find the subsample of microlensing events hosting stellar remnants, we first need to populate the $t_E-\pi_E$ diagram, which will allow us to select potentially massive lenses. First, as we only consider standard (point source, point lens; hereafter PSPL) events with parallax, we filter out all the events with clear caustic-crossing or finite-source signatures based on visual inspection. After this step, out of the ~950 Spitzer events, 720 remain. These events constitute our final sample that is used for the construction of the $t_E-\pi_E$ diagram.

### 2.1. Ground-based Data

During the modeling procedure (see Section 2.3), we are fitting the PSPL model with parallax to the joint set of OGLE +KMT+Spitzer data. In this step of unsupervised, automatic fitting we decide to omit MOA data for practical reasons, as OGLE+KMT sets are sufficient for creating an initial model. Later on, in a detailed analysis of individual events from the selected subsample, we use the full re-reduced data.

Triggering of targets in the Spitzer microlensing campaign was based on the OGLE EWS (A. Udalski et al. 1992, 2015a) and MOA (I. A. Bond et al. 2001) alerts. The KMTNet (S.-L. Kim et al. 2016) data were incorporated into the decision-making starting in 2016. The OGLE-IV data were collected with a large mosaic camera, consisting of 32 CCD chips, each with a resolution of 2048 × 4096 pixels, and a scale of 0.′′26 pixel$^{-1}$. The camera is mounted on the 1.3 m Warsaw Telescope, located in Las Campanas Observatory in Chile. The cadence for each event varied from ~4 hr$^{-1}$ up to ~0.5 day$^{-1}$, depending on the field, with the frequency of observations declining with the increasing (projected) distance from the Galactic center. The data were reduced with the OGLE-IV photometric pipeline (A. Udalski et al. 2015a), which implements an improved difference image analysis (DIA) procedure from P. R. Wozniak (2000). Full OGLE-IV light curves were used, with the data collected up until$^{35}$ HJD' ≈ 8920.$^{36}$

<footnote>$^{35}$ This marks the beginning of the 2020 bulge season and is a practical cutoff date for observations of 2019 (and earlier) events, given that OGLE paused its operations because of the COVID-19 outbreak at this time.</footnote>

<footnote>$^{36}$ HJD' ≡ HJD − 2450000.</footnote>
The KMTNet survey uses three 1.6 m telescopes, located at the Cerro Tololo Inter-American Observatory (CTIO, Chile), Siding Spring Observatory (SSO, Australia), and South African Astronomical Observatory (SAAO, South Africa). Each telescope has a mosaic camera with four 9 k × 9 k CCD chips mounted and a pixel scale of 0.4 pixel⁻¹. As the majority of events have good OGLE coverage, in this initial step of building an initial $t_E-\pi_E$ diagram it is sufficient to use publicly available KMTNet data from the automatic pipeline, which is only the part of the light curve from the discovery season for each event. In the detailed analysis of smaller subsamples we use all the available, re-reduced data.

## 2.2. Spitzer Data

The Spitzer Space Telescope microlensing campaign was conducted during the “warm” part of the mission, and so the data were collected using the IR, narrow bandwidth L, centered at the wavelength of 3.6 µm. While covering the peak of the event is the most beneficial for the microlensing parallax determination, the satellite observations usually only cover a part of the light curve, as there is a delay between the target selection and the actual observation time. Naturally, it is also not known a priori when exactly the event is going to peak from the Spitzer perspective.

Each year, targets could have been observed by Spitzer during either of the two, ~38-day-long windows (northern summer and winter), when the Galactic bulge was visible from the satellite’s location (S. Calchi Novati et al. 2015). In practice, most of the events were only observed within the summer window, as during the winter the Galactic center is not accessible for Earth-based instruments. There were a few exceptions, where the event was observed either in multiple years (mostly for the baseline information) or during the winter window. In fact, a few events from the selected “massive” subsample did get such additional winter observations, which significantly enhanced the microlensing parallax measurements (see Section 3). Nonetheless, in the full, 950-event sample of Spitzer microlensing events, most of the targets have only one patch of Spitzer data, taken within the summer window, with ~1-day cadence, over a period of ~2–5 weeks.

The rough description and numbers quoted above, while not detailing exact information for every specific event, provide the needed understanding of the Spitzer light curves’ coverage that is sufficient for the goals of our work. The specific details of the Spitzer microlensing campaign, its observing strategy, and its target selection criteria are beyond the scope of this paper. More information can be found in, e.g., J. C. Yee et al. (2015) and A. Udalski et al. (2015b). The detailed analysis of the whole Spitzer microlensing sample, which will be used for more general studies of the population of the Galactic planets and stellar remnants, will be presented in a separate paper.

## 2.3. Light-curve Modeling

To construct the initial $t_E-\pi_E$ plane and identify the events with potentially massive lenses, we fit a standard PSPL microlensing model with parallax to the 720 events that were not classified as “clearly binary” during the initial, visual inspection. In our light-curve analysis we use procedures from the `pyLIMA` package (E. Bachelet et al. 2017), which employs the `VBBinaryLensing` code for light-curve computation (V. Bozza et al. 2018).

First, we rescale error bars in all the data sets, which is necessary to obtain meaningful parameter uncertainties and to compare different models. We apply the rescaling procedure in steps to consecutive groups of data sets, fixing the error bars in the groups that have already been modified. We start with fitting to the OGLE data, as it has the most stable photometry and a long baseline. We first fit a PSPL model, which is then used as a seed for the two ($u_0 > 0$ and $u_0 < 0$) fits of PSPL + parallax models. Out of these two fits we pick the one with the better $\chi^2$ to be a reference. We then apply a standard formula for the new error bars (e.g., J. C. Yee et al. 2012):

$$ \sigma_{\text{new}} = \sqrt{(\gamma\sigma_{\text{old}})^2 + \epsilon^2}, \quad (3) $$

where we fix the value of the error floor $\epsilon = 0.005$ mag for all data sets but Spitzer, for which we do not set the floor for the error. We find the rescaling factor $\gamma$ by requiring $\chi^2/\text{dof} = 1$. This procedure is repeated multiple times, which iteratively removes outliers. After the whole procedure is finished for the OGLE data, we add all the other ground-based sets and, finally, the Spitzer data. In Table 1 we list the rescaling factors for all the events that were later included in the "massive" subsample.

Starting with a simple least-squares minimization, we fit a PSPL model to the rescaled ground-based data. We calculate the model using the standard formula for the magnification (e.g., B. Paczynski 1986):

$$ A(u) = \frac{u^2 + 2}{u\sqrt{u^2 + 4}}; \quad u(t) = \sqrt{u_0^2 + \left(\frac{t - t_0}{t_E}\right)^2}, \quad (4) $$

where $t_0$, $u_0$, and $t_E$ are the standard microlensing parameters: time of maximum, smallest projected separation in the units of Einstein radius, and Einstein timescale, respectively. The magnification enters the formula for total flux that is changing during the event:

$$ F_{\text{tot}}(t) = A(t)F_s + F_{\text{bl}}, \quad (5) $$

where $F_s$ is the flux from the source and $F_{\text{bl}}$ is the flux from the blend. During the modeling procedure, we include the blend flux through the blending parameter

$$ g = \frac{F_{\text{bl}}}{F_s}. \quad (6) $$

We use the derived Paczyński parameters ($t_0$, $u_0$, $t_E$) as a seed in the next fitting step. Then, before the final step, we fit a model with parallax to the ground-based data only, to assess the microlensing parallax signal resulting from Earth’s orbital motion. The model incorporates two additional parameters: northern and eastern components of the microlensing parallax vector $\pi_E = (\pi_{\text{EN}}, \pi_{\text{EE}})$ (see, e.g., A. Gould 2004 for details). Finally, the joint fit to all space- and ground-based data is performed, where both annual and space parallax effects have to be taken into account.

The space parallax also allows constraints on the ($\pi_{\text{EN}}, \pi_{\text{EE}}$) vector. Having a satellite at projected distance $D_\perp$ from Earth allows calculation of the microlensing parallax based on the difference between the $t_0$ and $u_0$ parameters measured from the ground and from space (e.g., S. Refsdal 1966; A. Gould 1994):

$$ \pi_E = \frac{au}{D_\perp}(\Delta\tau, \Delta u_0), \quad (7) $$
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Initial $t_E-\pi_E$ plane, populated with the results from the parallax PSPL fit to all Spitzer events (gray points). Colored crosses and circles mark all the events that were also observed by the Gaia mission and published either through GSA (blue) or in Ł. Wyrzykowski et al. (2023; green). Red diamonds mark those events that were selected for detailed analysis, as they host potentially massive lenses. For a given value of relative lens-source proper motion, one can draw a line of constant lens mass. Assuming a typical value of $\mu_{rel} = 4$ mas yr$^{-1}$, all the events lying below the gray line would be caused by lenses heavier than 1.4 $M_\odot$, while those below the black line are caused by lenses heavier than 2.5 $M_\odot$. To visualize the impact of proper motion, we also plot the dashed-dotted line for $\mu_{rel} = 2$ mas yr$^{-1}$ and $M_L = 2.5 M_\odot$. We note that there are two groups of events with faulty fits: one with high values of $t_E$ (~500 days), and the other with $\pi_E$ (~2); see Section 2.4 for details.</figcaption></figure>

where

$$ \Delta\tau = \frac{t_{0,sat} - t_0}{t_E}, \quad \Delta u_0 = u_{0,sat} - u_0, \qquad (8) $$

and subscript “sat” refers to parameters measured from the perspective of the satellite. The space-based measurement of the microlensing parallax is independent of the ground-based measurement, which allows for an additional cross-check between the two.

In both ground-based-only and full parallax fitting we explore $u_0 > 0$ and $u_0 < 0$ regimes separately, which usually results in two, similarly well-fitted models. While in principle we do expect up to four solutions (the so-called fourfold degeneracy; see, e.g., A. Gould 2004), at this stage of analysis we limit ourselves only to search for the $u_0+ \leftrightarrow u_0-$ degeneracy. Because during the fitting procedure we start from $\pi_E = (0, 0)$, we expect the solution to converge to the smaller values first, and thus we might not be finding solutions with large $\pi_E$ values. While for complete analysis of each event the full grid search has to be done, we accept this limitation in this work, as our goal is to find potentially massive lenses, i.e., those with smaller parallaxes. Then, for the selected nine candidates we explore the full $\pi_{EN}-\pi_{EE}$ space (see Section 3 for detailed analysis and selection criteria).

## 2.4. Initial $t_E-\pi_E$ Plane

Having PSPL models with parallax fitted to all the nonbinary light curves from the Spitzer sample of microlensing events, we were able to build an initial $t_E-\pi_E$ plane (Figure 1). All the events have two solutions, as for each event the $u_0-$ and

$u_0+$ planes were separately explored (see Section 2.3). Since in further analysis we are specifically interested in massive lens candidates, in Figure 1 we are plotting the solution that results in a higher expected mass of the lens (see Equation (1)). We mark the events that were observed by the Gaia mission (blue crosses and green circles). From among them we choose the ones that host potentially massive lens (red diamonds)—those are analyzed in detail as good candidates for dark remnants (see Section 3).

The majority of the events lie approximately in the middle of the presented space, but one might identify two distinct, smaller groups: one with large microlensing parallax values, clumped at $\pi_E \approx 2$, and one with long timescales of $t_E \approx 500$ days, which is the internal, upper limit on this parameter value in the fitting algorithm. Events in both groups (in total $\approx 40$) suffer from faulty fits. After visual inspection of the light curves and fitting results, we conclude that the reason is low signal-to-noise ratio of ground-based data, significant systematic trends (likely due to proper motion of the source/lens, which impacts the DIA reduction procedure), or a low number of Spitzer data points. All of these factors, and particularly their combination, lead to inaccurate microlensing parallax measurements, especially if the effect is small. As a result, the fitting algorithm converges to parameter values that are inaccurate or wrong, or does not converge at all. While most of these issues can be addressed by careful remodeling, detrending, or simply filtering out impacted light curves, such analysis is beyond the scope of this paper. There are no events observed by Gaia in these groups, and so they would not enter the final sample, although we still plot them here for the sake of completeness.
## 3. Massive Lens Candidates

In this section we describe the analysis of the subsample consisting of candidates for potentially massive lenses.

### 3.1. Sample Selection

There were three general criteria used during the selection procedure: long Einstein timescale $t_E$, low value of microlensing parallax $\pi_E$, and presence of Gaia data. The first two together make higher values of lens mass $M_L$ more likely, although we note that it also depends on the lens-source relative proper motion (Equation (1)). The last condition makes the mass of the lens (potentially) verifiable in the near future, thanks to the Gaia astrometry—see Section 4 for the analysis of this aspect. Construction of the initial $t_E-\pi_E$ diagram for all the PSPL Spitzer events allows easy selection of candidates according to the criteria described above. We refrain from applying strict cuts on the parameters, as there is no clear reason for such. Our subsample will not be representative for the stellar remnant population, and the selection procedure will be arbitrary anyway. We pick mostly events lying below the gray dashed line in Figure 1, i.e., those for which the more massive solution results in $M_L > 1.4 M_\odot$ (as a lower limit on the mass of neutron stars), under the assumption of $\mu_{rel} = 4$ mas yr$^{-1}$. We supplement this sample with event OGLE-2015-BLG-0145, which, although it lies above the gray line, has a long timescale of $t_E \gtrsim 100$ days and, with its relatively small $\pi_E$ value, might still be considered as a good candidate.

This results in a subsample of nine candidates (red diamonds on Figure 1) for stellar remnant lenses, for which we perform a more detailed analysis. With the additional source color analysis described below, we can measure the parallax signal more accurately, which in turn leads to more reliable predictions of the physical parameters of the lens, presented in Section 3.16.

### 3.2. Gaia Photometric Data

Currently there are two groups of microlensing events observed by Gaia that have photometric time series publicly available. The first group consists of all events published through Gaia Science Alerts (GSA; S. T. Hodgkin et al. 2013, 2021). Once an alert is announced through this channel, all the Gaia photometry for the target is published. The second group is the Gaia DR3 microlensing catalog (Ł. Wyrzykowski et al. 2023). Because it was constructed using DR3 data, it only contains measurements collected up to ~2017 May. In the last column of Table 3 we mark which events from our sample were alerted by GSA and thus have a full Gaia light curve available.

Most of the photometric observations from the Gaia satellite are taken in G band, and only these are used in this work during the construction of the light-curve model. The photometric measurements for the events detected through GSA do not have uncertainties reported, and so we estimate them based on the Gaia DR2 photometric content and validation paper (D. W. Evans et al. 2018).

On average, Gaia comes back to the same field every ~30 days, and every visit it takes two measurements separated by ~6 hr (each for one of the mirrors), although this number can vary depending on the gradient of the scanning angle (Gaia Collaboration et al. 2016). Although such frequency of

observations is insufficient to properly cover a microlensing light curve, Gaia photometry is only playing a supplementary role in the process of light-curve characterization. In particular, we do not detect any meaningful signal of space parallax in the Gaia data (which is expected, as Gaia is in orbit around Lagrange point L2, much closer to Earth than the Spitzer satellite). On the other hand, the astrometric data it will provide not only will be crucial for the light centroid shift detection (see Section 4) but also do not require as high cadence as photometry to be useful, due to the much longer effective timescale of the astrometric microlensing signal compared to its photometric counterpart (e.g., M. Dominik & K. C. Sahu 2000).

### 3.3. Additional Ground-based Data

In the analysis of the OGLE-2016-BLG-0293 event we also included data from the MOA collaboration (I. A. Bond et al. 2001). The photometric measurements were taken with the 1.8 m telescope located at Mt. John, New Zealand, using their standard broad *R* filter.

To determine the color of the source (see Section 3.5) for the OGLE-2015-BLG-0211 event, we used the *H*-band measurements from the ANDICAM instrument (D. L. DePoy et al. 2003), mounted on the 1.3 m SMARTS telescope located in the CTIO observatory in Chile.

### 3.4. Source Stars

Generally the information about the color of the source and the position of the source and the blend on the color–magnitude diagram (CMD) is useful to investigate potential inconsistencies in the model. It is also necessary to support the space parallax calculations based on the Spitzer data (see the next subsection). To perform a source color analysis, observations of the event in a second band are required. Having regular coverage in two bands that span a reasonably large range of magnifications throughout the event, we are able to derive the source color independently of the model. Linear regression is used to fit the relationship between the total fluxes in both bands, as the slope of the linear function is equivalent to the source flux ratio in the two bands, which directly translates to the source color (see Equation (5)). Whenever it is possible, the color analysis is made based on the OGLE data set, as it is the one with the longest baseline and most reliable photometry. Excluding OGLE-2015-BLG-0211, all the events in our sample have the OGLE V band available, and so the $(V-I)_{OGLE}$ color of the source is derived. Most of the sources reside in the red clump (see Figure 2), which suggests that they are red giants that belong to the Galactic bulge population. For the case of OGLE-2015-BLG-0211 we use the follow-up data in *H* band to construct the CMD and conclude that the source is also part of the red clump.

### 3.5. Spitzer Color Constraints

The Spitzer data often do not allow measurements of the *L*-band baseline flux, which makes the $\pi_E$ analysis less constrained. It is possible to enhance the microlensing parallax information inferred from Spitzer observations by finding the color of the source ($I-L)_{src}$ and using it as a constraint on the Spitzer source brightness *L* during the modeling process.<fnref n="37" /> To

<footnote n="37">³⁷ *L* denotes Spitzer’s 3.6 µm band.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 2.** CMDs for all the massive lens candidate events, constructed from the indicated catalog of stars (mostly OGLE, in all events but OGLE-2015-BLG-0211) in a 2' × 2' box around the magnified object. The large red circles mark positions of the source stars, while the cyan circle marks the color and brightness of the baseline (i.e., including blend). They overlap whenever the blending is low, which is the case for almost all events in the subsample, excluding OGLE-2016-BLG-0689.</figcaption></figure>

do that, we match stars detected in both OGLE and Spitzer frames and construct a color–color diagram, using *I*, *V*, and *L* bands (so-called VIL relation; see S. Calchi Novati et al. 2015,

for details). For this purpose we use stars from the red clump, and so we expect such a color–color relation to be linear for most of the sources. Using linear regression, we can find its
functional form and, by interpolating (or extrapolating) to the known values of the source color ($V - I$)$_{src}$, derive ($I - L$)$_{src}$.

Such a constraint is expected to have the highest impact on the light-curve fitting for the cases where Spitzer data only constrain the local slope of the light curve (which is often the case). We perform the modeling with and without the constraint and note significant improvement introduced by the color constraint for almost all events. The exception is the case of OGLE-2018-BLG-0662, where Spitzer covered part of the peak, already constraining the shape of the light curve sufficiently.

### 3.6. Light Curves and Modeling

To make sure there are no other solutions due to the fourfold degeneracy (A. Gould 2004), and to correctly probe the ($\pi_{EN}$, $\pi_{EE}$) space, we perform a dense grid search over these two parameters. We do the search separately for three fits: ground-based-only, Spitzer-"only" (A. Gould et al. 2020), and joint (using all data sets). The Spitzer-"only" approach is fitting the model to the ground-based data first and then, after fixing $t_0$, $u_0$, and $t_E$, fitting ($\pi_{EN}$, $\pi_{EE}$) and source/blend fluxes using only Spitzer data. Although the accuracy of this approach can be limited, as it does not involve simultaneous fitting of all the parameters, it is useful to gain better insight into the constraints on the parallax introduced by the ground- and space-based data. In theory, given the constraint on the source flux, each Spitzer measurement provides a circular constraint on the ($\pi_{EN}$, $\pi_{EE}$) plane (A. Gould 2019). In practice, the data are taken in different epochs and have a nonzero dispersion. As a consequence, the Spitzer-"only" parallax contours for real events form elongated arcs, which might be reduced further and provide constraints in both dimensions, for the cases where the Spitzer measurements cover larger parts of the light curve or are closer to the peak.

In Figures 3 and 4 we gathered relevant parts of the parallax contours for all three modeling approaches. Light curves for the most preferred solutions and a table with the fit parameters are presented in Figure 5 and Table 2, respectively. We present a more detailed discussion and comment on the analysis of each event separately in the following subsections.

### 3.7. OGLE-2015-BLG-0145

The $u_0$+ solution shows nonphysical negative blending level with $g_{OGLE} = -0.7 \pm 0.1$. The $u_0$- solution also yields negative blending, but consistent with zero at $g_{OGLE} = -0.06 \pm 0.12$. The $u_0$- model is also preferred in terms of the goodness of fit,<fnref n="38" /> as $\Delta\chi^2 \approx 18$. In addition, the Spitzer-"only" parallax measurement is slightly more compatible with the ground-based-only measurement for the $u_0$- case, although for the $u_0$+ solution the two fits also remain in reasonably good agreement.

### 3.8. OGLE-2015-BLG-0149

While OGLE covered the full event, KMT only observed this event starting in 2016, so there is only the tail of the declining part of the light curve available from this survey. The event lasts for ~100 days, and Spitzer measured only a part of the (rising) slope, so most of the parallax information comes from the ground-based light curve, even after including the VIL

<table><caption>**Table 1**
Rescaling Factors $\gamma$ for Each Telescope, Along with the Number of Photometric Data Points Used</caption><thead><tr><th>Event<sup>a</sup></th><th>Telescope</th><th>γ</th><th>N points</th></tr></thead><tbody><tr><td rowspan="3">OB150145</td><td>OGLE</td><td>1.4</td><td>554</td></tr><tr><td>Spitzer</td><td>1.6</td><td>48</td></tr><tr><td>Gaia</td><td>1.7</td><td>31</td></tr><tr><td rowspan="3">OB150149</td><td>OGLE</td><td>1.4</td><td>1653</td></tr><tr><td>Spitzer</td><td>4.5</td><td>50</td></tr><tr><td>Gaia</td><td>5.0</td><td>41</td></tr><tr><td rowspan="3">OB150211</td><td>OGLE</td><td>1.2</td><td>937</td></tr><tr><td>Spitzer</td><td>4.7</td><td>109</td></tr><tr><td>Gaia</td><td>1.5</td><td>28</td></tr><tr><td rowspan="3">OB160293</td><td>OGLE</td><td>1.6</td><td>2718</td></tr><tr><td>Spitzer</td><td>3.1</td><td>24</td></tr><tr><td>Gaia</td><td>1.9</td><td>32</td></tr><tr><td rowspan="6">OB160689</td><td>KMTS</td><td>1.0</td><td>1404</td></tr><tr><td>MOA</td><td>2.2</td><td>1798</td></tr><tr><td>OGLE</td><td>1.4</td><td>1162</td></tr><tr><td>Spitzer</td><td>2.7</td><td>16</td></tr><tr><td>Gaia</td><td>1.7</td><td>21</td></tr><tr><td>KMTA</td><td>1.1</td><td>442</td></tr><tr><td rowspan="3">OB180410</td><td>KMTC</td><td>1.0</td><td>726</td></tr><tr><td>KMTS</td><td>1.1</td><td>662</td></tr><tr><td>OGLE</td><td>1.7</td><td>2660</td></tr><tr><td rowspan="3">OB180483</td><td>Spitzer</td><td>1.9</td><td>29</td></tr><tr><td>Gaia</td><td>1.1</td><td>80</td></tr><tr><td>KMTA</td><td>3.1</td><td>190</td></tr><tr><td rowspan="3">OB180662</td><td>KMTC</td><td>1.6</td><td>369</td></tr><tr><td>KMTS</td><td>1.9</td><td>179</td></tr><tr><td>OGLE</td><td>1.4</td><td>1170</td></tr><tr><td rowspan="3">OB190169</td><td>Spitzer</td><td>1.9</td><td>27</td></tr><tr><td>Gaia</td><td>0.6</td><td>39</td></tr><tr><td>OGLE</td><td>2.3</td><td>952</td></tr><tr><td rowspan="3">OB190169</td><td>Spitzer</td><td>2.4</td><td>44</td></tr><tr><td>Gaia</td><td>0.5</td><td>59</td></tr><tr><td>KMTA</td><td>1.4</td><td>395</td></tr><tr><td rowspan="3">OB190169</td><td>KMTC</td><td>1.1</td><td>821</td></tr><tr><td>KMTS</td><td>1.3</td><td>359</td></tr><tr><td>OGLE</td><td>5.0</td><td>955</td></tr><tr><td rowspan="3">OB190169</td><td>Spitzer</td><td>2.5</td><td>33</td></tr><tr><td>Gaia</td><td>1.6</td><td>52</td></tr><tr><td>KMTS</td><td>2.1</td><td>3610</td></tr></tbody></table>

**Note.**

color constraint. As one can see from panels (c) and (d) of Figure 3, Spitzer-"only" and ground contours are consistent. Blending behaves well and is consistent with zero for both $u_0$+ and $u_0$- solutions. The $u_0$- is somewhat preferred according to the photometric models, with $\Delta\chi^2 \approx 19$.

### 3.9. OGLE-2015-BLG-0211

The event lies in a highly extincted field, and so there are no measurements in the $V_{OGLE}$ band that could be used to estimate the source color (from GDR3 we can have a very rough estimate of the $V - I$ color as $BP - RP = 4.52$ mag). Instead, we used additional observations in H band. The data points were distributed over a large part of the light curve, which allowed for the source color determination.

Spitzer data with a color constraint allow refinement of the microlensing parallax constraints obtained based on the ground-based data (panels (e) and (f) of Figure 3). Both

<footnote>a For simplicity, shortened versions of the event names are used in tables and figures throughout the paper: OBXXYYYY ≡ OGLE-20XX-BLG-YYYY.</footnote>

<footnote n="38">³⁸ The photometry error bars are renormalized, and so $\chi^2/\text{dof} \approx 1$ for all the events (see Table 2). Thus, when comparing different models, we use $\Delta\chi^2$ without quoting the number of dof throughout the text.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Microlensing parallax contours for $u_0+$ and $u_0^-$ solutions for the first four events in the selected subsample. Green and red outlines represent the $5\sigma$ level for ground-based and Spitzer-"only" (with color constraint) fits, respectively. The colored, filled contours changing from dark red to blue represent $1\sigma-5\sigma$ levels of the final, joint fit. For the events where additional ground-based data sets beyond OGLE were used, we also provide a ground parallax contour based only on the OGLE data (gray), to track potential systematics.</figcaption></figure>

$u_0+$ and $u_0^-$ solutions have similar $\chi^2 (\Delta\chi^2 \approx 2)$, and so one cannot assess which one is preferred based only on the photometric data.

### 3.10. OGLE-2016-BLG-0293

In the OGLE-I and MOA data, short- and long-timescale systematics are visible. That being said, the source is relatively bright, and so some level of baseline variability is expected. In the OGLE data it is comparable to the scatter and does not seem to have a large impact on the overall fit quality, but in the MOA data the long-timescale baseline systematics are much more prominent. We address that issue by checking for any color changes in the OGLE I-band and V-band data and potential correlations between the baseline color and magnitude changes, which would imply astrophysical origin. The additional analyses do not show such color changes or correlations.

In addition, trends visible in the MOA data are not compatible with the smaller trends in the OGLE data, and so we conclude that larger systematics visible in MOA should not be taken into account. Thus, we discard the part of the MOA data affected the most and only include the measurements collected at $HJD' > 7000$.

For the case of OGLE-2016-BLG-0293 we are able to discriminate between the two solutions, as there are multiple lines of evidence against the $u_0+$ one. First of all, there is a small difference $\Delta\chi^2 \approx 6$, which favors the negative solution. In addition, there is an evident disagreement between the ground-based-only and Spitzer-"only" parallax fits for the $u_0+$ case. Finally, the blending value obtained in the $u_0+$ solution is negative and not consistent with zero ($4\sigma$ away; see Table 2).

Including additional photometry beyond the OGLE data set causes increasing inconsistency between the Spitzer-"only"
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Same as Figure 3, but for the remaining five events. For the case of OGLE-2018-BLG-0483 the Spitzer data allowed only for estimation of the upper limit on the flux, which excludes part of the $\pi_{EN}-\pi_{EE}$ space (marked with red hatched area in panels (m) and (n)). The magenta contours in panels (k) and (l) represent ground-based-only fits with fixed blending parameter $g = 0$.</figcaption></figure>

parallax contours and the ground-based-only contours (see panel (h) of Figure 3). However, as we mentioned before, the Spitzer-"only" fit is only a diagnostic tool, and a general

agreement between the two contours is sufficient to claim that the two sources of the microlensing parallax measurement are consistent for the $u_0$-case.
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Light curves of all the events from our Spitzer–Gaia subsample, with all the data sets used during the modeling. Solid lines represent the best fit of the PSPL +parallax model, showing Earth/Gaia (black line) and Spitzer (red line) perspectives. The ground-based data from all the surveys are plotted together (blue circles). All models and data points are displayed with respect to *I*<sub>OGLE</sub> baseline brightness. For the case of OGLE-2018-BLG-0483 we only derive limits on the Spitzer flux, which are represented by downward-pointing triangles. The inset in the bottom left panel shows residuals of the region with a small deviation that might have been caused by a low-mass companion to the lens (see details in the text).</figcaption></figure>

### 3.11. OGLE-2016-BLG-0689

For this event, in addition to the “regular” Spitzer observations, the data were also taken in the December window, which resulted in covering both the rising and declining parts of the light curve.

Thanks to that, it was possible to determine the microlensing parallax relatively well.

While the blending for both of the solutions is consistent with zero, the posterior distribution for this parameter is also very wide, meaning that the blending parameter is not well
<table><caption>**Table 2**
The Final (Joint Ground+Space+Source Color Constraint) Light-curve Model Fit Results for All the Events</caption><thead><tr><th>Event</th><th>t<sub>0</sub> - 2,450,000<br>(days)</th><th>u<sub>0</sub></th><th>t<sub>E</sub><br>(days)</th><th>π<sub>EN</sub></th><th>π<sub>EE</sub></th><th>I<sub>0,OGLE</sub><br>(mag)</th><th>g<sub>OGLE</sub></th><th>χ<sup>2</sup>/dof</th></tr></thead><tbody><tr><td>OB150145</td><td>7297.8<sup>+1.2</sup><sub>-1.1</sub><br>7296.9<sup>+1.1</sup><sub>-1.0</sub></td><td>1.13<sup>+0.26</sup><sub>-0.17</sub><br>-0.622<sup>+0.048</sup><sub>-0.058</sub></td><td>76.6<sup>+8.0</sup><sub>-9.5</sub><br>116.9<sup>+5.2</sup><sub>-5.5</sub></td><td>0.108<sup>+0.014</sup><sub>-0.010</sub><br>-0.1064<sup>+0.0069</sup><sub>-0.0075</sub></td><td>0.111<sup>+0.020</sup><sub>-0.014</sub><br>0.104<sup>+0.010</sup><sub>-0.009</sub></td><td>15.683</td><td>-0.70<sup>+0.13</sup><sub>-0.12</sub><br>-0.06<sup>+0.13</sup><sub>-0.13</sub></td><td>668/622<br>650/622</td></tr><tr><td>OB150149</td><td>7277.8550<sup>+0.0097</sup><sub>-0.0095</sub><br>7277.8564<sup>+0.0095</sup><sub>-0.0094</sub></td><td>0.05349<sup>+0.00034</sup><sub>-0.00034</sub><br>-0.05382<sup>+0.00029</sup><sub>-0.00029</sub></td><td>101.82<sup>+0.44</sup><sub>-0.43</sub><br>100.77<sup>+0.44</sup><sub>-0.43</sub></td><td>0.0130<sup>+0.0079</sup><sub>-0.0083</sub><br>0.0512<sup>+0.0099</sup><sub>-0.0104</sub></td><td>-0.0446<sup>+0.0011</sup><sub>-0.0011</sub><br>-0.0458<sup>+0.0011</sup><sub>-0.0010</sub></td><td>15.922</td><td>0.0204<sup>+0.0062</sup><sub>-0.0060</sub><br>0.0140<sup>+0.0051</sup><sub>-0.0051</sub></td><td>1790/1733<br>1771/1733</td></tr><tr><td>OB150211</td><td>7209.85<sup>+0.37</sup><sub>-0.37</sub><br>7209.54<sup>+0.37</sup><sub>-0.36</sub></td><td>0.647<sup>+0.090</sup><sub>-0.070</sub><br>-0.636<sup>+0.069</sup><sub>-0.086</sub></td><td>95.8<sup>+7.6</sup><sub>-7.9</sub><br>96.0<sup>+7.9</sup><sub>-8.0</sub></td><td>0.035<sup>+0.011</sup><sub>-0.011</sub><br>-0.032<sup>+0.011</sup><sub>-0.012</sub></td><td>-0.046<sup>+0.011</sup><sub>-0.013</sub><br>-0.060<sup>+0.015</sup><sub>-0.018</sub></td><td>17.291</td><td>-0.24<sup>+0.16</sup><sub>-0.15</sub><br>-0.21<sup>+0.16</sup><sub>-0.16</sub></td><td>1076/1063<br>1078/1063</td></tr><tr><td>OB160293</td><td>7669.66<sup>+0.22</sup><sub>-0.21</sub><br>7669.56<sup>+0.22</sup><sub>-0.22</sub></td><td>0.576<sup>+0.011</sup><sub>-0.011</sub><br>-0.5352<sup>+0.0098</sup><sub>-0.0102</sub></td><td>107.3<sup>+1.2</sup><sub>-1.2</sub><br>114.3<sup>+1.4</sup><sub>-1.4</sub></td><td>0.0305<sup>+0.0027</sup><sub>-0.0027</sub><br>-0.0371<sup>+0.0030</sup><sub>-0.0030</sub></td><td>-0.0013<sup>+0.0014</sup><sub>-0.0014</sub><br>-0.0062<sup>+0.0014</sup><sub>-0.0014</sub></td><td>16.223</td><td>-0.113<sup>+0.028</sup><sub>-0.027</sub><br>-0.001<sup>+0.030</sup><sub>-0.030</sub></td><td>6096/5959<br>6090/5959</td></tr><tr><td>OB160689</td><td>7618.72<sup>+0.27</sup><sub>-0.27</sub><br>7618.54<sup>+0.28</sup><sub>-0.28</sub></td><td>1.08<sup>+0.19</sup><sub>-0.13</sub><br>-0.869<sup>+0.075</sup><sub>-0.095</sub></td><td>66.0<sup>+5.5</sup><sub>-6.5</sub><br>77.5<sup>+4.7</sup><sub>-4.9</sub></td><td>0.142<sup>+0.034</sup><sub>-0.033</sub><br>-0.131<sup>+0.033</sup><sub>-0.030</sub></td><td>-0.048<sup>+0.011</sup><sub>-0.012</sub><br>-0.063<sup>+0.014</sup><sub>-0.014</sub></td><td>16.954</td><td>-0.31<sup>+0.22</sup><sub>-0.21</sub><br>0.09<sup>+0.20</sup><sub>-0.20</sub></td><td>3059/3012<br>3063/3012</td></tr><tr><td>OB180410</td><td>8291.7666<sup>+0.0072</sup><sub>-0.0073</sub><br>8291.7566<sup>+0.0060</sup><sub>-0.0060</sub></td><td>0.11069<sup>+0.00069</sup><sub>-0.00069</sub><br>-0.11192<sup>+0.00073</sup><sub>-0.00073</sub></td><td>53.46<sup>+0.26</sup><sub>-0.26</sub><br>53.01<sup>+0.27</sup><sub>-0.27</sub></td><td>0.052<sup>+0.021</sup><sub>-0.022</sub><br>-0.136<sup>+0.022</sup><sub>-0.019</sub></td><td>0.1073<sup>+0.0099</sup><sub>-0.0109</sub><br>0.067<sup>+0.014</sup><sub>-0.014</sub></td><td>16.567</td><td>-0.0188<sup>+0.0067</sup><sub>-0.0067</sub><br>-0.0302<sup>+0.0069</sup><sub>-0.0067</sub></td><td>3599/3490<br>3580/3490</td></tr><tr><td>OB180483</td><td>8532.64<sup>+0.51</sup><sub>-0.50</sub><br>8533.34<sup>+0.56</sup><sub>-0.54</sub></td><td>0.1586<sup>+0.0076</sup><sub>-0.0071</sub><br>-0.1342<sup>+0.0052</sup><sub>-0.0056</sub></td><td>272.6<sup>+10.4</sup><sub>-9.9</sub><br>328.1<sup>+10.7</sup><sub>-10.4</sub></td><td>0.015<sup>+0.019</sup><sub>-0.020</sub><br>-0.064<sup>+0.029</sup><sub>-0.023</sub></td><td>0.1217<sup>+0.0068</sup><sub>-0.0075</sub><br>0.1033<sup>+0.0099</sup><sub>-0.0103</sub></td><td>16.567</td><td>0 (fixed)<br>0 (fixed)</td><td>3607/3490<br>3600/3490</td></tr><tr><td>OB180662</td><td>8303.484<sup>+0.021</sup><sub>-0.021</sub><br>8303.437<sup>+0.019</sup><sub>-0.018</sub></td><td>0.2624<sup>+0.044</sup><sub>-0.044</sub><br>-0.2612<sup>+0.0045</sup><sub>-0.0045</sub></td><td>64.99<sup>+0.86</sup><sub>-0.84</sub><br>65.20<sup>+0.88</sup><sub>-0.86</sub></td><td>0.0388<sup>+0.0068</sup><sub>-0.0065</sub><br>-0.0118<sup>+0.028</sup><sub>-0.0084</sub></td><td>0.0396<sup>+0.0028</sup><sub>-0.0027</sub><br>0.0431<sup>+0.0058</sup><sub>-0.0054</sub></td><td>17.600</td><td>-0.010<sup>+0.021</sup><sub>-0.021</sub><br>-0.006<sup>+0.022</sup><sub>-0.021</sub></td><td>2711/2613<br>2740/2613</td></tr><tr><td>OB190169</td><td>8657.902<sup>+0.073</sup><sub>-0.072</sub><br>8657.967<sup>+0.059</sup><sub>-0.058</sub></td><td>0.478<sup>+0.013</sup><sub>-0.012</sub><br>-0.478<sup>+0.016</sup><sub>-0.016</sub></td><td>79.3<sup>+1.5</sup><sub>-1.5</sub><br>78.9<sup>+2.0</sup><sub>-1.9</sub></td><td>0.0046<sup>+0.0062</sup><sub>-0.0061</sub><br>-0.085<sup>+0.012</sup><sub>-0.012</sub></td><td>-0.0447<sup>+0.0060</sup><sub>-0.0059</sub><br>-0.054<sup>+0.011</sup><sub>-0.011</sub></td><td>15.096</td><td>0.135<sup>+0.046</sup><sub>-0.045</sub><br>0.135<sup>+0.059</sup><sub>-0.057</sub></td><td>4799/4637<br>4768/4637</td></tr></tbody></table>

Note. Most of them have two ($u_0$ + and $u_0$−) solutions with comparable $\chi^2$, but for some it was possible to discriminate between them; see the text for more details. Blending parameter $g$ is the ratio of the blend flux $F_{bl}$ to the source flux $F_s$. The baseline brightness $I_0$ is recovered with 1–2 mmag precision in all cases. The photometry error bars are renormalized, and so $\chi^2/\text{dof} \approx 1$ for all the events. Thus, when comparing different models, we use $\Delta\chi^2$ without quoting the number of dof throughout the text.

constrained from the light curve. Indeed, the blending seems to be present, as the source color and baseline color are different (see Figure 2). Additionally, there are only two points in the OGLE V band on the magnified part of the light curve, which makes the source color determination more uncertain. Thus, we treat the Spitzer source color constraint with caution and conservatively use a wider prior on the Spitzer source color during the fitting ($3\sigma$ instead of $1\sigma$ resulting from the linear regression procedure described in Section 3.5). The $u_0$+ and $u_0$− solutions are similar in terms of the goodness of fit ($\Delta\chi^2 \approx 4$). The $u_0$+ has significant negative blending but is consistent with zero at $g_{OGLE} = -0.31 \pm 0.21$.

### 3.12. OGLE-2018-BLG-0410

After investigation of the photometric data, we noted that MOA data show large deviations that are not present in the remaining two surveys (especially at the wings of the event), and so we decide not to use it in the analysis. Spitzer data only cover part of the decline of the light curve and by themselves provide rather sparse limits on the microlensing parallax measurements. Nonetheless, performing Spitzer-"only" analysis and including the VIL color constraint shows that space-based parallax gives two relatively well-defined $\chi^2$ minima in $\pi_{EN}-\pi_{EE}$ space for each of the $u_0$+ and $u_0$− solutions. The $u_0$− is preferred with $\Delta\chi^2 \approx 19$ and is also much more convincing in terms of consistency between the Spitzer-"only" and ground-based-only parallax measurements, which can be seen from panels (k) and (l) of Figure 4. In addition, the blending is negative in both ground-based-only solutions and is not consistent with zero ($3\sigma-4\sigma$ away). To address this problem, we decide to redo the fitting procedure with blending

fixed to zero. Although the fit is clearly worse in terms of goodness (with the $\chi^2$ difference of 8 and 20; see Table 2), the tension between the ground-based-only and Spitzer-"only" parallax disappears, at least for the $u_0$− solution (see magenta contours in panels (k) and (l)).

### 3.13. OGLE-2018-BLG-0483

The event was observed by Spitzer in two seasons, separated by about a year, but Spitzer gives no significant constraints on the parallax. Due to the low signal-to-noise ratio and crowding of the sky region around the event, we were only able to set limits on the Spitzer flux (see the zoom-in in Figure 5). As a result, it was only possible to exclude part of the central region in the $\pi_{EN}-\pi_{EE}$ space (see Figure 4, panels (m) and (n)). Still, the microlensing parallax can be accurately determined from the ground, which is not surprising since the event is extremely long—the Einstein timescale is either ~275 or ~330 days, depending on the solution. The $u_0$+ solution has negative blending only marginally consistent with zero ($g_{OGLE} = -0.15 \pm 0.05$), which might suggest that the $u_0$− solution is the real one. More importantly, the $u_0$+ seems to be excluded by the Spitzer limits (Figure 4, panels (m) and (n)).

### 3.14. OGLE-2018-BLG-0662

Before fitting the final light-curve model, the OGLE data were corrected for a small linear trend ($\approx 0.04$ mag over 10 yr), likely caused by the change of the source position compared to the position on the reference image. The correction had a minor effect on the final results of the fit.
The microlensing parallax for the OGLE-2018-BLG-0662 event is weakly constrained from the ground, but with the addition of Spitzer data, a high-precision $\pi_E$ measurement was achieved, as the satellite covered part of the peak. In both the $u_0$+ and $u_0$− solutions, the blending is consistent with zero with small error bars, as in both cases $g = -0.01 \pm 0.02$, which makes this event a good candidate to host a dark lens (see Section 3.16). Comparison of the $\chi^2$ suggests that the $u_0$+ solution is preferred, as $\Delta\chi^2 \approx 28$. In addition, the Spitzer-"only" and ground-based-only fits are more compatible for that case (see panels (o) and (p) of Figure 4).

### 3.15. OGLE-2019-BLG-0169

The Spitzer measurements were only taken at the decline, but close to the peak, which helped to constrain the microlensing parallax from space. They played a significant role in the microlensing parallax determination—see panels (r) and (s) of Figure 4. There is a clear offset between the Spitzer-"only" and all ground-based parallax contours, while the OGLE-only solution is consistent with Spitzer-"only." There is a low-level (amplitude of ~20 mmag), irregular variability visible in the light curve, which might be the reason for systematic errors in ground-based microlensing parallax measurements. Nonetheless, the final result obtained from the joint fit is driven mostly by the Spitzer data and so is not strongly affected by the systematics in the ground-based data.

In the light curve of OGLE-2019-BLG-0169, at HJD' ≈ 8720, there is a clear deviation from the light-curve model lasting ~10 days (see inset in the bottom left panel of Figure 5). It is not an instrumental effect, as it appears in both OGLE and KMTNet data. While the scenario of the deviation being a planetary anomaly should not be completely excluded, the amplitude of the variability mentioned in the previous paragraph is comparable to the "anomaly." Thus, we conduct the analysis using the single lens model and attribute the feature to the variability of the source.

We note that one of the Gaia epochs was taken during the "anomaly." In an unlikely scenario of the deviation being due to a planet, it might be an interesting point of the analysis of the Gaia astrometric data. The planetary scenario will be investigated elsewhere, as it is beyond the scope of this paper.

### 3.16. Physical Parameters—Methodology

To assess the probability that the lens is dark, and to provide reasonable predictions of astrometric signal expected from the Gaia mission (see Section 4), we need to estimate the physical properties of the lens, namely its mass, distance, and brightness. All the events in the subsample analyzed here are standard events with microlensing parallax signal, which means that the light curve does not contain enough information to directly measure these properties. Nonetheless, the microlensing model provides some constraints that, coupled with the assumptions about the Milky Way kinematics and structure, can be used to evaluate Bayesian probabilities on the physical properties of the lens.

### 3.17. Mass and Distance

We adopt a similar approach to the one presented in K. Kruszynska et al. (2022) and K. Howil et al. (2024), which in turn is based on the procedure used by Ł. Wyrzykowski et al. (2016) and later refined by P. Mróz & Ł. Wyrzykowski (2021).

The technical details regarding the Milky Way model are summarized in P. Mróz & Ł. Wyrzykowski (2021) and in Appendix A of K. Howil et al. (2024). Below we provide a qualitative description of the analysis and comment on some aspects relevant for our use case.

The starting point of the procedure are posterior distributions of the light-curve parameters, obtained from the Markov Chain Monte Carlo (MCMC) modeling. To calculate the mass of the lens, the photometric model has to be supplemented with the source-lens relative proper motion $\mu_{rel} = |\mu_{rel}|$. Initially it is drawn from a wide, flat distribution [0, 30] mas yr$^{-1}$, and later it is weighted according to the Galactic model. Similarly, distance to the source $D_s$ is drawn from a flat distribution [0, 15] kpc. Then, for each link of the MCMC chains resulting from the light-curve modeling, we can calculate the mass of the lens (see Equation (2)) and also its distance as

$$D_L = \left( \theta_E \pi_E + \frac{1}{D_s} \right)^{-1}, \quad \theta_E = \mu_{rel} t_E. \qquad (9)$$

We then apply the "galactic prior" by weighing the resulting mass and distance using weights of the form (V. Batista et al. 2011)

$$w_{\text{Gal}} = \frac{4 D_L^4 \mu_{\text{rel}}^4 t_E}{\text{au} \pi_E} \nu_d f_{\mu} f_M M_L. \qquad (10)$$

The above expression combines three priors: the mass function $f_M$, relative proper motion prior $f_{\mu}$, and stellar density distribution $\nu_d$. Remaining quantities result from the transition between the physical parameters and the microlensing variables. The stellar density distribution $\nu_d$ consists of two separate expressions, with "double exponential" disk and barred bulge profiles (see V. Batista et al. 2011 and C. Han & A. Gould 2003). After applying the weights $w_{\text{Gal}}$ to lens mass and distance, we obtain the posterior distributions presented in Figure 7.

For the mass function we assume a power law $f_M \sim M^{\alpha}$, and for each event we address the impact of the assumed mass prior on the final distribution by comparing a "flat" prior $\alpha = -1$ with the Kroupa mass function (P. Kroupa 2001), where the slope is $\alpha = -2.35$ for the more massive ($M > 0.5 M_\odot$) tail of the distribution (see Table 4). While this is a simplification, as we disregard the different slopes for masses in the range $M_L < 0.5 M_\odot$, this approach is compatible with the selection process of our subsample, which favors more massive lenses. In addition, it is a conservative assumption in the context of claims of whether a lens is a dark remnant or not, because the slope of the power law at the massive end is much steeper.

The relative proper-motion prior $f_{\mu}$ is constructed based on the lens and source proper-motion assumptions. For the lens we consider two distinct cases: the lens lying in the Galactic disk and in the bulge. In the first scenario we assume a normal distribution of velocities in the galactic (l, b) coordinates: $V_l \approx N(220, 30)\text{km s}^{-1}$ and $V_b \approx N(0, 20)\text{km s}^{-1}$. For the bulge we assume $V_l = V_b \approx N(0, 100)\text{km s}^{-1}$ (C. Han & A. Gould 1995; V. Batista et al. 2011). These distributions are corrected for the motion of the Sun with respect to the local standard of rest (R. Schönrich et al. 2010). The mean of the disk velocity distribution also can vary, as it depends on the distance to the lens. A more detailed description is outlined in K. Howil et al. (2024). The physical velocities of a lens are then transformed to proper motions with
<table><caption>**Table 3**
Positions, Proper Motions, and RUWE (See Footnote 39) Parameters for All the Events in the Analyzed Sample, Taken from the GDR3 Catalog</caption><thead><tr><th>Event</th><th>R.A.<br>(J2016)</th><th>Decl.<br>(J2016)</th><th>μ<sub>R.A.*</sub><br>(mas yr<sup>-1</sup>)</th><th>μ<sub>Decl.</sub><br>(mas yr<sup>-1</sup>)</th><th>RUWE</th><th>GSA Name</th><th>GDR3 Source_id</th></tr></thead><tbody><tr><td>OB150145</td><td>270.17782</td><td>-35.15408</td><td>-1.29 ± 0.07</td><td>-6.51 ± 0.05</td><td>1.03</td><td>...</td><td>4041998223399082752</td></tr><tr><td>OB150149</td><td>270.28819</td><td>-32.55773</td><td>-2.01 ± 0.21</td><td>-5.03 ± 0.13</td><td>2.35</td><td>...</td><td>4042928139682133120</td></tr><tr><td>OB150211</td><td>262.35909</td><td>-30.98178</td><td>-3.50 ± 0.29</td><td>-7.22 ± 0.21</td><td>1.05</td><td>...</td><td>4058004814930630912</td></tr><tr><td>OB160293</td><td>268.16140</td><td>-32.48960</td><td>-1.49 ± 0.09</td><td>-7.48 ± 0.07</td><td>0.92</td><td>...</td><td>4043504794840743040</td></tr><tr><td>OB160689</td><td>261.14944</td><td>-30.13245</td><td>-1.80 ± 0.36</td><td>-7.25 ± 0.23</td><td>1.38</td><td>...</td><td>4059051309459806208</td></tr><tr><td>OB180410</td><td>272.22770</td><td>-27.16974</td><td>-4.51 ± 0.16</td><td>-6.84 ± 0.11</td><td>1.14</td><td>Gaia18cho</td><td>4063011505688647296</td></tr><tr><td>OB180483</td><td>262.64177</td><td>-27.49183</td><td>4.85 ± 1.23</td><td>-2.13 ± 0.74</td><td>1.34</td><td>Gaia18ayh</td><td>4061439448723558016</td></tr><tr><td>OB180662</td><td>266.87754</td><td>-32.52442</td><td>-5.74 ± 0.59</td><td>-6.09 ± 0.34</td><td>1.13</td><td>Gaia18cej</td><td>4054012488194100096</td></tr><tr><td>OB190169</td><td>265.98559</td><td>-32.87095</td><td>-0.87 ± 0.11</td><td>-4.54 ± 0.06</td><td>0.97</td><td>Gaia19drv</td><td>4054032245075925760</td></tr></tbody></table>

Note. In the last two columns we also provide the alert name, for those detected through GSA, and the GDR3 source identifier.

$\mu_L = 4.74V_L[\text{km s}^{-1}]/D_L[\text{kpc}]$ and, after subtraction of the source proper motion (see next paragraph), can be used as a final prior $f_\mu$.

The relative proper motion partially depends on the source proper motion, which is available from the Gaia DR3 catalog (Gaia Collaboration et al. 2023), and we incorporate this information in the prior. Values of GDR3 proper motions for each event are presented in Table 3. We also quote the RUWE<fnref n="39" /> parameter to quantify the credibility of the Gaia measurements. It is important to note that we can use quoted proper motions only for the events where all the light (or at least most of it) actually comes from the source. Otherwise, the proper motion measured by Gaia is a combination of that of the source and blend(s). The majority of events in our subsample do not exhibit large blending, and so we decide to include Gaia proper motions in the prior for eight out of nine events. In the remaining case of OGLE-2016-BLG-0689, where the blending appears to be more significant, we use a different approach; see Section 3.19.5 for more details.

It is also worth noting that the direction of the relative lens —source proper motion can be constrained from the light-curve model thanks to the microlensing parallax measurements:

$$ \hat{\mu}_{\text{rel}} = \frac{\mu_{\text{rel}}}{\mu_{\text{rel}}} = \frac{\pi_{\text{E}}}{\pi_{\text{E}}}, \qquad (11) $$

which has to be taken into account.

In principle, we could assume fixed distance to the source in our calculations—all the events analyzed here lie toward the Galactic center, and so it is expected that the sources belong to the bulge population, particularly when their location on the CMD coincides with the red clump (see Figure 2). Nonetheless, we decide to weight them with the stellar density distribution, which gives more realistic results. We also decide not to use parallax measurements from GDR3, or the distance estimates based on them (C. A. L. Bailer-Jones et al. 2021). The reason is that the sources mostly lie in the bulge, and so Gaia parallax is not measured accurately enough (with the signal of the order of ~0.1 mas), especially in such crowded fields. Indeed, after inspecting the parallax values quoted in GDR3, we found that the measured parallax is either negative or parallax_over_error ≲ 1, meaning that these measurements do not carry useful information.

It is worth noting that in principle mass and velocity distributions for stellar remnants are different from those of stars. Nonetheless, in this experiment we do not know a priori if the lenses belong to the stellar remnant population, and so we cannot assume that from the beginning and use priors for black holes/neutron stars. Instead, we use a Galactic model and the mass function based on stars, and after assessing the amount of light to the lens, we examine the scenario of it being a “regular” star. We note that if the velocity of a lens is higher, as expected for neutron stars (e.g., G. Hobbs et al. 2005) and some black holes (e.g., S. Repetto et al. 2017), then the derived mass can be regarded as a lower limit.<fnref n="40" />

### 3.18. Lens Light

To claim that a lens is a stellar remnant, not only the mass but also the blend flux has to be investigated. After constructing the posterior distributions for the lens mass and distance, one can estimate the brightness $I_{\text{lens}}$ that the object of these properties should emit, under the assumption that it is a main-sequence (MS) star. For this purpose we use the empirical mass–luminosity relation found by M. J. Pecaut & E. E. Mamajek (2013).<fnref n="41" /> If we denote this relation by $L(M_L)$, we can define $I_{\text{lens}}$ as

$$ I_{\text{lens}} = -2.5 \log \frac{L(M_L)}{4\pi D_L^2} + A_{\text{lens}}, \qquad (12) $$

where $A_{\text{lens}}$ is extinction to the lens.

To estimate $A_{\text{lens}}$, we use extinction maps from D. M. Nataf et al. (2013) for the four events lying in the OGLE-III fields covered in their analysis. For the remaining five events we derive the extinction value based on the red clump position on the CMD, using a similar procedure to that of Nataf et al., using their dereddened red clump brightness. The $A_{\text{lens}}$ parameter is the extinction integrated along the whole distance to the source, so we treat it as an upper limit on the lens extinction. Again, for the sake of determining whether the lens is a stellar remnant or not, this is a conservative assumption. In addition, this simplification is justified by the fact that most of the disk dust between the observer and the source resides within the first few kiloparsecs from Earth, especially for the events with larger galactic latitudes.

<footnote n="39">³⁹ Renormalized Unit Weight Error; see, e.g., Gaia Collaboration et al. (2021). Generally, for a well behaved model RUWE ≈ 1. As a rule of thumb, RUWE ≳ 1.4 suggests that the astrometric model might not be reliable.</footnote>

<footnote n="40">⁴⁰ It is important to note that in microlensing we access only transverse velocity—the natal-kick velocities of such objects can in principle have any direction, and so, due to projection, they will not necessarily be seen as high-velocity lenses.</footnote>

<footnote n="41">⁴¹ http://www.pas.rochester.edu/~emamajek</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Visualization of the $p_{\text{dark}}$ calculation procedure for the example event OGLE-2016-BLG-0293. Top panel: posterior distribution of the blend flux (in arbitrary units) resulting from the light-curve modeling. The red line shows the blend flux limit we assume for the $p_{\text{dark}}$ calculations. Bottom panel: posterior distribution $p(I_{\text{lens}})$ of the lens brightness $I_{\text{lens}}$ (see details in the text). The red line shows the upper limit on the blend brightness derived from the top histogram. All the samples for which the lens is brighter than the limit (to the left from the red line) contribute to the dark lens scenario.</figcaption></figure>

Having a distribution of $I_{\text{lens}}$, one can compare it with the brightness of the blend, which is one of the products of the light-curve modeling. In the case where the total blend brightness is higher than $I_{\text{lens}}$, we get more light from the blend than is expected from a lens of given mass at a given distance. Such a scenario is very common and easy to explain, as the excess can be attributed to any other sources lying on the same line of sight and not participating in microlensing—the lens is not the only light source contributing to the blending light, which is common for the Galactic bulge direction. On the other hand, if the blend brightness calculated from the photometric model is lower than $I_{\text{lens}}$, the situation is opposite—there is not enough light emitted by the blend to explain an MS star of given mass at a given distance. In other words, such an MS star would be too bright compared to the expected blend light.<fnref n="42" /> This suggests that the lensing object is in fact not luminous. Then, if it is massive enough, it is considered a candidate for a stellar remnant.

To perform the comparison between $I_{\text{lens}}$ and blend brightness, the latter has to be estimated based on the blend flux distribution $F_{\text{bl}}$, which is the product of the light-curve modeling. To be conservative, we choose a $3\sigma$ upper limit on the blend flux (red dashed line in the top panel of Figure 6). To compare it with $I_{\text{lens}}$, we translate this limit to magnitudes and call it $I_{\text{blend}}$ (red dashed line in the bottom panel of Figure 6).

As mentioned above, the lens is expected to be dark for each set of parameters (each link of MCMC chains) resulting in $I_{\text{lens}} < I_{\text{blend}}$. To formally assess the probability $p_{\text{dark}}$ that the lens is dark, we integrate all the samples for which $I_{\text{lens}} < I_{\text{blend}}$:

$$P_{\text{dark}} = \frac{\int_{0}^{I_{\text{blend}}} p(I_{\text{lens}}) dI_{\text{lens}}}{\int_{-\infty}^{\infty} p(I_{\text{lens}}) dI_{\text{lens}}}, \quad (13)$$

where $p(I_{\text{lens}})$ is the posterior distribution of the $I_{\text{lens}}$ brightness. We illustrate this procedure in Figure 6: all the samples for which lens brightness can be explained by the blend are grayed out in the bottom panel. For the remaining part of the distribution, the lens with given mass at a given distance would be too bright to explain it with an MS star, and thus this region corresponds to the dark lens scenario.

### 3.19. Physical Parameters—Results

Below we provide results of the Bayesian analysis for each event separately.

#### 3.19.1. OGLE-2015-BLG-0145

In Figure 7 we present the Bayesian analysis results for both the $u_0-$ and $u_0+$ solutions, although the $u_0+$ one is virtually excluded by the light-curve analysis (see Section 3.7). The $M_L$ versus $D_L$ distribution shows a bimodality, as the lens is either in the disk or in the bulge. This overall structure is visible in most of the events in the subsample, although usually a bulge lens is the preferred solution, as expected. Indeed, for the $u_0-$ solution here, the bulge scenario is much more preferred and yields a low lens mass $M_L = 0.30_{-0.12}^{+0.17}$ at $D_L = 6.12_{-0.75}^{+0.87}$, with a very low chance of being dark (see Table 4).

#### 3.19.2. OGLE-2015-BLG-0149

In the $M_L-D_L$ plane, the two solutions behave in a similar way, with only one maximum each. The preferred $u_0-$ solution places the lens at around 6.3 kpc with mass in the range 0.7–2 $M_{\odot}$. The 79% probability of being dark suggests that the lens might be a massive white dwarf or a neutron star. This is the only event in the sample that has an RUWE parameter value substantially higher than 1 (see Table 3), and so we decided to perform the Bayesian analysis with the source proper-motion prior (taken from GDR3 and thus, given high RUWE, unreliable) loosened. Nominally we use the $1\sigma$ value from GDR3 as the width of the prior. Here we are broadening the prior to 10 standard deviations, as the Gaia values should not be fully trusted. We do not observe significant change in the $M_L-D_L$ plane, apart from the anticipated broadening of the posterior distributions.

#### 3.19.3. OGLE-2015-BLG-0211

The $M_L-D_L$ plane is again similar for the two solutions, with only the bulge lens scenario being viable. In both distributions one can see an additional, sharper structure at around 8 kpc.

<footnote n="42">⁴² We treat MS stars as a reference in this analysis, but more evolved stars of a given mass would be even brighter.</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Posterior distributions for lens mass and distance, under the assumption of the Kroupa mass function (slope $\alpha = -2.35$)—see text for the details. The orange, yellow, cyan, light-blue, and dark-blue contours mark $1\sigma$, $2\sigma$, $3\sigma$, $4\sigma$, and $5\sigma$ confidence levels, respectively. Note that we show the distributions for all the possible solutions (both $u_0^+$ and $u_0^-$ models for each event), even though some of them can be rejected based on the light-curve analysis. We list the results for preferred solutions in Table 4.</figcaption></figure>

Loosening the Gaia prior makes it merge with the wider bulge distribution, which means that it is caused by the source proper-motion prior. Similar structure is visible in some of the other events in the subsample.

### 3.19.4. OGLE-2016-BLG-0293

The lens in this event is one of the best candidates for a massive remnant in the analyzed subsample. Multiple factors in the light-curve modeling suggest that the $u_0^-$ model is correct (see Section 3.6), and so we present values yielded by this model in Table 4.

Although $\theta_E$ does not have an extreme value, thanks to a small microlensing parallax, the resulting lens mass is somewhat large and is expected to lie in the range 1.7–4.9 $M_\odot$ for bulge lenses. This results in a 94% probability for the lens being dark. Hence, it is an excellent candidate for a neutron star or a stellar-mass black hole. There is a possibility that the lens lies in the disk, which yields even higher masses, but it is

strongly disfavored in our analysis (see relevant panel in Figure 7).

### 3.19.5. OGLE-2016-BLG-0689

The blending for this event is not well constrained from the photometric model, but it is likely nonnegligible (see Section 3.11). As a result, the proper motion detected by Gaia is a combination of source and lens proper motion, which means that we cannot use it directly in the prior as source proper motion. Instead, when calculating the prior on relative proper motion, we assume the source proper motion in galactic coordinates $(l, b)$ to be $(\mu_l, \mu_b) = (-6.12, -0.19) \pm 2.64 \text{ mas yr}^{-1}$ (R. Schönrich et al. 2010). This corresponds to a typical motion of the Galactic center relative to the Sun and is a reasonable assumption, as the source most likely resides in the bulge. The same procedure was applied in P. Mroz et al. (2021). Naturally, the resulting $M_L$ versus $D_L$ distribution is much wider compared to the one resulting from the approach with the Gaia value. Because of the low expected lens
<table><caption>**Table 4**

Results of Bayesian Analysis in the Form of Median Values of Marginalized Distributions for Lens Masses, Distances, and Einstein Radii, Along with the Probabilities That the Lens Is Dark</caption><thead><tr><th>Event</th><th>Model</th><th>Mass Prior</th><th>Lens Mass M<sub>L</sub><br>(M<sub>☉</sub>)</th><th>M<sub>L,2σ min</sub><br>(M<sub>☉</sub>)</th><th>Lens Distance D<sub>L</sub><br>(kpc)</th><th>θ<sub>E</sub><br>(mas)</th><th>p<sub>dark</sub></th></tr></thead><tbody><tr><td rowspan="2">OB150145</td><td>u<sub>0</sub>-</td><td>Flat</td><td>0.83<sup>+3.49</sup><sub>-0.53</sub></td><td>0.17</td><td>2.1<sup>+4.2</sup><sub>-1.4</sub></td><td>1.00<sup>+4.26</sup><sub>-0.63</sub></td><td>49.9%</td></tr><tr><td>u<sub>0</sub>-</td><td>Kroupa MF</td><td>0.30<sup>+0.17</sup><sub>-0.12</sub></td><td>0.09</td><td>6.15<sup>+0.83</sup><sub>-0.77</sub></td><td>0.36<sup>+0.19</sup><sub>-0.14</sub></td><td>5.8%</td></tr><tr><td rowspan="2">OB150149</td><td>u<sub>0</sub>-</td><td>Flat</td><td>1.80<sup>+0.83</sup><sub>-0.64</sub></td><td>0.70</td><td>6.00<sup>+0.93</sup><sub>-1.14</sub></td><td>0.97<sup>+0.42</sup><sub>-0.32</sub></td><td>92.5%</td></tr><tr><td>u<sub>0</sub>-</td><td>Kroupa MF</td><td>1.40<sup>+0.75</sup><sub>-0.56</sub></td><td>0.45</td><td>6.35<sup>+0.95</sup><sub>-0.95</sub></td><td>0.77<sup>+0.38</sup><sub>-0.29</sub></td><td>79.0%</td></tr><tr><td rowspan="2">OB150211</td><td>u<sub>0</sub>+</td><td>Flat</td><td>2.7<sup>+3.4</sup><sub>-1.7</sub></td><td>0.56</td><td>6.8<sup>+1.5</sup><sub>-1.1</sub></td><td>1.05<sup>+0.72</sup><sub>-0.61</sub></td><td>28.4%</td></tr><tr><td>u<sub>0</sub>+</td><td>Kroupa MF</td><td>1.04<sup>+1.41</sup><sub>-0.49</sub></td><td>0.30</td><td>8.20<sup>+0.17</sup><sub>-1.75</sub></td><td>0.43<sup>+0.57</sup><sub>-0.16</sub></td><td>4.7%</td></tr><tr><td rowspan="2">OB160293</td><td>u<sub>0</sub>-</td><td>Flat</td><td>4.0<sup>+1.3</sup><sub>-1.5</sub></td><td>1.46</td><td>6.23<sup>+0.86</sup><sub>-0.70</sub></td><td>1.21<sup>+0.54</sup><sub>-0.44</sub></td><td>98.9%</td></tr><tr><td>u<sub>0</sub>-</td><td>Kroupa MF</td><td>3.0<sup>+1.8</sup><sub>-1.3</sub></td><td>0.85</td><td>6.58<sup>+0.91</sup><sub>-0.80</sub></td><td>0.91<sup>+0.52</sup><sub>-0.39</sub></td><td>93.6%</td></tr><tr><td rowspan="2">OB160689</td><td>u<sub>0</sub>+</td><td>Flat</td><td>1.10<sup>+0.79</sup><sub>-0.48</sub></td><td>0.31</td><td>4.0<sup>+1.6</sup><sub>-1.0</sub></td><td>1.11<sup>+0.50</sup><sub>-0.42</sub></td><td>42.2%</td></tr><tr><td>u<sub>0</sub>+</td><td>Kroupa MF</td><td>0.68<sup>+0.54</sup><sub>-0.36</sub></td><td>0.13</td><td>4.5<sup>+2.1</sup><sub>-1.2</sub></td><td>0.79<sup>+0.48</sup><sub>-0.40</sub></td><td>16.2%</td></tr><tr><td rowspan="2">OB180410</td><td>u<sub>0</sub>-</td><td>Flat</td><td>0.58<sup>+0.24</sup><sub>-0.21</sub></td><td>0.16</td><td>5.62<sup>+0.89</sup><sub>-2.53</sub></td><td>0.59<sup>+0.45</sup><sub>-0.22</sub></td><td>...</td></tr><tr><td>u<sub>0</sub>-</td><td>Kroupa MF</td><td>0.42<sup>+0.24</sup><sub>-0.26</sub></td><td>0.08</td><td>6.09<sup>+2.23</sup><sub>-0.91</sub></td><td>0.42<sup>+0.30</sup><sub>-0.25</sub></td><td>...</td></tr><tr><td rowspan="2">OB180483</td><td>u<sub>0</sub>-</td><td>Flat</td><td>6.8<sup>+5.6</sup><sub>-2.6</sub></td><td>2.55</td><td>1.55<sup>+0.71</sup><sub>-0.62</sub></td><td>5.4<sup>+2.3</sup><sub>-2.0</sub></td><td>99.6%</td></tr><tr><td>u<sub>0</sub>-</td><td>Kroupa MF</td><td>4.7<sup>+3.1</sup><sub>-2.1</sub></td><td>0.47</td><td>2.07<sup>+1.03</sup><sub>-0.69</sub></td><td>3.7<sup>+2.4</sup><sub>-1.6</sub></td><td>92.3%</td></tr><tr><td rowspan="2">OB180662</td><td>u<sub>0</sub>+</td><td>Flat</td><td>3.33<sup>+0.65</sup><sub>-0.60</sub></td><td>2.02</td><td>5.28<sup>+0.80</sup><sub>-0.66</sub></td><td>1.48<sup>+0.25</sup><sub>-0.23</sub></td><td>99.9%</td></tr><tr><td>u<sub>0</sub>+</td><td>Kroupa MF</td><td>3.15<sup>+0.66</sup><sub>-0.64</sub></td><td>1.54</td><td>5.35<sup>+0.95</sup><sub>-0.67</sub></td><td>1.42<sup>+0.25</sup><sub>-0.26</sub></td><td>99.3%</td></tr><tr><td rowspan="2">OB190169</td><td>u<sub>0</sub>-</td><td>Flat</td><td>0.71<sup>+0.41</sup><sub>-0.28</sub></td><td>0.25</td><td>6.63<sup>+1.26</sup><sub>-0.91</sub></td><td>0.54<sup>+0.25</sup><sub>-0.20</sub></td><td>0.7%</td></tr><tr><td>u<sub>0</sub>-</td><td>Kroupa MF</td><td>0.50<sup>+0.32</sup><sub>-0.22</sub></td><td>0.15</td><td>7.1<sup>+1.0</sup><sub>-1.0</sub></td><td>0.40<sup>+0.22</sup><sub>-0.17</sub></td><td>0.1%</td></tr></tbody></table>

Note. The error values reflect the 68% confidence level intervals. We also show the 2σ lower limit on the lens mass $M_{L,2\sigma \min}$. For each event we show the preferred solution (either $u_0+$ or $u_0^-$) and results for the two priors on the mass function (see text for more information).

mass and presumably strong blending, it is very unlikely that this event hosts a remnant lens.

### 3.19.6. OGLE-2018-BLG-0410

The $u_0^-$ solution is preferred in terms of goodness of fit in the light-curve analysis ($\Delta\chi^2 \approx 19$). Redoing the photometric fit with the blending fixed to zero helps to resolve part of the tension between the microlensing parallax solutions (see Section 3.12), and so we use this model to estimate the physical parameters. With the blending parameter fixed to zero, we cannot assess the probability of the lens to be dark, as it is assumed to be dark in the first place. Nonetheless, the mass of the lens yielded by our analysis is very low, and so it is most likely an ordinary dwarf star lying in the bulge.

### 3.19.7. OGLE-2018-BLG-0483

There is a clear bimodality in the $M_L$ versus $D_L$ distribution, due to the duality in the possible lens populations. Each of the $u_0^+$ and $u_0^-$ solutions has the disk lens solution, in which the lens is relatively heavy, in the range 5–10 $M_☉$, and the bulge lens solution where $M_L \lesssim 1 M_☉$. The $u_0^+$ case is excluded with the flux Spitzer limits, so we know that $u_0^-$ is the correct one. From the posterior distribution of the lens mass and distance for the negative solution, it seems that the scenario of heavier lens located in the disk is more preferred, although there is a small region of similar probability for the bulge lens case. As the blending is very low for this event, the lens is a good remnant candidate with estimated mass $M_L = 4.7_{-2.1}^{+3.1}$ and $p_{\text{dark}} \approx 92\%$.

### 3.19.8. OGLE-2018-BLG-0662

The positive solution is favored in terms of both goodness of fit and compatibility of ground-based-only and Spitzer-"only" parallax solutions (see Section 3.14). The expected mass for

this event remains somewhat large for a range of lens distances (see Figure 7). Additionally, the blending level is very low, which makes this event one of the best candidates in our subsample to host a dark remnant with $M_L = 3.15 \pm 0.65$ and $p_{\text{dark}} \approx 99\%$.

### 3.19.9. OGLE-2019-BLG-0169

In the $M_L-D_L$ plane, in both of the solutions we can see the impact of the GDR3 source proper-motion prior, similarly to OGLE-2015-BLG-0211—the narrow structure gradually disappears with the increase of the width of the prior. None of the solutions show the prospect for the remnant lens, and the preferred $u_0^-$ case suggests particularly low probability for that with mass of the lens $M_L \approx 0.5 M_☉$ and $p_{\text{dark}} < 1\%$.

## 4. Gaia Predictions

All the events analyzed here were chosen under the condition that they had been observed by the Gaia mission and so there will be high-precision astrometry available for them in the Gaia Data Release 4 (GDR4, to be published ~2026). Knowing the epochs of astrometric measurements and the expected uncertainties and assuming the $\theta_E$ values derived in Section 3.16, we perform realistic simulations of the Gaia astrometry and investigate the prospects for measuring $\theta_E$ and consequently the masses of the lenses in the studied events.

The Gaia mission will provide 2D astrometry only for the brightest objects ($G > 13$ mag). For the remaining ones, it will provide only 1D astrometry, measured along the instantaneous scanning direction $\tilde{AL}$, which we will denote with angle $\psi$ measured from the north direction eastward. In the simulations, we use actual values of scanning angles and epochs calculated based on the Gaia scanning law, taken from GOST$^{43}$ (Gaia

<footnote>$^{43}$ https://gaia.esac.esa.int/gost/</footnote>
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Visualization of all the known information about the Gaia measurements, under the assumption of the photometric microlensing model and $\theta_E = 0.91$ mas (see Table 4), for the case of OGLE-2016-BLG-0293. Top panel: photometric model with marked Gaia measurements. Middle panel: absolute astrometric shift calculated based on the photometric model (gray line). The black arrows indicate how much the astrometric signal is “reduced” owing to the projection on the $\hat{AL}$ direction. Bottom panel: circles represent scanning angles of Gaia mission for this field. Dark-colored dashed lines show the direction of the astrometric shift vector $\delta$ ±180°, while light-colored dashed lines mark perpendicular directions. Gaia only provides measurements along the instantaneous scanning direction $\hat{AL}$. As a result, whenever the astrometric shift vector $\delta$ and $\hat{AL}$ are aligned or counteraligned ($|\cos \angle (\delta, \hat{AL})| \approx 1$), the signal measured by Gaia is larger. This occurs when the points in the bottom panel are darker and situated closer to one of the dark lines. Indeed, one can see that darker points in the bottom panel correspond to shorter arrows in the middle panel, meaning that more astrometric signal is “available” for Gaia.</figcaption></figure>

Observation and Forecast Tool). We treat each Gaia visit as a single epoch, without dividing into submeasurements from each astrometric field (AF).<fnref n="44" /> To estimate the error bars, we follow the conservative approach of K. A. Rybicki et al. (2018), where the centroiding errors from J. H. J. de Bruijne et al. (2014) are increased by 50% to account for potential systematics. They also take into consideration the fact that there are measurements from multiple AFs within one epoch, which scales down the error bars by a factor of $\sqrt{9} = 3$.

### 4.1. Centroid Trajectory Model

The shift of the centroid position from the position of the source can be expressed as (e.g., M. Dominik & K. C. Sahu 2000)

$$ \delta = \frac{u}{u^2 + 2} \theta_E. \qquad (14) $$

Because we are investigating the Gaia potential to detect the astrometric microlensing effect, for now we disregard the source proper motion and parallax. In simulated data they will be easily distinguished from the microlensing and thus, to first order, should not affect how well the microlensing signal can be recovered. On the other hand, while dealing with the real data, potentially contaminated by unknown systematics, a full model with parallax and proper motions will have to be applied.

Based on Equation (14), one can see that information about the astrometric shift comes from three somewhat distinct instances: the relative separation $u(t)$, the direction of the lens-source relative motion, and the Einstein radius. The first is almost always available, as the relative separation is a function of parameters easily derivable from the standard light curve, namely ($t_0$, $u_0$, $t_E$). This means that, for all standard events that are relatively well covered, one can already predict the astrometric shift of the centroid for each solution (which traces out an ellipse; see, e.g., M. Dominik & K. C. Sahu 2000, for details), although without orienting it on the sky or scaling it to

<footnote n="44">⁴⁴ Gaia constantly rotates, and so whenever an object is observed, it transits the focal plane, passing through (nominally) nine AFs.</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 9.** The cumulative plot of $\Delta\chi^2$ of the astrometric 1D fit for the event OGLE-2016-BLG-0293 for an example realization of the Gaia data and $\theta_E = 0.91$ mas (see Table 4). Top panel: mock Gaia astrometry (gray points) of the centroid displacement during microlensing (proper motion and parallax not included here). The model used to generate the data is marked by the dashed line. Middle panel: residuals of the data and the model shown in the top panel. Bottom panel: cumulative $\Delta\chi^2$, which is the difference between the null model and the astrometric microlensing model (see details in the text).</figcaption></figure>

the absolute angular units. For events considered in this analysis, we also have very strong constraints on the microlensing parallax vector, which in turn gives the direction of the lens-source relative proper motion (see Equation (11)). Consequently, we are able to align the centroid shift trajectory on the sky only using the light-curve model. As a result, for each photometric solution, the only information that needs to be recovered from the astrometric data is $\theta_E$—the scaling factor of the microlensing ellipse, whose shape and orientation are already known. This leads to the conclusion that even astrometric data of seemingly insufficient precision might provide enough information to determine the Einstein radius.

To construct the model of the astrometric microlensing shift as seen by Gaia, one has to project the vector $\delta = (\delta_N, \delta_E)$ onto the scan direction $\widehat{AL}$, and so the 1D astrometric microlensing signal observed by the satellite can be written as

$$ \delta_{\psi} = \delta \cdot \widehat{AL} = \delta_N \cos \psi + \delta_E \sin \psi. \qquad (15) $$

In Figure 8, we visualize the “transition” from the 2D signal to the 1D Gaia data. The absolute astrometric microlensing shift presented in the middle panel (gray solid line) is “reduced” by the projection on the $\widehat{AL}$ direction, which is denoted by arrows. The level of this reduction depends on the

angle between the astrometric shift and the scanning direction, which changes in time, as shown in the bottom panel of the figure.

We use posterior distributions on the physical parameters derived in Section 3.16 to calculate $\theta_E$ and generate 1D Gaia astrometry for all the events in the subsample.

## 4.2. Detectability

To quantify the detectability of astrometric microlensing, we first simulate the Gaia data, following steps and assumptions from the previous subsection. Then, we calculate the $\chi^2$ statistic for the “null model,” which is the case where astrometric microlensing is not present (effectively $\theta_E = 0$):

$$ \chi^2_{\text{null}} = \sum_i \frac{(\delta_{\psi,i} - \delta_{\psi,0,i})^2}{\sigma_i^2}, \qquad (16) $$

where the sum is evaluated over all Gaia measurements and $\sigma_i$ are their respective uncertainties. Because we only consider the astrometric signal from microlensing, our null model is simply the baseline level and thus $\delta_{\psi,0} = 0$. Finally, calculating the difference $\Delta\chi^2 = \chi^2_{\text{null}} - \chi^2$ between the $\chi^2$ of the null model and the correct model with astrometric microlensing allows us
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Detectability of the astrometric microlensing in Gaia data for the selected subsample of long-timescale Spitzer events, plotted as a function of Einstein radius $\theta_E$. The $\Delta\chi^2$ parameter is the difference between the null model and the best-fit astrometric microlensing model (see details in the text). Points represent expected $\theta_E$ values for each event that were calculated in Section 3.17 and listed in Table 4. Thickened lines mark $1\sigma$ errors on $\theta_E$. We also simulated Gaia astrometric data for a wider range of Einstein radii—the detectability for $\theta_E$ in the range from 0.1 to 5 mas is plotted with dashed lines.</figcaption></figure>

to evaluate the confidence level for detecting the astrometric signal present in the Gaia data.

In the top panel of Figure 9, we show an example of 1D astrometric Gaia data for OGLE-2016-BLG-0293, simulated with realistic scatter. As mentioned before, the only information that we need to extract from the astrometry is $\theta_E$—a scaling factor of the model whose shape is already known for a given set of photometric parameters. Thus, even though the visual inspection of the top panel of Figure 9 does not reveal any obvious signal, the bottom panel with the cumulative $\Delta\chi^2$ plot shows that (for event OGLE-2016-BLG-0293) the astrometric signal should be strong enough to detect it.

### 4.3. Results

For all the events in the subsample, we assess the detectability of the astrometric microlensing signatures by comparing the null model to the model generated with the microlensing signal, as explained in the previous section. Two events require some caution in interpreting the results, though. The first one is OGLE-2016-BLG-0689, which is likely blended. Then, the signal measured by Gaia is a combination of the flux from the source and lens/blend(s). As a consequence, the astrometric microlensing signal expected here will be much weaker and the measurement much less robust. The second event is OGLE-2015-BLG-0149. In that case the RUWE parameter is notably larger than unity, which suggests that the Gaia five-parameter astrometric solution can

be unreliable. As a result, the real Gaia data for this event might look different than what we simulate here, and so the measurement of astrometric microlensing signal may be more difficult.

We simulate the Gaia data and calculate $\Delta\chi^2$ for the nine events, using $\theta_E$ values spanning from 0.1 to 5 milliarcseconds, to construct a relation between $\theta_E$ and $\Delta\chi^2$. Having statistical estimates of the lens masses and Einstein radii, presented in Section 3.17, we are able to use this relation to determine for which events we expect to detect the astrometric microlensing signal in the Gaia data. The results are presented in Figure 10, where we show the relations between $\Delta\chi^2$ and Einstein radius, also marking the $\theta_E$ values that are expected from the Bayesian analysis (see Table 4). Out of the nine analyzed events, two of them show promising results in terms of Gaia astrometry detection capabilities: OGLE-2016-BLG-0293 (magenta marker in Figure 10) and OGLE-2018-BLG-0483 (black). The differences between the astrometric null model and the microlensing model for the two events are $\Delta\chi^2_{OB160293} \approx 35$ and $\Delta\chi^2_{OB180483} \approx 24$. Taking into account the one-parameter difference between the null model and microlensing model (1 dof difference), this translates to about a $6\sigma$ and $5\sigma$ expected detection level, respectively.

Generally, there are multiple factors impacting the detectability of the astrometric microlensing signal in the Gaia data. The most important ones are the Einstein radius, the brightness of the event, Gaia sampling, and the scanning angle configuration. It is not a surprise that the two events mentioned above are preferred. In the case of OGLE-2016-BLG-0293, all of the above requirements are met—in addition to the relatively large predicted $\theta_E = 0.91$ mas and several points covering the amplified part of the light curve, the event is bright, having $I_{0,OGLE} = 16.2$ mag. The latter is particularly important because Gaia astrometric accuracy declines steeply with decrease in brightness. On the other hand, OGLE-2018-BLG-0483 is faint, with baseline at $I_{0,OGLE} \approx 18.6$ mag (and even fainter in Gaia, as $I_{0,Gaia} \approx 20.0$ mag), but we expect it to be detected in the Gaia data thanks to its extreme timescale $t_E \approx 330$ days, which translates into a large Einstein radius. Another favorable consequence of the long timescale is the fact that there are more data points from Gaia throughout the (significantly) amplified part of the light curve, which further enhanced detectability.

The second group consists of four events that are less likely to be detected, with the expected astrometric signal on the level of $1.5\sigma-3\sigma$: OGLE-2015-BLG-0145 (red marker in Figure 10), OGLE-2015-BLG-0149 (dark-blue marker), OGLE-2019-BLG-0169 (cyan marker), and OGLE-2016-BLG-0689 (light-green marker). The events from this group are bright and relatively long (see Table 2), but according to our mass/distance analysis, the most probable Einstein radius values for favored solutions are likely too small to be robustly detected by Gaia, even for such bright targets. The OGLE-2015-BLG-0149 event stands out here, but, as already mentioned, it is more challenging to make predictions for owing to the higher RUWE value.

In the case of OGLE-2015-BLG-0211 (gray) and OGLE-2018-BLG-0410 (green) the expected detectability is even smaller, with the astrometric microlensing signal at or below the $1\sigma$ level. The reason is that the Einstein radii are rather small ($\theta_E \lesssim 0.5$ mas) and these events are fainter than the three targets mentioned before.
Finally, the mass measurement of the lens in the OGLE-2018-BLG-0662 (yellow) event is also rather unlikely through the detection of astrometric microlensing in the Gaia data, even though the expected lens mass and Einstein radius are large (see Table 4). In this case, the main factors affecting the detectability prospects are the relatively low brightness of the event ($I_{0,OGLE} = 17.6$; additionally the source is very red, meaning that the brightness in G band is significantly lower) and duration of the event—even though the relatively short timescale $t_E = 66.6$ days translates to a somewhat large Einstein radius $\theta_E = 1.47$ mas, it results in fewer Gaia data points covering the event. Moreover, the scanning angle configuration is such that the 1D astrometric signal is close to zero for a large fraction of the few important points. As a consequence, one of the best black hole candidates in our sample, with an expected dark lens of mass $2.5-4.2 M_\odot$, is unlikely to be confirmed by the Gaia astrometry.

## 5. Summary and Conclusions

In this work, we analyzed a subsample of nine events chosen from the whole population of microlensing events observed by the Spitzer Space Telescope in the years 2014–2019. The subsample was chosen from the events that had also been observed by the Gaia mission, with long timescales and small microlensing parallaxes. The last two requirements were used to identify candidates for massive lenses.

Based on this small subsample, we demonstrate the procedures that will be applied to the whole Spitzer sample of microlensing events, which will allow us to populate the $t_E-\pi_E$ diagram and conduct statistical analysis of microlensing parallax. Such analysis is necessary to complete the studies on the planet frequency in the galaxy, which is based on microlensing planets detected in the Spitzer campaign.

The detailed analysis of the selected subsample of nine events allowed us to identify candidates for black holes and neutron stars that can be later confirmed by the Gaia time-series astrometry, which is expected to be released in early 2026. Based on the Bayesian analysis incorporating the Galactic model and proper-motion information from GDR3 (Section 3.17), we found four candidates for dark remnant lenses: OGLE-2016-BLG-0293, OGLE-2018-BLG-0483, OGLE-2018-BLG-0662, and OGLE-2015-BLG-0149 (see Table 4 for estimated masses). The masses of the lenses of the four candidate events lie somewhat on the edge of the known distinction between the neutron stars and stellar black holes. In the case of OGLE-2015-BLG-0149, it is most likely a neutron star. As for the remaining three candidates, the median mass suggests black holes, although the error bars are large, and so the heavy neutron star scenario cannot be excluded. The cases of OGLE-2016-BLG-0293 and OGLE-2018-BLG-0483 are expected to have an astrometric microlensing signal detectable by Gaia, which in turn will allow us to confirm the masses of the lenses and their remnant nature.

We have found individual events that are likely caused by remnant lenses, but one cannot infer more general information about the dark remnant population based on this analysis. To do that, one would have to estimate the detection efficiency, which is a function of Einstein timescale, microlensing parallax, and blend flux. In addition, the whole population of Spitzer events would have to be thoroughly modeled and investigated in the same manner as the subsample of nine events analyzed here—only then would one have a robust picture of $t_E-\pi_E$ space populated by Spitzer events. Finally, the

selection bias related to Spitzer observations has to be taken into account. This last task is not trivial, as the selection procedure for Spitzer microlensing targets was somewhat complicated (although consistent and well-defined; see J. C. Yee et al. 2015). Taking into account these aspects requires a separate analysis that is beyond the scope of this paper and will be covered elsewhere.

## Acknowledgments

E.O.O. is grateful for the support of grants from the Willner Family Leadership Institute, André Deloro Institute, Paul and Tina Gardner, The Norman E. Alexander Family M Foundation ULTRASAT Data Center Fund, Israel Science Foundation, Israeli Ministry of Science, Minerva, NSF-BSF, Israel Council for Higher Education (VATAT), Sagol Weizmann-MIT, Yeda-Sela, and the Rosa and Emilio Segre Research Award. This research was supported by the Institute for Environmental Sustainability (IES) and the André Deloro Institute for Space and Optics Research at the Weizmann Institute of Science.

This research has made use of the KMTNet system operated by the Korea Astronomy and Space Science Institute (KASI) at three host sites of CTIO in Chile, SAAO in South Africa, and SSO in Australia. Data transfer from the host site to KASI was supported by the Korea Research Environment Open NETwork (KREONET). This research was supported by KASI under the R&D program (project No. 2024-1-832-01) supervised by the Ministry of Science and ICT. W.Z., H.Y., S.M., R.K., J.Z., and W.Z. acknowledge support by the National Natural Science Foundation of China (grant No. 12133005). W.Z. acknowledges the support from the Harvard-Smithsonian Center for Astrophysics through the CfA Fellowship. J.C.Y. and I.-G.S. acknowledge support from US NSF grant No. AST-2108414. Y.S. acknowledges support from BSF grant No. 2020740. Work by C.H. was supported by the grants of the National Research Foundation of Korea (2019R1A2C2085965 and 2020R1A4A2002885).

The MOA project is supported by JSPS KAKENHI grant Nos. JP24253004, JP26247023, JP16H06287, and JP22H00153.

## ORCID iDs

Krzysztof A. Rybicki https://orcid.org/0000-0002-9326-9329
Yossi Shvartzvald https://orcid.org/0000-0003-1525-5041
Jennifer C. Yee https://orcid.org/0000-0001-9481-7123
Calen Henderson https://orcid.org/0000-0001-8877-9060
Michael M. Fausnaugh https://orcid.org/0000-0002-9113-7162
Yoon-Hyun Ryu https://orcid.org/0000-0001-9823-2907
Michael D. Albrow https://orcid.org/0000-0003-3316-4012
Sun-Ju Chung https://orcid.org/0000-0001-6285-4528
Cheongho Han https://orcid.org/0000-0002-2641-9964
Kyu-Ha Hwang https://orcid.org/0000-0002-9241-4117
Youn Kil Jung https://orcid.org/0000-0002-0314-6000
In-Gu Shin https://orcid.org/0000-0002-4355-9838
Hongjing Yang https://orcid.org/0000-0003-0626-8465
Weicheng Zang https://orcid.org/0000-0001-6000-3463
Sang-Mok Cha https://orcid.org/0000-0002-7511-2950
Seung-Lee Kim https://orcid.org/0000-0003-0562-5643
Chung-Uk Lee https://orcid.org/0000-0003-0043-3925
Dong-Joo Lee https://orcid.org/0009-0000-5737-0908
Yongseok Lee https://orcid.org/0000-0001-7594-8072
Byeong-Gon Park https://orcid.org/0000-0002-6982-7722
Richard W. Pogge https://orcid.org/0000-0003-1435-3053
David P. Bennett https://orcid.org/0000-0001-8043-8413

Yuki K. Satoh https://orcid.org/0000-0002-1228-4122

## References

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016, PhRvL, 116, 061102

Bachelet, E., Norbury, M., Bozza, V., & Street, R. 2017, AJ, 154, 203

Bachelet, E., Zielinski, P., Gromadzki, M., et al. 2022, A&A, 657, A17

Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M., Demleitner, M., & Andrae, R. 2021, AJ, 161, 147

Batista, V., Gould, A., Dieters, S., et al. 2011, A&A, 529, A102

Bond, I. A., Abe, F., Dodd, R. J., et al. 2001, MNRAS, 327, 868

Bozza, V., Bachelet, E., Bartolić, F., et al. 2018, MNRAS, 479, 5157

Calchi Novati, S., Gould, A., Udalski, A., et al. 2015, ApJ, 804, 20

Cassan, A., & Ranc, C. 2016, MNRAS, 458, 2074

Cassan, A., Ranc, C., Absil, O., et al. 2022, NatAs, 6, 121

Dalal, N., & Lane, B. F. 2003, ApJ, 589, 199

de Bruijne, J. H. J., Rygl, K. L. J., & Antoja, T. 2014, The Milky Way Unravelled by Gaia: GREAT Science from the Gaia Data Releases (Les Ulis: EDP Sciences), 23

Delplancke, F., Górski, K. M., & Richichi, A. 2001, A&A, 375, 701

DePoy, D. L., Atwood, B., Belville, S. R., et al. 2003, Proc. SPIE, 4841, 827

Dominik, M., & Sahu, K. C. 2000, ApJ, 534, 213

Dong, S., Mérand, A., Delplancke-Ströbele, F., et al. 2019, ApJ, 871, 70

Dong, S., Udalski, A., Gould, A., et al. 2007, ApJ, 664, 862

Evans, D. W., Riello, M., De Angeli, F., et al. 2018, A&A, 616, A4

Fardeen, J., McGill, P., Perkins, S. E., et al. 2024, ApJ, 965, 138

Fukui, A., Suzuki, D., Koshimoto, N., et al. 2019, AJ, 158, 206

Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 649, A1

Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1

Gaia Collaboration, Vallenari, A., Brown, A. G. A., et al. 2023, A&A, 674, A1

Gould, A. 1994, ApJL, 421, L75

Gould, A. 2004, ApJ, 606, 319

Gould, A. 2019, JKAS, 52, 121

Gould, A. 2023, arXiv:2310.19164

Gould, A., Ryu, Y.-H., Calchi Novati, S., et al. 2020, JKAS, 53, 9

Gravity+ Collaboration, Abuter, R., Alarcon, P., et al. 2022, Msngr, 189, 17

Han, C., & Gould, A. 1995, ApJ, 449, 521

Han, C., & Gould, A. 2003, ApJ, 592, 172

Hobbs, G., Lorimer, D. R., Lyne, A. G., & Kramer, M. 2005, MNRAS, 360, 974

Hodgkin, S. T., Harrison, D. L., Breedt, E., et al. 2021, A&A, 652, A76

Hodgkin, S. T., Wyrzykowski, L., Blagorodnova, N., & Koposov, S. 2013, RSPTA, 371, 20120239

Howil, K., Wyrzykowski, Ł., Kruszyńska, K., et al. 2024, arXiv:2403.09006

Kim, S. -L., Lee, C.-U., Park, B.-G., et al. 2016, JKAS, 49, 37

Klüter, J., Bastian, U., & Wambsganss, J. 2020, A&A, 640, A83

Kroupa, P. 2001, MNRAS, 322, 231

Kruszyńska, K., Wyrzykowski, Ł., Rybicki, K. A., et al. 2022, A&A, 662, A59

Lam, C. Y., & Lu, J. R. 2023, ApJ, 955, 116

Lam, C. Y., Lu, J. R., Udalski, A., et al. 2022, ApJL, 933, L23

Lu, J. R., Sinukoff, E., Ofek, E. O., Udalski, A., & Kozlowski, S. 2016, ApJ, 830, 41

McGill, P., Anderson, J., Casertano, S., et al. 2023, MNRAS, 520, 259

Mróz, P., Udalski, A., & Gould, A. 2022, ApJL, 937, L24

Mroz, P., Udalski, A., Wyrzykowski, L., et al. 2021, arXiv:2107.13697

Mróz, P., & Wyrzykowski, Ł. 2021, AcA, 71, 89

Nataf, D. M., Gould, A., Fouqué, P., et al. 2013, ApJ, 769, 88

Nucita, A. A., Licchelli, D., De Paolis, F., et al. 2018, MNRAS, 476, 2962

Paczynski, B. 1986, ApJ, 304, 1

Pecaut, M. J., & Mamajek, E. E. 2013, ApJS, 208, 9

Rattenbury, N. J., & Mao, S. 2006, MNRAS, 365, 792

Refsdal, S. 1966, MNRAS, 134, 315

Repetto, S., Igoshev, A. P., & Nelemans, G. 2017, MNRAS, 467, 298

Rybicki, K. A., Wyrzykowski, Ł., Bachelet, E., et al. 2022, A&A, 657, A18

Rybicki, K. A., Wyrzykowski, Ł., Klencki, J., et al. 2018, MNRAS, 476, 2013

Sahu, K. C., Anderson, J., Casertano, S., et al. 2017, Sci, 356, 1046

Sahu, K. C., Anderson, J., Casertano, S., et al. 2022, ApJ, 933, 83

Sajadian, S., & Sahu, K. C. 2023, AJ, 165, 96

Schönrich, R., Binney, J., & Dehnen, W. 2010, MNRAS, 403, 1829

Udalski, A., Szymanski, M., Kaluzny, J., Kubiak, M., & Mateo, M. 1992, AcA, 42, 253

Udalski, A., Szymański, M. K., & Szymański, G. 2015a, AcA, 65, 1

Udalski, A., Yee, J. C., Gould, A., et al. 2015b, ApJ, 799, 237

Wozniak, P. R. 2000, AcA, 50, 421

Wyrzykowski, Ł., Kostrzewa-Rutkowska, Z., Skowron, J., et al. 2016, MNRAS, 458, 3012

Wyrzykowski, Ł., Kruszyńska, K., Rybicki, K. A., et al. 2023, A&A, 674, A23

Wyrzykowski, Ł., & Mandel, I. 2020, A&A, 636, A20

Yee, J. C., Shvartzvald, Y., Gal-Yam, A., et al. 2012, ApJ, 755, 102

Yee, J. C., Udalski, A., Calchi Novati, S., et al. 2015, ApJ, 802, 76

Zang, W., Dong, S., Gould, A., et al. 2020, ApJ, 897, 180

Zurlo, A., Gratton, R., Mesa, D., et al. 2018, MNRAS, 480, 236