PopSyCLE: A New Population Synthesis Code for Compact Object Microlensing Events

Galaxia is a synthetic stellar survey of the Milky Way (Sharma et al. 2011). Given a survey area and location, Galaxia will generate a stellar model of the Milky Way for a circular field of sky, which corresponds to a 3D conical volume. Galaxia is a resolved simulation, that is, for every single star in the field, it will return a position, velocity, age, mass, photometry in several filters, 3D extinction, metallicity, surface gravity, and more.⁵ Galaxia implements the Besançon analytic model for the Milky Way (Robin et al. 2003), with a modified version of the disk kinematics that adjusts the velocity in the azimuthal direction (Shu 1969). For a detailed description of the implementation and Galactic model parameters, see Sharma et al. (2011) and references therein; for a summary of the model, see Tables 1–3 in Sharma et al. (2011). We also perform additional modifications to Galaxia's bulge kinematics, as they are known to produce some inconsistencies with observations. This is described in more detail in Appendix A. In the rest of the main text, whenever reference is made to Galaxia, we specifically mean Galaxia with the modified bulge, unless otherwise mentioned. Additionally, there is also an option to implement the N-body stellar halo model of Bullock & Johnston (2005) instead of the smooth analytic stellar halo model from the Besançon model. We choose to use the smooth analytic model instead of the N-body model, as we are primarily considering lensing in the disk and bulge; stellar halo substructure is not a significant concern.

Galaxia includes a 3D extinction map; however, the built-in extinction law from Schlegel et al. (1998) is not well suited to the high-extinction regions toward the inner Galactic bulge. Unfortunately, the newest 3D extinction maps (e.g., Green et al. 2018) do not cover significant fractions of the Galactic bulge region. Similarly, the extinction maps of Nataf et al. (2013), which were constructed using Optical Gravitational Lensing Experiment (OGLE), VISTA Variables in the Via Lactea (VVV), and Two Micron All Sky Survey (2MASS) data, only cover the OGLE-III bulge fields. Instead, we modify the Galaxia output in order to adopt the reddening law of Damineli et al. (2016) as described in Appendix B.

PyPopStar is a software package that generates single-age, single-metallicity populations (i.e., star clusters) using adjustable initial mass functions (IMFs), multiplicity distributions, stellar evolution and atmosphere grids, and extinction laws; for a more detailed description, see Section 4 of Hosek et al. 2019.⁶ In summary, IMFs are parameterized as piecewise functions consisting of multiple power laws and generated according to the prescription from Pflamm-Altenburg & Kroupa (2006). We adopt the MESA Isochrones and Stellar Tracks (MIST) v1.2 stellar evolution models (Choi et al. 2016; Dotter 2016), which are calculated using the Modules for Experiments in Stellar Astrophysics (MESA) code (Paxton et al. 2011, 2013, 2015). These models provide theoretical isochrones with stellar temperatures and surface gravities as a function of mass at a given stellar age. Model atmosphere grids are then used to assign an atmosphere to each star. The stellar temperature determines which grid is used; an ATLAS9 grid is used for stars with T_eff > 5500 K (Castelli & Kurucz 2004), a PHOENIX grid is used for 3800 K < T_eff < 5000 K (version 16; Husser et al. 2013), and a BTSettl grid is used for 1200 K < T_eff < 3200 K (Baraffe et al. 2015). For temperatures at a transition region between grids (e.g., 5000–5500 K), an average atmosphere between the two model grids is used. For all models, solar metallicity is assumed. This assumption will not have a large effect on the BH population synthesis; see Section 8.4.2 for justification. Although the MIST stellar evolution models do not evolve stars all the way down to NSs and BHs, it does have some partial support for WDs; a full discussion of stellar evolution with WDs is presented in Section 2.2.2.

The  PyPopStar software as described in Section 4 of Hosek et al. (2019) would discard stars that had evolved into compact objects. To be able to perform population synthesis, we have updated  PyPopStar with an IFMR that will properly account for compact objects. Given the zero-age main-sequence (ZAMS) mass of a star, the IFMR describes the type and mass of compact object formed. The IFMR is an active area of research, as the ZAMS mass does not solely determine the evolutionary path of a star. Factors such as rotation (which induces mixing in the star) and mass loss due to stellar winds (which are a function of metallicity) are important in the final stages of stellar evolution. In particular, the pre-supernova core structure of the star strongly determines whether there is an explosion or not, which in turn affects the type of compact remnant formed (Sukhbold et al. (2018) and references therein). We modify  PyPopStar to implement a combination of IFMRs taken from recent simulations to build a population of WDs, BHs, and NSs as described in more detail below.

Note that  PyPopStar has support for stellar multiplicity. However, as Galaxia only has a treatment of single stars, this functionality was not used. A discussion of this is presented in Section 8.4.1.

2.2.1. Neutron Stars and Black Holes

2.2.2. White Dwarfs

One complicating factor in WD population synthesis is that the MIST models do not include the late stages of WD cooling. Thus, not all WDs are generated, in particular, the oldest and faintest ones. This means that a WD IFMR is also required. For all WDs not included in a set of evolutionary models, we use the empirically determined IFMR of Kalirai et al. (2008, hereafter K08), given by

$\begin{eqnarray}&&{M}_{\mathrm{WD}}=(0.109{M}_{\mathrm{ZAMS}}+0.394){M}_{\odot }.\end{eqnarray} \tag{ 1 }$

The data in K08 have initial mass in the range of 1.16 M_⊙ < M_ZAMS < 6.5 M_⊙. However, to have continuous coverage in M_ZAMS for the IFMR, we extend the domain to 0.5 M_⊙ < M_ZAMS < 9 M_⊙, where the lower mass limit was chosen to match the MIST models' lower mass limit, and the upper mass limit was chosen to match where the IFMR for NSs and BHs begins. The low-mass extrapolation covers all possible WDs formed within the age of the universe. The high-mass extrapolation ensures that there is an IFMR for all stellar masses, and the maximum WD mass would then be

$\begin{eqnarray*}&&{M}_{\mathrm{WD}}({M}_{\mathrm{ZAMS}}=9\,{M}_{\odot })=1.375\,{M}_{\odot }\lesssim {M}_{\mathrm{ch}},\end{eqnarray*}$

where ${M}_{\mathrm{ch}}=1.4{M}_{\odot }$ is the Chandrasehkar mass.

In particular for the MIST models, WDs are tracked only so far down the cooling curve, and the coolest/oldest WDs are usually dropped. For these objects, the WD IFMR (Equation (1)) is used to produce dark WDs, and their luminosity is neglected. This is justified, as these old WDs have absolute magnitudes of M ∼ 18 in red-optical filters and contribute negligible flux at typical distances for lenses and sources (e.g., Campos et al. 2016). Thus, a PopSyCLE population will contain both luminous WDs from the MIST models and dark WDs from the K08 IFMR.

The scale of the microlensing event is set by the the angular Einstein radius

$\begin{eqnarray}&&{\theta }_{{\rm{E}}}=\sqrt{\displaystyle \frac{4{GM}}{{c}^{2}}\left(\displaystyle \frac{1}{{d}_{L}}-\displaystyle \frac{1}{{d}_{S}}\right)},\end{eqnarray} \tag{ 3 }$

