Binary microlensing by high eccentric stellar-mass black hole binaries

Microlensing is a powerful tool that allows us to probe dark or faint celestial objects, such as planets, dwarfs, neutron stars, and black holes (BHs), as well as massive compact halo objects by magnifying the luminosity of background sources (Paczynski, 1986; Mao & Paczynski, 1991; Gould & Loeb, 1992). Thanks to several photometric surveys conducted hitherto—such as Expérience pour la Recherche d’Objets Sombres (EROS), MAssive Compact Halo Objects (MACHO) project, Microlensing Observations in Astrophysics (MOA), Optical Gravitational Lensing Experiment (OGLE), and Korea Microlensing Telescope Network (KMTNet)—tens of thousands of microlensing events within our own and nearby galaxies have been observed (Aubourg et al., 1993; Alcock et al., 2000; Sumi et al., 2011; Udalski et al., 2015; Yock & Muraki, 2018; Byun et al., 2022).

Stellar-mass BHs are one of the possible lens objects of microlensing events. They have been identified through X-ray observations  (Liu et al., 2007; Tetarenko et al., 2016), but microlensing provides a unique opportunity to detect isolated stellar-mass BHs as well. Based on the stellar mass function, stellar-mass BHs are expected to account for approximately 1% of microlensing events  (Gould, 2000). However, only about a dozen stellar-mass BH candidates have been found from microlensing events to date, even though they are anticipated to be overrepresented in the observation due to the longer event duration resulting from their substantial masses  (Agol et al., 2002; Bennett et al., 2002; Mao et al., 2002; Poindexter et al., 2005; Dong et al., 2007; Shvartzvald et al., 2015; Wyrzykowski et al., 2016; Mróz & Wyrzykowski, 2021; Sahu et al., 2022).

As for stellar-mass BH binaries, they have not been conclusively discovered through microlensing events, which contrasts sharply with the fact that nearly 100 BH binaries have been identified through gravitational-wave (GW) observations  (Abbott et al., 2019, 2021a, 2021b, 2021c). However, we should note that this situation arises from differences in the observable spatial volume and time, rather than suggesting a scarcity of wide-orbit BH binaries that produce microlensing events. This lack of detection could be attributed to the small microlensing parallax of the BH binary due to their large mass as well as the lens’s orbital motion that complicates the evaluation of the parallax, which makes it difficult to accurately determine the lens mass (Karolinski & Zhu, 2020; Ma et al., 2022).

It is known that BH binaries can form through the evolution of massive stellar binaries (Marchant & Bodensteiner, 2023; Langer et al., 2020). Also, they can form via dynamical interactions at the center of star clusters (O’Leary et al., 2009; Banerjee et al., 2010; Bae et al., 2014; Hong & Lee, 2015; Park et al., 2017). In particular, in highly dense regions like the center of nuclear star clusters, many BHs can be gathered in a small volume through the dynamical friction and mass segregation (Hailey et al., 2018), so the number of interactions of BHs increases, which triggers the dynamical formation of BH binaries.

In this study, we consider BH binary lenses with a range of eccentricities, from low to high, including those in unbound orbits with eccentricities equal to or exceeding unity. Generally, the formation of BH binaries from massive stellar binaries tends to favor low eccentricity, while dynamical formation processes markedly favor the production of BH binaries with high eccentricities (Bae et al., 2014; Hong & Lee, 2015). Furthermore, gravitational perturbations from the supermassive BH at the center of the galaxy can amplify the eccentricity of BH binaries, known as the Kozai-Lidov mechanism (Hoang et al., 2018). Additionally, in multi-body interactions—such as binary-single or binary-binary encounters—highly eccentric encounters occur repeatedly  (Samsing, 2018; Zevin et al., 2019), commonly resulting in unbound orbits. Thus, the occurrence rate of high-eccentricity encounters, especially in cluster centers, is significant and cannot be overlooked.

In the determination of the characteristics of binary microlensing from a light curve of a luminous celestial source, well-aligned configurations between the source, binary lens, and observer ought to be satisfied, and several fundamental parameters, e.g., the distances to the source and lens, the mass of the binary lens, pericenter distance or separation between the components of a binary, the trajectory and relative velocity of the source, and so on, should be taken into account. Assuming negligible parallax for BH binary lens due to its small value as previously mentioned, in this study, we focus on the effects of orbital properties of BH binary lenses, specifically including unbound BH trajectories, that can appear on the light curve of microlensing events. We explore the signature of binary lensing depending on the eccentricity of a BH binary and compare the cases of elliptic and hyperbolic orbits with varying other parameters characterising a binary lens system, such as the pericentre distance and mass ratio between two component BHs, and the inclination angle between the orbital axis and the line of sight. Additionally, we test whether we can restore the true value of the parameters of a BH binary lens by using a model-based fitting method, the Nelder-Mead fitting (Nelder & Mead, 1965), and Bayesian inference using the Markov chain Monte Carlo sampling method (Goodman & Weare, 2010).

We find from this study that the employed parameter estimation methods can retrieve true values of the binary microlensing parameters within the $1\sigma$ uncertainty of inferred values. For some samples of which the Nelder-Mead fitting failed in the uncertainty estimation, Bayesian inference based on the Markov chain Monte Carlo method can succeed in estimating the uncertainty. On the other hand, if the median value of the posterior distribution obtained by the Bayesian inference deviates from the true value more than that from the initial Nelder-Mead fitting, it turns out that using the uncertainty estimated from the Bayesian inference helps reduce overall deviations of all parameters.

