© The Authors 2025

1. Introduction

Microlensing is a powerful technique, allowing for the detection of isolated objects regardless of their luminosity. This makes it particularly valuable for studying dark, isolated stellar remnants in our Galaxy. The ability to identify more isolated black holes, neutron stars, and white dwarfs provides a unique opportunity to explore the late stages of stellar evolution and the nature of supernovae.

Despite its potential, microlensing remains a rare phenomenon. Detecting such events requires large-scale sky surveys that monitor millions of stars, such as Gaia (Gaia Collaboration 2016). Operating until January 2025, Gaia scanned the entire sky approximately once per month, leading to the discovery of numerous microlensing events (Wyrzykowski et al. 2023; Hodgkin et al. 2021). These observations have already resulted in the identification of several stellar remnant candidates (e.g. Jabłońska et al. 2022; Kruszyńska et al. 2022, 2024; Howil et al. 2025).

Determining the mass of a microlensing object requires two key parameters: the microlens parallax, π_E, and the angular Einstein radius, θ_E. Once these parameters are obtained, the lens mass, M_L, can be calculated using the relation (Gould 2000)

$$\begin{array}{c}\hfill M_{\mathrm{L}}=\frac{{\theta }_{\mathrm{E}}}{\kappa {\pi }_{\mathrm{E}}},\end{array}$$(1)

where $\kappa =8.144\frac{\mathrm{mas}}{M_{\odot }}$. If the distance to the source is known, the distance to the lens, D_L, can be determined from the equation

$$\begin{array}{c}\hfill D_{\mathrm{L}}=(\frac{{\theta }_{\mathrm{E}}{\pi }_{\mathrm{E}}}{1\mathrm{AU}}+\frac{1}{D_{\mathrm{S}}}{)}^{-1}.\end{array}$$(2)

Microlensing events involving stellar remnants typically have long timescales, allowing for the measurement of the annual microlens parallax from the light curve. This effect arises due to Earth’s motion around the Sun and appears as characteristic deviations, such as asymmetries, in the light curve (Gould 2000). This technique has been successfully applied to numerous microlensing events (e.g. Kruszyńska et al. 2022; Howil et al. 2025).

Measuring the Einstein radius is generally more challenging. One theoretically universal method applicable to any microlensing event involves detecting the astrometric counterpart – an effect known as astrometric microlensing. The astrometric deviation, δ, is directly proportional to the angular Einstein radius and reaches its maximum value at: δ_max ≈ 0.354θ_E (Dominik & Sahu 2000). Since these deviations are typically of the order of a milliarcsecond or smaller, they are challenging to measure. This has been achieved twice using the Hubble Space Telescope (HST). The lens in one of the events, OGLE-2011-BLG-0462, is an isolated black hole with a mass of ∼8 M_⊙ (Sahu et al. 2022; Lam et al. 2022; Mróz et al. 2022; Lam & Lu 2023), while in the other event it is a well-known white dwarf LAWD 37 (McGill et al. 2023). Apart from HST, Gaia has sufficient astrometric precision to detect astrometric microlensing. While these astrometric data are not yet publicly available, they are expected to provide numerous lens mass measurements in the future (Rybicki et al. 2018).

Another approach is the interferometric microlensing. During the microlensing event, where the lens is an isolated object, two source images with different brightness are produced. In this technique, the images are directly detected using interferometry. However, this is challenging in particular due to the small angular separation between the images. Therefore, this method has been successfully implemented only three times: in Dong et al. (2019), Cassan et al. (2022), and Mróz et al. (2025).

A more straightforward way to measure the angular Einstein radius is by determining the source-lens relative proper motion, μ_rel. Then it can be obtained as θ_E = μ_relt_E, where t_E is a timescale of the event. This has been done several times; for example, in Bennett et al. (2024), Rektsini et al. (2024), Bhattacharya et al. (2021), and Vandorou et al. (2020). Nevertheless, measuring the relative proper motion still requires advanced techniques and is only feasible if the lens is luminous, making it unsuitable for dark stellar remnants.

A less common method involves detecting the xallarap effect, which occurs when a source star is a binary system and its orbital motion influences the light curve. Despite the abundance of binary stars, this effect is expected to be present in just ∼20% of the microlensing events, as a source system must have an optimal period with respect to the timescale of the event (Poindexter et al. 2005; Hu et al. 2024). Einstein radius measurement or estimations with xallarap were performed in, for example, Hu et al. (2024) and Smith et al. (2002).

Finally, the finite source (FS) effect, the primary focus of this study, requires specific conditions to be detected, but has been extensively used to measure the Einstein radius. It occurs when the projected separation between a source and a caustic is comparable to the source size (Gould 1994). The magnitude of the FS effect is determined by the angular source radius in units of the angular Einstein radius of the event, typically denoted by ρ. By measuring this parameter and the angular source radius, the angular Einstein radius can be determined.

The FS effect is frequently detected in binary lens events. Most often such lenses are planetary systems (e.g. Wu et al. 2024; Shin et al. 2022; Yee et al. 2021) or binary stars (e.g. Rota et al. 2024; Jung et al. 2015). Masses of single lenses are also commonly measured. Most of them are single main-sequence stars (e.g. Rybicki et al. 2022; Zang et al. 2020; Yee et al. 2015). However, free-floating planets (e.g. Koshimoto et al. 2023; Mróz et al. 2020) and isolated brown dwarfs (e.g. Shvartzvald et al. 2019) have also been detected.

For the FS effect to be robust, there should be a moment during the event when the source-lens separation is of the order of the source size, and the event has to be actually observed during this period (Gould 1994). In single-lens events, this typically requires a precise source-lens alignment, resulting in moderate-to-high magnification events. However, some events with significantly magnified curves do not exhibit the FS effect. Two primary reasons for this include a small source size and a large Einstein radius. The latter possibility is particularly promising in the context of searching for massive lenses such as black holes and other stellar remnants, as it may imply high mass. In cases in which there is no clear detection of the FS effect, it is not possible to measure ρ directly, but it is possible to obtain its upper limit, ρ_lim, and consequently a lower limit of the angular Einstein radius, θ_E, lim, which may be sufficient to infer the lens nature if combined with a lens light analysis.

Similar approaches have been adopted in several studies. Smith et al. (2002) determined ρ_lim and θ_E, lim, but it is not a strong constraint in this case, as the angular Einstein radius limit is quite small, likely attributable to the relatively low magnification of the light curve (A_max ∼ 10). Shvartzvald et al. (2014), Zang et al. (2018), and Bachelet et al. (2022) focused on binary lens events with relatively distant source-caustic passages and also derived ρ_lim using similar methods. The first two additionally estimated the angular Einstein radius limits (both < 100 μas) and used Bayesian analysis for lens properties, while Bachelet et al. (2022) employed Markov chain Monte Carlo (MCMC) methods to refine constraints. In these cases, determining ρ_lim effectively extracts additional information from the light curve, leading to more meaningful constraints on lens properties. Our goal is to apply this technique more systematically.

In this study, we analyse four medium-to-high magnification microlensing events from Gaia that show no apparent FS effect and assess whether their lenses could be stellar remnant candidates. We constrain the angular Einstein radius and the lens mass using the distribution of the ρ parameter, quantifying the FS effect. This constraint is based solely on photometric data, making it a simple yet effective method that can be applied to a broader range of events.

The structure of this paper is as follows. In Section 2, we present the dataset and describe the cleaning process applied to ensure data quality. Section 3 outlines our methodology for analysing the microlensing events. The results for each event are detailed in Section 4. We provide a broader discussion of our findings in Section 5 and summarise our conclusions in Section 6.

2. Data and event selection

We focus on four microlensing events identified in Gaia data: Gaia21efs (ZTF21abqhbuh), Gaia21azb (ASASSN-21ht), GaiaDR3-ULENS-067 (OGLE-2015-BLG-0149), and Gaia21dpb (ZTF21abtulxp). In this section, we describe the photometric data collected for these events (Gaia photometry in Section 2.1 and ground-based observations in Section 2.2), as well as the process of data cleaning (Section 2.3). Additionally, we describe spectroscopic follow-up for Gaia21efs in Section 2.4.

The events were selected based on the following criteria. First, the amplitude of the event’s light curve had to be higher than 3^m. Second, the light curve had to be well covered, particularly at the peak. Third, the initial microlensing model (point source point lens model with parallax, more details in Section 3.1) had to satisfy the conditions u₀ ≲ 0.1 and ρ ≲ u₀, where u₀ is the impact parameter at the closest approach of the source and the lens. The method of constraining the angular Einstein radius, described in Section 3, applies to such events.

