New limits from microlensing on Galactic Black Holes in the mass range $10M_{\odot}<M<1000M_{\odot}$

Introduction

Observation of gravitational waves due to coalescence of massive objects ([1 b]) have demonstrated the existence of merging black holes heavier than $10M_{\odot}$, and has renewed interest in Black Holes as dark matter, especially primordial Black Holes (PBH) ([7]; [14]; [32]). Microlensing surveys toward the LMC, that allow to probe the content of the Galactic halo in massive compact objects, have shown that objects lighter than $10M_{\odot}$ do not significantly contribute to the hidden mass of the Galactic spherical halo of our galaxy ([36]; [3]; [39]). Since the typical duration of microlensing events increases with the lens mass, the detection of heavier objects such as those responsible for gravitational wave emissions needs 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. To explore the dark matter halo beyond this limit by searching for events with longer timescales, we have combined the databases of EROS-2 and MACHO, thus starting the program described in Mihosseeini & Moniez (2018).

In Sect. 2 we describe the microlensing effect, focusing on the LMC searches. Section 3 introduces the EROS-2 and MACHO surveys and their light curve datasets, summarizes how we associated objects between the two catalogs, and presents 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 fitting the observed light curves with theoretical microlensing curves. In Sect. 5 we explain the calculation of the efficiency for detection of microlensing events. We quantify the effect of blending on the detection efficiency by using 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 achieved with respect to previous results, and propose some perspectives for further data combinations.

Send offprint requests to: M. Moniez, e-mail: [email protected]

Microlensing toward 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 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 Large Magellanic Cloud (LMC) stars as a technique to search for massive compact objects in the Galactic halo was first described in Paczynski [27]. Reviews of the formalism can be found in Schneider et al. [33] and Rahvar [30]. 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\,\mathrm{kpc}$ [29], the total magnification of the source luminosity at a given time $t$ is the sum of the contributions of two images, given by:

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

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 r_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:

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 which 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 time scales.