This paper is organized as follows: In Sec. 2, we summarize how to compute the orbits of two black holes and how to simulate binary microlensing phenomena that can occur by the two black holes acting like a binary lens. We investigate the parameter dependency in the binary microlensing signature and show the resulting caustics and light curve in Sec. 3. We provide the result of parameter estimation for the light curve of a simulated binary microlensing event in Sec. 4. In Sec. 5, we discuss our result and the future prospect of observing binary microlensing events by stellar-mass binary black holes within our and/or nearby galaxies.

In this work, we assume observing a binary lensing event is made with a constant time interval for a total observation time. Following the assumption, we interpolate the trajectory of orbit obtained by solving the equation of motion for the given observation time so that the time steps of the trajectory are equally spaced. Consequently, we can simulate a binary lensing light curve with $\sim 1\leavevmode\nobreak\ \textrm{day}$ cadence. In addition, in order to ensure a sufficient observation time for identifying the signature of binary lensing, i.e., two sharp spikes connected with magnified ‘U’-shape light curve, we take twice of the Einstein cross time, calculable for the Einstein radius of the total mass of a BH binary, as the total observation time.

Regarding the consideration discussed thus far and denoting a model parameter of this work as $\lambda$, we configure a set ${\Lambda}$ as follows

where $D_{S}$ and $D_{L}$ are the distance to source and BH binary lens, $v_{r}$ is the relative velocity of source against to the lens, $\theta$ is the angle between the real axis $z_{x}$ on the lens plane and the source trajectory projected onto the lens plane, $m_{1}$ and $m_{2}$ are the mass of component BHs of a binary lens, $s$ is the pericenter distance between two BHs, $e$ is the eccentricity of the orbital motion of BH binary, and $\iota$ is the inclination angle defined by the angle between the orbital axis and the line of sight. In order to focus on the original and projected orbital motion of binary BHs, we fix $D_{S}=8\leavevmode\nobreak\ \textrm{kpc}$, $D_{L}=5\leavevmode\nobreak\ \textrm{kpc}$, $v_{r}=200\leavevmode\nobreak\ \textrm{km/s}$, $\theta=120^{\circ}$, and alter $m_{1}$, $m_{2}$, $s$, $e$, and $\iota$.

We demonstrate a fitting-based parameter estimation of binary microlensing events. In practice, a successful fitting of a given light curve data with a binary microlensing model can be a primary approach for estimating the most promising values of microlensing parameters from the value of model parameters used in the fitting. We adopt the Nelder-Mead (NM) method (Nelder & Mead, 1965) as the fitting method of this work.

We perform the fitting with the prior range tabulated in Table 1 and summarize the result obtained using the method in a Python package, lmfit (Newville et al., 2016), in Table 2 for simulated samples prepared with the combination of different values of selected parameters $q$, $s$, $e$, and $\iota$ including the samples studied in the previous section. Here, we test following parameter values, $q=\{0.3,1.0\}$, $s=\{10\leavevmode\nobreak\ \textrm{AU},20\leavevmode\nobreak\ \textrm{AU}\}$, $e=\{1.0,1.5,2.0,2.5,3.0\}$, and $\iota=\{5^{\circ},15^{\circ},30^{\circ}\}$.

As tabulated in Table 2, in order to measure the goodness-of-fitting, we compute the reduced $\chi^{2}$, i.e., $\chi^{2}_{r}$, which is defined as

where $\rho_{i}$ is the residual of magnitude subtracting the inferred magnitude $d({\Lambda})_{i}$ from the true magnitude $d_{T,i}$ and scaled by the data uncertainty $\sigma_{i}$, i.e.,

at each observation point $i$. In this work, we adopt the same form of $\sigma_{i}$ used in Ma et al. (2022) which is given by

where $\sigma_{\textrm{sys}}=0.004$ is the systematic floor and $\sigma_{0}=0.04$ is the uncertainty at the reference magnitude $d_{0}$ set as 18. $N$ and $N_{\textrm{vars}}$ denote the number of observation data points and the number of variables that are fixed as $400$ and $9$ in this work, respectively. Note that we take $\sum_{i}\rho_{i}$ as the cost function of the NM method for seeking the best fitting parameters. Along with $\chi_{r}^{2}$, for the figure-of-merit of the fitting-based parameter estimation, we calculate the averaged relative difference $\left<\delta_{\lambda}\right>$ as follows

where $\delta_{\lambda}$ is the relative difference between the true value $\lambda_{T}$ and inferred value $\lambda_{I}$ of each parameter $\lambda$, i.e., $\delta_{\lambda}=|\lambda_{T}-\lambda_{I}|/\lambda_{T}$, and summarize the value for each sample in Table 2.

Looking at Table 2, we can see that the 9 parameters of binary microlensing are mostly well recovered by the NM-based fitting. Also, the goodness-of-fit, represented by the value of $\chi^{2}_{r}$, of all tested samples results in around $1$ which means the extent of residuals is in accord with given degrees of freedom, i.e., $(N-N_{\textrm{var}})$, in the fitting. When we compare the range of $\left<\delta_{\lambda}\right>$ of samples in each case, the range of $\left<\delta_{\lambda}\right>$ of Case I, resulted in $\sim\!10^{-3}$–$10^{-2}$, are relatively larger than the ranges of Case II and III, $\sim\!10^{-4}$–$10^{-3}$ and $\sim\!10^{-4}$–$10^{-2}$, respectively. This result tells us that relatively smaller/narrower caustics of shorter $s$ such as $s=10\leavevmode\nobreak\ \textrm{AU}$ hinders precise estimations than longer $s$ such as $s=20\leavevmode\nobreak\ \textrm{AU}$.