where d_S is the observer–source distance, d_L is the observer–lens distance, and M is the mass of the lens. Note that in microlensing θ_E is unresolved, unlike in strong lensing. For a "typical" microlensing event consisting of a stellar bulge source at d_S = 8 kpc and a stellar disk lens at d_L = 4 kpc with a mass M = 0.5 M_⊙, the Einstein radius is θ_E = 0.71 mas. For context, the highest-resolution images from the largest ground-based or space-based telescopes (with a filled aperture) are 50 mas in the optical or infrared. However, in some cases Einstein radii of this scale are interferometrically resolvable (Dong et al. 2019), moving from the regime of microlensing to strong lensing.

The Einstein crossing time t_E, the time it takes for the lens to traverse the Einstein radius, can be inferred from the photometric light curve. The magnitude of the relative source–lens proper motion μ_rel relates the Einstein crossing time and Einstein radius via

$\begin{eqnarray}&&{\theta }_{{\rm{E}}}={\mu }_{\mathrm{rel}}{t}_{{\rm{E}}}.\end{eqnarray} \tag{ 4 }$

The dimensionless source–lens separation vector normalized in units of the Einstein radius is

$\begin{eqnarray}&&{\boldsymbol{u}}=\displaystyle \frac{{{\boldsymbol{\theta }}}_{S}-{{\boldsymbol{\theta }}}_{L}}{{\theta }_{{\rm{E}}}},\end{eqnarray} \tag{ 5 }$

where ${{\boldsymbol{\theta }}}_{L}$ and ${{\boldsymbol{\theta }}}_{S}$ denote the angular position of the lens and source on the sky, respectively. The magnitude of the separation $| {\boldsymbol{u}}| =u$ as a function of time, neglecting higher-order effects such as parallax, can be written as

$\begin{eqnarray}&&u(t)=\sqrt{{u}_{0}^{2}+{\left(\displaystyle \frac{t-{t}_{0}}{{t}_{{\rm{E}}}}\right)}^{2}},\end{eqnarray} \tag{ 6 }$

where t₀ is the time of closest approach and u₀ is the separation at t₀; u₀ is also known as the impact parameter.

The photometric amplification of the source is given by

$\begin{eqnarray}&&A(u)=\displaystyle \frac{{u}^{2}+2}{u\sqrt{{u}^{2}+4}}\end{eqnarray} \tag{ 7 }$

and is maximized at u = u₀. It can be seen from this formula that a smaller u₀ produces a larger amplification. Most microlensing surveys define a "microlensing event" as an event that has u₀ less than 1.

The source flux fraction, which is sometimes called the blend parameter, or, confusingly, the blend fraction, is the fraction of source flux over total flux in the telescope's photometric extraction aperture

$\begin{eqnarray}&&{b}_{\mathrm{SFF}}=\displaystyle \frac{{F}_{S}}{{F}_{S}+{F}_{L}+{F}_{N}},\end{eqnarray} \tag{ 8 }$

where F_S, F_L, and F_N are the fluxes due to the (unlensed) source, lens, and any neighboring stars. From this definition, a highly blended event will have a source flux fraction of b_SFF ≈ 0 and an event with no blending (i.e., only flux from the source) will have b_SFF = 1.

The photometric extraction aperture is typically proportional to the spatial resolution of the images; we adopt an aperture diameter of the FWHM, which is set by seeing for ground-based surveys like OGLE and the diffraction limit for space-based surveys such as WFIRST (these surveys are described in Section  5).   To give concrete numbers, the median seeing for the OGLE survey is roughly 13, which we take to be the FWHM. WFIRST, a planned 2.4 m space telescope observing at 1630 nm (H band), will have FWHM ∼ 017. These FWHM values for OGLE and WFIRST would correspond to aperture radii of roughly  065 and 009, respectively.

The bump magnitude, Δm, is the difference between the baseline magnitude m_base and the peak magnitude m_peak

$\begin{eqnarray}&&{\rm{\Delta }}m=-2.5{\mathrm{log}}_{10}\left(\displaystyle \frac{A({u}_{0}){F}_{S}+{F}_{L}+{F}_{N}}{{F}_{S}+{F}_{L}+{F}_{N}}\right).\end{eqnarray} \tag{ 9 }$

 Note that although the magnification of the source flux is achromatic, the bump magnitude and the source flux fraction are dependent on wavelength and the photometric extraction aperture (Di Stefano & Esin 1995). Thus, when events are selected on Δm or b_SFF, we explicitly note the filter and aperture used.

In addition to the selection parameters described in Section 4.1, there are a number of other measurable microlensing quantities. In particular, the astrometric and microlensing parallax signals are useful for identifying possible BH microlensing events.

The astrometric shift ${{\boldsymbol{\delta }}}_{c}$ of the source image centroid, assuming no blending, is given by

$\begin{eqnarray}&&{{\boldsymbol{\delta }}}_{c}=\displaystyle \frac{{\theta }_{{\rm{E}}}}{{u}^{2}+2}{\boldsymbol{u}}.\end{eqnarray} \tag{ 10 }$

In contrast to the photometric amplification, the maximum astrometric shift occurs at $u=\pm \sqrt{2}$. Thus, if ${u}_{0}\lt \sqrt{2}$, the maximum astrometric shift will occur before and after the maximum photometric amplification. However, if ${u}_{0}\geqslant \sqrt{2}$, the time at which maximum astrometric shift and photometric amplification occur will coincide. It should be noted that δ_c is the maximum possible astrometric shift for a given geometry (i.e., for given M, d_L, d_S), as it does not include blending.

Blending due to a luminous lens modifies the astrometric signal to

$\begin{eqnarray}&&{{\boldsymbol{\delta }}}_{c,{LL}}=\displaystyle \frac{{\theta }_{{\rm{E}}}}{1+g}\displaystyle \frac{1+g({u}^{2}-u\sqrt{{u}^{2}+4}+3)}{{u}^{2}+2+{gu}\sqrt{{u}^{2}+4}}{\boldsymbol{u}},\end{eqnarray} \tag{ 11 }$

where g = F_L/F_S (Dominik & Sahu 2000). Blending due to the lens and other neighboring stars dilutes the astrometric signal further to

$\begin{eqnarray}&&{{\boldsymbol{\delta }}}_{c,{LLN}}=\displaystyle \frac{\tilde{A}(u){F}_{S}{{\boldsymbol{\theta }}}_{S}+{F}_{N}{{\boldsymbol{\theta }}}_{N}}{A(u){F}_{S}+{F}_{L}+{F}_{N}}-\displaystyle \frac{{F}_{S}{{\boldsymbol{\theta }}}_{S}+{F}_{N}{{\boldsymbol{\theta }}}_{N}}{{F}_{S}+{F}_{L}+{F}_{N}},\end{eqnarray} \tag{ 12 }$

where ${{\boldsymbol{\theta }}}_{S}$ and ${{\boldsymbol{\theta }}}_{N}$ are the angular positions of the source and the centroid of all the on-sky neighbors contributing to blending, relative to the lens position; A(u) is given by Equation (7), and $\tilde{A}(u)$ is defined as

$\begin{eqnarray}&&\tilde{A}(u)=\displaystyle \frac{{u}^{2}+3}{u\sqrt{{u}^{2}+4}}.\end{eqnarray} \tag{ 13 }$

For long-duration events (t_E ≳ 3 months), the motion of Earth orbiting the Sun will modify the otherwise symmetric light curve. This signal is called the microlensing parallax π_E and is given by

$\begin{eqnarray}&&{\pi }_{{\rm{E}}}=\displaystyle \frac{{\pi }_{\mathrm{rel}}}{{\theta }_{{\rm{E}}}},\end{eqnarray} \tag{ 14 }$