Additionally, given our primary interest in identifying dark lenses, we required f_s ≥ 0.8 across most filters. Here, f_s is the blending parameter, which accounts for the potential influence of lens light on the light curve. Section 3.1.2 provides a more detailed description of this parameter.

The astrometric parameters for the sources in the selected events are presented in Table 1. They come from Gaia Data Release 3 (GDR3) (Gaia Collaboration 2023b).

2.1. Gaia photometry

For all events considered, Gaia jointly collected 654 data points in the G filter (Jordi et al. 2010). Three of these events – Gaia21efs, Gaia21azb, and Gaia21dpb – were discovered through the Gaia Science Alerts system (GSA) (Wyrzykowski et al. 2012; Hodgkin et al. 2013, 2021). GaiaDR3-ULENS-067 is from the first Gaia catalogue of candidate microlensing events, based on the GDR3 data (Wyrzykowski et al. 2023; Gaia Collaboration 2023b).

For the events from GSA, the magnitude in the G filter varies from 19.800^m (baseline magnitude for Gaia21dpb) to 12.180^m (maximum magnitude for Gaia21efs). While magnitude errors are available in GDR3, they are not provided for GSA data. To estimate the errors, we used the formula described in Kruszyńska et al. (2022). Accordingly, the error bar for G = 19.800^m is 0.077^m, and for G = 12.180^m it is 0.003^m. The maximum and minimum magnitudes for GaiaDR3-ULENS-067 event from GDR3 are G = 17.053^m ± 0.007^m and G = 14.557^m ± 0.002^m, respectively.

2.2. Ground-based photometric observations

Gaia data alone are not sufficient to accurately constrain the FS effect in the microlensing events, as it observed each object roughly once a month (Gaia Collaboration 2016). To achieve a proper fit, the light curve must be covered as densely as possible. Fortunately, other photometric surveys and numerous ground-based observatories have collected data for the events considered in this study. We briefly describe them in this section.

The Zwicky Transient Facility (ZTF) is one of the major contributors to the data collection for our study (Bellm et al. 2019). This optical time-domain survey uses the 1.2 m telescope at Palomar Observatory, scanning the northern sky every two days. The survey collects data in the ZTF g, r, and i filters. For our study, ZTF provided 7079 points for the events Gaia21efs, Gaia21azb, and Gaia21dpb.

The Asteroid Terrestrial-impact Last Alert System (ATLAS) also contributed data to our study. This survey comprises four 50 cm telescopes located in Hawaii, Chile, and South Africa, allowing it to scan the sky several times per night. ATLAS utilises a variety of filters, including c (cyan), o (orange), g, r, i, H, and OIII (Tonry et al. 2018; Smith et al. 2020). The data from the survey were downloaded from the ATLAS Forced Photometry website (Shingles et al. 2021). For our analysis, ATLAS provided a total of 6168 data points in the o, c, i, g, r, and H filters for the events Gaia21efs, Gaia21azb, and GaiaDR3-ULENS-067.

Data from the Optical Gravitational Lensing Experiment (OGLE) are also available for GaiaDR3-ULENS-067. This survey is focused on dense regions such as the Galactic Bulge and the Magellanic Clouds. It uses a 1.3 m telescope located at Las Campanas Observatory in Chile (Udalski et al. 2015). Microlensing and other transient events are detected using OGLE Early Warning System (EWS) (Udalski 2003). For GaiaDR3-ULENS-067, OGLE provided 3169 data points in OGLE I filter.

Additional surveys that observed the sources include the Panoramic Survey Telescope and Rapid Response System 1 (Pan-STARRS1, PS1) and the DECam Plane Survey (DECAPS), providing high-quality photometry (Chambers et al. 2016; Schlafly et al. 2018). PS1 data from 2011–2012 collected in the PS1 g, r, i, and z filters are available for Gaia21efs, Gaia21azb, and Gaia21dpb. Although these points were not used in modelling the microlensing curve, they were crucial for constraining the baseline magnitudes in filters where the baseline was not directly observed, and for estimating the angular radius of the sources, as described in Sections 3.2 and 3.3. DECAPS also provided precise photometry for GaiaDR3-ULENS-067 from 2017. These data were not used in our analysis, as a different method of source radius estimation was employed for this event (Section 4.3).

Data from infrared space telescopes, including the Two Micron All Sky Survey (2MASS) and the Near-Earth Object Wide-Field Infrared Survey Explorer (WISE/NEOWISE) (Skrutskie et al. 2006; Mainzer et al. 2011), were also available. However, they were not utilised in the analysis due to their sparse coverage.

Lastly, the most numerous contribution came from the follow-up observations of the events from GSA – 27 ground-based telescopes worldwide contributed a total of 8043 data points in UBVRI and ugriz filters. All of them are listed in Table A.1. The telescope apertures vary from 30 cm to 2 m. To process and manage this data, the Black Hole TOM (BHTOM¹) was employed. Images from the ground-based telescopes corrected for bias, dark current, and flat field were uploaded to this platform. It facilitated the extraction of instrumental photometry for the target object as well as other objects in the field. The photometry was then automatically standardised to Gaia Synthetic Photometry (GaiaSP) (Gaia Collaboration 2023a) as is described in Zieliński et al. (2019, 2020).

2.3. Data cleaning

Given the diverse data we have, data cleaning is an essential step to obtain accurate microlensing models. In this section, we give an overview of the data cleaning process for each event. All the photometric data available are listed in Tables B.1–B.4. Note that Gaia G data are available for every event and do not require pre-processing. Also, we do not mention data from infrared telescopes, or from PS1 and DECAPS, as they were not used in light curve modelling. Before any processing, we discarded possible duplicates. After cleaning, we scaled the data so that $\frac{{\chi }^2}{\mathrm{d}.\mathrm{o}.\mathrm{f}.}=1$ for each filter.

Gaia21efs. For this event, data from ZTF, ATLAS, and 20 ground-based telescopes are available. ZTF observations include data in all three filters: ZTF g, r, and i. We used only ZTF g and r data, as ZTF i data cover only the baseline. The majority of data points in the selected filters exhibit a characteristic curve on a magnitude-error plot, stemming from the Poissonian statistics of light. To ensure data uniformity, we removed points that were visually identifiable as obvious outliers in this plot. These outliers constitute < 15% of the r filter data points and < 5% of the g filter data points; thus, their removal does not significantly impact the curve coverage.

ATLAS data are available in ATLAS o, c, and H filters. We did not use the ATLAS H filter due to the small number of points. Some baseline data points in ATLAS o and c filters have significant error bars. Hence, we truncated the error distribution on the 75th percentile for each filter. Given the substantial datasets (> 500 points for the c filter and > 1800 points for the o filter), this still left us with a sufficient number of data points in the baseline.

The remaining data are available in filters U, B, V, R, I, u, g, r, i, and z. We excluded U, u, and z filters from our analysis because of the limited number of points. We also did not use data from TRT-SRO-0.7_Andor-934. Despite the abundance of data points from this telescope, they have higher error bars compared to other observations at the same magnitudes. The remaining data, even after cleaning, adequately covers the curve. For the B filter, we truncated the error distribution at the 95th percentile to eliminate outliers. We applied the same procedure to the r filter. For the g and i filters, we exclusively utilised LCO data due to their smaller error bars and predominant representation in the dataset.

Data in V, R, and I filters come from a large set of data sources, so the cleaning technique we applied to them was rather unconventional. Firstly, we manually removed several apparent outliers. Then, for each filter, we divided the data roughly at the midpoint of the light curve. Depending on the extent of the distribution tails, we truncated them at either the 75th or 90th percentile. This approach ensures good-quality data near the peak, while maintaining sufficient data in the lower parts of the curve and the baseline.

Gaia21azb. Data for this event include observations from ZTF, ATLAS, and 12 ground-based observatories. We used ZTF g and r data without additional pre-processing and excluded ATLAS data due to high error bars at the given brightness. The remaining data are in the U, B, V, R, I, u, g, r, and i filters. We excluded the U and u filters due to the small number of points, and the B filter because of its large error bars and sparse coverage. The V filter was also excluded, as it appears to exhibit a systematic bias that adversely affected the model fit.