• The magnification is a known function of time, depending on only 3 parameters ($u_0$, $t_0$, $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.

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 order $v_{\mathrm{rot}}^2/c^2$, where $v_{\mathrm{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 [15] (hereafter called S-model), with the mass density distribution

where $R_0=8.5\ \mathrm{kpc}$ is the Galactocentric distance to the Sun, $R_c=5\ \mathrm{kpc}$ the halo “core radius” and $r$ the Galactocentric radius. We use this model with the most recent values of the LMC distance $D_{LMC}=49.5\ \mathrm{kpc}$ [29].

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

with $v_0=120\ \mathrm{km\,s^{-1}}$; the velocity of the Sun is $\mathbf{v}_\odot=(11.1,251,7.3)\ \mathrm{km\,s^{-1}}$ [10] in galactocentric Cartesian coordinates (X pointing from the Sun to the Galactic center, Z pointing north), and the proper velocity of the LMC is $\mathbf{v}_{LMC}=(-57,-226,221)\ \mathrm{km\,s^{-1}}$ [19]. This model gives an optical depth to the LMC $\tau_{LMC}\sim4.7\times10^{-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$ [15], where $\langle t_E\rangle$ is the mean $t_E$ of all events¹. For lenses of mass $M_L$, the mean event duration for this benchmark halo model is $\langle t_E\rangle\sim63\ \mathrm{day}\sqrt{M_L/M_\odot}$.

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 of a 1 meter (F/5) diameter telescope, equipped with a dichroic beam splitter and two cameras, each with 8 2Kx2K CCD’s, covering 1. deg$^2$. MACHO used a slightly larger telescope (1.27m, F/3.9), equipped with two cameras, each with 4 2Kx2K CCD’s, covering a smaller field of view (0.5 deg$^2$).

	EROS-2	MACHO
telescope	1 m	1.27 m
pixel size	0.62″	0.63″
blue passband	[420, 720] nm	[450, 590] nm
red passband	[620, 920] nm	[590, 780] nm
median image quality	2″	2.1″

EROS-2 surveyed 88 fields (1. deg$^2$) toward LMC, and MACHO surveyed 82 (0.5 deg$^2$) LMC fields. MACHO adopted the same exposure time (300s) for all fields, whereas EROS-2 has chosen to adapt its exposure times according to the surface ends_mid=1 brightness of the fields, varying from 180s (for the central fields) to 900s (for the external fields). From 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² (Alcock et al. [4]). The catalog contains 22.3 $\times10^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. [3]).

The EROS-2 catalog for LMC was produced by Tisserand [35] for the final EROS LMC publication (Tisserand et al. [36]). 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 EDR3 astrometry (Gaia Collaboration et al. [13]), 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 $\sim1''$), 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 (Figure 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 $\times10^6$ objects cross-matched in the EROS-2 and MACHO catalogs; they benefit from a total of 10.6 years of luminosity measurements, of which 3.8 years overlap, during which the surveys used 4 different passbands. The complementary set, consisting of objects reconstructed in EROS-2 only or in MACHO only, comprising 22.7 $\times10^6$ objects monitored by a single survey, is also included in our analysis, although the objects are monitored in only 2 passbands for shorter times.

	EROS-2 only	MACHO only	common
Dates (month/yr)	7/96–2/03	7/92–1/00	7/92–2/03
$T_{\mathrm{obs}}$ (year)	6.7	7.7	10.6
$N_{\mathrm{objects}}$ ($\times10^6$)	15.8	6.9$^{(a)}$	14.1
central fields deg$^2$
stars/arcmin$^2$	$\sim70$	$\sim100$	$\sim70$
mag. lim. $V_{\mathrm{Cousins}}$	$\sim20.5$	$\sim20.5$	$\sim20.5$
# measurements B	500	1400	1900
# measurements R	600	1550	2150
outer fields deg$^2$
stars/arcmin$^2$	$\sim30$	$\sim20$	$\sim20$
mag. lim. $V_{\mathrm{Cousins}}$	$\sim22.5$	$\sim21.5$	$\sim21.5$
# measurements B	250	200	450
# measurements R	300	250	550

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 more than 5 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 [8]).

We further discarded measurements that deviated by more than $5\sigma$ from the median flux in a sliding window of 5 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 $\sim3\%$ ($\sim1\%$) of the measurements from the MACHO (EROS-2) data.

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 in general underestimated by an average factor of $\sim0.73$ while those of EROS-2 were incorrectly estimated by a factor varying from $\sim0.75$ (for stars brighter than $I = 19$) to $\sim1.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 neighbors:

Here, $\phi(t_i)$ is the flux measured at time $t_i$, $\phi_i^{\mathrm{int}}$ is the flux interpolated from measurements $i-1$ and $i+1$ at $t_i$, and $\sigma_{\mathrm{int}}$ is the uncertainty on $(\phi(t_i)-\phi_i^{\mathrm{int}})$ deduced from the cataloged uncertainties. The sum is over points with a precision better than 0.55 mag. Renormalizing squared uncertainties by dividing them by $X$ generally made the $\chi^2$ per DOF near unity for fits assuming slow flux variations over time. As the correction factors to be applied to photometric uncertainties vary with magnitudes, this renormalization procedure is not effective for light curves with large variations, specially for low fluxes with large uncertainties ($> 0.55$ mag) in EROS-2 data.

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 criteria. More precisely, for each object we performed a simultaneous 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$, $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 section 5, in order to determine the analysis efficiency.

The impact of the Earth’s rotation (parallax) around the Sun on the detection efficiency of multi-year events was discussed in [9]. The effect on a tolerant search algorithm such as the one we describe below was found to be negligible. 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.

Selection of candidates

The light curves of the $36.8 \times10^6$ cataloged objects underwent the following selection process to identify long-duration microlensing candidates. This process was tuned with the simulation described in Section 5.

• We eliminate objects with fewer than 200 total measurement points, all passbands included, and with fewer than 50 measurements in at least both EROS-2 or both MACHO passbands. This selection leaves $36 \times10^6$ objects.

• We require 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 require

where the function

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 accept 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 Section 5).

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

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

• We require that $t_E$ be in the range $100\ days < t_E < (t_{\max}-t_{\min})/2$. This criterion rejects most supernovae and short-lived fluctuations, and ensures that the fitted variations are sufficiently within the observing period. 34 light curves remain 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 $1M_\odot$ to 62% for $1000M_\odot$).

• We eliminate so-called “blue bumpers” by making a stricter selection on events in their region of the color-magnitude diagram. 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 [17], can extend over years. We reject these artifacts by requiring that the fitted value of $u_0$ be less than 0.9 for objects in the following CMD domain: $I_{\mathrm{Cousins}} < 19$ and $(V_{\mathrm{Johnson}}-I_{\mathrm{Cousins}}) < 0.5$.

28 candidate events remain after this selection.

Remaining candidates: rejection of known artefacts

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 eliminate three types of objects.

• Objects outside the region of the CMD containing 99% of the stars: We reject these outliers to discard in particular the rarest objects, likely to show variability. 24 objects remain after this CMD-based rejection, at the cost of 1% of the expected events.