where ${\pi }_{\mathrm{rel}}$ is the relative parallax

$\begin{eqnarray}&&{\pi }_{\mathrm{rel}}=1\,\mathrm{au}\left(\displaystyle \frac{1}{{d}_{L}}-\displaystyle \frac{1}{{d}_{S}}\right).\end{eqnarray} \tag{ 15 }$

Combining Equations (3) and (15) and solving for the lens mass M yields

$\begin{eqnarray}&&M=\displaystyle \frac{{\theta }_{{\rm{E}}}}{\kappa {\pi }_{{\rm{E}}}},\end{eqnarray} \tag{ 16 }$

where $\kappa \equiv \tfrac{4G}{1\,\mathrm{au}\,{c}^{2}}$. If both the photometric magnification A(u) and the astrometric shift δ_c(u, θ_E) can be measured, the Einstein radius θ_E can be deduced. Along with a measurement of the microlensing parallax π_E, the mass of the lens M can then be derived using Equation (16). Thus, by having both astrometric and photometric measurements, the degeneracies between lens mass, lens distance, and source distance can be broken.⁷ The lens mass can then be constrained without having to resort to assumptions about the lens and source distances. Many microlensing studies only measure the photometric signal and use a Galactic model to weight different source distances in order to derive the lens mass; however, this is a significant assumption. By combining both the photometric and astrometric microlensing signal, the lens mass can be determined without having to resort to modeling distances.

The method just described can be used to measure the masses of BHs. Since long-duration microlensing events are more likely to be BHs, probable BH microlensing events can be selected from ground-based photometry and then followed up astrometrically, since the astrometric signal peaks after the photometric signal when ${u}_{0}\lt \sqrt{2}$ (Lu et al. 2016).

With population synthesis complete, we can specify a duration and sampling cadence to perform a synthetic microlensing survey. Typical microlensing survey lengths are on the order of years, over which celestial motions are linear for nearly all lenses and sources. Thus, with an initial position and velocity, the position of any object at some later time in the survey can be calculated.

First, a grid of latitude and longitude bins is overlaid across the entire survey area. Next, four adjacent latitude and longitude bins of objects generated from population synthesis are read in and combined to make one large bin over which microlensing events are searched for at the sampled times within the survey window. This combination ensures that microlensing events that occur across the edges of two smaller bins are not missed. The algorithm for finding event candidates has four "cuts" as follows:

• 1.  
  Find the nearest on-sky neighbor for each object and calculate the separation Δθ. Pairs with separations of Δθ > θ_cut are rejected, where θ_cut is a user-specified separation. Analogous to the bin size selection, choosing too small a value for θ_cut will mean that some events are missed. Choosing too large a value for θ_cut is less problematic; it will merely slow down subsequent calculations by considering pairs that will not produce a detectable microlensing signal.
• 2.  
  Calculate the Einstein radius θ_E of the lenses from these initial microlensing candidate pairs. Only keep the pairs where Δθ/θ_E < u_cut, where u_cut is another user-specified value that sets the maximum impact parameter. For this further refined list of candidates, record all the stars that fall within the seeing disk radius θ_blend, which will be used to calculate blending later on. The seeing disk radius is also a user-specified parameter.
• 3.  
  From the new list of microlensing candidate pairs, calculate the time of closest approach t₀ when u = u₀. Events are kept if t₀ ∈ [t_i, t_f], where t_i and t_f are the start and end times for the survey, respectively.
• 4.  
  Remove unphysical "events" where the source is a dark object.

This yields the list of microlensing event candidates. The procedure is then repeated for all the bins. This algorithm will generate many repeats, due to the multiple overlaps in the bins and also the multiple evaluation times. Only unique events are counted, that is, a lens–source pair is only counted once; the event parameters that get recorded are the ones corresponding to when the source and lens are closest to each other. A cartoon schematic of this entire process is shown in Figure 2.

Zoom In Zoom Out Reset image size

With complete information about the 6D kinematics, masses, and photometry, any microlensing parameters of interest can be calculated for all event candidates.  It should also be noted that the choice of sampling cadence t_obs, defined as the time between consecutive observations, will change the number of event candidates. For example, t_obs = 1 day corresponds to a very dense sampling, while t_obs = 100 days corresponds to a very sparse sampling.  There will be a loss of sensitivity to events with t_E < t_obs since shorter events will fall between observations. It should also be emphasized that the sampling cadence of PopSyCLE has no relation to the observational survey cadence of real microlensing surveys. For example, as shown in Figure 3, PopSyCLE run with a sampling cadence of 100 days finds nearly all events with t_E ≳ 30 days since objects only have to approach each other within u_cut θ_E to be detected.

Zoom In Zoom Out Reset image size

Many of the microlensing events reported by PopSyCLE as described in Section 4.3 would be undetectable by a real survey. For example, there will be events too faint to be observed. Thus, we consider the list of events generated in Section 4.3 to be event candidates. The user can select from this list of event candidates those that satisfy some observational criteria to generate a final list of microlensing events.

We first use our population synthesis model to calculate the number of BHs and NSs in the Milky Way. First, Galaxia generates a random fraction of the entire Milky Way. We then perform a simplified version of the compact object population synthesis described in Section 3 by ignoring the spatial distribution and kinematics of the compact objects, instead only returning the masses of the compact objects. We then scale the number of compact objects in accordance to the fraction of stars generated to obtain the actual number.

With this method, we estimate that there are 2.2 × 10⁸ BHs and  4.4 × 10⁸ NSs in the entire Milky Way. Our findings compare well with previous theoretical estimates that predict around 10⁸–10⁹ stellar-mass BHs in the Milky Way, with a similar estimate for NSs, using methods based on the number of microlensing events toward the Galactic bulge (Agol & Kamionkowski 2002) and supernova explosion and Galactic chemical evolution models (Timmes et al. 1996).

Of the 2.2 × 10⁸ BHs produced by PopSyCLE,  8.5 × 10⁷ have masses greater than 10 M_⊙. This is comparable to the findings of Elbert et al. (2018), where they estimate that there should be around 10⁸ BHs with M > 10 M_⊙ in a galaxy with a stellar mass equal to that of the Milky Way (${M}_{* ,\mathrm{MW}}\approx 6\,\times {10}^{10}{M}_{\odot }$). Note that the models of Elbert et al. (2018) include metallicity and BH–BH mergers, which PopSyCLE currently does not include. The mass distribution of BHs is shown in Figure 5, and a discussion of the distribution is presented in Section 7.1.

We next compare stellar densities, event rates, and Einstein crossing times from PopSyCLE simulations with those from M19 and SP16. A summary is presented in Table 5.

Table 5.  Comparing PopSyCLE to Surveys

Field	ns (106 deg−2)		Γ (10−6 star−1 yr−1)		$\langle {t}_{{\rm{E}}}\rangle$ (days)
	Obs.	Mock	Obs.	Mock	Obs.	Mock