For the g, r, and i filters, we excluded the data points with the highest errors, truncating the error distribution of each filter at either 90th or 95th percentile, depending on its shape. The error bars of the V, R, and I filters are more significant than for the rest of the data. To make these datasets more uniform, we took the data only from the bigger telescopes (ASV1.4_Andor, TJO_MEIA2, TRT-SBO-0.7, pt5m). The error distributions of the selected data were additionally truncated at the 90th or 95th percentile, depending on their shapes.

GaiaDR3-ULENS-067. For this event, only data from surveys, namely from OGLE and ATLAS, are available. We used only OGLE data points without additional pre-processing, as they have small error bars, no significant outliers, and adequately cover both the light curve and the baseline.

Gaia21dpb. This event was observed by ZTF, LOIANO1.52_BFOSC, HAO68_G2-1600, and PIWNICE90_C4-16000EC. We used ZTF g and r data without additional pre-processing. The remaining data mostly covered the baseline and could not be used to model the microlensing curve.

2.4. Spectroscopic follow-up of Gaia21efs

In addition to photometry, Gaia21efs was bright enough – especially during the peak magnification – to obtain high-resolution spectroscopic data by using the Potsdam Echelle Polarimetric and Spectroscopic Instrument (PEPSI, Strassmeier et al. (2015)) mounted at the 2 × 8.4-m Large Binocular Telescope (LBT)² which is located on Mt. Graham, Arizona, US. The event was observed on November 3, 2021. The PEPSI configuration of a 300 μm fiber diameter and two cross-dispersers, CD2 (blue spectral arm) and CD5 (red spectral arm), were used simultaneously, which gave us a high signal-to-noise ratio (218) and a resolving power of R ≈ 40 000. After the standard calibration process (SDS4PEPSI, Ilyin 2000, PhD thesis, Univ. Oulu), we had obtained spectroscopic data divided into two parts: blue, covering the wavelength range between 422 and 478 nm, and red, covering the wavelength range between 624 and 743 nm.

3. Methodology

In this section, we describe the methodology we applied to measure the lower limit of the angular Einstein radius in every event studied. It comprises several steps.

First, we used the MCMC technique to explore the parameter space of the combined model described in Section 3.1 and derived posterior distributions for its parameters. The modelling process is outlined in Section 3.2.

Next, we calculated the angular source radius, θ_*, using colour-angular diameter relations, as is described in Section 3.3. Using this, we transformed a posterior distribution of ρ (the angular size of the source relative to the angular Einstein radius) into the Einstein radius distribution, from which we estimated the lower limit of this parameter. This forms the basis for further analysis and interpretation.

Then, we employed the DarkLensCode (DLC) software to estimate the lens mass, the lens distance, as well as the probability of the lens being a dark object (Section 3.4). Finally, using the angular Einstein radius limit and the DLC results, we aimed to draw conclusions regarding the possible nature of the lensing object.

3.1. Photometric model

We used a finite source point lens (FSPL) model in the vicinity of the peak and the point source point lens model (PSPL) in the remaining parts of the curve for all events. The specific time range for the FSPL model is given in Section 3.2. Both models incorporate the annual microlens parallax. We used the PSPL model with a parallax as the initial model at the stage of event selection (Section 2).

The expression for the photometric curve magnification, A(u), in the case of the point source is as follows (Paczynski 1986):

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

Here, u is the angular separation between the lens and the source in the units of angular Einstein radius in a given moment. It is usually represented as the combination of components τ and β:

$$\begin{array}{c}\hfill u(t)=\sqrt{\tau {(t)}^2+\beta {(t)}^2}.\end{array}$$(4)

For the PSPL model with the annual parallax, the components are given by

$$\begin{array}{c}\hfill \tau (t)=\frac{(t-t_0)}{t_{\mathrm{E}}}+\delta \tau \end{array}$$(5)

and

$$\begin{array}{c}\hfill \beta (t)=u_0+\delta \beta .\end{array}$$(6)

In these expressions, u₀ represents the impact parameter, defined as the minimum separation between the source and the lens, occurring at the time of the light curve maximum, t₀, and t_E denotes the timescale of the event. For each event, we find two very similar models corresponding to u₀ > 0 and u₀ < 0, as it is a common degeneracy. δτ and δβ are the corrections necessary to take into account the annual parallax. They are given by the following expression (Gould 2004):

$$\begin{array}{c}\hfill (\delta \tau ,\delta \beta )=({\mathit{\pi }}_{\mathrm{E}}·\mathrm{\Delta }\mathit{s},{\mathit{\pi }}_{\mathrm{E}}\times \mathrm{\Delta }\mathit{s}),\end{array}$$(7)

where π_E is the microlens parallax vector. Its northern and eastern components, denoted by π_EN and π_EE, respectively, are included in the model. Δs denotes the positional offset of the Sun in the geocentric frame of reference.

We selected events with medium-to-long timescales, where parallax effects are typically detectable. To ensure that the parallax signal is not an artefact introduced by a single dataset, we compared the χ² values for models with and without a parallax separately for each filter. In most cases, the inclusion of parallax led to a significant improvement in fit quality, with Δχ² > 5 in the majority of filters, and Δχ² > 10 in several – indicating medium to strong evidence of the presence of a parallax across multiple datasets.

3.1.1. Finite source effect

To incorporate the finite source (FS) effect, we integrated Formula 3 over the area of the source (Witt & Mao 1994):

$$\begin{array}{c}\hfill A_{\mathrm{FS}}(u)=\frac{1}{\pi {\rho }^2}{\displaystyle {\int }_{S_{\mathrm{source}}}}A(u)d{S}_{\mathrm{source}}.\end{array}$$(8)

The parameter ρ is included in the model and is given by

$$\begin{array}{c}\hfill \rho =\frac{{\theta }_{\ast }}{{\theta }_{\mathrm{E}}}.\end{array}$$(9)

In this equation, θ_* represents the angular radius of the source. Typically, the FS effect becomes significant when the source size and the source-lens separation are comparable, i.e. when ρ ≈ u₀. In our events, however, we find that ρ < u₀, meaning that ρ cannot be directly measured but only constrained.

We employed the high-magnification approximation method outlined in Gould (1994). It simplifies the magnification for small u to $A(u)\approx \frac{1}{u}$, facilitating easier integration. This approach is computationally efficient and remains sufficiently accurate for our analysis. Comparing its performance near the peak with the more precise FSPL method outlined in Lee et al. (2009), we found that the most pronounced difference between the models is ∼0.001%, with a root mean square error of ∼1.5 ⋅ 10⁻³ in proximity to the peak.

We do not include limb darkening in our model, as we focus on events with ρ < 0.1. In these instances, it would not significantly alter the results, but it would introduce unnecessary complexity.

3.1.2. Blending

Microlensing light curves can be affected by additional light that does not originate from the microlensed source. The most common causes include unresolved nearby sources in dense fields or the luminous lens. Therefore, the total flux, F(t), is given by the expression

$$\begin{array}{c}\hfill F(t)=A(t)F_{\mathrm{S}}+F_{\mathrm{B}},\end{array}$$(10)

where F_S represents the source flux and F_B denotes the blend flux, originating from sources other than the microlensed one. We calculated them separately for every filter and for each set of model parameters using the MulensModel Python package described below.

To quantify the contribution of the source flux to the total flux, the blending parameter, f_s, is typically introduced. It is defined as the fraction of the source flux in the total detected flux (Woźniak & Paczyński 1997):

$$\begin{array}{c}\hfill f_{\mathrm{s}}=\frac{F_{\mathrm{S}}}{F_{\mathrm{S}}+F_{\mathrm{B}}}.\end{array}$$(11)

For events with negligible blend flux (i.e. with f_s ≈ 1), we assume that all the light detected comes from the microlensed sources. Therefore, when the blending is close to unity in most filters for a particular event, we consider that the lens there is most likely dark.

With the source and blend fluxes, we calculated the apparent magnitude at baseline, denoted as m₀, using the formula

$$\begin{array}{c}\hfill m_0=22.0^m-2.5log(F_{\mathrm{S}}+F_{\mathrm{B}}),\end{array}$$(12)

where 22.0^m is assumed to be a zero point value.

3.2. Modelling

In summary, the model we used includes the following parameters: (t₀, u₀, t_E, π_EN, π_EE, ρ) and pairs (F_S, F_B) for each filter. We applied the MCMC method using the emcee Python package, which utilises an ensemble sampler, allowing us to explore parameter space in parallel with multiple chains (‘walkers’) and achieve faster convergence and better sampling accuracy (Foreman-Mackey et al. 2013).