• Echoes from SN1987A ([34]): We discard 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 remain, while almost 100% of the simulated events selected so far are retained.

• Variable objects identified in external catalogs: At the end of our selection process, we reject one object at $(\alpha,\delta)_{J2000}=(80.6256,-71.7500)$, identified on external catalogs as a

QSO [20] (Kozłowski et al. 2012), one at $(76.7367, -71.4573)$ as a Young Stellar Object [38], and one at $(87.9522, -74.5234)$ as a clear supernova, associated with a cataloged host-galaxy (LEDA database, Paturel et al. [28]).

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

	candidate ID	lm0690k17399	lm0073m17729
RA (J2000)		78.3674	93.3520
DEC (J2000)		-71.9644	-69.5183
R		$20.89 \pm0.04$	$20.89 \pm0.02$
B		$21.28 \pm0.06$	$21.48 \pm0.01$
$u_0$		$0.194 \pm0.010$	$0.413 \pm0.007$
$t_E$ (day)		$106.3^{+9.6}_{-8.0}$	$183.1^{+8.3}_{-7.8}$
$t_0$ (MJD)		$51129.7^{+1.6}_{-1.9}$	$51567.2 \pm2.0$
$\Delta\chi^2$		$106.16$	$85.02$

Candidate $lm0690k17399$ is one of the 2% faintest objects in the EROS-2 catalog, and $\approx20\%$ of its measurements had 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 supernovae [37]. Differential photometry indicates that the object’s position during the event was offset from the position outside the event by $\approx0.6''$, again being consistent with being a supernova in a background galaxy.

Candidate $lm0073m17729$ shows hints of variability outside the main event, and the post-event baseline in the blue appears to be $\approx0.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 will 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 $20M_\odot$) to the halo.

Efficiency and expected detection rate

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

where $\mathrm{d}\Gamma/\mathrm{d}t_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_{\mathrm{obs}}$ and $u_0 < 1$. Because of the “blending” effects, there is a non-trivial 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_{\mathrm{exp}}$ depends on the details of the observations, on the blending of stars, and on the halo model.

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, takes 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 generate microlensing events by modifying the light curves of a random sub-sample of objects as follows. First, for each object values of $t_E$, $u_0$, $t_0$ and parallax parameters are drawn randomly from the appropriate distributions. In particular, $u_0$ is drawn uniformly in the interval $0 < u_0 < 1.5$. The time of event maximum, $t_0$, is generated randomly over a time period that extends sufficiently beyond the first and last observations, since photometric fluctuations may make such events appear to be within the observing period. The $t_E$ distribution is derived from the mass, distance and velocity distributions of the lenses. The parallax depends in addition on the random orientation of the transverse velocity vector of the lens. For each of the $N_{\mathrm{star}}$ stars associated with the object, modified object light curves are 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 are modified to account for the increased flux. Each of the $N_{\mathrm{star}}$ set of light curves (or event) is 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)$ can 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 Figure 4.

Blending

The impact of blending was quantified in the simulation as a multiplicative factor of the number of amplifiable objects. We statistically assign to each object a list of HST contributing stars with their positions, and calculate the contribution of each star similarly to Tisserand et al. [36] and Wyrzykowski et al. [39], which 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 brighter than I=19.5). In the sparse 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 the 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; but a statistical analysis of binary systems closer than 600 pc extracted from the Gaia database (Blaineau [8]) showed us that such systems, located at the distance of the LMC (49.5 kpc), increase the chances that a star has a neighbor not resolved by HST. We therefore decided to examine further this particular regime of blend. In Blaineau [8] 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\,\mathrm{AU}$ when $M_L > 100M_\odot$). We indeed find that less than 7% of the LMC objects should be binaries unresolved by the HST with a projected orbital distance $> 50\,\mathrm{AU}$. Because of this, we ignore the effect in the calculation of the efficiency.