F00M19	2.62	2.52	1.3 ± 0.8	5.97	32.6	20.4
F01M19	4.54	2.92	5.5 ± 0.9	8.09	39.5	25.0
F03S	3.64	3.97	${14.0}_{-2.4}^{+2.9}$	22.74	25.5	18.5
F04S	1.11	1.12	${3.5}_{-1.7}^{+2.6}$	7.66	17.0	28.5
F06S	2.80	5.77	${36.6}_{-4.4}^{+4.9}$	50.67	17.4	18.4
F07M19	1.75	2.09	5.6 ± 1.4	4.11	51.8	39.5
F09M19	4.84	4.75	23.9 ± 2.0	43.69	18.8	17.4
F11M19	4.95	5.7	16.2 ± 1.3	32.2	21.8	17.0
F12M19	3.26	2.07	3.4 ± 1.1	7.25	36.7	16.7
F13M19	4.51	3.44	5.2 ± 1.1	11.25	30.8	21.9

Note. Fields with ^M19 indicate that the observed values come from M19, while those with ^S are from SP16.

Download table as:  ASCIITypeset image

For each field listed in Table 3, we apply the observational cuts from Tables 2 and 4 to generate a final list of detectable events as described in Section 5. In order to convert to an event rate in units of events star⁻¹ yr⁻¹, we also calculate the total number of detectable stars within the magnitude limits for the corresponding survey. The microlensing event rate, Γ, varies between different observational surveys and different fields. The rates for  Mock Sumi and Mock Mróz19 are presented in Table 6.

Table 6.  Mock Simulations

Field	Sim	ns	Γ	$\langle {t}_{{\rm{E}}}\rangle$	Med(tE)
		(×106 deg−2)	(×10−6 star−1 yr−1)	(days)	(days)
F00	S	1.46	5.87	17.3	11.2
	M19	2.52	5.97	20.4	14.4
F01	S	1.52	6.37	25.1	14.8
	M19	2.92	8.09	25.0	15.5
F02	S	2.70	19.09	20.2	14.2
	M19	6.11	25.66	19.5	13.9
F03	S	3.97	22.74	18.5	13.2
	M19	7.92	30.80	19.7	14.5
F04	S	1.12	7.66	28.5	21.6
	M19	2.37	9.96	19.6	17.3
F05	S	5.38	38.57	17.5	13.4
	M19	11.05	43.65	19.5	13.9
F06	S	5.77	50.67	18.4	12.4
	M19	12.53	59.68	18.9	12.8
F07	S	0.99	4.33	40.9	28.1
	M19	2.09	4.11	39.5	22.9
F08	S	2.07	33.68	20.1	16.0
	M19	5.02	40.04	18.1	13.1
F09	S	2.06	33.84	20.5	16.0
	M19	4.75	43.69	17.4	12.7
F10	S	0.33	19.69	13.4	13.0
	M19	0.64	33.32	24.5	15.1
F11	S	2.05	21.53	16.6	10.1
	M19	5.70	32.20	17.0	11.4
F12	S	1.11	5.80	18.1	10.1
	M19	2.07	7.25	16.7	9.1
F13	S	1.68	6.40	12.3	7.9
	M19	3.44	11.25	21.9	13.1

Note. For each field simulation (S = Mock Sumi, M19 = Mock Mróz19) we list the density of stars n_s brighter than the survey magnitude limit (no correction for blending or confusion), the event rate Γ, average Einstein crossing time $\langle {t}_{{\rm{E}}}\rangle$, and the median Einstein crossing time Med(t_E).

Download table as:  ASCIITypeset image

  The event rates from SP16 for fields F03, F04, and F06 (corresponding to MOA-II-gb14, gb21, and gb5) are 14.0, 3.5, and 36.6 events star⁻¹ yr⁻¹, respectively. For those same fields, the Mock Sumi PopSyCLE simulation gives event rates of  22.74, 7.66, and 50.67 events star⁻¹ yr⁻¹. The event rates from PopSyCLE are roughly  between 1.4 and 2.2 times higher than the ones reported in SP16.

The event rates from M19 for fields F00, F01, F07, F09, F11, F12, and F13 (corresponding to OGLE-IV-BLG547, 527, 672, 500, 611, 629, and 617) are 1.3, 5.5, 5.6, 23.9, 16.2, 3.4, and 5.2 events star⁻¹ yr⁻¹, respectively. For those same fields, the Mock Mróz19 PopSyCLE simulation gives event rates of 5.97, 8.09, 4.11, 43.69, 32.2, 7.25, and 11.25 events star⁻¹ yr⁻¹. The event rates from PopSyCLE are generally a factor of 2 times higher than the ones reported in M19, although for one field it is nearly 5 times higher and in one field it is only 0.7 of the observed rate.

Our event rate estimates from PopSyCLE use a total number of stars that is 100% complete down to the selected magnitude, as we have not imposed any confusion due to the finite beam size of a real telescope. We assume that the reported event rates in Sumi & Penny (2016) and Mróz et al. (2017, 2019) use completeness-corrected stellar number counts.

The average Einstein crossing time $\langle {t}_{{\rm{E}}}\rangle$ is a commonly reported parameter for microlensing distributions. Most authors refer to both the raw $\langle {t}_{{\rm{E}}}\rangle$ derived from the uncorrected distribution obtained directly after making all observational cuts and an efficiency-corrected $\langle {t}_{{\rm{E}}}{\rangle }_{\varepsilon }$, obtained after correcting for the detection efficiency. In this section, we compare the reported $\langle {t}_{{\rm{E}}}{\rangle }_{\varepsilon }$ from SP16 and M19 to those obtained from PopSyCLE simulations. For brevity, whenever the notation $\langle {t}_{{\rm{E}}}\rangle$ is used, it is understood to be efficiency corrected, unless otherwise stated.

The $\langle {t}_{{\rm{E}}}\rangle$ from SP16 for fields F03, F04, and F06 are 25.5, 17.0, and 17.4 days, respectively. For these same fields, the Mock Sumi PopSyCLE simulation gives $\langle {t}_{{\rm{E}}}\rangle$ of 18.5, 28.5, and 18.4 days; these values range from factors of 0.7 to 1.7 times the observed timescales.

The $\langle {t}_{{\rm{E}}}\rangle$ from M19 for fields F00, F01, F07, F09, F11, F12, and F13 are 32.6, 39.5, 51.8, 18.8, 21.8, 36.7, and 30.8 days, respectively. For these same fields, the Mock Mróz19 simulation gives $\langle {t}_{{\rm{E}}}\rangle$ of 20.4, 25.0, 39.5, 17.4, 17.0, 16.7, and 21.9 days; these values are all shorter than the observed timescales by factors of 0.5–0.9.

Values for $\langle {t}_{{\rm{E}}}\rangle$ for the events in M17 cannot be calculated, as they do not report or present individual event parameters (only the values binned into the histogram are given).   Thus, we compare the entire t_E distribution to PopSyCLE, and in particular the "peak" of the t_E distribution (${t}_{{\rm{E}},\mathrm{peak}}$), defined as the location of the maximum of the t_E distribution when the individual timescales are binned logarithmically. The peak of a logarithmically binned distribution is not the best quantity to compare, and we advocate performing a fit to the t_E distribution instead. However, for the moment we use ${t}_{{\rm{E}},\mathrm{peak}}$ as a proxy for some measure of central tendency (note that ${t}_{{\rm{E}},\mathrm{peak}}$ does not correspond to the mean, median, or mode of the t_E values).

From the efficiency-corrected t_E distributions presented in several papers, ${t}_{{\rm{E}},\mathrm{peak}}$ is located around  15–20 days (see Figure 13 of Sumi et al. 2013, Figure 17 of  Wyrzykowski et al. 2015, and Figure 2 of M17). As can be seen in Figure 4, ${t}_{{\rm{E}},\mathrm{peak}}$ from PopSyCLE falls in the lower end of that range.