We also used the MulensModel package to streamline the handling and analysis of microlensing data (Poleski & Yee 2019). It is helpful in particular in fitting (F_S, F_B) pairs for every model sampled. Additionally, in MulensModel, various methods for calculating magnification are available for both point and finite source scenarios. We employed the finite_source_uniform_Gould94 method in the vicinity of the peak, specifically during the time range corresponding to u ≤ 3u₀ (Gould 1994). This is the implementation of the Gould approximation described in Section 3.1.1. For the remainder of the curve, we employed the point_lens method, which is the default option.

In Gaia21efs and Gaia21azb, many filters lack sufficient baseline coverage, which leads to blending values exceeding one – unphysical in the standard microlensing formalism. To mitigate this, we applied Gaussian priors on the baseline magnitudes for selected filters:

$$\begin{array}{c}\hfill \ln (\mathrm{Prior})=-\frac{1}{2}{\left(\frac{m_0-m_{\exp }}{\sigma }\right)}^2.\end{array}$$(13)

Here, m₀ is the baseline magnitude inferred from the model, m_exp is the expected baseline magnitude based on independent photometry, and σ represents the uncertainty set to 1–5% of the expected magnitude. The expected baseline magnitudes for the BVRI and gri bands were derived from PS1 photometry, using the colour transformations provided by Tonry et al. (2012). Additionally, we applied a Gaussian prior to the blend flux:

$$\begin{array}{c}\hfill \ln (\mathrm{Prior})=-\frac{1}{2}{\left(\frac{F_{\mathrm{B}}}{0.05}\right)}^2.\end{array}$$(14)

This prior favours low blend flux values without imposing a hard constraint of zero. More details on using these priors can be found in Sections 4.1 and 4.2.

The filters lacking baseline coverage include observations from both large and small telescopes. Therefore, anomalous blending values in these bands could also be a result of unresolved nearby sources or photometric systematics in the lower-resolution data from smaller telescopes. We tested models excluding data from small telescopes, but blending did not improve noticeably without them. Therefore, we retained these datasets to preserve good coverage of the curve.

The posterior distributions of all parameters except ρ typically exhibit Gaussian behaviour in our analysis. We characterised these distributions by determining the median as the parameter value. Subsequently, their uncertainties were quantified as the differences between the medians and the 16th or 84th percentiles. However, in cases in which the distribution displays noticeable skewness, we employed an alternative procedure. Specifically, the parameter value is a mode of distribution, and its uncertainties are the differences between the mode and the boundaries of the highest density interval, encompassing 68% of the samples.

The posterior distribution of ρ, in turn, tends to be concentrated near 0. We characterised this distribution by the upper limit, further denoted as ρ_lim. We defined it to be at the 95th percentile of its posterior distribution, corresponding to a one-sided Bayesian credible interval. The 95% threshold is a common and statistically sound choice, excluding the least probable tail of the distribution without being overly conservative. It is worth noting that this limit may be sensitive to the starting value of ρ in the MCMC process. To prevent it from being underestimated, we set the initial value of ρ to be approximately equal to u₀. This is the highest reasonable value, as if the FS effect is not present, the ρ distribution is expected to be smaller than u₀.

3.3. Source stars

We used parallaxes and proper motions of the sources published in GDR3. These parameters are listed in Table 1. Our analysis also requires their angular radius, θ_*, distance, D_S, and extinction in the G filter, A_G. We now describe how we obtained these parameters.

3.3.1. Angular radius

We computed the angular radius using empirical colour-angular diameter relations provided in Adams et al. (2018). These relations are constructed for giants or dwarfs and subgiants. To discern the type of microlensed source, we constructed colour-magnitude diagrams (CMDs) using objects from the GDR3 within 5′ circles around the source. Based on them, we treated the sources in Gaia21efs, Gaia21azb, and GaiaDR3-ULENS-067 as giants and the source in Gaia21dpb as a main-sequence star. For Gaia21efs, this classification is consistent with the spectroscopic analysis. Then, we employed the appropriate relations to calculate the source sizes.

We calculated colours for Gaia21dpb, Gaia21azb, and Gaia21efs using PS1 photometry. We corrected these for extinction using values for PS1 filters available in Schlafly & Finkbeiner (2011). Next, we converted PS1 filters to the UBVRI system. We did this by applying relations from Tonry et al. (2012). For the source of GaiaDR3-ULENS-067, a different set of filters and procedures were employed, as it is a red clump (RC) star from the Bulge. This is detailed in Section 4.3.

Given the posterior distribution of ρ and an estimate of the source angular radius, θ_*, we derived a lower limit on the Einstein radius, θ_E, lim. To avoid unphysically large θ_E values caused by the few extremely small ρ samples, we discarded the lowest ρ values. The cut-off was placed where the number of samples begins to rise rapidly, indicating the start of the well-constrained portion of the posterior. We then generated a Gaussian distribution of θ_* with the same number of samples as ρ, accounting for its measurement uncertainty. The resulting θ_E distribution was computed via sample-wise division (transformation of Formula 9). θ_E, lim was then taken as the fifth percentile of this distribution, corresponding to a one-sided 95% lower credible interval, analogous to ρ_lim estimation.

3.3.2. Distance

The distances to the sources were determined based on various considerations specific to each event. When the RUWE of the Gaia astrometric solution approximately equals unity, we can reasonably rely on Gaia parallaxes. This applies to the sources of Gaia21efs and Gaia21azb (RUWE values are given in Table 1). Accordingly, for Gaia21azb, we adopted the distance value from Bailer-Jones et al. (2018, 2021). In contrast, for Gaia21efs, we used a spectroscopic distance derived from a high-resolution spectrum. For GaiaDR3-ULENS-067 and Gaia21dpb, we used other approaches, described in Sections 4.3 and 4.4, respectively. Below, we describe how the atmospheric parameters and spectroscopic distance for Gaia21efs were obtained.

We analysed a high-resolution spectrum from LBT/PEPSI using the iSpec³ framework (Blanco-Cuaresma et al. 2014; Blanco-Cuaresma 2019). The analysis employed the SPECTRUM⁴ radiative transfer code, MARCS model atmospheres (Gustafsson et al. 2008), and solar abundances from Grevesse et al. (2007).

Figure 1 presents the observed LBT/PEPSI spectrum together with the best-matching spectrum synthesised for the following atmospheric parameters: effective temperature T_eff = 4600 ± 22 K, surface gravity log g = 1.70 ± 0.07, metallicity [M/H] = −0.19 ± 0.03 dex, and microturbulence velocity v_t = 1.26 ± 0.02 km s⁻¹. Based on this solution, we classify the source star in Gaia21efs event as G8 II spectral type and estimate its absolute magnitude to be M_V = 0.9^m ± 0.3 (Straižys 1992). Assuming an apparent magnitude of V = 16.546^m (from the microlensing model) as well as a line-of-sight extinction of A_V = 2.235^m (Schlafly & Finkbeiner 2011), we were able to estimate the spectroscopic distance to the source star as D_S = 4.81 ± 0.91 kpc. The distances presented in Bailer-Jones et al. (2018, 2021) are in agreement with this value within 1σ uncertainty. Therefore, in the case of Gaia21efs, the above spectroscopic distance was adopted in order to constrain the distance and mass to the lens in this specific microlensing event.

Fig. 1. Spectrum of the Gaia21efs obtained with LBT/PEPSI on November 3, 2021 (blue) and the best-matching fit (red) synthesised for the specific parameters. The Fe I regions in blue (top) and red spectral arms (bottom) are presented.

3.3.3. Extinction

Extinction in the G filter is provided in the GDR3 as the ag_gspphot parameter. It is available only for one of the sources considered in this work: Gaia21azb. For Gaia21efs and Gaia21dpb, we adopted a similar procedure as is described in Kruszyńska et al. (2022). Specifically, we calculated the average ag_gspphot extinction from approximately four nearest sources within 20″ circle around the source. Finally, we computed extinction for GaiaDR3-ULENS-067 using the CMD (Section 4.3).

3.4. DarkLensCode

The DarkLensCode (DLC) is a software used to obtain the posterior distribution of the lens mass and distance, as well as to estimate its probability of being dark. It utilises the posterior distributions of the photometric model parameters and a Galactic model. DLC is based on the approach described by Wyrzykowski et al. (2016), Mróz & Wyrzykowski (2021) and is detailed in Howil et al. (2025). Below, we briefly summarise the main steps of the DLC procedure.

