# New limits from microlensing on Galactic black holes in the mass range 10 M < M < 1000 M

T. Blaineau, M. Moniez, C. Afonso, J.-N. Albert, R. Ansari, E. Aubourg, C.

Coutures, J.-F. Glicenstein, B. Goldman, C. Hamadache, et al.

► To cite this version:

T. Blaineau, M. Moniez, C. Afonso, J.-N. Albert, R. Ansari, et al. New limits from microlensing on Galactic black holes in the mass range 10 M < M < 1000 M. Astronomy and Astrophysics - A&A, 2022, 664, pp.A106. 10.1051/0004-6361/202243430. hal-03602697

HAL Id: hal-03602697

https://hal.science/hal-03602697

Submitted on 13 Apr 2023

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
# New limits from microlensing on Galactic black holes in the mass range $10 M_{\odot} < M < 1000 M_{\odot}$

T. Blaineau<fnref n="1" />, M. Moniez<fnref n="1" />●, C. Afonso<fnref n="2" />, J.-N. Albert<fnref n="1" />, R. Ansari<fnref n="1" />, E. Aubourg<fnref n="2" />,<fnref n="3" />, C. Coutures<fnref n="2" />, J.-F. Glicenstein<fnref n="2" />, B. Goldman<fnref n="4" />,<fnref n="5" />, C. Hamadache<fnref n="1" />, T. Lasserre<fnref n="2" />, L. Le Guillou<fnref n="6" />, E. Lesquoy<fnref n="2" />, C. Magneville<fnref n="2" />, J.-B. Marquette<fnref n="7" />, N. Palanque-Delabrouille<fnref n="2" />,<fnref n="8" />, O. Perdereau<fnref n="1" />, J. Rich<fnref n="2" />, M. Spiro<fnref n="2" />, and P. Tisserand<fnref n="9" />

Received 27 February 2022 / Accepted 3 June 2022

## ABSTRACT

We searched for long-duration microlensing events originating from intermediate-mass black holes (BH) in the halo of the Milky Way, using archival data from the EROS-2 and MACHO photometric surveys towards the Large Magellanic Cloud (LMC). We combined data from these two surveys to create a common database of light curves for 14.1 million objects in the LMC, covering a total duration of 10.6 years, with flux series measured in four wide passbands. We carried out a microlensing search on these light curves, complemented by the light curves of 22.7 million objects, observed only by EROS-2 or only by MACHO, over about 7 years, with flux series measured in only two passbands. A likelihood analysis, taking into account the LMC self-lensing and Milky Way disk contributions, allows us to conclude that compact objects with masses in the range $10–100 M_{\odot}$ cannot make up more than ~15% of a standard halo total mass (at a 95% confidence level). Our analysis sensitivity weakens for heavier objects, although we still rule out the possibility of ~50% of the halo being made of ~1000 $M_{\odot}$ BHs. Combined with previous EROS results, an upper limit of ~15% of the total halo mass can be obtained for the contribution of compact halo objects in the mass range $10^{-6}–10^2 M_{\odot}$.

**Key words.** gravitational lensing: micro – surveys – stars: black holes – Galaxy: halo – Galaxy: kinematics and dynamics – dark matter

## 1. Introduction

