© T. Blaineau et al. 2022

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_⊙, 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_⊙ 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_⊙ 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 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:

$$\begin{array}{c}\hfill A(t)=\frac{u{(t)}^2+2}{u(t)\sqrt{u{(t)}^2+4}},\end{array}$$(1)

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

$$\begin{array}{c}\hfill r_{\mathrm{E}}=\sqrt{\frac{4G{M}_{\mathrm{L}}}{c^2}{D}_{\mathrm{LMC}}x(1-x)}\simeq 10.0\mathrm{AU}\times {\left[\frac{M_{\mathrm{L}}}{M_{\odot }}\right]}^{\frac{1}{2}}\frac{{[x(1-x)]}^{\frac{1}{2}}}{0.5}.\end{array}$$(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 r_{\text{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:

$$\begin{array}{c}\hfill u(t)=\sqrt{u_0^2+{(t-t_0)}^2/t_{\mathrm{E}}^2},\end{array}$$(3)

where t_E = r_E/v_T is the Einstein radius crossing time, and u₀ is the minimum distance to the undeflected line of sight at time t₀.

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₀, t₀, 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₀ and of the impact parameters u₀ 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 the order of $v_{\text{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

$$\begin{array}{c}\hfill {\rho }_H(r)=0.0078\frac{R_0^2+R_c^2}{r^2+R_c^2}{M}_{\odot }{\mathrm{pc}}^{-3},\end{array}$$(4)

where R₀ = 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

$$\begin{array}{c}\hfill p(v_{\mathrm{L}})={\left(\frac{1}{2\pi v_0^2}\right)}^{3/2}4\pi v_{\mathrm{L}}^2{\mathrm{e}}^{-v_{\mathrm{L}}^2/2v_0^2},\end{array}$$(5)

with v₀ = 120 km s⁻¹; the velocity of the Sun is v_⊙ = (11.1, 251, 7.3) km s⁻¹ (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) km s⁻¹ (Kallivayalil et al. 2013).

This model gives an optical depth to the LMC τ_LMC ∼ 4.7 × 10⁻⁷ if the halo is completely made of compact objects. The total rate (events per star per unit time) for u₀ < 1 is Γ = (2/π)τ_LMC/⟨t_E⟩ (Griest 1991), where ⟨t_E⟩ is the mean t_E of all events¹. For lenses of mass M_L, the mean event duration for this benchmark halo model is $⟨t_{\mathrm{E}}⟩\sim 63\mathrm{day}\sqrt{M_{\mathrm{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 × 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 × 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 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² (Alcock et al. 1999). The catalog contains 22.3 × 10⁶ 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.

Fig. 1. Distribution of the EROS-2 and MACHO magnitude difference Δm = mBE − mBM as a function of the EROS-2 blue magnitude mBE for the 14.1 million associated objects; the values of mBM 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 Δ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 mBE = 18.8). White dots and bars show the average and standard deviations of Δm for each mBE slice. The color scale is expressed in millions of objects per squared magnitude.

Table 2 summarizes the characteristics of the two surveys and the cross-matched catalog. The common set consists of the 14.1 × 10⁶ 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 four different passbands. The complementary set, consisting of objects reconstructed only in EROS-2 or only in MACHO, comprising 22.7 × 10⁶ 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σ) 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σ 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 χ², 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:

$$\begin{array}{c}\hfill X=\frac{1}{N-2}{\displaystyle \sum _{i=2}^{N-1}{(\varphi (t_i)-{\varphi }_i^{\mathrm{int}})}^2/{\sigma }_{\mathrm{int}}^2.}\end{array}$$(6)

Here, ϕ(t_i) is the flux measured at time t_i, ${\varphi }_i^{\text{int}}$ is the flux interpolated from measurements i − 1 and i + 1 at t_i, and σ_int is the uncertainty on $(\varphi (t_i)-{\varphi }_i^{\text{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 χ² 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 magnitudes, 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 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₀, u₀, 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 determine 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 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.

4.1. Selection of candidates

The light curves of the 36.8 × 10⁶ cataloged objects underwent the following selection process in order to identify long-duration microlensing candidates. This process was tuned with the simulation described in Sect. 5.

We eliminated 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 left 36 × 10⁶ objects.

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

$$\begin{array}{c}\hfill \mathrm{\Delta }{\chi }^2=\frac{{\chi }_{\text{const.}}^2-{\chi }_{\mathrm{ML}}^2}{{\chi }_{\mathrm{ML}}^2/N_{\text{d.o.f.}}}\frac{1}{\sqrt{2N_{\text{d.o.f.}}}}>80,\mathrm{and}\frac{{\chi }_{\mathrm{ML}}^2}{N_{\text{d.o.f.}}}<f(\mathrm{\Delta }{\chi }^2),\end{array}$$(7)

where the function

$$\begin{array}{c}\hfill f(\mathrm{\Delta }{\chi }^2)=1.44+0.26{log}_{10}\left(\frac{\mathrm{\Delta }{\chi }^2}{80}\right)+0.23{log}_{10}^2\left(\frac{\mathrm{\Delta }{\chi }^2}{80}\right)\end{array}$$(8)

is tuned to accept 92.5% of the simulated light curves for each given interval of Δχ², and is in the range 1.4 < f(Δχ²) < 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₀ < 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 Δχ² > −0.1, where Δχ² 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_min + 200 days < t₀ < t_max, where t_min and t_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_max − t_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_⊙ to 62% for 1000 M_⊙).

We eliminated 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 (Hubert et al. 2007), can extend over years. We rejected these artifacts by requiring that the fitted value of u₀ be less than 0.9 for objects in the following CMD domain: I_Cousins < 19 and (V_Johnson − I_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° < α < 84.064° and −69.34° < δ < −69.20°. After this exclusion, 5 objects remained, while almost 100% of the simulated 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 (α, δ)_J2000 = (80.6256, −71.7500), identified 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 cataloged host galaxy (LEDA database, Paturel et al. 1995).

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.

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.

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 geniune microlensing event.

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_⊙) to the halo.

5. Efficiency and expected detection rate

For N_objects objects observed over a time T_obs, the expected number of detected events is

$$\begin{array}{c}\hfill N_{\exp }=N_{\mathrm{objects}}{T}_{\mathrm{obs}}{\displaystyle \int }\mathrm{d}{t}_{\mathrm{E}}ϵ(t_{\mathrm{E}})\frac{\mathrm{d}\mathrm{\Gamma }}{\mathrm{d}{t}_{\mathrm{E}}}\end{array}$$(9)

where dΓ/dt_E is the rate per t_E interval and ϵ(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₀ is within the observation time T_obs and u₀ < 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 ϵ(t_E). As explained in the following subsections, the number N_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, ϵ(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₀, t₀, and parallax parameters were drawn randomly from the appropriate distributions. In particular, u₀ was drawn uniformly in the interval 0 < u₀ < 1.5. The time of event maximum, t₀, 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_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_star set of light curves (or event) was then subjected to the selection process described in the previous section. The efficiency ϵ(t_E, u₀) 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, ϵ(t_E, u₀) could be greater than unity, essentially near u₀ = 0.

The calculation of the expected number of events requires ϵ(t_E), defined as the number of simulated events passing the cuts divided by the number of objects whose event was simulated with u₀ < 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₀. 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.

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 tE > 200 days. The gray histograms (not normalized) are the expected tE distributions for a halo made of 10, 100, and 1000 M⊙ compact objects (from left to right).

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 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_⊙). 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 }^{(\mathrm{h})}(M)$ expected from lenses in the standard isothermal spherical halo (S-model) with τ_LMC ∼ 4.7 × 10⁻⁷, 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_objects, we subtract from the 36.8 × 10⁶ 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 }^{(\mathrm{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 }^{(\mathrm{h})}(M)$ curve, weighted by the deflector mass distribution function.

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 tE > 200 days in the analysis. Bottom: 95% CL upper limits on the fraction of the halo mass in the form of compact objects f = τLMC/4.7 × 10−7. Limits obtained in this analysis are shown in red, and in black if we require that tE > 200 days. The gray curves correspond to the latest limits published by MACHO, EROS-2 and OGLE-III.

The expected number of detected events (Eq. (9)) varies with the adopted halo model through the optical depth, τ_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 details. Isothermal models consistent with a flat rotation curve (e.g., model S of Alcock et al. 1996), all yield values of $N_{\exp }^{(\mathrm{h})}$ similar to that of our baseline model, while $N_{\exp }^{(\mathrm{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 τ_LMC ≈ 3.0 × 10⁻⁷, and $⟨t_{\mathrm{E}}⟩\approx 120\mathrm{day}\sqrt{M/M_{\odot }}$, $N_{\exp }^{(\mathrm{h})}$ is decreased by a factor ≈2 for M ≈ 10 M_⊙ and the maximum mass for which one would expect three events is reduced by a factor ≈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_{\exp }^{(\mathrm{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 ≤ τ_LMC ≤ 4.7 × 10⁻⁶ (i.e., up to ten times the expected maximum optical depth). For N_obs observed events with durations t_i, i = 1…N_obs, the likelihood for the fractional optical depth f = τ_LMC/4.7 × 10⁻⁷ is defined for a given lens mass M as:

$$\begin{array}{cc}& \mathcal{L}(f;M|\mathrm{obs})=\frac{{\mathrm{e}}^{-N_{\exp }}}{N_{\mathrm{obs}}}{\displaystyle \prod _{i=1}^{N_{\mathrm{obs}}}[\frac{\mathrm{d}{\mathrm{\Gamma }}_{\exp }^{(\mathrm{s})}}{\mathrm{d}{t}_{\mathrm{E}}}(t_i)+f\times \frac{\mathrm{d}{\mathrm{\Gamma }}_{\exp }^{(\mathrm{h})}(M)}{\mathrm{d}{t}_{\mathrm{E}}}(t_i)]}\hfill \\ \hfill & {\displaystyle \int }\mathrm{d}{\mathrm{\Gamma }}_{\exp }^{(\mathrm{s})}=N_{\exp }^{(\mathrm{s})},{\displaystyle \int }\mathrm{d}{\mathrm{\Gamma }}_{\exp }^{(\mathrm{h})}=N_{\exp }^{(\mathrm{h})}\hfill \\ \hfill & N_{\exp }=N_{\exp }^{(\mathrm{s})}+f\times N_{\exp }^{(\mathrm{h})}\hfill \end{array}$$(10)

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

$$\begin{array}{c}\hfill p(\mathrm{obs}|f)=\frac{\mathcal{L}(f;M|\mathrm{obs})}{{\int }_0^{10}\mathcal{L}(f;M|\mathrm{obs})\mathrm{d}f}·\end{array}$$(11)

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_{\exp }^{(\mathrm{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_{\exp }^{(\mathrm{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_⊙, we can use the results of EROS-1 (Renault et al. 1997) and EROS-2 reported in Tisserand et al. (2007). They expected ≳40 events over a wide range 10⁻⁶ < M < 1 M_⊙ peaking at 120 events for 10⁻² M_⊙, but no events were seen. Combining the results of the analysis reported here with those of Tisserand et al. (2007) we would expect ≈40 events for a flat distribution of log M extending over the range 10⁻⁷ < M < 10³ M_⊙. 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⁻⁶ < M < 10² M_⊙ cannot contribute more than f ∼ 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⁻¹⁰ < M < 10² M_⊙. At higher masses, 10³ < M < 10¹⁵ M_⊙, a variety of techniques have provided stronger limits, as low as 10⁻⁴. At lower masses, PBHs with M < 10⁻¹⁷ M_⊙ 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_⊙ (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_⊙ 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_⊙ 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 ≳3000 M_⊙. 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).

Acknowledgments

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

All Tables

All Figures

Fig. 1. Distribution of the EROS-2 and MACHO magnitude difference Δm = mBE − mBM as a function of the EROS-2 blue magnitude mBE for the 14.1 million associated objects; the values of mBM 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 Δ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 mBE = 18.8). White dots and bars show the average and standard deviations of Δm for each mBE slice. The color scale is expressed in millions of objects per squared magnitude.

In the text

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.

In the text

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 geniune microlensing event.

In the text

	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 tE > 200 days. The gray histograms (not normalized) are the expected tE distributions for a halo made of 10, 100, and 1000 M⊙ compact objects (from left to right).
In the text

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 tE > 200 days in the analysis. Bottom: 95% CL upper limits on the fraction of the halo mass in the form of compact objects f = τLMC/4.7 × 10−7. Limits obtained in this analysis are shown in red, and in black if we require that tE > 200 days. The gray curves correspond to the latest limits published by MACHO, EROS-2 and OGLE-III.

In the text