We provide the event co-ordinates, the adopted mass function, and the posterior distributions of the parameters: t₀, u₀, t_E, π_EN, π_EE, m_0, G, and f_S, G (where the subscript G denotes quantities related to the Gaia G filter). Additionally, the source distance, D_S, components of the source proper motion, μ_S, and extinction, A_G, are included as inputs.

The procedure begins by drawing a random sample from the posterior distribution. The relative proper motion of the lens and source, μ_rel, is sampled from a uniform distribution in the range (0, 30) mas yr⁻¹, while D_S is drawn from either a uniform or Gaussian distribution, depending on the configuration parameter ds_weight. Using these values, the angular Einstein radius, θ_E, lens mass, M_L, and lens distance, D_L, are computed for each sample. This process is repeated for the desired number of iterations. Each sample is subsequently weighted based on the Galactic priors.

These priors include the relative proper motion, the lens distance, and the lens mass. The distance prior comes from Han & Gould (2003) and Batista et al. (2011). The prior for relative proper motion is assumed to be Gaussian, derived from separate Gaussian distributions for the source and lens proper motions. The source proper motion is characterised by an input mean value, μ_S, and variance, σ_S (see Table 1). The lens proper motion mean, μ_L, is computed using standard kinematic relations from Reid et al. (2009). The variance of the lens proper motion is determined by its position in the Galaxy: (σ_l, σ_b) = (100, 100) km s⁻¹ for Bulge lenses and (σ_l, σ_b) = (30, 20) km s⁻¹ for disc lenses (Howil et al. 2025).

The lens mass prior, i.e. the mass function, is another input that we define. Since the true mass function for dark stellar remnants remains uncertain, we explored several plausible assumptions. Specifically, we tested the standard Kroupa mass function for stars (Kroupa 2001), an approximate mass function for stellar remnants from Mroz et al. (2021), and a simple power-law form, f(M)∝M⁻¹. Among these, the Kroupa mass function yields the lowest inferred lens masses and the most conservative estimates of dark lens probabilities. We therefore adopted it as our default in the analysis. However, if the actual mass function more closely resembles that proposed by Mroz et al. (2021), the evidence supporting dark lens candidates would be even stronger.

The dark lens probability was computed by comparing the inferred brightness of the lens with the expected brightness of a main-sequence star of the same mass. The actual lens brightness distribution was derived from the blending and baseline magnitude posterior distributions provided as input. For each sample, the expected main-sequence brightness was calculated based on its estimated mass and distance, using empirical relations from Pecaut & Mamajek (2013). DLC computes the probability for zero extinction and for the extinction provided as an input, which is considered the maximum possible extinction. Consequently, we obtained a range for the dark lens probability for zero and maximum extinction.

As we have determined the lower limit for the angular Einstein radius, we included only the samples for which θ_E ≥ θ_E, lim in the resulting posterior distributions of mass, distance, and lens brightness. From these filtered samples, we computed the median mass and distance, as well as the dark lens probability.

It is also useful to calculate the transverse velocity of the lens v_⊥, L, which is given by

$$\begin{array}{c}\hfill v_{\perp }={\mu }_{\mathrm{L}}{D}_{\mathrm{L}},\end{array}$$(15)

where the proper motion of the lens μ_L can be expressed as the vector sum μ_L = μ_rel + μ_S (Gould 2000). DLC provides us with the posterior distributions of μ_rel, μ_S, and D_L. Since the microlensing parallax vector π_E and the relative proper motion μ_rel are aligned, we can determine the angle between μ_S and μ_rel using the dot product, allowing us to compute the full vector, μ_L(Gould 2014). Combining all components, we derived the posterior distribution of the lens transverse velocity, which is weighted as the standard DLC posteriors. As before, we used the θ_E, lim constraint for this distribution and computed its median.

4. Results

In this section, we present the results of our analysis for each microlensing event. The key parameters, including the Einstein radius limit, lens mass, and distance, along with the dark lens probability from DLC, are summarised in Table 2. Additionally, we list the lower mass limit and upper distance limit of the lens, calculated using θ_E, lim with Formula 1 and Formula 2, respectively.

We present light curves of the events and two types of plots generated by the DLC. The first is the mass-distance plot. The second is the blend-lens plot representing the blend light compared to the lens light, assuming it is a main-sequence star. The dashed line in the blend-lens plot indicates the extinction’s lower limit, and the solid line shows its upper limit. The colours indicate the probability density on a logarithmic scale (base 10).

The bright samples in the DLC plots meet the condition θ_E ≥ θ_E, lim. They are overlaid on the shaded samples, which represent the original distribution without considering this limit. All plots representing the light curve or the DLC results show only the model with u₀ > 0, as the models with u₀ < 0 yielded very similar plots. In the figures in which the light curves are presented, both models FSPL and PSPL include the annual parallax.

4.1. Gaia21efs

Initial modelling using all available datasets for this event resulted in unphysical blending values (f_S > 1) in many filters, including Gaia G, likely due to the lack of baseline coverage in some of them. To investigate this, we initially modelled the event using only the datasets that provided good coverage of both the magnification peak and the baseline. Specifically, we began with the following filters: Gaia G, ZTF g, ZTF r, ATLAS o, and ATLAS c. This combination yielded a model with blending parameters close to unity and consistent with physical expectations. We refer to this as the base model, which serves as a reference point in our analysis; the inclusion of additional filters or priors is expected to preserve the core parameters of this model.

Fig. 2. Light curve of Gaia21efs with the data from different surveys as well as the follow-up data collected with BHTOM. The solid line shows the parallax FS model, while the dashed line shows a parallax model for the point source. The right figure shows the zoom-in on the peak of the light curve.

Fig. 3. DLC results for Gaia21efs. Top: Mass-distance plot. Bottom: Blend-lens plot.

We then added the g- and i-band data from LCO, which, despite lacking a baseline, are homogeneous in quality and sample the peak well. With this subset (Gaia G, ZTF g, ZTFr, ATLAS o, ATLAS c, g, and i), physically plausible blending ratios were recovered by applying Gaussian priors on the blend flux and baseline magnitude in the ATLAS o, ATLAS c, g, and i bands, as is described in Section 3.2. Although ATLAS data provide baseline coverage, their blending estimates remained sensitive to the addition of filters without baseline, and thus also required priors. The baseline magnitude prior means for the ATLAS filters were set to the values from the base model.

Finally, we included the remaining filters while maintaining the same priors and fixing blend flux in the r filter to zero. This allowed us to incorporate all available data without introducing unphysical blending. This constitutes the final model presented in this work. It differs from the base model at ∼2σ level, particularly in the microlens parallax, which is sensitive to the additional filters.

The light curve of Gaia21efs is shown in Figure 2. The lens mass limit for this event exceeds 0.7 M_⊙, and the dark lens probability is over 85%. The top panel of Figure 3, which displays the mass-distance plot, highlights the significant influence of this limit on the DLC results. Without taking it into account, the median mass from the DLC analysis would be around ∼0.4 M_⊙. In contrast, when considering this limit, the inferred mass is ∼1.0 M_⊙. This discrepancy suggests that such high-mass objects are relatively rare in this region of the Galaxy. This observation motivated us to estimate the transverse velocity of the lens, as described in Section 3.4, as it could have travelled to this region from a different place in a Galaxy. A high velocity, potentially caused by a supernova kick, would provide additional evidence supporting its compact nature.

When taking into account the Einstein radius limit, the median transverse velocity of the lens is ∼120 km s⁻¹, with the lower bound of ∼90 km s⁻¹. These values are comparable with the mean transverse velocity of recycled pulsars (87 ± 13 km s⁻¹) reported by Hobbs et al. (2005), but significantly lower than of non-recycled pulsars (246 ± 22 km s⁻¹). At the same time, Wegg & Phinney (2012) suggested that white dwarfs with masses greater than 0.75 M_⊙, as in our case, typically exhibit transverse velocities up to 60 km s⁻¹. Given that the lens velocity exceeds this range, it is consistent with either a relatively high-velocity white dwarf or a low-velocity neutron star.

4.2. Gaia21azb

Similarly to Gaia21efs, the initial model of Gaia21azb exhibited unphysical blending in some filters because of incomplete baseline coverage, so we applied a similar approach as for the previous event. In this case, the filters covering both the baseline and the curve are Gaia G, ZTF g, and ZTF r. The blending parameters for this combination were physically reasonable and close to one, forming the base model for this event.

