RESOLVING MICROLENSING EVENTS WITH TRIGGERED VLBI

To generate mock microlensing events we require the spatial and velocity distributions of the stellar targets and lenses.

For the structure of the Galaxy, we use the thin disk model and standard bulge model from Binney & Tremaine (1987). The density distribution in the disk is modeled in cylindrical coordinates by a double exponential function,

$$\begin{eqnarray}&&{\rho }_{D}(R,z)=\displaystyle \frac{{\rm{\Sigma }}}{2H}\exp \left(\displaystyle \frac{-(R-{R}_{\odot })}{{R}_{d}}\right)\exp \left(\displaystyle \frac{-| z| }{H}\right),\end{eqnarray} \tag{ 4 }$$

where Σ is the column density of the disk at the Sun's position, H the height scale and R_d the length scale of the disk. The distribution of the lens transverse velocity with respect to the line of sight is established from solar motion and the local lens velocity distributions (see Table 1).

Table 1.  Assumed Parameters of the Galactic Model; the Model is Partly Adopted from (Rahvar et al. 2003)

Structure	Parameter		Value
	${\rm{\Sigma }}\ ({M}_{\odot }\ {\mathrm{pc}}^{-2})$		50
	$H\ (\mathrm{kpc})$		0.325
Disk	${R}_{d}\ (\mathrm{kpc})$		3.5
	velocity	${\sigma }_{r}\,(\mathrm{km}\,{{\rm{s}}}^{-1})$	34.
	disper-	${\sigma }_{\theta }\,(\mathrm{km}\,{{\rm{s}}}^{-1})$	28.
	sions	${\sigma }_{z}\,(\mathrm{km}\,{{\rm{s}}}^{-1})$	20.
	${M}_{B}\ ({M}_{\odot })$		$1.7\times {10}^{10}$
	$a\ (\mathrm{kpc})$		1.49
Bulge	$b\ (\mathrm{kpc})$		0.58
	$c\ (\mathrm{kpc})$		0.40
	velocity	$\sigma \,(\mathrm{km}\,{{\rm{s}}}^{-1})$	110
	${\rho }_{0}\ ({M}_{\odot }\,{\mathrm{pc}}^{-3})$		$0.932\times {10}^{-5}$
Spheroid	${a}_{c}\ (\mathrm{kpc})$		0.5
	velocity	$\sigma \,(\mathrm{km}\,{{\rm{s}}}^{-1})$	120
	${\rho }_{H\odot }\ ({M}_{\odot }\,{\mathrm{pc}}^{-3})$		0.008
	${R}_{c}\ (\mathrm{kpc})$		5.0
Halo	$M\ \mathrm{in}\ 60\ \mathrm{kpc}\ ({10}^{10}{M}_{\odot })$		51
	velocity	$\sigma \,(\mathrm{km}\,{{\rm{s}}}^{-1})$	200

Download table as:  ASCIITypeset image

The bar is described in a Cartesian frame positioned at the Galactic center with the major axis x tilted by $\phi =45^\circ$ with respect to the Galactic center–Sun line. The bar density is given by

$$\begin{eqnarray}{\rho }_{B} & = & \displaystyle \frac{{M}_{B}}{6.57\pi {abc}}{e}^{-{r}_{s}^{2}/2}\\ {r}_{s}^{4} & = & {\left[{\left(\displaystyle \frac{x}{a}\right)}^{2}+{\left(\displaystyle \frac{y}{b}\right)}^{2}\right]}^{2}+\displaystyle \frac{{z}^{4}}{{c}^{4}},\end{eqnarray} \tag{ 5 }$$

where M_B is the bulge mass, and a, b and c are the scale length factors. The dispersion velocity of bulge stars is $110\,\mathrm{km}\,{{\rm{s}}}^{-1}$. For the density of the spheroid structure we take the following function (Robin et al. 2003),