Zoom In Zoom Out Reset image size

The BH PDMF encodes information about the BH IFMR and stellar IMF, and to a lesser degree the star formation history (SFH). It also provides information about BH binaries and their formation channel(s). With a sufficiently large sample of BH mass measurements from both LIGO (extragalactic) and microlensing (Galactic), the BH PDMF can be measured and the IFMR can be constrained.

We use PopSyCLE to generate the Milky Way BH PDMF, which is shown in Figure 5. The R18 IFMR clearly shows the "mass gap," as the lowest-mass BH is around 5 M_⊙, which is greater than the largest possible NS mass ${M}_{\mathrm{NS},\max }\sim 2-3\,{M}_{\odot }$ (Özel & Freire 2016; but also see Margalit & Metzger 2017, which suggests an upper limit closer to $\sim 2.2\,{M}_{\odot }$). As discussed in Section 2.2.1, there are no $\gtrsim 30\,{M}_{\odot }$ BHs, as the R18 IFMR assumes only solar-metallicity progenitor stars and PopSyCLE currently implements only single BHs. For the Milky Way, the SFH (according to the Besançon model) does not influence the BH PDMF very much, as most stars are over 10⁹ yr old. The minimum ZAMS mass for a star to form a BH is ∼15 M_⊙, and the corresponding main-sequence lifetime of such a star is ∼10⁷ yr. Thus, the vast majority of BHs produced when their progenitor star dies will have already been formed.

Zoom In Zoom Out Reset image size

The R18 IFMR also produces structures such as peaks and gaps in the BH PDMF. The spike around 6 M_⊙ is due to a combination of stars with ZAMS masses between 15–20 M_⊙ and 70–120 M_⊙. Although there are more stars of lower ZAMS masses owing to the IMF, only 34% of stars within $15\mbox{--}20\,{M}_{\odot }$ form BHs, while 60% of the >70 M_⊙ stars form BHs. The paucity around $10\,{M}_{\odot }$ is due to the fact that most stars between 25 and 27 M_⊙ form NS and not BHs. The general decrease in number of BHs greater than 8 M_⊙ is simply due to the IMF; high-mass stars are rarer. These trends are shown visually in Figure 6, which combines the IFMR with the IMF and SFH. We reemphasize that this structure is  specific to the R18 IFMR; a different IFMR will produce a different BH PDMF. However, with this assumption, the structure in the BH PDMF may be detectable with a sample of ∼100 BHs. This will become a possibility with WFIRST, as discussed in Section 7.5.

Zoom In Zoom Out Reset image size

Identifying a BH with microlensing requires showing that the lens mass exceeds 5 M_⊙ and is not luminous. Since 5 M_⊙ stars are extremely bright, it is relatively easy to rule out massive stellar lenses once a lens mass is in hand. However, as described in Section 4.2, photometric microlensing alone is not sufficient to determine the mass of the lens in a microlensing event; astrometry is also needed to break degeneracies between lens mass and source and lens distances. Thus, data from a photometric survey must be combined with astrometric follow-up of a smaller set of BH candidates. As an example, we investigate the strategy of Lu et al. (2016) of selecting astrometric follow-up candidates with long Einstein crossing times t_E ≳ 120 days, high magnification ${A}_{\max }\gtrsim 10\equiv {u}_{0}\lesssim 0.1$, and high source flux fraction b_SFF ≈ 1 from photometric surveys such as the OGLE EWS (described in Section  5).

7.2.1. Selecting Black Hole Candidates with OGLE

Out of all the Mock EWS events, ∼1% of them have a BH lens. This is consistent with the fraction of BH lensing events in the unfiltered PopSyCLE event candidate list as well. While BHs make up ≲0.1% of objects in the Milky Way by number, their lensing rate is enhanced by their ∼10× larger Einstein radius. Thus, the observational selection criteria of the OGLE survey do not dramatically change the probability of detecting BH lensing events. In Figure 7, the various types of lenses contributing to the t_E distribution are shown, which is comparable to the  efficiency-corrected t_E curves produced by surveys. Figure 8 shows roughly the number of candidates that can be followed up based on EWS photometry per season. Although there usually are nearly 100 events with t_E > 120 days, often only 10 or so of those have sufficiently high signal-to-noise ratio (S/N) data and good coverage of the event, which are necessary later when trying to fit the event and extract the microlensing parameters t_E and π_E.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Selecting  Mock EWS events with t_E > 120 days increases the probability of the lens being a BH to 40% (Figure 9). Note that a different Galactic model can increase this fraction by up to a factor of two (see Appendix A). We then explored two additional selection criteria in addition to long events.

Zoom In Zoom Out Reset image size

First, we tried selecting on events from PopSyCLE simulations with b_SFF ≈ 1. Intuitively, one would expect a source to be less blended with a BH lens than with a stellar lens, as the BH does not contribute any flux. However, for ground-based surveys the seeing disk is so large that many background stars fall within the disk, and whether the lens is luminous or not does not significantly change b_SFF. In principle, additional high-resolution imaging from space- or ground-based telescopes with adaptive optics provides a much smaller seeing disk, reducing contamination from neighboring stars. This would then circumvent the aforementioned problem, meaning that BH lensing events would truly have a high b_SFF. However, this strategy fails, as OGLE microlensing events are observationally biased toward brighter sources and the average lens is quite faint (I ∼ 26). Thus, the distributions of b_SFF for BH and stellar lenses are similar in a sample limited to I ≲ 21.

Similarly, selecting events with small u₀ does not increase the probability of finding a BH lens. Although the average BH is more massive than the average star and thus has a larger θ_E, it does not have a correspondingly smaller u₀, as u₀ is independent of mass.

7.2.2. Individual Candidate Follow-up

As our ultimate goal is to measure the number of BHs in the Milky Way, we now consider how many BHs can be expected to be detected via astrometric follow-up. We estimate that  ∼12 BH candidates (corresponding to five expected BH detections) are required to constrain the total number of BHs in the Milky Way to 50%, assuming Poisson statistics (Figure 10). With an astrometric follow-up program, due to practicalities such as limited telescope time on facilities able to perform such measurements, three to four candidates can be observed each year, and the probability of the candidate being a BH is ∼40%; thus, about one to two BHs per year can be expected to be detected. Over 5 yr, this would result in a total of 5–10 BHs. However, the BH PDMF would be difficult to constrain based on current sample sizes and sporadic astrometric follow-up, due to the inefficiency of the process. Dedicated astrometric surveys are needed to place useful constraints on the BH PDMF.

Zoom In Zoom Out Reset image size

BH candidates identified from photometric microlensing surveys can be followed up with high-precision astrometric measurements in order to measure the astrometric microlensing shift, δ_c,max. When combined with t_E, π_E, and u₀ from the photometric light curve, the astrometric shift yields a constraint on the lens mass. In Figures 11 and 12, microlensing events are plotted in ${\delta }_{c,\max }-{\pi }_{{\rm{E}}}$ space. In Figures 13 and 14 the microlensing events are plotted in ${\pi }_{{\rm{E}}}-{t}_{{\rm{E}}}$ and ${\delta }_{c,\max }-{t}_{{\rm{E}}}$ space, respectively. We note that blending due to the lens (if it is luminous) is incorporated into the astrometric shift calculations, assuming a ∼100 mas aperture size that is roughly equivalent to infrared imaging with JWST or with an 8–10 m telescope and an adaptive optics system. BHs occupy a very particular region in ${\delta }_{c,\max }-{\pi }_{{\rm{E}}}$ space (large ${\delta }_{c,\max }$ and small π_E), as they are massive. Although both a massive luminous stellar lens and equally massive BH will have equal Einstein radii, all else equal, the astrometric shift will be larger for the BH, as it does not blend with the source images (see Equation (10) vs. Equation (11)).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The very sharp delineation between the stars, WDs, NSs, and BHs in Figure 12 can be simply explained. Consider a lens of mass M. Using Equations (3), (14), and (15) and the definition of κ given in Section 4.2, the microlens parallax will simply be