Fig. 4. Light curve of Gaia21azb with the data from different surveys as well as the follow-up data collected with BHTOM. The solid line shows the parallax FS model, while the dashed line shows a parallax model for the point source. The right figure shows the zoom-in on the peak of the light curve.

Subsequently, we modelled the event using an extended set of filters: Gaia G, ZTF r, and ZTF g, g, and i. The photometry in the g and i bands comes from LCO and was treated with priors as is described in Section 3.2, leading to physically consistent blending values across all filters. Adding the remaining filters with the same treatment produced similar results, so we adopted this as our final model. Its parameters are consistent with the base model within 1σ, except for u₀ and t_E, which agree within 3σ.

The light curve of Gaia21azb is shown in Figure 4, and its results for the DLC analysis are in Figure 5. Its higher distance limit is ∼3.7 kpc, and the lens lower mass limit is ∼0.5 M_⊙.

Fig. 5. DLC results for Gaia21azb. Top: Mass-distance plot. Bottom: Blend-lens plot.

The median mass derived from the DLC analysis, considering the Einstein radius limit, is ∼1.0 M_⊙. The dark lens probability for Gaia21azb falls within the range of 86-95%. Additionally, we estimated the transverse velocity of the lens to exceed 150 km s⁻¹, which lies between the typical velocities of recycled and ordinary pulsars reported by Hobbs et al. (2005). Based on these findings, we conclude that the lens is most likely dark and may plausibly be either a high-mass white dwarf or a low-mass neutron star.

4.3. GaiaDR3-ULENS-067

Fig. 6. Light curve of GaiaDR3-ULENS-067 with the data from Gaia and OGLE. The solid line shows the parallax FS model, while the dashed line shows a parallax model for the point source. The right figure shows the zoom-in on the peak of the light curve.

Fig. 7. Colour-magnitude diagram (CMD) for GaiaDR3-ULENS-067. The red square represents the source, the yellow circle marks the RC centroid, and the black dots denote objects within a 5′ radius of the source.

Fig. 8. DLC results for GaiaDR3-ULENS-067. Top: Mass-distance plot. Bottom: Blend-lens plot.

Based on the CMD, the source in GaiaDR3-ULENS-067 is identified as a RC star from the Galactic Bulge (Figure 7). Consequently, we estimate the distance to it to be 8 ± 2 kpc (Zoccali & Valenti 2016). To derive the dereddened colour of the source, crucial for determining its angular radius, we employed the method outlined in Yoo et al. (2004). Specifically, we utilised the following formula:

$$\begin{array}{c}\hfill {((G_{\mathrm{BP}}-G_{\mathrm{RP}}),G)}_0={((G_{\mathrm{BP}}-G_{\mathrm{RP}}),G)}_{\mathrm{RC},0}+\mathrm{\Delta }((G_{\mathrm{BP}}-G_{\mathrm{RP}}),G).\end{array}$$(16)

In this equation, ((G_BP − G_RP),G)₀ and ((G_BP − G_RP),G)_RC, 0 represent the dereddened colours and apparent magnitudes of the source and the centroid of the RC, respectively. Δ((G_BP − G_RP),G) denotes the difference between the measured source and RC centroid parameters extracted from the CMD. For the source they are equal to ((G_BP − G_RP),G) = (2.32, 16.98 ± 0.02)⁵ and for the RC ((G_BP − G_RP),G)_RC = (2.10, 16.83 ± 0.03).

We obtained the values for the colour (G_BP − G_RP)_RC, 0 and absolute magnitude, M_G, RC, 0, of the RC centroid from Plevne et al. (2020). This paper provides G_RC, 0 two values for distances: calculated using the inverted Gaia parallax and taken from Bailer-Jones et al. (2018, 2021). Furthermore, it offers values for two RC populations: high-α and low-α. Since the RC in the Bulge belongs to the high-α population, we calculated G_RC, 0 using the average over two distances corresponding to this population. The resulting values are ((G_BP − G_RP), G)_RC, 0 = (1.21±0.07, 15.09±0.24).

Finally, combining all of the necessary values, we obtained ((G_BP − G_RP),G)₀ = (1.43 ± 0.072, 15.24 ± 0.24). Now we could use G₀ to calculate extinction in the G filter as follows: A_G = G − G₀ = 1.74 ± 0.24. Then we transformed the colour (G_BP − G_RP)₀ to ((V − I), I)₀ using relations from GDR3 documentation⁶ and calculated the source size with a relation for giants from Adams et al. (2018).

The light curve for GaiaDR3-ULENS-067 is shown in Figure 6 and the DLC results are in Figure 8. The estimated mass limit of 0.3 − 0.4 M_⊙, along with a dark lens probability exceeding 80%. While this mass limit does not significantly affect the DLC result, the median mass derived from the DLC for the u₀ > 0 model is ∼1.70 M_⊙. This value exceeds the Chandrasekhar limit, indicating that a dark remnant with this mass would be a neutron star. The median mass for the u₀ < 0 model is slightly lower at ∼1.18 M_⊙, indicating that a high-mass white dwarf or low-mass neutron star could act as a lens.

The transverse velocity estimate further supports the possibility that this lens is a neutron star. For u₀ > 0, the median transverse velocity of the lens is ∼320 km s⁻¹, and for u₀ < 0, it is ∼300 km s⁻¹. The median value is higher than the mean for the pulsar transverse velocity distribution reported by Hobbs et al. (2005).

Our results are in good agreement with Rybicki et al. (2024). They analyse this event (listed as OGLE-2015-BLG-0149) using a similar method and incorporating Spitzer data. They obtain a similar mass estimation and also point out that this lens as a strong stellar remnant candidate.

4.4. Gaia21dpb

Fig. 9. Light curve of Gaia21dpb with the data from Gaia and ZTF. The solid line shows the parallax FS model, while the dashed line shows a parallax model for the point source. The right figure shows the zoom-in on the peak of the light curve.

Fig. 10. DLC results for Gaia21dpb. Top: Mass-distance plot. Bottom: Blend-lens plot.

As the blending parameter was around 1.1 in the ZTF g filter, we fixed the blend flux in the ZTF g filter to 0 during the modelling. This made our model more physically plausible without significantly changing the results.

Based on the CMD, the source in Gaia21dpb is identified as a main-sequence star. Leveraging this information, we estimated the distance to the source. Initially, we converted the dereddened PS1 colours to UBVRI colours, specifically (V − I)₀ and (V − R)₀, as is described in Section 3.3. These values are (V − I)₀ = 0.40^m ± 0.02 and (V − R)₀ = 0.15^m ± 0.02. Next, to determine the star type, we compared these colours with the values provided by Pecaut & Mamajek (2013). According to this comparison, the source is classified as approximately F0V to F1V type. By averaging the absolute magnitudes of these types, we obtained a value of M_V = 2.665^m. Pecaut & Mamajek (2013) did not provide errors for the absolute magnitude, so we added a conservative error of 0.5^m, similar to Kruszyńska et al. (2022). We then calculated the distance using the well-known formula

$$\begin{array}{c}\hfill D_{\mathrm{S}}={10}^{\frac{V_0-M_{\mathrm{V}}+5}{5}}.\end{array}$$(17)

The V₀ value was also derived from dereddened PS1 photometry and is equal to 17.34^m ± 0.01. The calculated distance falls within the range of 6.84 to 10.83 kpc, with a central value of 8.61 kpc, which is subsequently utilised in the DLC procedure.

The light curve of Gaia21dpb is demonstrated in Figure 9, while the DLC results are in Figure 10. The lower mass limit for this event is below ∼0.014 M_⊙, providing limited information about its mass or nature. However, the DLC results still suggest the possibility of a low-mass white dwarf, since its mass is ∼0.4 M_⊙. Furthermore, its probability of being dark falls within the range of 80–90%, indicating a high likelihood of it being a dark object.

5. Discussion

In this work, we have analysed four unusual microlensing events: Gaia21efs, Gaia21azb, GaiaDR3-ULENS-067, and Gaia21dpb. Their light curves are moderately or highly magnified, but they do not exhibit any distinct FS effect. We constrained this effect and obtained the lower limits for the angular Einstein radius, and hence for the lens mass for each event. One of the possible explanations for the absence of the FS effect, which motivated us to pursue this work, is a large Einstein radius, as it is inversely proportional to the parameter, ρ, measuring the FS effect. This is promising in the context of searches for massive stellar remnants, as a large Einstein radius typically implies a large lens mass.