$$\begin{eqnarray}{\rho }_{\mathrm{spher}}=\left\{\begin{array}{ll}{\rho }_{0}{({a}_{c}/8.5\mathrm{kpc})}^{-2.44} & a\leqslant {a}_{c}\\ {\rho }_{0}{(a/8.5\mathrm{kpc})}^{-2.44} & a\geqslant {a}_{c}\ \end{array}\right.\end{eqnarray} \tag{ 6 }$$

${a}^{2}={R}^{2}+{z}^{2}/{(0.76)}^{2}$, a_c, ${\rho }_{c}$ and spheroid velocity dispersion are given in Table 1.

Combining the density distribution of the Galaxy, kinematics, mass function and stellar population of stars, we can generate microlensing events in terms of t_E, M_l, D_l, D_s, ${\theta }_{E}$ and π, where π is the parallax parameter described in the Appendix. We use the detection efficiency of OGLE survey in terms of Einstein crossing time (i.e., $\epsilon ({t}_{E})$) to select the observed microlensing events in our simulation (Wyrzykowski et al. 2015). A general overview about the simulation of microlensing events can be found in Rahvar (2015).

In addition to the known stellar populations, gravitational lens candidates include the Galactic population of compact objects. In principle, these may be comprised of the remnants of stellar evolution and primordial objects. We make the conservative assumption that the latter are absent, and consider only the remnants of massive stars.

The end point of stellar evolution depends primarily on the mass of the progenitor. We consider two classes of remnants, white dwarfs and black holes. We ignore intermediate neutron stars, which we expect to contribute marginally to the high-mass tail of the white dwarf population and otherwise be indistinguishable from them using microlensing observations alone. Hence in practice the remnants are effectively white dwarfs/neutron stars and black holes.

Encoded in the remnant population is the Galactic high-mass star formation history. Once the relationship between the zero-age main sequence progenitor mass (${M}_{\mathrm{MS}}$) and the final remnant mass (M_R, where R may be WD or BH) is specificed, the mass function⁷ ($\phi ({M}_{R})$) is given solely in terms of the star formation rate (${\rm{\Lambda }}(t)$) and initial mass function of main sequence stars (${\rm{\Phi }}({M}_{\mathrm{MS}})$):

$$\begin{eqnarray}&&\phi ({M}_{R})=\int {dt}\displaystyle \frac{{\rm{\Lambda }}(t-{T}_{\mathrm{MS}}){\rm{\Phi }}({M}_{\mathrm{MS}})}{d\,{\mathrm{log}}_{10}\,{M}_{R}/d\,{\mathrm{log}}_{10}\,{M}_{\mathrm{MS}}},\end{eqnarray} \tag{ 7 }$$

where ${T}_{\mathrm{MS}}\approx 10{({M}_{\mathrm{MS}}/{M}_{\odot })}^{-1.5}\,\mathrm{Gyr}$ is the lifetime of a main sequence star with mass ${M}_{\mathrm{MS}}$ associated with the remnant mass M_R, and the integration is over the entire Galactic history. Effectively, this is simply the number of remnants formed over the history of Galactic star formation, excluding those stars that remain on the main sequence. For this we adopt the star formation rate from Hartwick (2004), which peaks roughly 10 Gyr ago, and the "Disk and Young Clusters" stellar initial mass function ($\mathrm{IMF}(M)$) from Chabrier (2003), though employing the "universal" IMF from Kroupa (2001) makes no discernible difference.

We relate ${M}_{\mathrm{WD}}$ and ${M}_{\mathrm{MS}}$ via the piecewise continuous expression presented in Salaris et al. (2009), obtained for white dwarfs in open clusters by comparing the white dwarf ages inferred by cooling models and cluster age from isochrone fitting:

$$\begin{eqnarray}&&{M}_{\mathrm{WD}}/{M}_{\odot }\\ =\left\{\begin{array}{ll}0 & {m}_{\mathrm{MS}}\leqslant 0.5588\\ 0.01{m}_{\mathrm{MS}}+0.5418 & 0.5588\lt {m}_{\mathrm{MS}}\leqslant 1.7\\ 0.134{m}_{\mathrm{MS}}+0.331 & 1.7\lt {m}_{\mathrm{MS}}\leqslant 4\\ 0.047{m}_{\mathrm{MS}}+0.679 & 4\leqslant {m}_{\mathrm{MS}}\lt 8\\ 0 & 8\leqslant {m}_{\mathrm{MS}},\end{array}\right.\end{eqnarray} \tag{ 8 }$$

where ${m}_{\mathrm{MS}}\equiv {M}_{\mathrm{MS}}/{M}_{\odot }$.

The relationship between ${M}_{\mathrm{BH}}$ and ${M}_{\mathrm{MS}}$ is highly dependent upon the uncertain evolution of massive stars and the details of core-collapse supernovae. As a result, above $30{M}_{\odot }$ it is unclear what the slope of the ${M}_{\mathrm{BH}}$–${M}_{\mathrm{MS}}$ relationship is (see Figure 9 of Fryer et al. 2012). Here we make a moderately conservative assumption that the resulting black hole mass is independent of ${M}_{\mathrm{MS}}$, setting

$$\begin{eqnarray}&&{M}_{\mathrm{BH}}/{M}_{\odot }\\ =\,\left\{\begin{array}{ll}0 & {m}_{\mathrm{MS}}\leqslant 8\\ 0.3{m}_{\mathrm{MS}} & 8\lt {m}_{\mathrm{MS}}\leqslant 33\\ 10 & 33\lt {m}_{\mathrm{MS}},\end{array}\right.\end{eqnarray} \tag{ 9 }$$

which lies in the middle of the permitted range presented in Fryer et al. (2012) and is roughly consistent with narrow black hole mass range inferred from X-ray binaries by Özel et al. (2010).

The resulting stellar and remnant mass functions for the disk, bulge and spheroid are shown in Figures 1 and 2. The spikes in the low- and high-mass ends of the white dwarf and black hole mass functions are associated with regions in which the remnant mass is nearly independent of the progenitor mass, and contains a finite number of objects. For both the disk and the bulge there is a break in the stellar mass function where stars produced during the star formation peak are leaving the main sequence, roughly $1{M}_{\odot }$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Because neither white dwarfs nor black holes are expected to experience strong kicks at formation, we assume the same velocity dispersion for the remnants as the surrounding stellar population. This is notably not the case for neutron stars, which often experience sufficient kicks to be launched into the halo, diluting their density in the disk and therefore their representation in the lens sample for current and planned microlensing surveys. The observation of pulsars with ages less than $3\,\mathrm{Gyr}$, indicates the pulsar birth velocity of about $400\,\mathrm{km}\,{{\rm{s}}}^{-1}$ (Hobbs et al. 2005).

In order to simulate the microlensing events, we choose the source stars from a Hipparcos-like color–magnitude distribution (Perryman et al. 1997), distributed proportional to the density of matter. After selecting the source star with a given absolute color and magnitude, the position of the star in the CMD changes due to the distance modulus and interstellar extinction. Those stars brighter than the limiting magnitude of the microlensing survey (which is taken about 19 in the I-band in this simulation) are selected for the observation. It should be noted that the results are not sensitive to the limiting magnitude and could be applied to microlensing surveys with different limiting magnitudes.

We note that not all the microlensing events generated in this simulation are observable in reality. Depending on the duration of observation, cadence of photometric data, and non-observing nights due to bad weather and technical failures, only a fraction of events as a function of Einstein crossing time, t_E, can be observed. We adapt the detection efficiency function (i.e., $\epsilon ({t}_{E})$) for selecting the observed events from OGLE published function (Wyrzykowski et al. 2015). A detailed description of gravitational microlensing formalism can be found in Rahvar (2015).

Each event is constructed by first choosing a random source with the appropriate properties, subsequently a random lens with the appropriate properties, and finally applying the OGLE detection efficiency. That is,

• 1.  
  Source star selection:

  • (a)  
    A source star position is selected with a probability distribution given by the stellar density and geometric factors, i.e.,

    $$\begin{eqnarray}&&\displaystyle \frac{{dP}}{{{dD}}_{S}}\propto \rho ({D}_{S}){D}_{S}^{2},\end{eqnarray} \tag{ 10 }$$
  • (b)  
    The source star stellar type is randomly chosen according to the Hipparcos color–magnitude distribution.
  • (c)  
    The stellar component (bulge, disk, spheroid) is selected with weights equal to their density at D_S.
  • (d)  
    The source star velocity is the sum of the bulk Galactic motion for the chosen stellar component at the chosen source location and a random contribution pulled from a Gaussian distribution with the appropriate dispersion (listed in Table 1).
• 2.  
  Lens selection:

  • (a)  
    A lens position is selected from the stellar density, modified by the lensing cross section, i.e.,

    $$\begin{eqnarray}&&\displaystyle \frac{{dP}}{{{dD}}_{L}}\propto \rho ({D}_{L})\sqrt{\displaystyle \frac{{D}_{L}({D}_{S}-{D}_{L})}{{D}_{S}^{2}}}.\end{eqnarray} \tag{ 11 }$$
  • (b)  
    A lens component (bulge, disk, spheroid) is randomly chosen according to the local densities of each at D_L.
  • (c)  
    A lens component type (star, white dwarf, black hole) is randomly chosen according to the local fraction of each for the given component.
  • (d)  
    A lens mass is selected according to the total mass function of the chosen component modified by the approximate lensing cross section; for a fixed observing time the probability of detecting a lens is $\propto {\theta }_{E}\propto {M}^{1/2}$.
  • (e)  
    A lens velocity is the sum of the bulk Galactic motion for the chosen component and a random contribution pulled from a Gaussian distribution with the appropriate dispersion (listed in Table 1).
• 3.  
  Detection efficiency cut: Given the lens mass, source and lens velocities, and source and lens positions t_E is computed, and the event is accepted with a probability set by the OGLE detection efficiency.

The result is a set of lensing events with known intrinsic properties (source and lens masses, distances, velocities, and types) selected using a realistic microlensing survey biases.

Figure 3 shows the distribution in Einstein crossing time of simulated events. Since in the simulation we flag lenses that are either main sequence or remnants, we can classify lenses based on lens type. Here we compare the overall microlensing events with those produced by black holes specifically; the fraction of events with black hole lenses is 0.035. This fraction is comparable to the estimate of Gould (2000), which finds 0.01 assuming a 100% detection efficiency, and thus neglecting the higher detection efficiency of Mira events.

Zoom In Zoom Out Reset image size

Another important point is the average $\mathrm{log}({t}_{E})$ for the overall events and black hole events are 1.5 and 2.0, respectively. This is unsurprising since the duration of events scales with lens mass as ${M}^{1/2}$, and is an immediate consequence of the fact that the average black hole mass is roughly an order of magnitude larger than that of the overall sample. Hence, black hole lensing events are typically more than three times longer than the typical microlensing event, with important consequences for identification and parameter estimation.

Figure 4 shows the distribution of lens locations in our simulated survey for each of the Galactic components and lens types. Events are dominated by lenses within the disk (65%) with the remainder from the bulge (35%); the spheroid contributes negligibly. Of note is that the lens distribution within the disk is nearly evenly distributed between the bulge and a heliocentric distance of 2 kpc, implying that microlensing surveys necessarily probe the lens population throughout the intervening disk. This remains true independent of the lens type.

Zoom In Zoom Out Reset image size

Miras provide a natural class of optically bright sources with strong radio emission in the form of circumstellar SiO masers. In most microlensing surveys Miras have already been identified due to their variable nature. Where they have not been, they may be crudely identified in the optical CMD, shown for the OGLE III sample in Figure 5. In particular, Miras produce a bright, red cluster that only marginally overlaps with the remainder of the survey targets. From direct inspection of the OGLE III fields we estimate that as many as 5% of the objects near the center of the Mira cluster in Figure 5 are identified Miras.⁸ In Section 5 we show that there should be nearly six times more Miras in the bulge than have been detected thus far, suggesting that it might be possible to increase the fraction of Miras in this region up to 10%–15%, after accounting for the field of view of the OGLE III survey. Given their high radio brightness, confirmation that a given candidate source is a Mira can then be quickly obtained by direct radio observation with existing large radio telescopes. Equally important is that there is little chance for lensing-induced source confusion—blending of source and lens stars does not move objects appreciably within the region populated by Miras.

Zoom In Zoom Out Reset image size

Rapid intrinsic stellar variability provides an obvious impediment to the identification and characterization of microlensing events. For this reason variable stars have been generally identified and excluded from past and ongoing surveys. Thus, while there are microlensing events reported by OGLE III that penetrate the Mira cluster (see Figure 5), these are almost certainly not associated with Miras. Nevertheless, the presence of such events makes clear that apart from variability there are no intrinsic barriers to including stars lying in the Mira-cluster region of the CMD in microlensing surveys.

The variable nature of Miras can be mitigated in a number of ways. Typical Miras vary with periods of $\approx 1\,\mathrm{year}$, I-band variability amplitudes between $1.4\,\mathrm{mag}$ and $3\,\mathrm{mag}$ (Soszyński et al. 2013), and V-band variability amplitudes typically two–three times larger (Kholopov 1990). However, black hole microlensing events can be easily distinguished from the intrinsic variability as a result of three important differences. All of these may be exploited in part due to the high luminosity of Miras, and therefore the high photometric accuracy with which they may be monitored.

First, microlensing light curves have a well-understood structure set by the lens–source geometry and nature of gravitational lensing, characterized by a divergent rise as the source approaches the lens caustics (moderated for the largest magnification events by the finite source size; see Section 4.2 and Appendix A.1). This differs from the intrinsic variations of Miras. Thus careful light curve fitting should be able to rapidly distinguish between microlensing and intrinsic variability.

Second, the amplitude of the flux variation during microlensing events is typically larger than that due to intrinsic variability. Typical magnifications are of order 200%, while high-magnification events can exceed 1000%. These correspond to magnitude variations of 0.75 and $2.5\,\mathrm{mag}$, comparable the that caused by intrinsic variations in the I-band. Hence, the microlensing signal is an order unity modification of the underlying source variability in the I-band. At longer wavelengths the intrinsic variability decreases further, on average variations in the magnitude at 1.25 μm are 20% of those at optical wavelengths (Smith et al. 2002).⁹

Third, the flux variations due to microlensing are achromatic. Again this is in stark contrast to the intrinsic variability of Miras, which varies dramatically with wavelength. Thus, multi-band monitoring should be readily able to disentangle the two variability components.

We make no further attempt to assess the efficiency with which Mira microlensing events can be identified, i.e., we presume that the intrinsic variability of Miras does not result in a significant inefficiency. As a result, our rate estimates may be considered to be optimistic.

The OGLE III survey monitored approximately 150 million stars in the Galactic bulge from 2001 to 2009 in search for gravitational microlensing events and provided us with the largest catalog of microlensing events available. OGLE used the EWS to detect ongoing microlensing events (Udalski et al. 1994). The EWS uses difference image analysis photometry and, by analyzing all the stars in the field, flags the potential microlensing events (Alard & Lupton 1998). In doing so it currently identifies and filters out variable stars. The flagged stars go through several tests including a visual inspection of the light curve and after satisfying all the criteria they are announced as a microlensing candidate. Figure 5 shows the CMD for OGLE III microlensing events (Wyrzykowski et al. 2015). The lensing events have not gone through blending correction for consistency reasons since it cannot be done for non-lensed stars.

Miras are a few orders of magnitude brighter than M-dwarfs, increasing their representation in microlensing surveys. An advantage of using bright stars is that blending by background and lens stars does not affect the brightness of the source star significantly (Figure 5). Moreover, from the lens equation, and estimating the mass of the lens star, we can calculate the contribution of blending and correct the position of source star in the CMD (Miyake et al. 2011; Muraki et al. 2011).

Within the OGLE III catalog there are currently 6528 Miras (Soszyński et al. 2013). The larger number of Miras inferred in the bulge suggest that this can be readily enlarged by surveying the entire bulge region. A rough estimate of the number of Miras in the bulge can be made via the infrared surface brightness.

We follow the approach used in Matsunaga et al. (2005). The method is based on using the empirical relation between Mira number density and infrared surface brightness to infer the total number of Miras in the bulge.

We employ the infrared brightness measurements from the Diffuse Infrared Background Experiment (DIRBE) aboard the Cosmic Background Explorer (COBE) satellite. DIRBE made full-sky brightness maps in ten infrared bands ranging from 1.25 to $240\,\mu {\rm{m}}$. For our purpose we make use of the zodical-light-subtracted infrared mission-average map at $2.2\,\mu {\rm{m}}$. The maps are available online at http://lambda.gsfc.nasa.gov/product/cobe/dirbe_prod_table.cfm.

Before using the DIRBE brightness measurements they must be corrected for dust extinction, and the contribution from the Galactic disk removed. Following Matsunaga et al. (2005) we correct for dust using an empirical relationship between the K-band extinction A_K and infrared colour:

$$\begin{eqnarray}&&{A}_{K}=0.73\times (-2.5\mathrm{log}({I}_{1.25}/{I}_{2.2})+0.14),\end{eqnarray} \tag{ 12 }$$

where ${I}_{1.25}$ is the $1.25\,\mu {\rm{m}}$ zodiacal-light-subtracted mission-average DIRBE flux.

To remove the Galactic disk contribution we make use of fits to the infrared brightness maps outside the bulge. Matsunaga et al. (2005) adopted an exponential function in Galactic longitude to model the disk

$$\begin{eqnarray}&&I(l,b)=I(0,b)\exp (-| l| /{l}_{0}(b)),\end{eqnarray} \tag{ 13 }$$

where the scale-height in Galactic longitude ${l}_{0}(b)$ is a function of Galactic latitude. This function was fit to the $2.2\,\mu {\rm{m}}$ DIRBE maps at high Galactic latitudes, $10^\circ \lt | l| \lt 45^\circ$. Motivated by the exponential vertical structure of the Galactic disk, we also estimated the disk contribution by fitting an exponential function in Galactic latitude. The resulting bulge brightness estimates, and thus the number of bulge Miras, are insensitive to the fitting function used.

As shown in Figure 6, the number of Miras in a given field is strongly correlated with the corrected infrared surface brightness. Like Matsunaga et al. (2005) we find a high-quality linear relationship between the Mira number counts in OGLE II fields and the corrected DIRBE $2.2\,\mu {\rm{m}}$. This remains true when Mira number counts in OGLE III fields are added. Linear fits to the latter gives an updated ratio of Mira number density to $2.2\,\mu {\rm{m}}$ surface brightness of $4.3\,\mathrm{Sr}/\mathrm{MJy}$, similar to that found by Matsunaga et al. (2005). With a total $2.2\,\mu {\rm{m}}$ flux of $0.9\mathrm{MJy}$ in the bulge ($-10^\circ \lt l\lt 10^\circ$ and $-10^\circ \lt b\lt 10^\circ$), this implies that there are around 4 × 10⁴ bulge Miras.¹⁰

Zoom In Zoom Out Reset image size

At least two improved microlensing surveys are underway, exceeding OGLE III in three ways: higher cadence, deeper magnitude limits, and increased sky coverage. OGLE IV, in operation since 2010, has a bulge field of $\approx 330\,{\deg }^{2}$, 3.3 times larger than that of OGLE III. The resulting OGLE IV EWS event rate is three times larger than that of OGLE III; OGLE III EWS detected about 650 events in 2008 whereas OGLE IV EWS detected almost 2000 events in 2013 and over 2000 in 2014. The forthcoming KMTNet consists of three dedicated $1.6\,{\rm{m}}$ telescopes located in Australia, South Africa, and South America, each with a $4\,{\deg }^{2}$ field of view. Thus, KMTNet anticipates detecting $\approx 2300$ microlensing events per year with a constant cadence of almost 10 minutes (Henderson et al. 2014).

The high luminosity of Miras and the long timescales of black hole microlensing events imply that the first two, higher cadence and deeper magnitude limits, will at best produce modest improvements in the Mira microlensing rate. High cadences will vastly over-sample black hole microlensing light curves, and thus while important for planet searches, are unlikely to be helpful for the study of Galactic remnants.

Even within the extincted region near the Galactic center Miras remain visible due to their high intrinsic luminosity. Since optically dim stars do not produce strong radio emission (stellar photospheres of dwarf stars are essentially undetectable with current radio telescopes), expanding the survey sample to less luminous objects provides little value for the study of black holes. Nevertheless, the deeper magnitude limits of future surveys will permit the inclusion of regions of high extinction, i.e., dark areas of the bulge, in survey fields (Schultheis et al. 2009), marginally increasing the number of Miras that can be monitored. This may be ameliorated by operating microlensing surveys in the infrared as well.

Much more important is the expanded sky coverage. We have already assumed that the entire Galactic bulge will be monitored in our previous rate estimate. Further increasing the number of monitored Miras substantially in an optical survey necessarily requires expanding the surveys beyond the bulge. Moderate modifications to the observing strategies of existing optical surveys present an obvious way in which radio-bright microlensing event rates can be substantially increased. Because Miras are intrinsically luminous, they lie well above the detection threshold of all ongoing microlensing surveys located anywhere within the Galaxy. Therefore, the relative paucity of Miras does suggest an alternative strategy: broad and shallow instead of narrow and deep.

Including the Galactic disk increases the number of stellar targets by a factor of ≈6–7. Assuming the stellar population does not vary widely between the bulge and the disk this would produce a commensurate growth in the Mira sample. Hence, we expect $\approx 2\mbox{--}3\times {10}^{5}$ Miras in total within the Milky Way. Unfortunately, the microlensing optical depth in the disk is roughly a third of that of the bulge as a result of the former's smaller stellar density. Therefore, after taking into account the disparity in optical depth, the net increase in Mira microlensing events is reduced to a factor of 2–3. As a result, we estimate that a survey that monitors the entire Galactic plane to an I-band magnitude limit of 15 would capture radio-bright microlensing events at a rate of $2\,{\mathrm{yr}}^{-1}$, corresponding to a radio-bright black hole microlensing events at a rate of $0.1\,{\mathrm{yr}}^{-1}$.

Importantly, these survey strategies need not be exclusive. Since the objects of primary interest for radio imaging, black holes, result in long-duration microlensing events, even sparse time sampling is sufficient (on the order of once per week). The OGLE IV survey has already adopted such a strategy, monitoring the entire Galactic disk with a cadence of 1–2 days and the bulge with a cadence $\lesssim 1\,\mathrm{hr}$ (Udalski et al. 2015).

The Large Synoptic Survey Telescope (LSST) is an 8.4 m telescope equipped with a wide-field camera with 9.6 square degrees field of view. It will provide a deep survey of the southern sky (over 20,000 square degrees) giving 1000 exposures of each patch during 10 years of observations. The optimal cadence of the survey for different areas of the sky is yet to be determined. The single exposure depth in the r-filter will be $\approx 24.5\,\mathrm{mag}$ and it can yet go deeper ($\approx 26.5\,\mathrm{mag}$) by co-adding the images. Being a ground-based survey, the angular resolution is limited by seeing (07) and ultimately by its pixel size (02) (Ivezic et al. 2008).

The mean cadence of 3–4 days is much shorter than the typical timescales of the long-duration black hole lensing events. Furthermore the deep survey combined with the fact that Mira variables are bright stars guarantees that most of the Mira lensing events in the Galaxy, LMC, and SMC will be detected by the LSST. Thus, for finding black hole lensing events amenable to radio imaging, the LSST presents a similar capability to current generation microlensing surveys operated with modified observation strategy.

Were a microlensing survey performed in the radio directly, the need to identify optically bright radio counterparts would be eliminated altogether, resulting in the finding of events that are amenable to radio imaging with a near 100% efficiency. Moreover such a radio survey would permit monitoring the much more numerous compact continuum radio sources. With the advent of a number of rapid-survey radio telescopes (e.g., the Square Kilometre Array, Taylor 2007; Canadian Hydrogen Intensity Mapping Experiment, Bandura et al. 2014, etc.), enabling rapid nearly all-sky surveys, high-cadence radio transient searches will become common. A radio microlensing survey is a natural byproduct of these.

The stringent requirements on source size ($\lesssim 1\,\mathrm{mas}$) restrict the potential radio microlensing survey source targets to stars (e.g., Miras) or primarily extragalactic objects. The details of the lens parameter estimation are only weakly dependent on the distance to the radio source; placing sources at extragalactic distances decreases the Einstein angle by at most 30%.

More important is the overall number of sources, which directly translates into the predicted rate enhancement. These, in turn, depend on the flux limit and wavelength of interest, both set by the VLBA. For a $10\,{M}_{\odot }$ black hole lens multiple image components are easily visible by the VLBA for frequencies above 3 GHz, providing a natural upper limit of 10 cm on the observation wavelength. For 10 minute integration times with 500 MHz bandwidths the VLBA flux limit should be near 20 μJy and 40 μJy at 22 GHz and 43 GHz, respectively; where required we adopt a flux limits of 30 μJy.¹¹ Large-amplification events will reduce the effective flux limit further. As we will see below, either improvements in the sensitivity of the VLBA or new radio-bright source classes will be necessary to leverage a radio microlensing survey.

Radio-bright AGNs provide a natural class of sources, being both numerous and, more importantly, compact. Conveniently, they are also typically flat-spectrum radio sources, meaning that neither the source identification nor the subsequent microlensing survey need be performed at the same wavelengths at which the events are ultimately imaged. Thus, to estimate the number of potential sources above the VLBA flux limit we employ the NRAO VLA Sky Survey at 1.4 GHz reported in Condon et al. (1998) (though see also Condon 1984; Condon et al. 2012) finding roughly $1.6\times {10}^{7}$ and $1.5\times {10}^{8}$ radio-bright AGNs in ellipticals and spirals, respectively, with fluxes above above 30 μJy.

These are distributed isotropically, and thus a large fraction of AGNs will be visible at high Galactic latitude with correspondingly smaller lensing optical depth. As a result, the average lensing optical depth for an AGN survey to that for the optical surveys that monitor the Galactic bulge alone is approximately 0.01. Thus, despite having nearly 290 times the number of targets, a radio microlensing survey of objects above 30 μJy would produce a black hole microlensing event rate 1 ${\mathrm{yr}}^{-1}$, about the same as that from Mira-based optical surveys.

To obtain an event rate of 10 ${\mathrm{yr}}^{-1}$, doubling the number of known black holes in two years, requires a flux limit of 3 μJy. This is would require corresponding improvements in the VLBA. However, in principle this may be achieved via an extended integration time combined with an expanded bandwidth. Currently, 4 GHz bandwidths supported by 16 Gbps recorders are planned stations participating in millimeter-wavelength VLBI observations (Whitney et al. 2012), and could be deployed to VLBA stations for use at centimeter wavelengths. Combined with 2 hr integration times, these reach the 3 μJy flux limit needed. Below 3 μJy the number of sources scale approximately as ${N}_{\gt S}\propto {S}^{-0.5}$, and thus further growth in the number of sources is a slow function of flux limit.

SiO masers at $\approx 43\,\mathrm{GHz}$ ($\approx 7$ mm) have routinely been detected around late-type giants in the central parsec of the Galaxy (Reid et al. 2003). Located within the extended atmospheres of Miras, these are typically within 8 au of their parent, corresponding to an astrometric offset of 1 mas in the Galactic center. Note that this is comparable to the Einstein angle, ${\theta }_{E}$, and therefore the lensing of the maser emission is distinct from that of optical emission of the star due to the finite-size effect of the source star (Witt & Mao 1994).

The masing spots are resolved both because of their intrinsic structure and as a result of an interstellar scattering screen that scatter-broadens images of the Galactic center (Backer 1978; Bower et al. 2006; Gwinn et al. 2014). Because of the latter effect, even point sources exhibit an extended source structure with a full width at half maximum (FWHM) of $\approx 0.7$ mas, effectively limiting the longest baselines that can be employed by the VLBA. However, this scattering is highly localized, occurring only for sources in the immediate vicinity of the Galactic center, and unlikely to limit the resolution attainable by efforts to image Miras throughout the Galactic bulge.

Typical intrinsic sizes for nearby SiO individual maser spots are 1 au, corresponding to $\approx 0.1$ mas at $\approx 8\,\mathrm{kpc}$ distance, and therefore also unlikely to significantly limit the resolution of VLBA observations. We ignore the impact of intervening scatter broadening and model the intrinsic emission from an individual spot by a Gaussian intensity profile with a FWHM of 0.1 mas:

$$\begin{eqnarray}&&{I}_{\mathrm{int}}(\beta )={I}_{0}{e}^{-{| \beta | }^{2}/2{\sigma }_{\mathrm{spot}}^{2}},\end{eqnarray} \tag{ 14 }$$

where ${\sigma }_{\mathrm{spot}}=0.1\,\mathrm{mas}/\sqrt{8\mathrm{ln}2}=0.042\,\mathrm{mas}$.

In practice, SiO masers from Mira variables often form arc-like structures that are dominated by a handful of individual masing spots. For black hole lenses the angular sizes of these structures are typically smaller than ${\theta }_{E}$, and thus multiple spots are likely to be strongly lensed simultaneously. For less massive lenses the ring angular angular scale and ${\theta }_{E}$ may be more similar. Nevertheless, because the emission between maser spots remains incoherent.¹² the resulting lensed image is a linear superposition of a number of individual spots. Thus, here we consider the simpler problem of a single masing spot to assess the size of the constraints that can be placed on the lens in principle.

Net velocity offsets between the star and masing spots are systematically incorporated into the radio source velocities and thus do not present an additional systematic uncertainty. However, the masing spots can also evolve in size and luminosity. Both are unlikely to produce substantial complications since they do not impact the separation of the lense images (see Section 4.4). While the latter can complicate the determination of the radio light curve, this may be ameliorated via the optical light curve.

Time sequences of mock images are generated via a multi-step process: beginning with the specification of the positions of the source, lens, and observer, mapping of the source to the image plane via the thin-lens equation, and convolution with a realistic beam.

Over the duration of a lensing event, the positions of the source and lens are assumed to evolve with a fixed velocity:

$$\begin{eqnarray}&&{{\boldsymbol{x}}}_{S}={{\boldsymbol{x}}}_{S,0}+{{\boldsymbol{v}}}_{S}t\quad \mathrm{and}\quad {{\boldsymbol{x}}}_{L}={{\boldsymbol{x}}}_{L,0}+{{\boldsymbol{v}}}_{L}t,\end{eqnarray} \tag{ 15 }$$

respectively. In practice only the transverse motion is important. In contrast, the observer position, i.e., that of the Earth, is orbital, and thus includes the orbital acceleration:

$$\begin{eqnarray}&&{{\boldsymbol{x}}}_{\oplus }={{\boldsymbol{x}}}_{\oplus ,0}+\int {dt}\,{{\boldsymbol{\Omega }}}_{\oplus }\times {{\boldsymbol{r}}}_{\oplus }.\end{eqnarray} \tag{ 16 }$$

The lens geometry is then fully defined by the distances

$$\begin{eqnarray}&&{D}_{L}\equiv | {{\boldsymbol{x}}}_{L}-{{\boldsymbol{x}}}_{\oplus }| \quad \mathrm{and}\quad {D}_{S}\equiv | {{\boldsymbol{x}}}_{S}-{{\boldsymbol{x}}}_{\oplus }| ,\end{eqnarray} \tag{ 17 }$$

and the transverse angular displacement

$$\begin{eqnarray}&&\beta \equiv \displaystyle \frac{({{\boldsymbol{x}}}_{S}-{{\boldsymbol{x}}}_{\oplus })}{{D}_{S}}-\displaystyle \frac{({{\boldsymbol{x}}}_{S}-{{\boldsymbol{x}}}_{\oplus })}{{D}_{S}}\cdot \displaystyle \frac{({{\boldsymbol{x}}}_{L}-{{\boldsymbol{x}}}_{\oplus })}{{D}_{L}}\displaystyle \frac{({{\boldsymbol{x}}}_{L}-{{\boldsymbol{x}}}_{\oplus })}{{D}_{L}}.\end{eqnarray} \tag{ 18 }$$

Neglecting the acceleration in the Earth's motion and the line-of-sight motion of the source and lens permits a simplification of $\beta$ to a linear function of time, though here we will make use of the more general expression in Equation (18).

For a point-mass lens the observed position on the sky, ${\boldsymbol{\theta }}$ is given by the thin-lens equation. Usually this is expressed as a potentially multi-valued function of the projected source position. However, we need only identify the projected source position in terms of the observed position ${\boldsymbol{\theta }}$, given by

$$\begin{eqnarray}&&\beta ={\boldsymbol{\theta }}\left(1-\displaystyle \frac{{\theta }_{E}^{2}}{{| {\boldsymbol{\theta }}| }^{2}}\right).\end{eqnarray} \tag{ 19 }$$

Since surface brightness is conserved, this then immediately defines the lensed image in terms of the intrinsic

$$\begin{eqnarray}&&{I}_{\mathrm{lens}}({\boldsymbol{\theta }})={I}_{\mathrm{int}}[\beta ({\boldsymbol{\theta }})].\end{eqnarray} \tag{ 20 }$$

Finally, we convolve the lensed image with a realistic, if not pessimistic, beam. For this we assume a 43 GHz beam, $B({\boldsymbol{\theta }})$ is an anisotropic Gaussian, consistent with the beam in Lu et al. (2011): semiminor axis of 0.5 mas and semimajor axis of 1.4 mas, oriented 12° east of north. The large aspect ratio is a result of a combination of the low declination of the Galactic center and the North American location of the VLBA antennae. In practice, this resolution is lower than may be achieved by the VLBA at 43 GHz; in Reid et al. (2003) the interstellar scatter broadening limited the size of the array that could be effectively employed. Outside of the Galactic center it should be possible to increase the resolution by a factor of 2–3, though we adopt the Lu et al. (2011) beam. The resulting observed intensity distribution is then

$$\begin{eqnarray}&&{I}_{\mathrm{obs}}({\boldsymbol{\theta }})=\int {d}^{2}\theta ^{\prime} B({\boldsymbol{\theta }}-{\boldsymbol{\theta }}^{\prime} ){I}_{\mathrm{lens}}({\boldsymbol{\theta }}^{\prime} ),\end{eqnarray} \tag{ 21 }$$

where prior to the convolution the angular positions are converted into equatorial coordinates.

All that remains is to specify the relative positions, velocities, and mass of the source and lens. We do this for a handful of illustrative examples in the following sections.

As fiducial source parameters we assume ${D}_{S}=8\,\mathrm{kpc}$ and a proper velocity of $120\,\mathrm{km}\,{{\rm{s}}}^{-1}$, consistent with sources in the Galactic bulge. The fiducial lens parameters are ${D}_{L}=4\,\mathrm{kpc}$ and a proper velocity of $30\,\mathrm{km}\,{{\rm{s}}}^{-1}$, consistent with the velocity dispersion in the Galactic disk. The relative positions of the source and lens are chosen so that minimum impact parameter is 0.3 mas oriented in declination. These result in typical images that are illustrative of the general situation. All that remains is to specify the mass of the lens.

The first case we consider is that of a $10{M}_{\odot }$ black hole, for which ${\theta }_{E}=3.2$ mas. A time-sequence of images for our fiducial case is shown in Figure 7, covering nearly a year. Since the minimum impact parameter is small in comparison to ${\theta }_{E}$, multiple images are present. Because ${\theta }_{E}$ is also considerably larger than the intrinsic source and beam sizes, these are well resolved, suggesting that detecting and interpreting the features of strong lensing by black holes will be straightforward. In particular, a direct measurement of ${\theta }_{E}$ is possible simply by fitting the spot separations, though we defer to how well this may be done in practice to Section 4.4.

Zoom In Zoom Out Reset image size

The same qualitative conclusions hold for a $1{M}_{\odot }$ lens, indicative of a white dwarf, shown in Figure 8.¹³ Despite the smaller Einstein angle (${\theta }_{E}=1.0$ mas), multiple images are resolved, again implying that ${\theta }_{E}$ can be directly measured. Where these differ most dramatically is near peak magnification (center panel in Figures 7 and 8); the more massive lens produces correspondingly larger distortions in the image as a result of the better resolution of the Einstein ring. The microlensing event is also significantly shorter.

Zoom In Zoom Out Reset image size

The effect of the accelerated observer frame, due to the Earth's orbital motion, is evident in the asymmetric entrance and exist from the microlensing event. However, this effect may be confused with features arising from the impact of the asymmetric beam and arbitrary event orientation. Thus in Figures 7 and 8 we have also shown by the red contours the images with the Earth's orbital motion neglected. From these it is clear that asymmetry in the temporal evolution of the image structure at early and late times is a robust indicator of parallax.

The primary parameters of interest are those of the lens: mass, velocity, distance. Given only the magnification light curve these are degenerate, constrained only by the event timescale and peak magnification.¹⁴ As a consequence, all that can be recovered even in principle are the projected impact parameter in unites of ${\theta }_{E}$, t_E, and the time of maximum magnification. However, images like those in Figures 7 and 8 introduce at least two key additional observables: ${\theta }_{E}$ and a lensing parallax.

The peak magnification occurs at the projected closest approach, for which the projected impact parameter is ${\beta }_{\min }$ and the magnification is

$$\begin{eqnarray}&&{A}_{\max }=\displaystyle \frac{{\beta }_{\min }^{2}+2{\theta }_{E}^{2}}{{\beta }_{\min }\sqrt{{\beta }_{\min }^{2}+4{\theta }_{E}^{2}}}.\end{eqnarray} \tag{ 22 }$$

The degree to which ${A}_{\max }$ may be measured depends upon the quality of the photometry that can be performed, and thus is likely to be limited by systematic errors in the radio flux calibration (flux calibration accuracy for the VLBA and the VLA at 43 GHz would be roughly 10%). At this time the separation between the centroid of the primary and secondary images is

$$\begin{eqnarray}&&{\rm{\Delta }}\theta =\sqrt{{\beta }_{\min }^{2}+4{\theta }_{E}^{2}},\end{eqnarray} \tag{ 23 }$$

which may be measured directly from the images to a precision that exceeds the interferometric beam width by up to an order of magnitude for high signal-to-noise detections (Reid et al. 2003). Thus, together, these yield a high-precision estimate of ${\theta }_{E}$:

$$\begin{eqnarray}{\theta }_{E} & = & \displaystyle \frac{{\rm{\Delta }}\theta }{\sqrt{2}}{[{A}_{\max }\sqrt{{A}_{\max }^{2}-1}-({A}_{\max }^{2}-1)]}^{1/2}\\ & \approx & \displaystyle \frac{{\rm{\Delta }}\theta }{2}\left(1-\displaystyle \frac{1}{8{A}_{\max }^{2}}\right),\end{eqnarray} \tag{ 24 }$$

where the latter expression is for large ${A}_{\max }$, and a fractional precision of

$$\begin{eqnarray}&&\displaystyle \frac{{\sigma }_{{\theta }_{E}}^{2}}{{\theta }_{E}^{2}}\approx \displaystyle \frac{{\sigma }_{{\rm{\Delta }}\theta }^{2}}{{\rm{\Delta }}{\theta }^{2}}+\displaystyle \frac{{\sigma }_{{A}_{\max }}^{2}}{4{A}_{\max }^{6}}.\end{eqnarray} \tag{ 25 }$$

When the uncertainty is dominated by the astrometric uncertainty, this implies that for a $10{M}_{\odot }$ black hole ${\theta }_{E}$ can be measured to within $\approx 1 \%$.

Even without detailed event modeling, ${\theta }_{E}$ immediately enables the determination of the relative proper motion from Equation (2):

$$\begin{eqnarray}&&{{\boldsymbol{\mu }}}_{L}-{{\boldsymbol{\mu }}}_{S}=\displaystyle \frac{{\theta }_{E}}{{t}_{E}}\hat{u},\end{eqnarray} \tag{ 26 }$$

where $\hat{u}$ is the asymptotic direction of the vector between the first and second images. The source proper motions can be directly measured by follow-up observations; the angular velocities of individual maser spots have been measured with an accuracy better than $1\,\mathrm{mas}\,{\mathrm{yr}}^{-1}$. As with the source positions, this is limited by the electron scattering (Reid et al. 2003), and thus for sources beyond the central pc it should be possible to measure this nearly an order of magnitude better, giving a typical accuracy of $0.1\,\mathrm{mas}\,{\mathrm{yr}}^{-1}$, corresponding to a physical velocity of $4\,\mathrm{km}\,{{\rm{s}}}^{-1}$ at the distance of the Galactic center. With a typical lens apparent velocity of $30\,\mathrm{km}\,{{\rm{s}}}^{-1}$, comparable to the velocity dispersion within the disk, this induces to a roughly 10% uncertainty in the inferred ${{\boldsymbol{\omega }}}_{L}$.

The remaining degeneracy with D_L may be broken via the measurement of a parallax from the microlensing event itself. The sources of most interest also are expected to have the longest t_E, and thus encompass a substantial fraction of a year. As made explicit in Figures 7 and 8, long t_E permit measurements of the impact of the Earth's orbital motion, i.e., a parallax. This comes in two forms.

First, the light curve itself is asymmetric as a result of the modified evolution of the Earth–lens–source alignment (Smith et al. 2003). Figure 9 shows the effect of parallax on light curve of a $10\,{M}_{\odot }$ black hole lensing event. Since the Earth accelerates substantially during this period, the additional component is distinguishable from the otherwise unknown but essentially fixed velocities of the source and lens, resulting in an asymmetry in the magnification light curve. While the maser emission and optical stellar emission are not spatially coincident, the offset in the distance is much smaller than D_S, and thus D_L can be reconstructed from either the optical or infrared light curves. Typical uncertainties of $\% 10$ in radio flux calibration may preclude the use of the radio light curve for this purpose. Since Miras are among the brightest stars in the microlensing survey fields, they may be good candidates for obtaining D_L in this way, assuming the underlying variability can be adequately modeled.

Zoom In Zoom Out Reset image size

Second, as described in the previous section, the underlying impact on the lensed images themselves is directly measurable. Note that this is essentially the same effect; the asymmetry in the magnification light curve arises from the asymmetry in the evolving image structure. However, unlike the parallax effect in the light curve, the image structure is insensitive to the Mira variability. While it is possible to attempt a full fit to the sequence of images, most of the information is contained in the locations of the multiple images. Thus, shown in Figure 10 is the angular separation of image components as a function of time for the sequence of images in Figure 7.

Zoom In Zoom Out Reset image size

The mock data shown in Figure 10 provide a convenient way to estimate the accuracy with which the parallax parameter can be reconstructed. We compute the centroid positions of the multiple image components by performing a maximum-likelihood fit of the positions of beam-convolved point sources. The uncertainty in the centroid is constructed assuming a flux limit of 10% of the maximum brightness of an unlensed maser spot, and introduce corresponding Gaussian fluctuations in the mock data (see Figure 10). At late times the uncertainties grow as a result of the dimming of one of the two images far from peak magnification.

As a simple model we write the impact parameter as a linear combination of impact parameters with and without parallax effect included:

$$\begin{eqnarray}&&\beta =\lambda {\beta }_{p}+(1-\lambda ){\beta }_{{np}},\end{eqnarray} \tag{ 27 }$$

where ${\beta }_{{np}}$ and ${\beta }_{p}$ are the angular separations when the orbital acceleration of the Earth is ignored and assuming ${D}_{L}={D}_{S}/2$. As shown in Appendix A.3, the interpolation parameter is related to the lens distance by $\lambda =({D}_{S}-{D}_{L})/{D}_{L}$, and hence a measurement of λ corresponds to a measurement of the lens position. That is, keeping all angular measurements fixed, fitting the evolving image separation directly probes the lens distance.

A maximum-likelihood fit recovers $\lambda =0.84\pm 0.20$, where the 1σ errors are indicated. That is, for typical parameters the impact of parallax can be detected within the radio images at more that 4σ. For a given λ the corresponding estimates of the lens distance and its uncertainty are

$$\begin{eqnarray}&&{D}_{L}=\displaystyle \frac{{D}_{S}}{1+\lambda }\quad \mathrm{and}\quad \displaystyle \frac{{\sigma }_{{D}_{L}}^{2}}{{D}_{L}^{2}}=\displaystyle \frac{{\sigma }_{\lambda }^{2}}{{(1+\lambda )}^{2}},\end{eqnarray} \tag{ 28 }$$

and thus we recover ${D}_{L}=4.3\pm 0.5\,\mathrm{kpc}$. That is, typically, D_L can be reconstructed with an accuracy of roughly 10%.

With an estimate of the lens distance it is possible to reconstruct the lens mass and transverse velocity from ${\theta }_{E}$ and ${{\boldsymbol{\omega }}}_{L}$. The mass may be estimated given measurements of ${\theta }_{E}$ and λ via

$$\begin{eqnarray}&&M=\displaystyle \frac{{c}^{2}}{4G}{\theta }_{E}^{2}\displaystyle \frac{{D}_{S}}{\lambda }.\end{eqnarray} \tag{ 29 }$$

Assuming the uncertainty in ${\theta }_{E}$ is negligible and in D_S is of order 10%, the fractional uncertainty in M is given by that in λ, and thus the mass can be estimated with an accuracy of roughly 14%. The lens transverse velocity is given by

$$\begin{eqnarray}&&{{\boldsymbol{v}}}_{L}={D}_{L}{{\boldsymbol{\omega }}}_{L}.\end{eqnarray} \tag{ 30 }$$

Assuming the typical uncertainties of 10% in both quantities, the transverse velocity can be typically reconstructed to 14% as well.

Resolving the lensed images provides a natural way to distinguish between single and double lenses, i.e., isolated and binary black holes. The latter, black hole–black hole binary systems are of particular interest as a potential gravitational wave source (Abadie et al. 2012). At the same time they have the potential to both inform and leverage the observations of X-ray binaries (see, e.g., Remillard & McClintock 2006). Figure 11 shows a typical binary image, which may be immediately distinguished from Figure 7 by both the complicated lensed-image morphology and, more directly, the presence of a third image.

Zoom In Zoom Out Reset image size

The efficiency with which binaries can be detected depends strongly on their projected angular separation. Sufficiently compact binaries will appear as a single lens, while the components of sufficiently wide binaries will produce independent lensing events. After conducting a suite of numerical experiments with binaries in different orientations, peculiar velocities, and impact parameters, we have found that binaries with projected angular separations spanning an order of magnitude about the Einstein angle can be readily identified. The case shown in Figure 11 has a binary angular separation comparable to the individual Einstein radii of the black holes, 3.2 mas.

As with the isolated case, the motion of the Earth induces shifts in the image locations due to parallax, and thus the lens distance can in principle be extracted. Unlike the isolated case, typically the image morphology during maximum magnification is complicated. This is mitigated by the fact that the parallax signal is largest when the lens–source separation is largest, and thus where the binary will appear most point-like. Nevertheless, we leave binary parameter estimation for future work.

The large event rates enable the systematic construction of a black hole mass function over a decade. For radio microlensing surveys (see below) this rate can be an order of magnitude larger. Unlike X-ray binaries, for microlensing the observing biases are well understood and can be addressed via direct modeling, resulting in an accurate representation of the mass distribution of stellar mass black holes in the Milky Way.

The mass distribution of black holes would be immediately diagnostic of the dynamics of core-collapse supernovae, the cataclysmic events surrounding their formation. Typically, accretion of fallback material produces a mass gap between neutron stars and black holes, a region of relative paucity in the black hole mass function extending to masses well above $2{M}_{\odot }$ (Timmes et al. 1996; Woosley et al. 2002; Zhang et al. 2008). In all cases, key uncertainties in the stellar evolution (e.g., wind loss rates, supernovae energetics, etc.) impact the final mass distributions (Gallart et al. 2005). There have been studies on the effect of variations in neutrino mechanism for core-collapse supernovae on the theoretical black hole mass function (Pejcha & Thompson 2015) and using the observed black hole mass function there already have been constraints on core-collapse supernova (Kochanek 2015). For compact binary systems, the location and depth of the mass gap is further sensitive to the binary evolution history (Fryer et al. 2002). While tentative evidence for a mass gap has already been reported based on the inferred masses of black holes in X-ray binaries (Özel et al. 2010), this is necessarily subject to the variety of strong selection effects and substantial uncertainties mentioned above. Moreover, studies of X-ray binaries are fundamentally limited by the small number of binaries currently known.

In contrast, radio-imaged microlensing events provide a means to systematically accrue large, debiased samples of Galactic black holes, limited only by survey duration. Thus, it offers both a method to assess the poorly known biases inherent in X-ray binaries and ultimately to directly access key elements of massive star evolution.

Furthermore, radio-imaged microlensing events will produce estimates of the lens positions and tangential velocities. Thus, together with the mass function it is possible in principle to construct a Galactic black hole distribution function. The velocity distribution will necessarily be strongly impacted by the dynamics of core-collapse supernovae via supernova kicks (Willems et al. 2005), as well as the subsequent evolution of the Galactic remnant population.¹⁵

Galactic black holes also serve as a fossil record of massive star formation. Therefore, detailed comparisons of the local¹⁶ densities of stellar mass black holes, white dwarfs, and sub-solar stars yield an estimate of the potentially temporally and spatially varying stellar initial mass function. In particular, this provides a unique ability to constrain the history of the high-mass end of the initial mass function.

The measurement of the mass of a single isolated stellar-mass black hole would provide an important calibration for population synthesis computations, performed primarily to produce rate estimates for gravitational wave detectors (Belczynski et al. 2002). This is the consequence of the typical nature of lenses in contrast to the particular evolutionary history of X-ray binaries, currently used to calibrate the variety of uncertainties in the rate estimates.

More generally, the binarity rate and source distribution are key inputs into black hole binary gravitational wave source population estimates. High-mass X-ray binaries, the objects likely progenitors of the black hole–black hole binaries visible by experiments like LIGO (Abadie et al. 2010), are necessarily young by virtue of the large mass of the secondary, providing a strong bias toward populations associated with the recent star formation history of the Milky Way. In stark contrast, binary mergers are expected to occur very long times after their formation, and are therefore indicative of the integrated star formation history of the Milky Way. The detection of a black hole–black hole binary presents a means to study the binary black hole population over the entirety of the Galactic history. Unfortunately, the lower limit on the projected binary separation for which radio-imaged microlensing events can efficiently detect binarity (roughly 1 mas at 8 kpc) corresponds to orbital periods of years, and thus will not present candidate gravitational wave sources. Nevertheless, it would provide a significant calibration of binary formation rates, and therefore black hole–black hole binary populations.

Thus far we have focused primarily on black hole lenses. Nevertheless, it is possible to directly detect solar mass objects as well. While we have included a discussion of the white dwarf lensing rates we have largely ignored lensing events from isolated neutron stars. An estimate of the neutron star lensing rate is nevertheless possible using the black hole lensing rate. The two key assumptions are that the massive stellar progenitors of black holes and neutron stars follow a Salpeter initial mass function and that the neutron stars receive kicks during formation.

From the first we estimate the ratio of the total number of Galactic neutron stars to the total number of Galactic black holes:

$$\begin{eqnarray}&&\displaystyle \frac{{N}_{\mathrm{NS}}}{{N}_{\mathrm{BH}}}\approx {\left(\displaystyle \frac{{M}_{\mathrm{ZNS}}}{{M}_{\mathrm{ZBH}}}\right)}^{-1.35}\approx 4,\end{eqnarray} \tag{ 31 }$$

where ${M}_{\mathrm{ZNS}}\approx 8\,{M}_{\odot }$ and ${M}_{\mathrm{ZBH}}\approx 21\,{M}_{\odot }$ are the minimum zero-age main-sequence masses of stars that form neutron stars and black holes, respectively (Woosley et al. 2002).

The second assumption modifies the volume occupied by neutron stars. Stars formed at an initial radius ${r}_{\mathrm{init}}$ with a natal kick (${v}_{\mathrm{kick}}$) larger than the circular velocity (${v}_{\mathrm{circ}}$) will isotropize within a radius of

$$\begin{eqnarray}&&{r}_{\max }\approx {r}_{\mathrm{init}}{e}^{{v}_{\mathrm{kick}}^{2}/2{v}_{\mathrm{circ}}^{2}}\approx 7{r}_{\mathrm{init}},\end{eqnarray} \tag{ 32 }$$

where we have used typical values ${v}_{\mathrm{kick}}\approx 400\,\mathrm{km}\,{{\rm{s}}}^{-1}$ and ${v}_{\mathrm{circ}}\approx 200\,\mathrm{km}\,{{\rm{s}}}^{-1}$. Thus, neutron stars formed in the Galactic bulge will be distributed across a volume nearly 400 times larger than similarly formed black holes, and thus exhibit a number density 400 times smaller. The net result is that for typical numbers the rate of neutron star microlensing events is expected to be roughly 1% of that for the black holes, justifying our neglect of neutron stars.

However, we caution that this conclusion is extremely sensitive to the typical kick velocities. If the typical kick velocity is $300\,\mathrm{km}\,{{\rm{s}}}^{-1}$ the neutron star lensing rate rises to 10% of the black hole lensing rate. More importantly, if the neutron star kick velocity distribution contains a low-velocity tail or is bimodal, as suggested by the large number of known Galactic neutron stars (see, e.g., Fryer et al. 1998; Podsiadlowski et al. 2005), the bulge may retain a large fraction of neutron stars originally formed within it. In this cases the neutron star and black hole lensing rates can be comparable. Thus, the relative frequency of neutron star and black hole events provides an additional probe of the distribution of neutron star formation kicks.

If large numbers of neutron star lensing events are observed, radio-imaged microlensing provides a novel way in which to directly measure the neutron star mass function, independent of the biases inherent in the study of neutron star binaries. Unfortunately, given their large typical age ($\gtrsim 1\,\mathrm{Gyr}$) these are all likely to be exceedingly dim, with luminosities of ≈10⁻¹¹–${10}^{-9}{L}_{\odot }$, and therefore not amenable to direct size measurements, complicating any effort to directly constrain the high-density nuclear equation of state.

Nominally, these objects would also fall below the death line for typical pulsar magnetic field strengths (${10}^{12}\,{\rm{G}}$) and periods (1 s). However, due to magnetic field decay, implicated by the low surface fields in recycled millisecond pulsars, these objects could still exhibit observable magnetospheric emission; typical surface fields after a Gyr would then be expected to lie near ${10}^{11}\,{\rm{G}}$. Nevertheless, many would still live near the death line, and therefore would at best be transient radio sources. Hence, in this case the radio emission of isolated neutron stars found by radio-imaged microlensing would provide an unbiased probe of neutron star magnetization and its evolution.

The single lens equation is given by a geometrical relation between the angular position of source β, angular position of the image θ and Einstein angle ${\theta }_{E}$ as follows:

$$\begin{eqnarray}&&{\theta }^{2}-\theta \beta -{\theta }_{E}^{2}=0,\end{eqnarray} \tag{ 33 }$$

where β and θ are given with respect to the center of the coordinate system where the lens is located. In gravitational microlensing β changes with time on a straight line as

$$\begin{eqnarray}&&{\beta }^{2}={\beta }_{0}^{2}+{\theta }_{E}^{2}{\left(\displaystyle \frac{t-{t}_{0}}{{t}_{E}}\right)}^{2},\end{eqnarray} \tag{ 34 }$$

where ${\beta }_{0}$ is the minimum impact factor and t₀ is the time of closest angular approach of the lens and source. From Equation (33), we get two solutions for the position of two images at two sides of the lens as follows:

$$\begin{eqnarray}&&{\theta }^{\pm }=\displaystyle \frac{\beta \pm \sqrt{{\beta }^{2}+4{\theta }_{E}^{2}}}{2}.\end{eqnarray} \tag{ 35 }$$

Substituting Equation (34) in (35) results in the radial dynamics of the images as a function of time.

The angular distance between the two images (${\rm{\Delta }}\theta ={\theta }^{+}-{\theta }^{-}$), regardless of the position of the source, can be written in terms of the impact parameter and the Einstein angle:

$$\begin{eqnarray}&&{\rm{\Delta }}\theta =\sqrt{{\beta }^{2}+4{\theta }_{E}^{2}}\end{eqnarray} \tag{ 36 }$$

At the minimum impact parameter where $\beta ={\beta }_{0}$, the two images attain the nearest approach to each other and in this case the trajectory of the source is perpendicular to the line connecting these two images:

$$\begin{eqnarray}&&{\beta }_{0}=\sqrt{{\rm{\Delta }}{\theta }^{2}-4{\theta }_{E}^{2}}.\end{eqnarray} \tag{ 37 }$$

We note that on the right-hand side the term in the square-root should be positive which implies ${\rm{\Delta }}\theta \geqslant 2{\theta }_{E}$. So the separation between the two images is always larger than the diameter of the Einstein ring while one image is inside and the other one is outside the ring.

The other observable parameter which can be measured from the photometry is the magnification of the source during the lensing. Since in the gravitational lensing the flux of light is conserved, the apparent change of area for lensed images would lead to a change in the energy we receive from each image. The area of images compared to the source size is given by the inverse of determinant of a Jacobian. Denoting the areas for the two images by ${A}^{+}$ and ${A}^{-}$ the differential element of the area of an image to that of the source is given by:

$$\begin{eqnarray}&&{A}^{\pm }=\displaystyle \frac{{\theta }^{\pm }}{\beta }\displaystyle \frac{\partial {\theta }^{\pm }}{\partial \beta }.\end{eqnarray} \tag{ 38 }$$

Substituting Equation (35) in (38) the area for each image can be written as:

$$\begin{eqnarray}&&{A}^{\pm }=\displaystyle \frac{1}{2}\left(\displaystyle \frac{{\beta }^{2}+2{\theta }_{E}^{2}}{\beta \sqrt{{\beta }^{2}+4{\theta }_{E}^{2}}}\pm 1\right).\end{eqnarray} \tag{ 39 }$$

The overall magnification is the sum of these two terms,

$$\begin{eqnarray}&&A=| {A}^{+}| +| {A}^{-}| =\displaystyle \frac{{\beta }^{2}+2{\theta }_{E}}{\beta \sqrt{{\beta }^{2}+4{\theta }_{E}^{2}}}.\end{eqnarray} \tag{ 40 }$$

Substituting equation for (34) the magnification is independent of Einstein angle and only depends on t₀, t_E and ${u}_{0}={\beta }_{0}/{\theta }_{E}$.

At the time of minimum approach the magnification would have its maximum as a function of ${\beta }_{0}$ and ${\theta }_{E}$. Substituting ${\beta }_{0}$ from Equation (37) in (40), maximum magnification can be obtained as a function of Einstein angle and angular separation between the images as follows:

$$\begin{eqnarray}&&{A}_{\max }=\displaystyle \frac{{\rm{\Delta }}{\theta }^{2}-2{\theta }_{E}^{2}}{{\rm{\Delta }}\theta \sqrt{{\rm{\Delta }}{\theta }^{2}-4{\theta }_{E}^{2}}}.\end{eqnarray} \tag{ 41 }$$

As a result the Einstein angle can be directly extracted from ${A}_{\max }$ as described in Equation (24).

If the distance of the source from Earth is known (which is the case for microlensing events in the Galactic center) measuring ${\theta }_{E}$ can put constraints on the mass and the lens distance from the observer. On the other hand, parallax measurements can impose another constraint on the lens parameters and in theory can resolve them.

The microlensing parallax effect results from the accelerating motion of the Earth around the Sun. It manifests itself as the relative motion of the source with respect to lens deviating from a straight line. Assuming the Earth moves around the Sun in a circular orbit, the relative angular motion of the lens and source with respect to the observer in the case of a stationary source and lens is given by

$$\begin{eqnarray}&&{\pi }_{x}=\pi \cos \xi (t),\end{eqnarray} \tag{ 42 }$$

$$\begin{eqnarray}&&{\pi }_{y}=\pi \cos \beta \sin \xi (t),\end{eqnarray} \tag{ 43 }$$

where

$$\begin{eqnarray}&&\pi =\displaystyle \frac{1}{{\theta }_{E}}\left(\displaystyle \frac{1\,\mathrm{au}}{{D}_{l}}-\displaystyle \frac{1\mathrm{au}}{{D}_{s}}\right),\end{eqnarray} \tag{ 44 }$$

β is the angular deviation of the orbital plane with respect to our line of sight, $\xi (t)=\omega t+{\xi }_{0}$, ω is the orbital angular frequency and ${\xi }_{0}$ is the initial phase of the Earth orbit which can be set to zero at the equinox. On the other hand, the relative displacement of the source with respect to position of the Sun on the lens plane is given by

$$\begin{eqnarray}&&{\beta }_{x}^{S}=\beta (t)\cos \alpha -{\beta }_{0}(\cos \alpha +\sin \alpha ),\end{eqnarray} \tag{ 45 }$$

$$\begin{eqnarray}&&{\beta }_{y}^{S}=\beta (t)\sin \alpha +{\beta }_{0}(\cos \alpha -\sin \alpha ),\end{eqnarray} \tag{ 46 }$$

where the superscript represents relative motion of the source with respect to the Sun,  ${\beta }_{0}$ is the closest approach of the source to the Sun, and α is the angle between the trajectory of source and semimajor axis of the Earth orbit. Now combining these two motions the relative motions of the source with respect to the Earth is given by

$$\begin{eqnarray}&&{\boldsymbol{\beta }}(t)={{\boldsymbol{\beta }}}^{S}(t)-{\boldsymbol{\pi }}(t).\end{eqnarray} \tag{ 47 }$$

The second term corresponds to a perturbation compared to simple microlensing events. The photometric effect of parallax can be measured as an asymmetry of the microlensing light curve (Alcock et al. 1995). Also the astrometric effect of parallax is observable if the lensed images are resolved. By combining Equations (24) and (44) the distance of the lens from the observer (D_l) can be extracted. Substituting in the definition of ${\theta }_{E}$, the mass of the lens can be obtained. In fact the lens mass can be obtained only given ${\theta }_{E}$ and ${\pi }_{r}$ using the relation:

$$\begin{eqnarray}&&\displaystyle \frac{M}{{M}_{\odot }}=\displaystyle \frac{1}{8{\pi }_{r}}(\displaystyle \frac{{\theta }_{E}}{1\,\mathrm{mas}}).\end{eqnarray} \tag{ 48 }$$

From photometric observations and microlensing light curves the physical parameter that is obtained by fitting to the light curve is the Einstein crossing time. Having measured this parameter the transverse velocity of the lens can be derived as ${v}_{t}={D}_{l}{\theta }_{E}/{t}_{E}$.

The acceleration of the Earth impacts the image separations for the same reason it affects the light curve. Here we estimate the size of this effect, explicitly determining the way in which the lens distance enters.

We begin by defining averaged and perturbed lens and source positions,

$$\begin{eqnarray}&&{\boldsymbol{s}}\equiv {{\boldsymbol{x}}}_{S}-{{\boldsymbol{x}}}_{\oplus }=\bar{{\boldsymbol{s}}}-\delta \quad \mathrm{and}\quad {\ell }\equiv {{\boldsymbol{x}}}_{L}-{{\boldsymbol{x}}}_{\oplus }=\bar{{\ell }}-\delta ,\end{eqnarray} \tag{ 49 }$$

where $\bar{{\boldsymbol{s}}}$ and $\bar{{\ell }}$ include the linear lens and source motions and the linear component of the Earth's motion. The acceleration of the Earth is contained in $\delta$ which is by construction quadratic in time to lowest order. We will presume this to be small in comparison to $\bar{{\boldsymbol{s}}}$ and $\bar{{\ell }}$.

Inserting these into Equation (18), linearizing in $\delta$, and ignoring terms of order ${\beta }^{2}$ yields

$$\begin{eqnarray}\beta & \equiv & \left(1-\displaystyle \frac{{\ell }{\ell }}{{D}_{L}^{2}}\right)\cdot \displaystyle \frac{{\boldsymbol{s}}}{{D}_{S}}\\ & \approx & \left(1-\displaystyle \frac{\bar{{\ell }}\bar{{\ell }}}{{D}_{L}^{2}}\right)\cdot \displaystyle \frac{\bar{{\boldsymbol{s}}}}{{D}_{S}}\\ & & +\left(\displaystyle \frac{{D}_{S}}{{D}_{L}}-1\right)\left(1-\displaystyle \frac{\bar{{\ell }}\bar{{\ell }}}{{D}_{L}^{2}}\right)\cdot \displaystyle \frac{\delta }{{D}_{S}}\\ & & +\displaystyle \frac{\bar{{\ell }}}{{D}_{L}}\displaystyle \frac{\bar{{\boldsymbol{s}}}}{{D}_{L}}\cdot \left(1-\displaystyle \frac{\bar{{\ell }}\bar{{\ell }}}{{D}_{L}^{2}}\right)\cdot \displaystyle \frac{\delta }{{D}_{S}}.\end{eqnarray} \tag{ 50 }$$

Therefore, noting that $\beta \cdot {\ell }=0$,

$$\begin{eqnarray}&&\beta \approx {\beta }_{{np}}+\left(\displaystyle \frac{{D}_{S}-{D}_{L}}{{D}_{L}}\right)\displaystyle \frac{{\beta }_{{np}}}{{\beta }_{{np}}}\cdot \displaystyle \frac{\delta }{{D}_{S}},\end{eqnarray} \tag{ 51 }$$

where ${\beta }_{{np}}\equiv (1-\bar{{\ell }}\bar{{\ell }}/{D}_{L}^{2})\cdot \bar{{\boldsymbol{s}}}/{D}_{S}$ is the angular position of the lens ignoring the Earth's acceleration. The corresponding angular separation between the multiple images is then given by Equation (36)

$$\begin{eqnarray}&&{\rm{\Delta }}\theta \approx {\rm{\Delta }}{\theta }_{{np}}+\left(\displaystyle \frac{{D}_{S}-{D}_{L}}{{D}_{L}}\right)\displaystyle \frac{{\beta }_{{np}}}{{\rm{\Delta }}{\theta }_{{np}}}\cdot \displaystyle \frac{\delta }{{D}_{S}}.\end{eqnarray} \tag{ 52 }$$

Note that Equation (51) interpolates between ${\beta }_{{np}}$ and β as the lens distance changes. With the definition

$$\begin{eqnarray}&&{\beta }_{p}\equiv {\beta }_{{np}}+\displaystyle \frac{{\beta }_{{np}}}{{\beta }_{{np}}}\cdot \displaystyle \frac{\delta }{{D}_{S}},\end{eqnarray} \tag{ 53 }$$

i.e., when ${D}_{L}={D}_{S}/2$, we recover Equation (27),

$$\begin{eqnarray}&&\beta =\lambda {\beta }_{p}+(1-\lambda ){\beta }_{{np}},\end{eqnarray} \tag{ 54 }$$

where

$$\begin{eqnarray}&&\lambda \equiv \displaystyle \frac{{D}_{S}-{D}_{L}}{{D}_{L}}.\end{eqnarray} \tag{ 55 }$$