$\begin{eqnarray}&&{\pi }_{{\rm{E}}}=\sqrt{\displaystyle \frac{{\pi }_{\mathrm{rel}}}{\kappa M}}.\end{eqnarray} \tag{ 17 }$

Assuming an impact parameter of ${u}_{0}\leqslant \sqrt{2}$, the maximum astrometric shift, assuming no blending, can then be written as

$\begin{eqnarray}&&{\delta }_{c,\max }=\sqrt{\displaystyle \frac{\kappa M{\pi }_{\mathrm{rel}}}{8}}\end{eqnarray} \tag{ 18 }$

by combining Equations (3), (10), and 15. Thus, for a given mass M, both π_E and δ_c,max scale as $\sqrt{{\pi }_{\mathrm{rel}}}$. Hence, when plotting δ_c,max against ${\pi }_{{\rm{E}}}$ for a given mass lens, the slope is 1. Since PopSyCLE currently assumes that all NSs are of the same mass, they lie artificially on a straight line. Similarly, as there exists a minimum mass for BHs, WDs, and stars, there is a hard edge on the right diagonal side of those populations. There is some downward scatter in δ_c,max due to some events with $\sqrt{2}\lt {u}_{0}\lt 2$, as the maximum astrometric shift is less than what would be given in Equation (18). Additionally when blending is included, δ_c,max scatters lower for some stars, depending on the lens luminosity. The larger scatter from the slope of 1 in the BH, WD, and stellar populations is simply due to the fact that there are a range of masses.

Equations (17) and (18) can also be used to understand the t_E gradient stretching from the bottom right (large π_E, small δ_c,max) to the top left (small π_E, large δ_c,max) shown in Figure 11. Since t_E scales with the square root of the lens mass, heavier lenses will have smaller π_E, larger δ_c,max, and longer t_E, on average. Scatter in this relation is due to the fact that t_E is also degenerate with the relative source–lens proper motion μ_rel.

The microlensing parallax π_E, combined with a measurement of the Einstein crossing time t_E, appears to be a powerful means of picking out BHs (Figure 13). Unfortunately, π_E can only be observed significantly after the photometric peak and thus cannot be used as a BH-candidate selection criterion for astrometric follow-up. However, the combination of ${\pi }_{{\rm{E}}}$ and t_E can still be used to obtain a sample of microlensing events with a very high fraction of BH lenses, as compared to looking at t_E alone. By selecting events with both t_E > 120 days and a microlensing parallax of π_E < 0.08, the detection rate of BHs is  ∼ 85%; by using an even lower value of π_E, the minimum value of t_E could be shifted to lower values and still preserve the high fraction of BHs. This has the additional advantage that t_E and π_E are quantities that can be fit from photometry alone. Thus, with a set of photometric light curves, events can be sorted into BH and non-BH lenses with high statistical confidence.

This method is only useful if the π_E distribution does not vary dramatically as the Galactic model changes. It is not entirely clear how distinct this separation will be if different assumptions are made about the underlying distributions of M, d_L, d_S, and μ_rel. We evaluate the impact on the ${\pi }_{{\rm{E}}}-{t}_{{\rm{E}}}$ space by using the scaling relations

$\begin{eqnarray*}\begin{array}{rcl}{t}_{{\rm{E}}} & \propto & \displaystyle \frac{\sqrt{{\pi }_{\mathrm{rel}}M}}{{\mu }_{\mathrm{rel}}}\\ {\pi }_{{\rm{E}}} & \propto & \sqrt{\displaystyle \frac{{\pi }_{\mathrm{rel}}}{M}}\end{array}\end{eqnarray*}$

as illustrated by the arrows in Figure 13.

Of all the populations in PopSyCLE, the most uncertainty lies in the BH spatial positions, velocities, and masses. We consider the effects of changing these parameters. The following trends are illustrated with specific values in Figure 13.

• 1.  
  Our BH mass distribution ranges from 5 to 16 M_⊙; by including more high-mass BHs, such as the ones discovered by LIGO, in the distribution, the BH population would shift toward longer t_E and smaller π_E.
• 2.  
  The kick velocities of BHs at birth, if there are any, are unknown. Changes in kick velocity are manifested in changes in μ_rel. If μ_rel increases (holding d_L and d_S fixed), t_E becomes shorter; if μ_rel increases, t_E becomes  shorter.
• 3.  
  The spatial distribution of BHs is also not known. Two extremes could be considered: a centrally concentrated population of BHs in the bulge, or a distribution spread throughout the stellar or dark matter halo. Bulge lenses have smaller π_rel, while disk lenses have a larger π_rel. Thus, having only BH bulge lenses would cause the average π_rel to increase, causing both π_E and t_E to increase. On the other hand, with a noncentrally concentrated distribution of BHs, this would cause more lenses to fall closer to Earth, meaning that the average d_L and thus ${\pi }_{\mathrm{rel}}$ would decrease, causing π_E and t_E to decrease.

Although this method of selecting BH lenses does not allow for direct confirmation or mass measurement, it does allow for a statistical analysis of BH lensing events. For example, the number of BHs in the Milky Way could be constrained; their masses could also be estimated by invoking some type of Galactic model, as is commonly done. Its strength is that no additional data are necessary and the improved BH selection process can still lead to improved physical constraints.

The challenge associated with this method is constraining π_E accurately down to the level that is necessary for this measurement. However, with high-cadence observations and high photometric precision, this is possible.

Lastly, we consider future microlensing surveys and the prospects they hold for measuring BH masses. As described in Section  5, the WFIRST mission includes a microlensing survey designed for exoplanet detections. Some primary differences between the OGLE and WFIRST surveys will be their filters, sensitivity, and resolution (seeing limited I ≲ 22 vs. diffraction limited H ≲ 26, respectively). An additional important difference is that with WFIRST, the astrometry will be obtained at the same time as the photometry. Thus, it is not necessary to do individual follow-up of candidate BH events; all the information will be contained within the WFIRST survey data. Simultaneous astrometry also has the advantage that it is possible to go "back in time" to look at astrometric data from before the photometric peak (i.e., before the microlensing event is recognized photometrically).

P19 presents thorough simulations and a detailed analysis to calculate the expected yield of bound planet detections with WFIRST microlensing. We present here a complementary analysis with PopSyCLE to determine the number of BHs we may expect a WFIRST-like survey to find. Detailed simulations of different season distributions and understanding how the BH yield changes are beyond the scope of this paper and are left as future work. However, our continuous 1000-day survey simulation with Mock WFIRST parameters (Table 4) is a good order-of-magnitude estimation, considering the other uncertainties in the PopSyCLE simulation and the WFIRST survey design itself.