We find that the lenses in Gaia21efs, Gaia21azb, and GaiaDR3-ULENS-067 are the stellar remnant candidates, while Gaia21dpb presents a slightly weaker case. We make this statement based on several reasons. First, the lower mass limits for the first three events are > 0.2 M_⊙, which is roughly the lowest possible mass for a white dwarf. Second, the blending is close to unity in these events, implying a dark lens. Finally, the lens light distribution analysis from the DLC procedure favours the hypothesis that these are indeed dark lenses rather than faint low-mass main-sequence stars.

In terms of particular remnant types, the lens in Gaia21efs is most likely a massive white dwarf; however, the estimate for its transverse velocity is relatively high, and could indicate the presence of a neutron star. Gaia21azb presents a similar case. Although a white dwarf is also suggested by the DLC analysis for Gaia21dpb, the evidence is not as strong due to the low mass limit. In GaiaDR3-ULENS-067, the large estimated mass and high transverse velocity of the lens strongly suggest that it is a neutron star. However, the possibility that it is a massive white dwarf cannot be ruled out.

There are other works in which dark isolated stellar remnant candidates are identified using microlensing light curves from Gaia (Section 1). Along with these candidates, the four lenses found in this work may contribute to the sample of isolated remnants in the Galaxy, allowing us to study their mass and spatial distribution.

As is described in Section 1, various methods exist for measuring lens masses, but they often present technical challenges or are applicable only to specific, uncommon cases. In this study, we focus on constraining the lens mass, which, as we demonstrate, can provide valuable insights into the lens nature. Our methodology is straightforward, applicable to events with u₀ ≲ 0.1, and exceptionally high magnification is not required for obtaining meaningful results. Applying this methodology to a larger set of events could be a valuable direction for future research.

Similar approaches were mentioned in Section 1 (Smith et al. 2002; Shvartzvald et al. 2014; Zang et al. 2018; Bachelet et al. 2022). These studies determine ρ_lim, and the latter three use techniques analogous to the DLC procedure. Although they share similarities with our work, we apply our methodology to single-lens, single-source events as a systematic approach to identifying stellar remnant candidates.

A primary limitation of our approach is the difficulty in predicting whether the Einstein radius limit will be meaningful from initial observations. Analysing our four events along with those from the four cited papers, we find no clear distinguishing characteristics. While magnification can be a factor, it appears that the method can be applied to moderately magnified events. For instance, Gaia21efs, the event with the highest Einstein radius and mass limit, has a maximum magnification of ∼28. In contrast, Gaia21azb, the most highly magnified event with a maximum magnification of ∼80, does not exhibit a particularly high mass limit. However, data quality near the light curve peak is crucial for accurately measuring ρ. Therefore, while there may be numerous sufficiently magnified events, the method may be ineffective if data are limited or have high uncertainties.

When searching for massive objects, it is advisable to assess the source radius beforehand to ensure it is not excessively small. It is generally more advantageous to plan follow-up observations with a larger source radius. With a smaller Einstein radius, the FS effect will probably be observed, enabling the lens mass to be derived, although it might not be significant. Otherwise, there is a higher chance of a substantial mass limit. This strategy may create promising prospects for faint yet highly magnified microlensing events discovered by surveys such as LSST, which could then be followed up using smaller telescopes.

Additionally, other parameters such as source angular radius, distance, and extinction are required, and in the case of a lack of additional measurements, empirical methods may be necessary to estimate these values. Thus, future improvements to our results may stem from more precise measurements of the sources. Given the spectroscopic follow-up of numerous microlensing events, it is realistic to obtain distance (as for Gaia21efs) and/or extinction. While interferometry can measure the angular radius, this level of complexity would be excessive for our purposes, especially considering our inability to measure ρ and derive the lens mass.

It is worth noting that the choice of the 95% lower credible interval for defining θ_E, lim is not unique and represents a methodological choice rather than a strict standard. Different studies adopt different thresholds depending on their goals: for example, Smith et al. (2002) use 2σ, Shvartzvald et al. (2014) and Zhu et al. (2016) adopt 3σ, and Bachelet et al. (2022) use a more conservative 10σ limit. A stricter threshold may be justified when the outcome influences costly follow-up observations. In our case, however, all data have already been collected, and the main purpose is to prioritise the analysis of potentially interesting events. We therefore opt for the 95% credible interval as a practical compromise between strictness and informativeness. However, even when adopting a stricter 5σ threshold, the lower mass limits for Gaia21efs and Gaia21azb remain above 0.2 M_⊙. The limit for GaiaDR3-ULENS-067 is somewhat lower under this criterion, but its DLC results (high mass, dark lens probability, and transverse velocity) still strongly support its classification as a stellar remnant.

The compact nature of the candidate lenses can be most accurately validated through astrometric counterparts of the analysed photometric microlensing events. The astrometric series from Gaia Data Release 4 will encompass data from the first 5.5 years of the mission; therefore, only the GaiaDR3-ULENS-067 event (from 2015) can be analysed using these data. However, it will not be publicly available until mid-2026. The remaining events, occurring in 2021, can be properly analysed with Gaia Data Release 5, the full Gaia catalogue, estimated for release in 2030. However, some astrometric shifts might be detectable even in the Gaia Data Release 4 data. The accuracy of the mass measured using these data depends on factors such as the source brightness and lens mass. The most precise mass measurements are expected for bright sources lensed by massive objects, assuming average Gaia coverage. Therefore, Gaia21efs with a G-band magnitude range of 15.7^m to 12.1^m, has promising potential for accurate mass measurement.

Apart from astrometric microlensing, imaging can also be used to verify the nature of the lenses. For example, young neutron stars would be detectable in X-rays and young white dwarfs in ultraviolet, provided they are relatively nearby. Therefore, there is a chance of detecting lenses in Gaia21efs or Gaia21dpb, which are situated at a distance of ∼1.8 kpc. Conversely, since low-mass main-sequence stars emit primarily in the infrared, the lack of lens detection in this band can decisively rule out the presence of such a star. This absence would strengthen the case for the lens being a compact object, potentially with a high temperature and emitting mainly in ultraviolet or X-rays. A similar approach has previously been used, as is demonstrated in Blackman et al. (2021).

Future surveys could benefit from an adaptation of the methodology used in this study. For instance, the Rubin Observatory is anticipated to observe thousands of microlensing events with dense temporal coverage during its operation (Sajadian & Poleski 2019). If the analysis of light curves similar to those analysed in this work produces significant mass limits, it will provide a strong argument for conducting follow-up astrometric time-series observations of these events, such as with the Roman Space Telescope, thereby enabling a more efficient use of observational resources.

6. Conclusions

Gaia21efs, Gaia21azb, GaiaDR3-ULENS-067, and Gaia21dpb are microlensing events with moderately or highly magnified light curves. Despite this, they lack the FS effect, potentially due to the large Einstein radius, and hence also possibly due to the high mass of the lens. For these events, we derived lower limits for their angular Einstein radii and lens masses by constraining the ρ distribution.

Our methodology involved three key steps: modelling the light curve, determining source parameters, and conducting the DLC analysis. It is most effective for events with u₀ < 0.1 and sufficient observational data. Given the challenges associated with lens mass measurement methods, this approach provides a reliable estimate while remaining relatively simple. Therefore, it can be incorporated into follow-up optimisation efforts for large numbers of microlensing events to be detected by facilities such as Rubin and Roman.

Gaia21efs, Gaia21azb, and GaiaDR3-ULENS-067 are estimated to have lower mass limits in the range of 0.3 − 0.8 M_⊙. In contrast, the mass limit for Gaia21dpb is significantly lower, at ∼0.014 M_⊙, which minimally contributes to this analysis. Considering the mass limits and the DLC analysis, we conclude that the lenses in Gaia21efs and Gaia21azb are most likely white dwarfs with masses of ∼1 M_⊙, although their relatively high transverse velocities leave open the possibility that they could be neutron stars. Gaia21dpb is likely a white dwarf with a mass of ∼0.4 M_⊙. The lens in GaiaDR3-ULENS-067 can be a neutron star, with a mass estimated between ∼1.18 − 1.70 M_⊙ depending on the model, and a high transverse velocity of ∼300 − 320 km s⁻¹. The nature of lenses can be confirmed by measuring their masses using astrometric data from Gaia or using imaging techniques.

If confirmed, these lenses, together with stellar remnant candidates identified in other studies, could significantly improve our knowledge of the mass and spatial distribution of dark remnants. In turn, this could provide valuable insights into the evolution of massive stars and binary systems, shedding light on how supernovae have shaped the history of the Milky Way.