Figure 4 shows the detection efficiency for the analysis just described (labelled $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.

The expected number of detected events

Top panel of Figure 5 shows the number of detected microlensing events $N_{\rm exp}^{(h)}(M)$ expected from lenses in the standard isothermal spherical halo (S-model) with $\tau_{\rm LMC} \sim4.7 \times10^{-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 Figure 4 for the simulation with blending. For the number of source objects $N_{\rm objects}$, we subtract from the $36.8 \times10^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 et al. [13]). Finally, as in Tisserand et al. [36] and Wyrzykowski et al. [39], 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_{\rm 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_{\rm 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_{\rm LMC}$, and $t_E$ distributions. The range of plausible models was discussed in [5] and we refer to their Table 2 for model details. Isothermal models consistent with flat rotation curve, e.g. model S of Alcock et al. [5], all yield values of $N_{\rm exp}^{(h)}$ similar to that of our baseline model, while $N_{\rm 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. [6]). For the model C of Alcock et al. [5] with $\tau_{\rm LMC} \approx3.0 \times10^{-7}$, and $\langle t_E\rangle\approx120\,\mathrm{day}\sqrt{M/M_\odot}$, $N_{\rm exp}^{(h)}$ is decreased by a factor $\approx2$ for $M \approx10M_\odot$ and the maximum mass for which one would expect three events is reduced by a factor $\approx3$.

New limits on compact objects in the halo

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

where $\mathrm{d}\Gamma_{\rm exp}^{(s)}/\mathrm{d}t_E$ and $\mathrm{d}\Gamma_{\rm exp}^{(h)}(M)/\mathrm{d}t_E$ are the differential microlensing event rates due to the Milky Way disk and the LMC self-lensing ($\Gamma_{\rm exp}^{(s)}$), and due to a halo ($\Gamma_{\rm exp}^{(h)}$) completely made of compact objects with identical mass $M$; $N_{\rm exp}^{(s)}$ and $N_{\rm 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:

The halo fraction excluded with 95% CL through this Bayesian analysis, as a function of deflector mass $f(M)$, is shown on the bottom panel of Figure 5 as the red curve. Assuming $N_{\rm exp}^{(s)}=0$ has 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_{\rm exp}^{(s)}$ becomes negligible ($<0.05$ event) (Calchi Novati & Mancini [11]), and we are left with zero event in the data. The resulting exclusion limit, using Poisson statistics, shown as the black curve, is weaker at the lower mass end, but unchanged on the high mass side. We have also shown as gray curves the previously published excluded halo fractions by MACHO (Alcock et al. [3]), EROS-2 (Tisserand et al. [36]), and OGLE-III (Wyrzykowski et al. [39]).

To test halos with mass distribution extending below $1M_\odot$, we can use the results of EROS-1 (Renault et al. [31]) and EROS-2 reported in Tisserand et al. [36]. EROS-2 reported in Tisserand et al. [36] expected $\ge40$ events over a wide range $10^{-6}<M<1M_\odot$, peaking at 120 events for $10^{-2}M_\odot$, whereas no events were seen. Combining the results of the analysis reported here with those of Tisserand et al. [36] we would expect $\sim40$ events for a flat distribution of $\log M$ extending over the range $10^{-7}<M<10^3M_\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^2M_\odot$ can not contribute more than $f\sim15\%$ at 95% CL to the halo mass, assuming a standard spherical halo.

Constraints on the existence of primordial black holes are reviewed in Carr et al. [12]. Limits at the level of 10% on the halo fraction have been established for most of the range $10^{-10}<M<10^2M_\odot$. At higher masses, $10^3<M<10^{15}M_\odot$, a variety of techniques have provided stronger limits, as low of $10^{-4}$. At lower masses, primordial Black Holes with $M<10^{-17}M_\odot$ are ruled out because of their rapid evaporation via Hawking radiation.

Discussion, conclusions, and perspectives

In this search, we did not consider the 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 checking the stability of the object outside the microlensing effect.

• 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 [23]). Note also 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, contribute for $4.0/5.55 = 72\%$ of the expected events if all the lenses have a mass of $1000M_\odot$ (Figure 5-Top). This indicates that the detection efficiency on the jointly measured curves is increased by a factor of $\sim4$ 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 sept. 2001 to march 2009 and OGLE-IV started in march 2010. It therefore still has potential, which has been studied in Mirhosseini & Moniez [22]. 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 6) could certainly further improve the limits obtained here, or result in the detection of some very long duration events. Mirhosseini & Moniez [22] predicted that if the halo was completely made up of $1000M_\odot$ objects there should be at least 3 microlensing events towards LMC in the combined data from all previous microlensing search programs (EROS-2, MACHO, 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 analysed here, allows us to conclude that $1000M_\odot$ objects filling the halo should have produced $\sim6$ 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 millions 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 $\gtrsim3000M_\odot$. The Rubin/LSST survey (Ivezić et al. [18]), which will further extend the monitoring duration of the LMC objects with better photometric precision (also allowing for 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 the ones 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 the previous searches and analysis (Hamadache et al. [16]; Moniez et al. [24]; 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 Miltzstajn, Luciano Moscardi, and Cécile Renault.