To normalize the Mock WFIRST survey area to that of the actual WFIRST microlensing survey, the number of events is multiplied by the ratio of the areas (a factor of 1.93). Microlensing events from the Mock WFIRST simulation with BH lenses and an Einstein crossing time of 90 days < t_E < 300 days are defined to have "measurable" BH masses. Note that there will certainly be many BH microlensing events with t_E that fall above or below this range; however, we do not consider those to have measurable masses, for the following reason. As described in Section 4.2, in addition to the photometric brightening, it is necessary to measure the astrometric shift and microlensing parallax to constrain the lens mass. The choice for the lower bound on t_E takes into consideration that an event needs to have a somewhat long duration (${t}_{{\rm{E}}}\gtrsim 3$ months) to have a potentially measurable microlensing parallax π_E. The choice for the upper bound on t_E takes into consideration that the astrometric signal falls off much more slowly than the photometric signal; a rough guess is that to measure δ_c,max requires astrometric data for ∼5t_E–10t_E. Note that the amount of astrometric data required also is dependent on when the astrometric measurements occur. To properly determine this requires detailed simulation (e.g., Section 8 of Lu et al. 2016) and is beyond the scope of this paper.

With this definition, the PopSyCLE Mock WFIRST survey produces  ∼1000 BH lensing events with measurable masses, nearly 100 times more than with an individual astrometric follow-up program, as estimated in Section 7.2.2. With this number of mass measurements, we can constrain the present-day BH mass function (Figure 15) and the results of supernova simulations.

Zoom In Zoom Out Reset image size

A major source of uncertainty in this estimate is due to the Galactic model, in both the stellar and compact object components. As described in Appendix A, we consider two different angles (α = 111 and α = 28°) for the Galactic bar/bulge; this modification significantly changes the number of stars along a given line of sight. The more tilted bar with α = 111 produces 4–5 times less microlensing event candidates than the less tilted bar with α = 28°. By applying the Mock WFIRST criterion to the six fields in Table 7, the number of events with measurable BH masses is decreased by a factor of about 4. To take into account these uncertainties, we estimate that the number of BHs with measurable masses is ${ \mathcal O }(100\mbox{--}1000)$.

Table 7.  Comparing PopSyCLE to Surveys (Galaxia v3)

Field	${n}_{s}({10}^{6}/{\deg }^{2})$		${\rm{\Gamma }}({10}^{-6}\,{\mathrm{star}}^{-1}\,{\mathrm{yr}}^{-1})$		$\langle {t}_{{\rm{E}}}\rangle (\mathrm{days})$
	Obs.	Mock	Obs.	Mock	Obs.	Mock
F00M19	2.62	1.48	1.3 ± 0.8	4.35	32.6	14.8
F01M19	4.54	2.04	5.5 ±  0.9	3.69	39.5	46.2
F03S	3.64	2.47	${14.0}_{-2.4}^{+2.9}$	8.26	25.5	20.3
F11M19	4.95	3.66	16.2 ± 1.3	17.88	21.8	19.7
F12M19	3.26	1.49	3.4 ± 1.1	2.88	36.7	51.1
F13M19	4.51	2.19	5.2 ± 1.1	4.9	30.8	28.2

Note. Identical analysis to that presented in Table 5, but using Galaxia v3 (with the tilted shorter bar) instead of v2. Fields with ^M19 indicate that the observed values come from M19, while those with ^S are from SP16.

Download table as:  ASCIITypeset image

8.1.1. Microlensing Simulations

8.1.2. Compact Object Population Synthesis

We advocate for defining a microlensing event using quantities that are directly observable/measured rather than fit. This is described in Dominik (2009) as a reparametrization when fitting microlensing events. For example, in the survey papers discussed, cuts were made on events whose source magnitude, m_S, is fainter than some value, where m_S is generally set by the magnitude limit of the telescope and camera. However, the more easily observed quantity is the baseline magnitude m_base since m_S is generally a parameter derived from fitting and is subject to degeneracy with other parameters (Dominik 2009). When fitting light curves, blending is degenerate with the observed timescale and peak magnification of the event, and if blending is not taken into account correctly, the microlensing parameters derived will be incorrect; masses and timescales will be systematically underestimated (Di Stefano & Esin 1995; Woźniak & Paczyński 1997; Han 1999). In other words, events with small t_E, large u₀, and large b_SFF are degenerate with events with large t_E, small u₀, and small b_SFF (Sumi et al. 2011). In order to facilitate easier comparison between simulations and observations and explore this degeneracy, we recommend that future microlensing analyses adopt cuts using observable quantities, as advocated for in Dominik (2009). However, we also suggest that such observational cuts be applied for sample selection, as these are more easily reproduced and are less dependent on differences in fitting codes and prior assumptions.

 Additionally, although the mean Einstein crossing time $\langle {t}_{{\rm{E}}}\rangle$ is the most commonly reported parameter by surveys, the t_E distribution is asymmetric in linear t_E bins;¹⁰ thus, the mean is easily skewed and depends on the range of t_E used to estimate the mean. For future comparisons across observational surveys and simulations, the median is a better choice and is less impacted by particular cuts. It is also important to note the range of t_E over which the distribution is made (e.g., as in Sumi et al. 2013).

Over the years, many types of follow-up observations have been proposed and used to constrain lens properties. Agol & Kamionkowski (2002) proposed to constrain lens masses of BH microlensing candidates by searching for their X-ray emission due to accretion from the interstellar medium or stellar winds. Nucita et al. (2006) and Maeda et al. (2005) used XMM-Newton and Chandra observations, respectively, to search for X-rays from the extremely long-duration BH candidates MACHO-96-BLG-5 (Mao et al. 2002); however, no significant detections were made.

HST high-resolution imaging follow-up of BH candidates has also been used to measure the degree of blending (Bennett et al. 2002). As discussed in Section 7.2.1, this is not particularly useful for selecting photometric candidates for astrometric follow-up; however, it still may be a useful means of confirming that the lens is not a massive star. Poindexter et al. (2005) also suggested that HST observations could be used to measure the source proper motion; however, they note that many degeneracies remain even with this measurement. High-resolution images can be taken many years after an event, when it may be possible to resolve the source and lens (Batista et al. 2015; Bennett et al. 2015; Alcock et al. 2001; F. N. Abdurrahman et al. 2019, in preparation) and the absence of a lens may provide additional confirmation of a BH.

For nearby sources, it is sometimes also possible to measure the source distance via parallax, which would give complete event parameters. The combination of photometry and astrometry is powerful not only for BH lenses but also for obtaining precise mass measurements of any type of lens, as was shown with the first astrometric microlensing signal detected outside our solar system (Sahu et al. 2017). Another method of breaking degeneracies is to have the ability to resolve the images themselves, which is possible for stellar-mass lenses using interferometric techniques (Dong et al. 2019).

We have not yet considered the number of BHs detectable with the Gaia satellite, an astrometric space mission that has been operating since 2013 (Gaia Collaboration et al. 2016). Although Gaia has incredible astrometric precision for parallax measurement (down to ∼10 μas for some stars by the end of the mission), it does not perform as well for single-epoch astrometric measurements in the bulge owing to the decreased observing cadence, as well as significant stellar crowding (>1 mas; Rybicki et al. 2018). Additionally, Gaia observes at green optical wavelengths, which cannot probe dusty, high-extinction regions like the bulge. However, there are estimates of the number of BHs that Gaia will be able to detect over its lifetime. Mashian & Loeb (2017) estimated that ∼2 × 10⁵ BHs in astrometric binary systems can be detected over Gaia's 5 yr mission; however, this is a severe overestimate, as this calculation has neglected extinction and crowding effects. Rybicki et al. (2018) estimated that a few isolated stellar-mass BHs should be detectable astrometrically at the end of the Gaia mission. Wyrzykowski & Mandel (2019) used Gaia DR2 data to calculate distances and proper motions for sources in OGLE-III microlensing events from 2001 to 2009, providing additional information to perform a more careful reanalysis of the lens masses to determine whether they could be BHs.