References

Appendix A: List of BHTOM facilities

The following table lists BHTOM facilities which collected data for this study.

Appendix B: Photometric data

The following tables present all the photometric data available for each event, detailing the number of data points per facility before and after cleaning. The tables also include the minimum and maximum MJD for the cleaned data. Facilities whose data was used in the analysis are highlighted in bold, as are the filters applied. A superscript ^B indicates that data from a facility was processed using BHTOM.

Appendix C: Photometric curve parameters

The following tables present the fitted parameters.

Appendix D: Acknowledgements

The authors would like to thank Dr Radek Poleski for valuable discussions and comments during the development of this work.

This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

We acknowledge ESA Gaia, DPAC and the Photometric Science Alerts Team (http://gsaweb.ast.cam.ac.uk/alerts).

BHTOM.space is based on the open-source TOM Toolkit by LCO and has been developed with funding from the OPTICON-RadioNet Pilot (ORP) of the European Union’s Horizon 2020 research and innovation programme under grant agreement No 101004719 (2021-2025). This project has received funding from the European Union’s Horizon Europe Research and Innovation programme ACME under grant agreement No 101131928 (2024-2028).

Based on observations obtained with the Samuel Oschin Telescope 48-inch and the 60-inch Telescope at the Palomar Observatory as part of the Zwicky Transient Facility project. ZTF is supported by the National Science Foundation under Grants No. AST-1440341 and AST-2034437 and a collaboration including current partners Caltech, IPAC, the Oskar Klein Center at Stockholm University, the University of Maryland, University of California, Berkeley, the University of Wisconsin at Milwaukee, University of Warwick, Ruhr University, Cornell University, Northwestern University and Drexel University. Operations are conducted by COO, IPAC, and UW.

This work has made use of data from the Asteroid Terrestrial-impact Last Alert System (ATLAS) project. The Asteroid Terrestrial-impact Last Alert System (ATLAS) project is primarily funded to search for near earth asteroids through NASA grants NN12AR55G, 80NSSC18K0284, and 80NSSC18K1575; byproducts of the NEO search include images and catalogs from the survey area. This work was partially funded by Kepler/K2 grant J1944/80NSSC19K0112 and HST GO-15889, and STFC grants ST/T000198/1 and ST/S006109/1. The ATLAS science products have been made possible through the contributions of the University of Hawaii Institute for Astronomy, the Queen’s University Belfast, the Space Telescope Science Institute, the South African Astronomical Observatory, and The Millennium Institute of Astrophysics (MAS),Chile.

The Pan-STARRS1 Surveys (PS1) and the PS1 public science archive have been made possible through contributions by the Institute for Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max-Planck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg and the Max Planck Institute for Extraterrestrial Physics, Garching, The Johns Hopkins University, Durham University, the University of Edinburgh, the Queen’s University Belfast, the Harvard-Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics and Space Administration under Grant No. NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, the National Science Foundation Grant No. AST-1238877, the University of Maryland, Eotvos Lorand University (ELTE), the Los Alamos National Laboratory, and the Gordon and Betty MooreFoundation.

Support for The Omega Key Project is provided by ANID’s Millennium Science Initiative through grant ICN12_009, awarded to the Millennium Institute of Astrophysics (MAS), and by ANID’s Basal project FB210003.

YT acknowledges the support of DFG priority program SPP 1992 "Exploring the Diversity of Extrasolar Planets" (TS 356/3-1). RAS and EB gratefully acknowledge support from the NASA XRP Program, through grant number 80NSSC19K0291.

This work was authored by employees of Caltech/IPAC under Contract No. 80GSFC21R0032 with the National Aeronautics and Space Administration.

The Omega Key Project has received funding from the European Union’s Horizon 2020 research and innovation program under grant agreement No. 101004719 (OPTICON - RadioNet Pilot).

This work is supported by the Polish MNiSW grant DIR/WK/2018/12. HHE thanks Türkiye National Observatories for partial support in using the TUG-T100 telescope with project number 21AT100-1799.

Eda Sonbas thanks the Adiyaman University Astrophysics Application and Research Center for their support in the acquisition of data with the ADYU60 telescope.

Josep Manel Carrasco work was (partially) supported by the Spanish MICIN/AEI/10.13039/501100011033 and by ”ERDF A way of making Europe” by the ”European Union” through grant PID2021-122842OB-C21, and the Institute of Cosmos Sciences University of Barcelona (ICCUB, Unidad de Excelencia ”María de Maeztu”) through grant CEX2019-000918-M and the project 2021-SGR-00679 GRC de l’Agència de Gestió d’Ajuts Universitaris i de Recerca (Generalitat de Catalunya). The Joan Oró Telescope (TJO) of the Montsec Observatory (OdM) is owned by the Catalan Government and operated by the Institute for Space Studies of Catalonia (IEEC).

G.D., M.S. and M.D.J acknowledge support by the Astronomical Station Vidojevica and the Ministry of Science, Technological Development and Innovation of the Republic of Serbia (MSTDIRS) through contract no. 451-03-136/2025-03/200002 made with Astronomical Observatory (Belgrade), by the EC through project BELISSIMA (call FP7-REGPOT-2010-5, No. 256772), the observing and financial grant support from the Institute of Astronomy and Rozhen NAO BAS through the bilateral SANU-BAN joint research project "GAIA astrometry and fast variable astronomical objects", and support by the SANU project F-187. Also, this research was supported by the Science Fund of the Republic of Serbia, grant no. 6775, Urban Observatory of Belgrade - UrbObsBel.

T.G. and the IST60 Telescope has been supported by the Scientific Research Projects Coordination Unit of Istanbul University (Project No.: TSG-2023-40046 and FBG-2017-23943) Turkish Republic, Directorate of Presidential Strategy and Budget project, 2016K121370.

This work is partially supported by the Fundamental Fund of Thailand Science Research and Innovation (TSRI) through the National Astronomical Research Institute of Thailand (Public Organization) (FFB680072/0269).

J.Z. and E.P. acknowledge funding from the Research Council of Lithuania (LMTLT, grant No. S-LL-24-1).

Some of the observations have been obtained with the 90cm Schmidt-Cassegrain Telescope (TSC90) in Piwnice Observatory, Institute of Astronomy of the Nicolaus Copernicus University in Toruń (Poland).

We acknowledge the use of ChatGPT-3.5 from OpenAI as a tool to improve the clarity and flow of the text.

All Tables

All Figures

	Fig. 1. Spectrum of the Gaia21efs obtained with LBT/PEPSI on November 3, 2021 (blue) and the best-matching fit (red) synthesised for the specific parameters. The Fe I regions in blue (top) and red spectral arms (bottom) are presented.
In the text

	Fig. 2. Light curve of Gaia21efs with the data from different surveys as well as the follow-up data collected with BHTOM. The solid line shows the parallax FS model, while the dashed line shows a parallax model for the point source. The right figure shows the zoom-in on the peak of the light curve.
In the text

	Fig. 3. DLC results for Gaia21efs. Top: Mass-distance plot. Bottom: Blend-lens plot.
In the text

	Fig. 4. Light curve of Gaia21azb with the data from different surveys as well as the follow-up data collected with BHTOM. The solid line shows the parallax FS model, while the dashed line shows a parallax model for the point source. The right figure shows the zoom-in on the peak of the light curve.
In the text

	Fig. 5. DLC results for Gaia21azb. Top: Mass-distance plot. Bottom: Blend-lens plot.
In the text

	Fig. 6. Light curve of GaiaDR3-ULENS-067 with the data from Gaia and OGLE. The solid line shows the parallax FS model, while the dashed line shows a parallax model for the point source. The right figure shows the zoom-in on the peak of the light curve.
In the text

	Fig. 7. Colour-magnitude diagram (CMD) for GaiaDR3-ULENS-067. The red square represents the source, the yellow circle marks the RC centroid, and the black dots denote objects within a 5′ radius of the source.
In the text

	Fig. 8. DLC results for GaiaDR3-ULENS-067. Top: Mass-distance plot. Bottom: Blend-lens plot.
In the text

	Fig. 9. Light curve of Gaia21dpb with the data from Gaia and ZTF. The solid line shows the parallax FS model, while the dashed line shows a parallax model for the point source. The right figure shows the zoom-in on the peak of the light curve.
In the text

	Fig. 10. DLC results for Gaia21dpb. Top: Mass-distance plot. Bottom: Blend-lens plot.
In the text