Observations of gravitational waves due to the coalescence of massive objects (Abbott et al. 2016a,b) have demonstrated the existence of merging black holes (BHs) heavier than $10 M_{\odot}$, which has renewed interest in BHs as dark matter, especially primordial black holes (PBHs; Bird et al. 2016; Green & Kavanagh 2021; Sasaki et al. 2016). Microlensing surveys toward the Large Magellanic Cloud (LMC), which allow the content of the Galactic halo in massive compact objects to be probed, have shown that objects lighter than $10 M_{\odot}$ do not significantly contribute to the hidden mass of the Galactic spherical halo of our galaxy (Tisserand et al. 2007; Alcock et al. 2001; Wyrzykowski et al. 2011). Since the typical duration of microlensing events increases with the lens mass, in order to detect heavier objects (such as those responsible for gravitational wave emissions), time-series data longer than the durations of each of the EROS-2 (Expérience de Recherche d'Objets Sombres) and MACHO (MAssive Compact Halo Objects) surveys (two of the first microlensing surveys that operated in the years 90s and 2000s) are needed. To explore the dark

matter halo beyond this $10 M_{\odot}$ limit by searching for events with longer timescales, we have combined the databases of EROS-2 and MACHO, thus starting the program described in Mirhosseini & Moniez (2018).

In Sect. 2 we describe the microlensing effect, focusing on the LMC searches. In Sect. 3 we introduce the EROS-2 and MACHO surveys and their light curve data sets. We summarize how we associated objects in the two catalogs, and present the procedures that allowed us to remove defective images and measurements. In Sect. 4 we describe the selection of candidates for gravitational microlensing events, based on a fitting of the observed light curves with theoretical microlensing curves. In Sect. 5 we explain the calculation of the efficiency of the detection of microlensing events. We quantify the effect of blending on the detection efficiency by using *Hubble* Space Telescope (HST) data, and discuss the impact of the binary sources. In Sect. 6, we confront the number of selected candidates with the number of events expected from the dark matter halo, from the Galactic disk, and from the LMC itself. We then derive a new upper limit on the contribution of compact objects to the halo. Finally, in Sect. 7 we list the sources of improvement we

<footnote n="1">¹ Laboratoire de Physique des 2 infinis Irène Joliot-Curie, CNRS Université Paris-Saclay, Bât. 100, Faculté des Sciences, 91405 Orsay Cedex, France
e-mail: moniez@lal.in2p3.fr</footnote>

<footnote n="2">² IRFU, CEA, Université de Paris-Saclay, 91191 Gif-sur-Yvette, France</footnote>

<footnote n="3">³ Université de Paris, CNRS, Astroparticule et Cosmologie, 75013 Paris, France</footnote>

<footnote n="4">⁴ The International Space University, 1 Rue Jean-Dominique Cassini, 67400 Illkirch-Graffenstaden, Germany</footnote>

<footnote n="5">⁵ Université de Strasbourg, CNRS, Observatoire Astronomique, CNRS, UMR 7550, 67000 Strasbourg, France</footnote>

<footnote n="6">⁶ Sorbonne Université, CNRS/IN2P3, Laboratoire de Physique Nucléaire et de Hautes Énergies (LPNHE), 75005 Paris, France</footnote>

<footnote n="7">⁷ Laboratoire d'astrophysique de Bordeaux, Univ. Bordeaux, CNRS, B18N, Allée Geoffroy Saint-Hilaire, 33615 Pessac, France</footnote>

<footnote n="8">⁸ Lawrence Berkeley National Laboratory, One Cyclotron Road, Berkeley, CA 94720, USA</footnote>

<footnote n="9">⁹ Sorbonne Universités, UPMC Univ Paris 6 et CNRS, UMR 7095, Institut d'Astrophysique de Paris, IAP, 75014 Paris, France</footnote>
achieved with respect to previous results, and propose some perspectives for further data combinations.

## 2. Microlensing toward the LMC

A gravitational microlensing effect occurs when a massive compact object (called a lens or a deflector in the following) passes close enough to the line of sight of a star to produce gravitational images that are not intercepted by the lens. The size of the opaque part of the lens and the relative positions of the source, lens, and observer must be such that the rays of the two gravitational images are not occulted (no eclipse). For the typical lens-source configurations considered here, the angular separation of the two images is too small to be resolved in telescopes. The detection of the event is made possible by the relative motion of the lens in the observer-source frame, which produces a transient variation of the source brightness.

Microlensing of the LMC stars as a technique to search for massive compact objects in the Galactic halo was first described in Paczynski (1986). Reviews of the formalism can be found in Schneider et al. (2006) and Rahvar (2015). When a point-like object (lens) of mass $M_L$ located at distance $D_L$ passes close to the line of sight of a point source located at distance $D_{LMC} = 49.5$ kpc (Pietrzyński et al. 2019), the total magnification of the source luminosity at a given time $t$ is the sum of the contributions of two images, given by:

$$A(t) = \frac{u(t)^2 + 2}{u(t) \sqrt{u(t)^2 + 4}}, \qquad (1)$$

where $u(t)$ is the distance of the lens to the undeflected line of sight, divided by the Einstein radius $r_E$,

$$r_E = \sqrt{\frac{4GM_L}{c^2} D_{LMC} x(1-x)} \approx 10.0 \text{ AU} \times \left[ \frac{M_L}{M_\odot} \right]^{\frac{1}{2}} \frac{[x(1-x)]^{\frac{1}{2}}}{0.5}. \qquad (2)$$

Here $G$ is the Newtonian gravitational constant and $x$ the lens to source distance ratio $x = D_L/D_{LMC}$. The Einstein radius of the lens is such that a point source that is behind the Einstein disk (of surface $\pi v_E^2$), sees its apparent luminosity magnified by a factor greater than 1.34. Assuming that the lens has a constant relative transverse velocity $v_T$, $u(t)$ is given by:

$$u(t) = \sqrt{u_0^2 + (t - t_0)^2/t_E^2}, \qquad (3)$$

where $t_E = r_E/v_T$ is the Einstein radius crossing time, and $u_0$ is the minimum distance to the undeflected line of sight at time $t_0$.

In the approximation of a point lens acting on a point source, with a uniform relative motion with respect to the line of sight (hereafter called PSPL approximation), the microlensing effect has several characteristic features that allow one to discriminate it from any known intrinsic stellar variability: Given the low probability of alignment required for a measurable microlensing effect to occur, it is expected that such an event will not be repeated for a given source or lens on typical human timescales; the magnification is a known function of time, depending on only three parameters ($u_0$, $t_0$, and $t_E$), with a symmetrical shape, and independent of the passband; since the geometric configuration of the source-deflector system is random, the prior distributions of $t_0$ and of the impact parameters $u_0$ of the events must be uniform; the passive role of the lensed stars implies that they should be representative of the monitored sample.

<table><caption>**Table 1.** Characteristics of EROS-2 and MACHO setups.</caption><thead><tr><th></th><th>EROS-2</th><th>MACHO</th></tr></thead><tbody><tr><td>Telescope</td><td>1 m</td><td>1.27 m</td></tr><tr><td>Pixel size</td><td>0.62"</td><td>0.63"</td></tr><tr><td>Blue passband</td><td>[420, 720] nm</td><td>[450, 590] nm</td></tr><tr><td>Red passband</td><td>[620, 920] nm</td><td>[590, 780] nm</td></tr><tr><td>Median image quality</td><td>2"</td><td>2.1"</td></tr></tbody></table>

If the gravitational field of the Galaxy is entirely due to the lenses, the optical depth for lensing, (i.e., the probability that the line of sight to an LMC star is within one $r_E$ of a lens) is of the order of $v_{rot}^2/c^2$, where $v_{rot}$ is the orbital velocity around the Galaxy at the LMC position. We use as a benchmark model the isotropic and isothermal halo first studied by Griest (1991); hereafter called the S-model, with the mass density distribution

$$\rho_H(r) = 0.0078 \frac{R_0^2 + R_c^2}{r^2 + R_c^2} M_\odot \text{pc}^{-3}, \qquad (4)$$

where $R_0 = 8.5$ kpc is the Galactocentric distance to the Sun, $R_c = 5$ kpc is the halo “core radius”, and $r$ is the Galactocentric radius. We use this model with the most recent values of the LMC distance $D_{LMC} = 49.5$ kpc (Pietrzyński et al. 2019).

The velocity $v_L$ of the halo objects follows the Maxwell distribution

$$p(v_L) = \left(\frac{1}{2\pi v_0^2}\right)^{3/2} 4\pi v_L^2 e^{-v_L^2/2v_0^2}, \qquad (5)$$

with $v_0 = 120 \text{ km s}^{-1}$; the velocity of the Sun is $v_\odot = (11.1, 251, 7.3) \text{ km s}^{-1}$ (Brunthaler et al. 2011) in galactocentric Cartesian coordinates (X pointing from the Sun to the Galactic center, Z pointing north), and the proper velocity of the LMC is $v_{LMC} = (-57, -226, 221) \text{ km s}^{-1}$ (Kallivayalil et al. 2013).

This model gives an optical depth to the LMC $\tau_{LMC} \sim 4.7 \times 10^{-7}$ if the halo is completely made of compact objects. The total rate (events per star per unit time) for $u_0 < 1$ is $\Gamma = (2/\pi)\tau_{LMC}/\langle t_E \rangle$ (Griest 1991), where $\langle t_E \rangle$ is the mean $t_E$ of all events<fnref n="1" />. For lenses of mass $M_L$, the mean event duration for this benchmark halo model is $\langle t_E \rangle \sim 63 \text{ day} \sqrt{M_L/M_\odot}$.

## 3. Combining EROS-2 and MACHO data

The EROS-2 and MACHO surveys were performed with similar setups, respectively installed at the La Silla Observatory (ESO, Chile) and at the Mount Stromlo Observatory (Australia) (see Table 1). EROS-2 used a 1 m (F/5) diameter telescope, equipped with a dichroic beam splitter and two cameras, each with 8 $2K \times 2K$ CCD's, covering 1 deg². MACHO used a slightly larger telescope (1.27 m, F/3.9), equipped with two cameras, each with 4 $2K \times 2K$ CCD's, covering a smaller field of view (0.5 deg²).

EROS-2 surveyed 88 fields (1. deg²) toward the LMC, and MACHO surveyed 82 (0.5 deg²) LMC fields. MACHO adopted the same exposure time (300 s) for all fields, whereas EROS-2 adapted its exposure times according to the surface brightness

<footnote n="1">¹ Rahvar (2015) uses $\langle 1/t_E \rangle$ instead of $1/(\langle t_E \rangle)$ to estimate the event rate. The former refers to an average over all the events in progress (within an Einstein ring) at a given time, whereas the latter refers to an average over all the events occurring within a given time interval.</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Distribution of the EROS-2 and MACHO magnitude difference $\Delta m = m_{BE} - m_{BM}$ as a function of the EROS-2 blue magnitude $m_{BE}$ for the 14.1 million associated objects; the values of $m_{BM}$ are derived from the original MACHO magnitudes using a first-degree color equation, so as to match the EROS blue magnitudes on average. As a result, the mean value of $\Delta m$ varies by a few percent depending on the color of the sources, in particular between the main sequence and the red giants branch (around $m_{BE} = 18.8$). White dots and bars show the average and standard deviations of $\Delta m$ for each $m_{BE}$ slice. The color scale is expressed in millions of objects per squared magnitude.</figcaption></figure>

of the fields, varying from 180 s (for the central fields) to 900 s (for the external fields). After the end of 1999, the longest exposure times were reduced in EROS-2, in order to increase the overall sampling of the survey. The EROS-2 object catalog was produced after co-adding at least 10 images per field, and rejecting the faintest and brightest objects. This allowed it to partially compensate for the effect of a shorter exposure in the densest fields and a smaller telescope diameter compared to MACHO.

The MACHO light curves and images are publicly available<fnref n="2" /> (Alcock et al. 1999). The catalog contains $22.3 \times 10^6$ objects with magnitudes $V < 21.5$, of which 6% are duplicate objects due to overlapping fields. The light curves cover a duration of 7.7 years, longer than the duration analyzed in the last MACHO publication (5.7 years; Alcock et al. 2001).

The EROS-2 catalog for the LMC was produced by Tisserand (2004) for the final EROS LMC publication (Tisserand et al. 2007). It will be made public to allow future expansion of the work described in this paper.

To associate objects in the two catalogs, available MACHO and EROS-2 sky coordinates (RA,Dec) were first refined using Gaia Early Data Release 3 astrometry (Gaia Collaboration 2021), correcting for local shifts up to 2'' for MACHO and 0.5'' for EROS-2. After this correction, we could associate the objects of the two surveys with a precision of better than 0.1''. Given the typical spread of the light of the stars on the best images (FWHM ~ 1''), and the similar resolutions of the surveys, the associated reconstructed objects of each survey contain the same stars, with little variation of the blend components. This is confirmed by the good correlation observed between the EROS-2 and MACHO fluxes (Fig. 1), which is compatible with the photometric accuracy.

Table 2 summarizes the characteristics of the two surveys and the cross-matched catalog. The common set consists of the $14.1 \times 10^6$ objects cross-matched in the EROS-2 and MACHO catalogs; they benefit from a total of 10.6 years of luminosity

<table><caption>**Table 2.** Statistics of the EROS-2 and MACHO surveys.</caption><thead><tr><th></th><th>EROS-2 only</th><th>MACHO only</th><th>Common</th></tr></thead><tbody><tr><td>Dates (month/yr)</td><td>7/96-2/03</td><td>7/92-1/00</td><td>7/92-2/03</td></tr><tr><td>T<sub>obs</sub> (year)</td><td>6.7</td><td>7.7</td><td>10.6</td></tr><tr><td>N<sub>objects</sub>(×10<sup>6</sup>)</td><td>15.8</td><td>6.9 <sup>(a)</sup></td><td>14.1</td></tr><tr><td colspan="4">Central fields</td></tr><tr><td>deg<sup>2</sup></td><td>~10</td><td>~10</td><td>~10</td></tr><tr><td>stars/arcmin<sup>2</sup></td><td>~70</td><td>~100</td><td>~70</td></tr><tr><td>Mag. lim. V<sub>Cousins</sub></td><td>~20.5</td><td>~20.5</td><td>~20.5</td></tr><tr><td># measurements B</td><td>500</td><td>1400</td><td>1900</td></tr><tr><td># measurements R</td><td>600</td><td>1550</td><td>2150</td></tr><tr><td colspan="4">Outer fields</td></tr><tr><td>deg<sup>2</sup></td><td>~77</td><td>~39</td><td>~39</td></tr><tr><td>stars/arcmin<sup>2</sup></td><td>~30</td><td>~20</td><td>~20</td></tr><tr><td>mag. lim. V<sub>Cousins</sub></td><td>~22.5</td><td>~21.5</td><td>~21.5</td></tr><tr><td># measurements B</td><td>250</td><td>200</td><td>450</td></tr><tr><td># measurements R</td><td>300</td><td>250</td><td>550</td></tr></tbody></table>

**Notes.** Showing survey durations, the number of monitored sources, the median stellar densities, approximate limiting magnitudes, and median numbers of flux measurements per object after cleaning. (a) Number of unique objects.

measurements (of which 3.8 years overlap), during which the surveys used four different passbands. The complementary set, consisting of objects reconstructed only in EROS-2 or only in MACHO, comprising $22.7 \times 10^6$ objects monitored by a single survey, is also included in our analysis, although the objects are monitored in only two passbands for shorter times.

## 3.1. Removing problematic images and measurements

We found that bad images and/or measurements still polluted our light curves sample. MACHO images with more than 5% of measurements that were more than five standard deviations ($5\sigma$) away from the reference magnitude were found to be mostly faulty (visual inspection showed blurred images or with guiding or readout defect), and we then discarded them from our analysis. We also rejected low-quality EROS-2 images using similar criteria (Blaineau 2021). We further discarded measurements that deviated by more than $5\sigma$ from the median flux in a sliding window of five consecutive measurements, unless this occurred in more than 10% of the cases (so as not to discard curves with many large and rapid variations). This operation would be penalizing when searching for events of very short duration, but it has no impact on the efficiency for events lasting more than a few months. These measurement-quality cuts removed about ~3% (~1%) of the measurements from the MACHO (EROS-2) data.

## 3.2. Corrections on flux uncertainties

By comparing the mean photometric uncertainties with the point-to-point flux dispersions along the light curves, it appeared that the MACHO uncertainties were underestimated by an average factor of ~0.73, while those of EROS-2 were incorrectly estimated by a factor varying from ~0.75 (for stars brighter than $I = 19$) to ~1.45 (for $I > 21$). In order to balance the weights of each survey in the calculations of the goodness-of-fit statistics $\chi^2$, we renormalized the squared uncertainties for each light curve (one per passband) with a quantity $X$ that characterizes the average fluctuations around a global trend, estimated from the differences between the measurements and their nearest

<footnote n="2">² https://macho.nci.org.au/</footnote>
neighbors:

$$X = \frac{1}{N-2} \sum_{i=2}^{N-1} \frac{(\phi(t_i) - \phi_i^{\text{int}})^2}{\sigma_{\text{int}}^2} \quad (6)$$

Here, $\phi(t_i)$ is the flux measured at time $t_i$, $\phi_i^{\text{int}}$ is the flux inter- polated from measurements $i-1$ and $i+1$ at $t_i$, and $\sigma_{\text{int}}$ is the uncertainty on $(\phi(t_i)-\phi_i^{\text{int}})$ deduced from the cataloged uncertain- ties. The sum is over points with a precision better than 0.55 mag. Renormalizing squared uncertainties by dividing them by $X$ gen- erally made the $\chi^2$ per degree of freedom (d.o.f.) near unity for fits assuming slow flux variations over time. As the correction factors to be applied to photometric uncertainties vary with mag- nitudes, this renormalization procedure is not effective for light curves with large variations, especially for low fluxes with large uncertainties (>0.55 mag) in the EROS-2 data.

# 4. Search for long-timescale microlensing events

Thanks to the increase in computing power available compared to the 1990s and 2000s, we were able to conduct an analysis essentially based on the comparison of the fit of a microlensing effect with that of a constant light curve, without preselection cri- teria. More precisely, for each object we performed a simultane- ous PSPL microlensing fit to the available light curves, with one base flux line per passband and a set of microlensing parameters ($t_0$, $u_0$, and $t_E$) common to all passbands (fit without blending or parallax). The selections described below were also applied to the simulated events, as described in Sect. 5, in order to deter- mine the analysis efficiency.

The impact of the Earth's rotation (parallax) around the Sun on the detection efficiency of multiyear events was discussed in (Blaineau & Moniez 2020). The effect on a tolerant search algo- rithm, such as the one we describe below, was found to be negligi- gible. Nevertheless, simulated light curves included the parallax effect, even though for the search in the data (real and simulated) we only tried to fit a PSPL microlensing effect.

## 4.1. Selection of candidates

The light curves of the $36.8 \times 10^6$ cataloged objects underwent the following selection process in order to identify long-duration microlensing candidates. This process was tuned with the simu- lation described in Sect. 5.

We eliminated objects with fewer than 200 total measure- ment points, all passbands included, and with fewer than 50 mea- surements in at least both EROS-2 or both MACHO passbands. This selection left $36 \times 10^6$ objects.

We required that the light curves of the object simultaneously fit well a PSPL microlensing event, with a global $\chi^2$ (summed for all passbands) significantly smaller than that for constant curves. Specifically, we required that

$$\Delta\chi^2 = \frac{\chi_{\text{const.}}^2 - \chi_{\text{ML}}^2}{\chi_{\text{ML}}^2/N_{\text{d.o.f.}}} \frac{1}{\sqrt{2N_{\text{d.o.f.}}}} > 80, \quad \text{and} \quad \frac{\chi_{\text{ML}}^2}{N_{\text{d.o.f.}}} < f(\Delta\chi^2), \quad (7)$$

where the function

$$f(\Delta\chi^2) = 1.44 + 0.26 \log_{10} \left( \frac{\Delta\chi^2}{80} \right) + 0.23 \log_{10}^2 \left( \frac{\Delta\chi^2}{80} \right) \quad (8)$$

is tuned to accept 92.5% of the simulated light curves for each given interval of $\Delta\chi^2$, and is in the range $1.4 < f(\Delta\chi^2) < 3$ for our entire data set. These criteria accepted 352 light curves,

and between 28% and 48% of the events simulated within the observation duration with $u_0 < 1$ and $100 < t_E < 1000$ days (see Sect. 5).

We required that the global microlensing fit did not result in a significant degradation of $\chi_{\text{ML}}^2$ compared to $\chi_{\text{const.}}^2$ for any of the four passbands. Specifically, we required that $\Delta\chi^2 > -0.1$, where $\Delta\chi^2$ was restricted to each light curve. This requirement was satisfied by all the simulated events selected so far, and 226 light curves remained at this stage. All known events from the MACHO and EROS publications remained at this stage, since the aforementioned criteria did not select preferentially long- duration events.

We required that events had a maximum well within the observing period: $t_{\text{min}}+200$ days $< t_0 < t_{\text{max}}$, where $t_{\text{min}}$ and $t_{\text{max}}$ are the start and end dates of the light curves. This asymmetric requirement, which allowed sufficient information about the rise phase to be obtained, eliminated transients and some long-period variable stars that usually have a shorter rise time than fall time, leaving 148 light curves.

We required that $t_E$ be in the range $100$ days $< t_E < (t_{\text{max}} - t_{\text{min}})/2$. This criterion rejected most supernovae and short-lived fluctuations, and ensured that the fitted variations were sufficiently within the observing period.

34 light curves remained after these two last requirements. Since these criteria apply to the time parameters, they have an impact on the simulation that varies strongly with the mass of the lenses (from 17% relative acceptance for $1 M_\odot$ to 62% for $1000 M_\odot$).

We eliminated so-called “blue bumpers” by making a stricter selection on events in their region of the color-magnitude dia- gram. These stars are Be class stars located in the blue and bright zone of the Color-Magnitude Diagram (CMD). They sometimes present asymmetric bumps with a faster rise than fall, and their luminosity variations, probably related to the dynamics of the decretion disk (Hubert et al. 2007), can extend over years. We rejected these artifacts by requiring that the fitted value of $u_0$ be less than 0.9 for objects in the following CMD domain: $I_{\text{Cousins}} < 19$ and ($V_{\text{Johnson}} - I_{\text{Cousins}}) < 0.5$. 28 candidate events remained after this selection.

## 4.2. Remaining candidates: rejection of known artifacts

At this point, there are still some physical phenomena that cause changes in the objects' brightness, which can be mistaken for long-duration microlensing effects. We eliminated three types of objects:

We rejected the objects located outside the CMD region containing 99% of the stars, to discard in particular the rarest objects, likely to show variability. 24 objects remained after this CMD-based rejection, at the cost of 1% of the expected events.

To avoid echoes from SN1987A (Suntzeff et al. 1988) we discarded the objects located in a zone defined by $83.670^\circ < \alpha < 84.064^\circ$ and $-69.34^\circ < \delta < -69.20^\circ$. After this exclusion, 5 objects remained, while almost 100% of the simu- lated events selected so far were retained.

As a final step in our selection process, we queried external catalogs to help identify known variabilities. We rejected one object at $(\alpha, \delta)_{J2000} = (80.6256, -71.7500)$, iden- tified on external catalogs as a quasi-stellar object (Kim et al. 2012; Kozłowski et al. 2012), one at $(76.7367, -71.4573)$ as a young stellar object (Whitney et al. 2008), and one at $(87.9522, -74.5234)$ as a clear supernova, associated with a cat- aloged host galaxy (LEDA database, Paturel et al. 1995).
<table><caption>**Table 3.** Characteristics of the candidates, with fitted microlensing parameters.</caption><thead><tr><th>Candidate ID</th><th>lm0690k17399</th><th>lm0073m17729</th></tr></thead><tbody><tr><td>RA (J2000)</td><td>78.3674</td><td>93.3520</td></tr><tr><td>Dec (J2000)</td><td>-71.9644</td><td>-69.5183</td></tr><tr><td>R</td><td>20.89 ± 0.04</td><td>20.89 ± 0.02</td></tr><tr><td>B</td><td>21.28 ± 0.06</td><td>21.48 ± 0.01</td></tr><tr><td>u<sub>0</sub></td><td>0.194 ± 0.010</td><td>0.413 ± 0.007</td></tr><tr><td>t<sub>E</sub> (day)</td><td>106.3<sup>+9.6</sup><sub>-8.0</sub></td><td>183.1<sup>+8.3</sup><sub>-7.8</sub></td></tr><tr><td>t<sub>0</sub> (MJD)</td><td>51129.7<sup>+1.6</sup><sub>-1.9</sub></td><td>51567.2 ± 2.0</td></tr><tr><td>Δχ<sup>2</sup></td><td>106.16</td><td>85.02</td></tr></tbody></table>

After elimination of these objects, only two candidates remained. Both were observed by EROS-2 but are not within any MACHO field. Their characteristics and light curves are shown in Table 3 and Figs. 2 and 3.

Candidate lm0690k17399 (Fig. 2) is one of the faintest objects (top 2%) in the EROS-2 catalog, and ≈20% of its measurements has negative flux. The event itself is chromatic and has an asymmetric shape with a factor of flux increase >5 (probably even larger when blend is taken into account), characteristic of a type II-L or II-P supernova (Valenti et al. 2015). Differential photometry indicates that the object's position during the event was offset from the position outside the event by ≈0.6", again being consistent with being a supernova in a background galaxy.

Candidate lm0073m17729 (Fig. 3) shows hints of variability outside the main event, and the post-event baseline in the blue passband appears to be ≈0.4 mag fainter than the pre-event baseline.

The light curves of these candidates bear little resemblance to the light curves of actual microlensing effects. However, in our analysis, we cannot formally exclude these candidates without additional data or stricter selection. Therefore, we have chosen to consider them, but only to establish upper bounds, and not to derive a microlensing optical depth due to halo compact objects. We show in Sect. 5 that considering (or not) these two candidates has no consequence on the upper limit of the contribution of high-mass compact objects (heavier than 20 $M_\odot$) to the halo.

## 5. Efficiency and expected detection rate

For $N_{\text{objects}}$ objects observed over a time $T_{\text{obs}}$, the expected number of detected events is

$$N_{\text{exp}} = N_{\text{objects}} T_{\text{obs}} \int dt_{\text{E}} \epsilon(t_{\text{E}}) \frac{d\Gamma}{dt_{\text{E}}} \quad (9)$$

where $d\Gamma/dt_E$ is the rate per $t_E$ interval and $\epsilon(t_E)$ is the detection efficiency, corresponding to the fraction of microlensing events accepted by the selection cuts, relative to the events for which $t_0$ is within the observation time $T_{\text{obs}}$ and $u_0 < 1$. Because of the “blending” effects, there is a nontrivial relation between the number of objects and the number of stars, so Eq. (9) effectively defines $\epsilon(t_E)$. As explained in the following subsections, the number $N_{\text{exp}}$ depends on the details of the observations, on the blending of stars, and on the halo model.

### 5.1. Simulations and efficiency

The efficiency, $\epsilon(t_E)$ was determined by superimposing microlensing events on the light curves of a representative random subsample of the observed objects, and then subjecting these new simulated light curves to the standard analysis procedure. Two types of simulations were performed. The first assumed that each observed object corresponds to one (and only one) star and that variations of the associated light curves can be induced only by microlensing of that one star. The second, more realistic simulations, took into account “blending” by assigning a set of stars to each object in a way that is consistent with the observed density of LMC stars in HST images, as described in the next subsection.

We generated microlensing events by modifying the light curves of a random subsample of objects as follows. First, for each object values of $t_E$, $u_0$, $t_0$, and parallax parameters were drawn randomly from the appropriate distributions. In particular, $u_0$ was drawn uniformly in the interval $0 < u_0 < 1.5$. The time of event maximum, $t_0$, was generated uniformly over a time period that extended sufficiently beyond the first and last observations, since photometric fluctuations could make such events appear to be within the observing period. The $t_E$ distribution was derived from the mass, distance, and velocity distributions of the lenses. To generate the parallax, it was also necessary to take into account the (random) orientation of the transverse velocity vector of the lens. For each of the $N_{\text{star}}$ stars associated with the object, modified object light curves were constructed by modifying the star's flux according to the theoretical microlensing light curve, and then modifying the object light curves in the appropriate manner, taking into account the star's contribution to the object. The uncertainties associated with the flux measurements in the object light curves were modified to account for the increased flux. Each of the $N_{\text{star}}$ set of light curves (or event) was then subjected to the selection process described in the previous section. The efficiency $\epsilon(t_E, u_0)$ is then the number of simulated events passing the cuts divided by the number of objects. Because in the simulation with blending there is more than one star per object, $\epsilon(t_E, u_0)$ could be greater than unity, essentially near $u_0 = 0$.

The calculation of the expected number of events requires $\epsilon(t_E)$, defined as the number of simulated events passing the cuts divided by the number of objects whose event was simulated with $u_0 < 1$. Figure 4 shows the resulting efficiency on simulations with and without blending for the analysis described in Sect. 4, and also by imposing the additional constraint $t_E > 200$ days. The effects of blending are similar to the modest effects seen in earlier studies. Blending increases the true number of stars subject to lensing, but it also affects the shape of the light curve by decreasing the apparent duration $t_E$ and increasing the apparent impact parameter $u_0$. The decrease in effective $t_E$ lowers the efficiency for events at low $t_E$ but increases it at high $t_E$. This change in efficiency is seen in Fig. 4.

### 5.2. Blending

The impact of blending was quantified in the simulation as a multiplicative factor of the number of amplifiable objects. We statistically assigned to each object a list of HST contributing stars with their positions, and calculated the contribution of each star, similarly to Tisserand et al. (2007) and Wyrzykowski et al. (2011), who relied on artificial star additions. In dense fields near the LMC bar, we find that, on average, 2.07 HST stars contribute at least 10% of the flux of cataloged objects (1.58 for objects
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Light curves (magnitude vs. time) of candidate lm0690k17399 in the EROS blue passband (upper panel), and the EROS red passband (lower panel). The black solid lines show the best no-blend microlensing fit. No data is available in the MACHO catalog. The chromaticity, asymmetric shape, and large flux variation indicate that this is most likely a Type II-L or II-P supernova rather than a microlensing event.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Light curves (magnitude vs. time) of candidate lm0073m17729 in the EROS blue passband (upper panel), and the EROS red passband (lower panel). The black solid lines show the best no-blend microlensing fit. No data is available in the MACHO catalog. The hints of variability outside the main bump and the asymmetry, especially visible in the change of the blue baseline after the maximum, make it unlikely that this event is a genuine microlensing event.</figcaption></figure>

brighter than *I* = 19.5). In the sparsest fields, this number drops to 1.46 (1.1 for bright objects). Overall, we find that 90% of the events passing all cuts are due to microlensing of the brightest star associated with an object.

This procedure, based on HST image analysis, is sufficient as long as the spatial distribution of the blend constituents is random, as appears to be the case in the LMC HST images. However, a statistical analysis of binary systems closer than 600 pc extracted from the Gaia database (Blaineau 2021) shows us that such systems located at the distance of the LMC (49.5 kpc) have an increased chance of having a star whose neighbor is not resolved by HST. We therefore decided to examine this particular regime of blend further. In Blaineau (2021) and in a forthcoming paper, we show that the components of most of the LMC binaries are too close together for them to experience very different magnifications through a heavy lens (which almost always have $r_E > 50$ AU when $M_L > 100 M_\odot$). We indeed find that less than 7% of the LMC objects are binaries unresolved by the HST with a projected orbital distance >50 AU. Because of this, we neglect the impact of binarity in our efficiency calculation. Figure 4 shows the detection efficiency for the analysis just described (labeled $t_E > 100$ days), as well as the efficiency for the same analysis but with a stricter requirement on the event duration $t_E > 200$ days.

### 5.3. Expected number of detected events

The top panel of Fig. 5 shows the number of detected microlensing events $N_{\exp}^{(h)}(M)$ expected from lenses in the standard isothermal spherical halo (S-model) with $\tau_{LMC} \sim 4.7 \times 10^{-7}$, assumed to all have the same mass $M$. This number is calculated as a function of $M$ from Eq. (9), using the efficiency shown in Fig. 4 for the simulation with blending. For the number of source objects $N_{\text{objects}}$, we subtract from the $36.8 \times 10^6$ objects in our catalog the contamination by Galactic stars, which we estimated to be less than 5% by counting stars in the Gaia catalog of fields located at the same Galactic latitude, but away from the LMC (Gaia Collaboration 2021). Finally, as in Tisserand et al. (2007) and Wyrzykowski et al. (2011), we consider that 10% of microlensing events may escape detection due to lens binarity (Mróz et al. 2019), to conservatively infer the expected number of events $N_{\exp}^{(h)}(M)$ as a function of the deflector’s mass $M$. For a wide mass distribution, the number of expected events from halo lenses can simply be calculated by integrating the $N_{\exp}^{(h)}(M)$ curve, weighted by the deflector mass distribution function.

The expected number of detected events (Eq. (9)) varies with the adopted halo model through the optical depth, $\tau_{LMC}$, and $t_E$ distributions. The range of plausible models was discussed in Alcock et al. (1996) and we refer to their Table 2 for model
<figure><img src="image_5.png" /><figcaption>Fig. 4. Detection efficiency: Solid (dashed) lines show efficiencies, taking (or not) into account blending effects. The gray lines show efficiencies of the analysis adding the constraint $t_E > 200$ days. The gray histograms (not normalized) are the expected $t_E$ distributions for a halo made of 10, 100, and 1000 $M_{\odot}$ compact objects (from left to right).</figcaption></figure>

details. Isothermal models consistent with a flat rotation curve (e.g., model S of Alcock et al. 1996), all yield values of $N_{\text{exp}}^{(h)}$ similar to that of our baseline model, while $N_{\text{exp}}^{(h)}$ is significantly reduced for models with a falling rotation curve, as might be suggested by some Milky Way rotation velocity measurements (e.g., Bhattacharjee et al. 2014). For model C of Alcock et al. (1996) with $\tau_{\text{LMC}} \approx 3.0 \times 10^{-7}$, and $\langle t_E \rangle \approx 120$ day $\sqrt{M/M_{\odot}}$, $N_{\text{exp}}^{(h)}$ is decreased by a factor $\approx 2$ for $M \approx 10 M_{\odot}$ and the maximum mass for which one would expect three events is reduced by a factor $\approx 3$.

## 6. New limits on compact objects in the halo

We took into account the contribution of the Galactic disk lenses and the LMC self-lensing by adapting the results of Calchi Novati & Mancini (2011) to our efficiency estimate. We therefore expect $N_{\text{exp}}^{(s)} \sim 0.64$ disk and self-lensing events with $t_E > 100$ days to be compared with two observed candidates. In the same way as Calchi Novati & Mancini (2011), we performed Bayesian inference, assuming a weakly constraining uniform prior for the optical depth of the halo compact objects $0 \le \tau_{\text{LMC}} \le 4.7 \times 10^{-6}$ (i.e., up to ten times the expected maximum optical depth). For $N_{\text{obs}}$ observed events with durations $t_i$, $i = 1... N_{\text{obs}}$, the likelihood for the fractional optical depth $f = \tau_{\text{LMC}}/4.7 \times 10^{-7}$ is defined for a given lens mass $M$ as:

$$ \mathcal{L}(f; M|\text{obs}) = \frac{e^{-N_{\text{exp}}}}{N_{\text{obs}}} \prod_{i=1}^{N_{\text{obs}}} \left[ \frac{d\Gamma_{\text{exp}}^{(s)}}{dt_{\text{E}}}(t_i) + f \times \frac{d\Gamma_{\text{exp}}^{(h)}(M)}{dt_{\text{E}}}(t_i) \right] \quad (10) \\ \int d\Gamma_{\text{exp}}^{(s)} = N_{\text{exp}}^{(s)}, \quad \int d\Gamma_{\text{exp}}^{(h)} = N_{\text{exp}}^{(h)} \\ N_{\text{exp}} = N_{\text{exp}}^{(s)} + f \times N_{\text{exp}}^{(h)} $$

where $d\Gamma_{\text{exp}}^{(s)}/dt_{\text{E}}$ and $d\Gamma_{\text{exp}}^{(h)}(M)/dt_{\text{E}}$ are the differential microlensing event rates due to the Milky Way disk and the LMC self-lensing ($\Gamma_{\text{exp}}^{(s)}$), and due to a halo ($\Gamma_{\text{exp}}^{(h)}$) completely made of compact objects with identical mass $M$; $N_{\text{exp}}^{(s)}$ and $N_{\text{exp}}^{(h)}$ are the corresponding numbers of expected events. The probability of observing the two events (with $t_E = 106$ days and $t_E = 183$ days), for a given halo fraction $f$ is then:

$$ p(\text{obs}|f) = \frac{\mathcal{L}(f; M|\text{obs})}{\int_{0}^{10} \mathcal{L}(f; M|\text{obs}) \, df} \quad (11) $$

<figure><img src="image_6.png" /><figcaption>Fig. 5. Contribution of compact objects to the halo mass: Top: number of events expected from a halo S-model entirely composed of compact objects of mass $M$: blue (green) line, from source objects monitored only by MACHO (EROS-2); orange line, from source objects monitored by both surveys; full red line shows the total; black line shows the total adding the constraint $t_E > 200$ days in the analysis. Bottom: 95% CL upper limits on the fraction of the halo mass in the form of compact objects $f = \tau_{\text{LMC}}/4.7 \times 10^{-7}$. Limits obtained in this analysis are shown in red, and in black if we require that $t_E > 200$ days. The gray curves correspond to the latest limits published by MACHO, EROS-2 and OGLE-III.</figcaption></figure>

The halo fraction excluded with a 95% confidence level (CL) through this Bayesian analysis, as a function of deflector mass $f(M)$, is shown in the bottom panel of Fig. 5 as the red curve. Adopting $N_{\text{exp}}^{(s)} = 0$ (i.e., ignoring the disk and self-lensing contribution) would have a negligible impact on the exclusion curve, as the expected signal from the halo is dominant, as shown in the upper part of the figure.

We tested the robustness of this result by changing the $t_E > 100$ days requirement to the stricter $t_E > 200$ days cut. In doing so $N_{\text{exp}}^{(s)}$ becomes negligible (<0.05 event; Calchi Novati & Mancini 2011), and we are left with a zero-event in the data. The resulting exclusion limit, using Poisson statistics (shown as the black curve in Fig. 5), is weaker at the lower mass end, but unchanged on the high-mass side. We also show as gray curves the previously published excluded halo fractions from MACHO (Alcock et al. 2001), EROS-2 (Tisserand et al. 2007), and OGLE-III (Wyrzykowski et al. 2011).

To test halos with mass distribution extending below $1 M_{\odot}$, we can use the results of EROS-1 (Renault et al. 1997) and EROS-2 reported in Tisserand et al. (2007). They expected $\ge 40$ events over a wide range $10^{-6} < M < 1 M_{\odot}$ peaking at 120 events for $10^{-2} M_{\odot}$, but no events were seen. Combining the results of the analysis reported here with those of Tisserand et al. (2007) we would expect $\approx 40$ events for a flat distribution of $\log M$ extending over the range $10^{-7} < M < 10^3 M_{\odot}$. More generally, the combination of the limits obtained by EROS-1, EROS-2, and the present analysis shows that deflectors with any mass distribution in the range $10^{-6} < M < 10^2 M_{\odot}$ cannot contribute more than $f \sim 15\%$ at a 95% CL to the halo mass, assuming a standard spherical halo.
Constraints on the existence of PBHs are reviewed in Carr et al. (2021). Limits at the level of 10% on the halo fraction have been established for most of the range $10^{-10} < M < 10^2 M_\odot$. At higher masses, $10^3 < M < 10^{15} M_\odot$, a variety of techniques have provided stronger limits, as low as $10^{-4}$. At lower masses, PBHs with $M < 10^{-17} M_\odot$ are ruled out because of their rapid evaporation via Hawking radiation.

# 7. Discussion, conclusions, and perspectives

In this search, we did not consider the Small Magellanic Cloud (SMC) data because of their low statistical impact (5.2 million light curves). We attribute the significant improvement in sensitivity over previously reported results to several factors:

- Past analysis focused on finding shorter duration events and always imposed strong constraints on the maximum value of $t_E$ to allow the stability of the object outside the microlensing effect to be checked.

- The addition of the expected signals in EROS-2 and MACHO, taking care to count only once the objects followed simultaneously, already mechanically improves the sensitivity (Moniez 2010). It is important to also note that we were able to use more archived data (7.7 years) from MACHO than the data analyzed by this collaboration (5.7 years).

- A significant part of the improvement comes from combining the light curves (extended to 10.6 years) for the 14.1 million objects jointly observed by EROS-2 and MACHO. These objects, which represent 38% of the catalog, account for $4.0/5.55 = 72\%$ of the expected events if all the lenses have a mass of $1000 M_\odot$ (Fig. 5-Top). This indicates that the detection efficiency on the jointly measured curves is increased by a factor of ~4 compared to objects measured by only one survey.

The combined data set from EROS-2, MACHO, and OGLE surveys now spans 30 consecutive years, since OGLE-III observed from September 2001 to March 2009, and OGLE-IV started in March 2010. It therefore still has potential, which has been studied in Mirhosseini & Moniez (2018). A group that would analyze all these surveys on the search principle applied here of simultaneous fitting of all light curves in all available passbands (up to six) could certainly further improve the limits obtained here, or result in the detection of some very long-duration events. Mirhosseini & Moniez (2018) predicted that, if the halo was completely made up of $1000 M_\odot$ objects, there should be at least three microlensing events toward the LMC in the combined data from all previous microlensing search programs (EROS-2, MACHO, and OGLE through 2018). In fact, these predictions were based on a conservative extrapolation of the detection efficiency, and the combination of only EROS-2 and MACHO data, as analyzed here, allows us to conclude that $1000 M_\odot$ objects filling the halo should have produced ~6 detectable events. The limit of this catalog fusion technique will be reached when all existing catalogs (EROS-2+MACHO+OGLE completed) have been used. If an analysis of a combined catalog with 20 million objects monitored for 30 years has a mean microlensing detection efficiency of 0.20, and assuming no event is detected, then the (ultimate) limit on the mass of the compact halo objects as unique component of the standard halo could be pushed to $\ge 3000 M_\odot$. The Rubin/LSST survey (Ivezić et al. 2019), which will further extend the monitoring duration of the LMC objects with better photometric precision (also allowing for a better detection efficiency), should further enhance the value of the historical surveys whose data could also be aggregated, or at least used for verification purposes in case of detection.

In conclusion, it seems that BHs with masses up to a thousand solar masses, similar to those observed by LIGO and Virgo as binary BH mergers, do not make up a major fraction of the Milky Way dark matter, at least if assumed to be distributed as a standard spherical halo. Such BHs are more likely to be found in structures following the visible mass distribution, and could be searched for through microlensing toward the Galactic bulge and spiral arms, in the long $t_E$ tail of event duration distribution, extending previous searches and analyses (Hamadache et al. 2006; Moniez et al. 2017; Mróz et al. 2020).

**Acknowledgements.** We thank Lucasz Wyrzykowski and Christopher Stubbs for providing useful information. This work was supported by the Paris Ile-de-France Region. This paper uses public domain data obtained by the MACHO Project, jointly funded by the US Department of Energy through the University of California, Lawrence Livermore National Laboratory under contract No. W-7405-Eng-48, by the National Science Foundation through the Center for Particle Astrophysics of the University of California under cooperative agreement AST-8809616, and by the Mount Stromlo and Siding Spring Observatory, part of the Australian National University. This article is dedicated to the memory of our EROS collaborators Johannes Andersen, Pierre Bareyre, Florian Bauer, Sergio Char, Eric Maurice, Alain Milsztajn, Luciano Moscoso, and Cécile Renault.

# References

Abbott, B., Abbott, R., Abbott, T., et al. 2016a, Phys. Rev. Lett., 116, 241103

Abbott, B., Abbott, R., Abbott, T., et al. 2016b, Phys. Rev. Lett., 116, 061102

Alcock, C., Allsman, R. A., Axelrod, T. S., et al. 1996, ApJ, 461, 84

Alcock, C., Allsman, R. A., Alves, D. R., et al. 1999, PASP, 111, 1539

Alcock, C., Allsman, R. A., Alves, D. R., et al. 2001, ApJ, 550, L169

Bhattacharjee, P., Chaudhury, S., & Kundu, S. 2014, ApJ, 785, 63

Bird, S., Cholis, I., Muñoz, J. B., et al. 2016, Phys. Rev. Lett., 116, 201301

Blaineau, T. 2021, PhD Thesis, Université Paris-Saclay, France

Blaineau, T., & Moniez, M. 2020, A&A, 636, L9

Brunthaler, A., Reid, M. J., Menten, K. M., et al. 2011, Astron. Nachr., 332, 461

Calchi Novati, S., & Mancini, L. 2011, MNRAS, 416, 1292

Carr, B., Kohri, K., Sendouda, Y., & Yokoyama, J. 2021, Rep. Progr. Phys., 84, 116902

Gaia Collaboration (Brown, A. G. A.) 2021, A&A, 649, A1

Green, A. M., & Kavanagh, B. J. 2021, J. Phys. G Nucl. Phys., 48

Griest, K. 1991, ApJ, 366, 412

Hamadache, C., Le Guillou, L., Tisserand, P., et al. 2006, A&A, 454, 185

Hubert, A. M. 2007, in Active OB-Stars: Laboratories for Stellare and Circumstellar Physics, eds. A. T. Okazaki, S. P. Owocki, & S. Stefl, ASP Conf. Ser., 361, 27

Ivezić, Ž., Kahn, S. M., Tyson, J. A., et al. 2019, ApJ, 873, 111

Kallivayalil, N., van der Marel, R. P., Besla, G., Anderson, J., & Alcock, C. 2013, ApJ, 764, 161

Kim, D.-W., Protopapas, P., Trichas, M., et al. 2012, ApJ, 747, 107

Kozłowski, S., Kochanek, C. S., Jacyszyn, A. M., et al. 2012, ApJ, 746, 27

Mirhosseini, A., & Moniez, M. 2018, A&A, 618, L4

Moniez, M. 2010, Gen. Relativ. Grav., 42, 2047

Moniez, M., Sajadian, S., Karami, M., Rahvar, S., & Ansari, R. 2017, A&A, 604, A124

Mróz, P., Udalski, A., Skowron, J., et al. 2019, ApJS, 244, 29

Mróz, P., Udalski, A., Szymański, M. K., et al. 2020, ApJS, 249, 16

Paczynski, B. 1986, ApJ, 304, 1

Paturel, G., Vauglin, I., Andernach, H., et al. 1995, A&AS, 112, 525

Pietrzyński, G., Graczyk, D., Gallenne, A., et al. 2019, Nature, 567, 200

Rahvar, S. 2015, Int. J. Mod. Phys. D, 24, 1530020

Renault, C., Afonso, C., Aubourg, E., et al. 1997, A&A, 324, L69

Sasaki, M., Suyama, T., Tanaka, T., & Yokoyama, S. 2016, Phys. Rev. Lett., 117, 061101

Schneider, P., Kochanek, C., & Wambsganss, J. 2006, Gravitational Lensing: Strong, Weak and Micro (Berlin: Springer)

Suntzeff, N. B., Heathcote, S., Weller, W. G., et al. 1988, Nature, 334, 135

Tisserand, P. 2004, PhD Thesis, Université Nice-Sophia Antipolis, France

Tisserand, P., Le Guillou, L., Afonso, C., et al. 2007, A&A, 469, 387

Valenti, S., Sand, D., Stritzinger, M., et al. 2015, MNRAS, 448, 2608

Whitney, B. A., Sewilo, M., Indebetouw, R., et al. 2008, AJ, 136, 18

Wyrzykowski, L., Skowron, J., Kozłowski, S., et al. 2011, MNRAS, 416, 2949