In Section 7, two different methods to better understand BHs in the Milky Way are discussed. The first involves using a combination of astrometry and photometry to measure t_E, ${\pi }_{{\rm{E}}}$, and δ_c,max, which allows lens masses to be measured for individual microlensing events. The second involves photometry only to measure t_E and π_E, enabling statistical constraints on a population of BHs. Although the BH nature of individual lenses cannot be confirmed, ensemble information can still be gleaned. To obtain masses with the second method would involve assuming some type of Galactic model or spatial distribution. It is interesting to note that searching for a small/undetectable parallax signal to identify BHs is the opposite of previous approaches in the literature (Poindexter et al. 2005; Wyrzykowski et al. 2016; Drlica-Wagner et al. 2019).

In the next version of PopSyCLE, we plan to add support for stellar and compact object binary systems, compact object mergers, metallicity-dependent IFMRs, different compact object spatial distributions, and primordial BHs. This will allow for a more realistic simulation and a larger exploration of parameter space. The lack of these features at the present put some caveats on this work, which we discuss.

8.4.1. Binarity and Mergers

8.4.2. Metallicity

8.4.3. Compact Object Distributions and Kinematics

There are three specific properties of the bulge we consider modifying:

• 1.  
  Pattern speed: Galaxia implements a pattern speed of Ω = 71.62 km s⁻¹ kpc⁻¹ (see Section 3.3 in Sharma et al. 2011). However, this is nearly twice as fast as the most recent values reported in the literature, which range from around 36 to 44 km s⁻¹ kpc⁻¹ (Bovy et al. 2019; Clarke et al. 2019; Sanders et al. 2019), determined using combinations of Gaia DR2, VVV, and APOGEE.
• 2.  
  Velocity dispersion: Galaxia implements velocity dispersions of σ_R = σ_ϕ = 110 km s⁻¹. However, this produces microlensing events with timescales that are too short, which suggests that smaller velocity dispersions might be more appropriate.¹¹ In reality, σ_R and σ_ϕ should vary with latitude and longitude (Howard et al. 2008, 2009), but to implement this in Galaxia would require significant rewriting of the code, which is far beyond the scope of this current work.
• 3.  
  Bar angle and length: Galaxia implements a bar angle of α = 111, where α is the angle from the Sun–Galactic center line of sight, and a major-axis scale length of x₀ = 1.59 kpc. It is suggested that the angle should be closer to α = 28° (Wegg & Gerhard 2013; Wegg et al. 2015), with a shorter scale length of 0.7 kpc. In fact, Portail et al. (2017) performed dynamical modeling using α = 28° to find a bar pattern speed of around 40 km s⁻¹ kpc⁻¹. However, there is still considerable debate in the literature about the value of α and the scale length, with values for α ranging from 10° to 45° (see Robin et al. 2012 for a summary and references).

We create two new variations of the original Galactic model implemented in Galaxia. We dub "v2" to be the version of Galaxia with a pattern speed Ω = 40 km s⁻¹ kpc⁻¹ and bulge velocity dispersions σ_R = σ_ϕ = 100 km s⁻¹. We dub "v3" to be the version where the bar is short and tilted, with ${\rm{\Omega }}=40$ km s⁻¹ kpc⁻¹, bulge velocity dispersions σ_R = σ_ϕ = 100 km s⁻¹, bar angle α = 28°, and major-axis scale length x₀ = 0.7 kpc.

Figure 16 compares the Einstein crossing time distribution of Galaxia assuming either the original Galactic model or the modified v2 model. Tables 5 and 7 compare PopSyCLE stellar densities, event rates, and Einstein crossing times to results from SP16 and M19 for v2 and v3, respectively. Figure 17 summarizes the results of the two tables together. The stellar densities of v2 match reasonably well with the observed number counts; v3 is consistently too low. Note that although in projection the length of the bar is the same, as $\sin (11\buildrel{\circ}\over{.} 1)\cdot 1.59\,\mathrm{kpc}\approx \sin (28^\circ )\,\cdot 0.7\,\mathrm{kpc}$ (P19), the number counts are quite different, due to extinction. The event rates for v3 match reasonably well with the observed rates; v2 is consistently too high. With respect to the average Einstein crossing time, it is not as clear whether v2 or v3 matches the observed values better.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Based on this analysis, we chose to use v2 instead of v3 for the analysis in the main text of this paper. Stellar density is the more fundamental observational quantity; hence, we prefer that the simulation reproduce this aspect accurately. The event rate is microlensing specific and dependent on many more factors (e.g., the detection efficiency correction). However, it is curious in and of itself that the tilted bar creates agreement in one regime but not the other. Additional modifications (such as the 3D $E(B-V)$ map) may be necessary to bring the models into better agreement with observation. We do not explore this further here, and we leave detailed Galactic modeling to future work. However, it is worth noting how this uncertainty affects some of the results of this paper. In Section 7.2.1, the fraction of BH events at long times can be up to a factor of 2 higher for v3 than for v2. In Section 7.5 we consider the number of BH masses WFIRST can measure; the number is about a factor of 4 higher for v2 than for v3.

The rest of the analysis and validation performed in this appendix also uses v2.

Following Section 6.2.2 of P19, we compare the results of our simulation to observational studies of bulge kinematics. In Clarkson et al. (2008), a study of proper motions in the Galactic bulge is performed using observations of the HST SWEEPS field (centered at $(l,b)=(1.25,-2.65)$ covering an area of 11 arcmin²) in the HST F606W and F814W bands. To compare, we created a synthetic survey in the same direction of the same area using Galaxia; we use I and R bands, as these have the closest effective wavelength to the HST F606W and F814W bands. First, we select stars in Galaxia photometrically (Figure 18), similarly to Clarkson et al. (2008), to obtain a red population (bulge proxy) and a blue population (disk proxy). We obtain 347 blue stars and 699 red stars (compare with P19's 37 blue stars and 105 red stars, drawing from an area of 1.44 arcmin²).

Zoom In Zoom Out Reset image size

We then compare our results to those of Clarkson et al. (2008) and P19, who use a different version of the Besançon model, named BGM1106 for short. The results are summarized in Table 8 and illustrated in Figures 19 and 20. In particular, the match with Clarkson et al. (2008) is improved in ${\sigma }_{l,* }$ for the red population.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 8.  Clarkson et al. (2018), BGM1106, and Galaxia Proper-motion Comparison

Model/Data	${\rm{\Delta }}{\mu }_{l* }$	${\rm{\Delta }}{\mu }_{b}$	${\sigma }_{l* }$, Blue	σb, Blue	σl*, Red	σb, Red
Clarkson et al. (2008)	3.24 ± 0.15	−0.81 ± 0.12	2.2	1.3	3.0	2.8
BGM1106 (P19)	3.53 ± 0.65	−0.12 ± 0.32	2.47 ± 0.29	1.11 ± 0.13	5.19 ± 0.36	2.64 ± 0.18
Galaxia	4.51	−0.42	2.48	1.42	3.74	2.32

Note. All units are mas yr⁻¹. Δ is defined as blue − red.

Download table as:  ASCIITypeset image