The initial-final mass relation (IFMR) maps the initial mass of stars on the main sequence to the mass of their compact remnants. The form of the IFMR is not well determined due to limitations in both theory and observation, and thus represents a very active area of research in stellar physics (Heger et al. 2003;Ertl et al. 2016;Lu et al. 2019;Costa et al. 2021). The final mass of a compact object depends not only on its zero-age main sequence (ZAMS) mass, but also on factors such as rotation (which causes mixing in a star), metallicity (which governs mass-loss rates due to stellar winds, particularly during the post-main sequence), multiplicity (stars in close binaries can evolve differently), and core structure just prior to explosion (which can determine what type of compact object will be formed; see Sukhbold et al. 2018, and references therein).
While the IFMR for low-mass stars is measured using direct observation of white dwarfs (WDs; Kalirai et al. 2008;Cummings et al. 2018) the IFMR for high-mass stars is an open question, due to difficulty in obtaining mass measurements of massive compact remnants. While it is possible to measure the mass of some isolated, young neutron stars (NSs) using pulsar timing (Lorimer 2008) or measure the mass of NSs in X-ray binaries (Steiner et al. 2013), obtaining mass measurements of a large number of isolated stellar-mass black holes (BHs) presents an even greater challenge.
One method to increase the number of mass measurements for dark, isolated compact objects is to use gravitational microlensing. When a massive object passes in front of a luminous background source, a transient brightening and positional shift of the background star occurs. While photometric observations are sufficient to detect microlensing events, precise astrometric observations of the centroid shift during the microlensing event are required to break degeneracies and measure the lens mass directly (Lu et al. 2016;Sahu et al. 2017). The first detection of an isolated dark object using photometric and astrometric microlensing was recently reported; the analysis of Lam et al. (2022) suggests the object is a NS or low-mass BH, while the analyses of Sahu et al. (2022) and Mróz et al. (2022) suggest a BH is the only possibility.
To date, thousands of photometric microlensing events have been detected by dedicated microlensing surveys such as the Optical Gravitational Lensing Experiment (OGLE; Udalski et al. 1992) and the Microlensing Observations in Astrophysics (Muraki et al. 1999), as well as other astronomical surveys such as the Zwicky Transient Facility (Bellm et al. 2019) and Gaia (Wyrzykowski et al. 2022). However, only a few tens of these photometric events have also been observed astrometrically, as the follow-up process is extremely resource intensive and few facilities have the requisite precision. In the absence of mass measurements due to lack of astrometric microlensing measurements, statistical constraints placed on the IFMR from photometric microlensing alone are the best way to compare our theoretical models for latestage stellar evolution and stellar death with observation.
This work explores whether or not photometric microlensing surveys will be able to place meaningful constraints on the IFMR. In Section 2 we discuss the Population Synthesis for Compact object Lensing Events (PopSyCLE) microlensing simulation and the modifications that we have made to it for this work. In Section 3 we discuss the implementation and characteristics of different IFMRs added to the Stellar Population Interface for Stellar Evolution and Atmospheres (SPISEA) simple stellar population synthesis code. In Section 4 we present our findings on the effect of the different IFMRs on the BH microlensing population (Section 4.3), and whether or not these differences are detectable with OGLE (Section 4.4) or the upcoming Roman Space Telescope's microlensing survey (Section 4.5). In Section 5 we compare the Galactic BH distribution predicted by the different IFMRs to the extragalactic BH distribution detected via gravitational-wave mergers, as well as discuss further enhancements to PopSyCLE. We finish in Section 6 with a summary of our main conclusions.
PopSyCLE is a microlensing population synthesis tool for the Milky Way (Lam et al. 2020). Given a survey location and area as well as other parameters like the length of the survey, number of observations, reddening law, and filter, PopSyCLE will return a list of observable microlensing events and the parameters associated with them (e.g., Einstein crossing time, microlensing parallax, magnitude of brightening). In this work we make several modifications to PopSyCLE in order to explore the effect of different IFMRs on the Milky Way BH microlensing population, which are described in detail below. For a full description of PopSyCLE, see Lam et al. (2020).
Galaxia is a resolved stellar simulation of the Milky Way (Sharma et al. 2011), which serves as the foundation for stellar population synthesis in PopSyCLE. Given a survey area and location Galaxia will return all the stars located in the conical volume of the projected circular survey area centered on the specified coordinates. Compact objects are not included in the output of Galaxia. Each star returned by Galaxia has a position, velocity, mass, age, metallicity, among other parameters.
The Galaxia stellar simulation is based on 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 summary of the relevant distributions from which the stellar parameters are drawn for various populations of stars and a brief description of each population, see Tables 123in Sharma et al. (2011).
In the PopSyCLE simulations presented in this paper, we use the "v3" Galactic model, described in Appendix A of Lam et al. (2020), which differs from the default Galactic model of Galaxia. The most salient change is to the Galactic bar. The angle of the line connecting the Sun and Galactic center, α, and the major axis scale length of the Galactic bar, x 0 , are changed from α = 11.1°to α = 28°and x 0 = 1.59 kpc to x 0 = 0.7 kpc (Wegg & Gerhard 2013;Wegg et al. 2015). All simulations discussed in the main text of this paper are run using these parameters.
Galaxia is a stellar survey, which means it does not include compact objects. In order to include compact objects in our microlensing simulation we must inject them from another source. SPISEAfoot_0 (Hosek et al. 2020) is a software package which generates a single-age, single-metallicity stellar population (i.e., a star cluster) based on user-controlled parameters such as the initial mass function, evolution models, atmosphere models, extinction maps, and multiplicity distributions. The MESA Isochrones and Stellar Tracks (MIST) stellar evolution models (Choi et al. 2016) are used to evolve the SPISEA clusters to determine which progenitors have left the post-main sequence and become compact objects. The IFMR is then used to assign a remnant type (WD, NS, BH) and mass to the compact objects. For details on the SPISEA input parameters used in the PopSyCLE simulation, see Section 2.2 of Lam et al. (2020).
In the SPISEA code the user can choose the IFMR used. Prior to this work, SPISEA contained a single IFMR object, hereafter referred to as Raithel18, based on Raithel et al. (2018), for BHs and NSs, and Kalirai et al. (2008) for WDs.
This work adds two additional IFMR objects to SPISEA.foot_1 One, called Spera15, is based on the Stellar EVolution N-body (SEVN) codefoot_2 (Spera et al. 2015), which is described in Section 3.3. The Spera15 IFMR is a function of progenitor metallicity as well as ZAMS mass, but does not take into account models of explosion physics like Raithel et al. (2018). The other new IFMR, based on simulations by Sukhbold et al. (2016) and Sukhbold & Woosley (2014), is called Sukh-boldN20 and is described in Section 3.4. The SukhboldN20 IFMR includes metallicity dependence and the explosion physics of Raithel et al. (2018). Both new IFMR objects also use the Kalirai et al. (2008) WD IFMR described in Section 3.1 for low-ZAMS-mass stars.
In the original PopSyCLE simulation, to perform the population synthesis of compact objects using SPISEA, the stars from Galaxia are binned according to population (i.e., thin disk, thick disk, bulge) and age (see Section 3 of Lam et al. 2020 for more details). For this work, each age bin is further divided by metallicity. This is necessary as SPISEA can create only single-age, single-metallicity populations. For PopSyCLE runs using the Raithel18 IFMR, which has no metallicity dependence, the only metallicity sub-bin is at solar metallicity and encompasses the full range of metallicities produced by Galaxia.
In order to determine the most appropriate bins for simulations using the Spera15 and SukhboldN20 IFMRs, we first look at the metallicity distribution at two different pointings in the Milky Way using Galaxia (Figure 1). We chose metallicity bins such that the number of distinct metallicity isochrones required by SPISEA (which must be generated and stored in order to perform population synthesis) is minimized in order to save on disk space and computational time, and such that the distribution of BH masses given by the Spera15 IFMR did not have arbitrary mass gaps introduced by the metallicity binning.
Based on these constraints, we have chosen four metallicity bins for use with the Spera15 and SukhboldN20 IFMRs. The first bin contains stars with [Fe/H] < -1. In future work, modifications made to the SPISEA population synthesis code will allow for interpolation between isochrones, which PopSyCLE will be able to take advantage of with finer metallicity binning.
One other change that has been made to the original PopSyCLE simulation is the addition of a more realistic distribution of birth kick velocities for NSs and BHs. Birth kick velocities are an additional velocity resulting from asymmetries in supernovae explosions, where excess mass loss in one direction or anisotropic neutron emission will result in an additional velocity in a random direction for the compact object left behind by the explosion (Janka & Mueller 1994;Kusenko & Segre 1996;Tamborra et al. 2014).
In the original version of PopSyCLE, each BH or NS was assigned a constant birth kick velocity (350 km s -1 for NSs and 100 km s -1 for BHs) in a random direction, which was then added to the existing "stellar" velocity assigned based on the distribution of stellar velocities coming from Galaxia (see Section 5 of Lam et al. 2020). The value chosen for the NS kick velocity was the average of a Maxwellian distribution reported based on observations of pulsar proper motions in Hobbs et al. (2005). In the current version of PopSyCLE, we now implement a more realistic Maxwellian kick distribution for the NS and BH populations, instead of applying a single-valued kick as before.
Because NS birth kicks follow a Maxwellian distribution, we might also expect BH birth kick velocities to be Maxwellian. The Maxwellian birth kick velocity distribution used in the simulations for this paper have averages which match the original PopSyCLE values for NS and BH birth kick velocities, respectively (i.e., the average of the randomly assigned values for NS kick velocity drawn from the Maxwellian distribution was chosen to be 350 km s -1 , while the average of the kick velocities assigned to the BHs was chosen to be 100 km s -1 ).
The WD IFMR used in this work is based on Kalirai et al. (2008). This IFMR is used in all of the SPISEA IFMR objects to derive WD masses for low-ZAMS-mass progenitors, except those young, luminous WDs already included in the MIST models. The Kalirai et al. (2008) WD IFMR is empirically determined based on observational data in the initial mass range 1.16 M e < M ZAMS < 6.5 M e , and is given by
We extend the range of this IFMR to 0.5 M e < M ZAMS < 9 M e for the Raithel18 and SukhboldN20 IFMR objects and to 0.5 M e < M ZAMS < 7 M e for the Spera15 IFMR object, where the upper mass range is chosen such that the IFMR is defined for all ZAMS masses below the lower limit of the NS/ BH IFMRs. The lower mass range is chosen to match the lower mass limit of the MIST model objects. See Section 2.2.2 of Lam et al. (2020) for more details.
The original PopSyCLE simulation as described in Lam et al. (2020) implemented a BH/NS IFMR based on Raithel et al. (2018) as well as the WD IFMR from Kalirai et al. (2008) described above in Section 3.1. The model includes a stochastic process to determine whether or not a BH or NS is ultimately formed from each progenitor, but does not include metallicity dependence. Because the simulations used to produce the Raithel18 IFMR assumed solar metallicity, all compact objects are assumed to be have solar-metallicity progenitor stars in the original PopSyCLE simulation. For a full set of the equations and a description of how the Raithel18 IFMR was implemented, see Appendix C of Lam et al. (2020).
As the Raithel18 IFMR is based entirely on progenitors stars with solar metallicity, it is missing the most massive compact objects formed from the low-metallicity population of stars in the Milky Way (Meng et al. 2008). The masses of BHs formed using the Raithel18 IFMR range from around 5 to 16 M e .
In the original PopSyCLE simulation every NS produced by the Raithel18 IFMR was assumed to have a mass of 1.6 M e . In the version of PopSyCLE used for this paper, the mass distribution of NSs produced by the Raithel18 IFMR is instead drawn from a Gaussian distribution with an average 1.36 M e and a standard deviation 0.09 M e , based on a compilation of NS masses from several observational studies (Appendix A).
The IFMR is dependent on the metallicity of progenitor stars as well as on their ZAMS mass (Heger et al. 2003;Meng et al. 2008). The SEVN code is a a software package that models latestage stellar evolution and supernovae physics including the effects of metallicity-dependent mass loss (Spera et al. 2015). This software, in addition to the stellar evolution models of the PAdova and TRieste Stellar Evolution Code Vertical black lines demarcate the chosen metallicity binning for the SPISEA clusters run with a metallicity-dependent IFMR, with the metallicity assigned to each bin labeled in red. Note that for the Raithel18 IFMR all SPISEA clusters are solar metallicity only. (Bressan et al. 2012(Bressan et al. , 2013;;Chen et al. 2014;Tang et al. 2014;Chen et al. 2015) were used by Spera et al. (2015) to create an analytical formula for stellar remnant mass as a function of ZAMS mass and the mass fraction of metals Z. This analytical formula can be found in Appendix C of Spera et al. (2015).
Note that the Spera15 IFMR takes as its argument the mass fraction in metals Z rather than [Fe/H] as returned by Galaxia. To convert between Z and [Fe/H] we use the equation
Fe H
where Z e = 0.014 (Ekström et al. 2012).
We used this analytical formula to create an IFMR object in SPISEA, hereafter referred to as Spera15. The function is defined for M ZAMS 7 M e and all metallicities. It works by first calculating the core mass based on progenitor ZAMS mass and metallicity, and then uses the core mass to calculate the final remnant mass. Only once the remnant mass has been calculated is the object determined to be a BH, NS, or WD. We use the Chandrasekhar mass of 1.4 M e as the limit for WD masses, i.e., all Spera15 remnants with masses less than 1.4 M e are assigned to be WDs. All Spera15 remnants with masses between 1.4 and 3 M e are assigned to be NSs, where the upper limit of a NS's mass is roughly estimated from Özel & Freire (2016). All Spera15 remnants with masses greater than 3 M e are assigned to be BHs. One issue with this method of assigning remnant types is that the outcome is continuous and deterministic based on initial mass. Several simulations, including those by Sukhbold et al. (2016), suggest that there is no mass above which stars become BHs and below which they become NSs. In addition, this method does not produce NSs with masses less than the Chandrasekhar mass, when in reality a significant fraction of NSs are less than 1.4 M e . The result is that the average mass of the compact objects generated is too high.
Because the Spera15 IFMR includes low-metallicity remnant populations, it will return more massive BHs than the Raithel18 IFMR. As compared to the fairly narrow mass distribution of BHs allowed by Raithel18 (5-16 M e ), the Spera15 IFMR allows for BHs as massive as 90 M e to form (Figure 2).
The SukhboldN20 IFMR is shown in Figure 3 and includes zero-metallicity models from Sukhbold & Woosley (2014), solar-metallicity models from the N20 set of Sukhbold et al. (2016), and pulsational-pair instability models from Woosley (2017) and Woosley et al. (2020). We collectively refer to these as the SukhboldN20 simulations.
Any object with remnant mass M rem < 3 M e is a NS, and those with M rem > 3 M e are BHs.
To obtain the BH IFMR, we use least-squares minimization to find the best-fit line through BH masses M BH for BHs in the zero-metallicity N20 models, where 15 M e < M ZAMS < 70 M e . These best-fit lines are shown in Figure 3, labeled as Z = 0.0 Z e and Z = 1.0 Z e , respectively. The zero-metallicity models are fit by
remnant mass is independent of metallicity, and Equation (3) is applicable for all stars with M ZAMS 40 M e . However, above M ZAMS ∼ 40 M e , the relationship between M ZAMS and M rem becomes metallicity dependent. We fit a line using least-squares minimization to the BH points between 40 M e < M ZAMS < 70 M e in the solar-metallicity N20 models. The solar-metallicity N20 models above M ZAMS ∼ 40 M e are fit by
To obtain the remnant mass as a function of metallicity for M ZAMS 40 M e , we linearly interpolate between the zero-and solar-metallicity models (Equations ( 3) and (4)). The choice of linear interpolation is arbitrary, as the uncertainties in massive stellar evolution are so large the proper interpolation scheme is unknown. However, the main trend is captured: for zero metallicity, the mass of the remnant BH always increases, while for solar metallicity, mass loss eventually catches up for the high-mass stars and the remnant mass decreases, with intermediate metallicities having behavior in between those two cases.
In addition, we include pulsational-pair instability supernovae (PPISN) by extrapolating Equations (3) and (4) out to 120 M e . For high-mass stars, 100 M e M ZAMS 140 M e , electron-positron pair production robs the star of energy, causing it to eject a substantial fraction of its mass, resulting in a core-collapse supernova and remnant BHs of ∼35-50 M e . The M ZAMS corresponding at which a star reaches the PPISN region increases with increasing metallicity. At sufficiently high metallicities, there is no M ZAMS that reaches the PPISN region.
Putting all the above together, the N20 BH IFMR is given by the following piecewise function:
where
The condition on M up (Equations ( 7) and ( 8)) ensures M BH 3 M e , while also restricting the maximum M ZAMS to 120 M e . We also assume supersolar-metallicity stars have the same behavior as solar-metallicity stars (Equation ( 6)).
Note also that the BH masses from Equation (5) are lower limits on the BH masses from the simulation, as they correspond to the helium (He) core mass of the star at the time of implosion. It is uncertain whether the envelope is entirely ejected or whether some of it falls back onto the core and contributes to the BH mass. The difference between the He core mass and the total pre-supernova mass provides the limits of the remnant BH mass (Sukhbold et al. 2016;Raithel et al. 2018). Compared to observational data of BHs, Raithel et al. (2018) found a high ejection fraction of the envelope and that the remnant BH mass was quite similar to the He core mass for solar-metallicity stars; however, it is not known whether this result is metallicity dependent.
Simulations indicate there is no M ZAMS above which BHs are always formed and below which NSs are formed (Sukhbold et al. 2016). Other factors such as the metallicity of the star or core structure directly prior to explosion also determine what type of remnant is left behind. To include this stochasticity in the SukhboldN20 IFMR, we assign the different outcomes probabilistically based on M ZAMS . The probabilities for NS versus BH formation are taken from the SukhboldN20 simulations. We follow the approach of Raithel et al. (2018), by choosing the fewest number of bins possible to capture the different probability regions where no/some/only BHs are formed. Below M ZAMS ∼ 60 M e , the probability of NS or BH formation are independent of metallicity. Similarly to the M BH -M ZAMS relationship, to determine the probability of a BH or NS remnant above M ZAMS ∼ 60 M e , we interpolate linearly between the zero-and solar-metallicity models (Table 1). For M up M ZAMS 120 M e , the probability of BH formation is zero. Although we take the NS formation probabilities from the SukhboldN20 simulations, we do not use the masses. Instead, we draw from the same NS mass distribution as for the Raithel18 IFMR (Appendix A).
Using the updated version of PopSyCLE, we simulated a total of nine 0.34 deg 2 fields, with each field centered at the location of an OGLE-IV bulge field (Table 2). Each simulation was 1000 days long, with a sampling cadence that detected all events with t E 3 days. 8 The chosen fields span Galactic longitudes 0°-8°and Galactic latitudes -4°-5°.
The different IFMRs presented in this paper result in different underlying mass distributions for the Milky Way compact object populations. Figure 4
Table 1 Compact Object Formation Probabilities Mass Range (M e ) P WD P NS P BH 0.5 M ZAMS < 9.0 1.00 0.00 0.00 9.0 M ZAMS < 15.0 0.00 1.00 0.00 15.0 M ZAMS < 21.8 0.00 0.75 0.25 21.8 M ZAMS < 25.2 0.00 0.00 1.00 25.2 M ZAMS < 27.4 0.00 1.00 0.00 27.4 M ZAMS < 60.0 0.00 0.00 1.00 60.0 M ZAMS < M up 0.00 0.8f Z 1-0.8f Z M up M ZAMS 120.0 0.00 1.00 0.00 Note. Probabilities of forming a white dwarf P WD , neutron star P NS , or black hole P BH as a function of ZAMS mass, M ZAMS , for the SukhboldN20 IFMR. f Z is defined in Equation (6) and M up is defined in Equation (7).
object populations produced by the different IFMRs in the simulated OGLE-IV-BLG611 field.
The WD population does not change significantly since all the IFMRs for low-mass stars are based on Kalirai et al. (2008). The average mass of a WD in field OGLE-IV-BLG611 is 0.65 M e , 0.66 M e , and 0.64 M e , for the Raithel18, Spera15, and SukhboldN20 IFMRs, respectively.
The NS population produced by the Spera15 IFMR tends to be about 30% more massive than those produced by the Raithel18 and SukhboldN20 IFMRs. Although many NSs have masses below 1.4 M e (Appendix A), the Spera15 IFMR labels all compact objects below 1.4 M e as WDs, only allowing NSs in the 1.4-3 M e mass range (Section 3.3) This causes the NSs from the Spera15 IFMR to be more massive on average than the NS populations of Raithel18 and SukhboldN20, both of which have their masses drawn from a Gaussian distribution with a mean of 1.36 M e (Appendix A). The average mass of a NS in field OGLE-IV-BLG611 is 1.36 M e , 1.75 M e , and 1.36 M e , for the Raithel18, Spera15, and SukhboldN20 IFMRs, respectively.
The most significant differences between the IFMRs are found in the BH mass distribution. The average mass of a BH in field OGLE-IV-BLG611 is 9.32 M e , 14.74 M e , and 9.99 M e , for the Raithel18, Spera15, and SukhboldN20 IFMRs, respectively. The similarity in average BH mass between the Raithel18 and SukhboldN20 IFMRs is indicative of the fact that most stars in the line of sight toward the Bulge are solar or supersolar metallicity, and the Raithel18 and SukhboldN20 IFMRs are similar in this metallicity regime (Figure 2). In addition, the Spera15 IFMR tends to produce much more massive BHs than either the Raithel18 or SukhboldN20 IFMRs, which is why its average BH mass is higher (Figure 2).
In order to validate the results from PopSyCLE, we compare the stellar density and event rates for the OGLE-like simulated survey to the efficiency-corrected results presented in Table 7 of Mróz et al. (2019). To replicate the observing conditions of the OGLE survey, in the PopSyCLE simulation we use a seeing-limited blending radius of 0 65 and make observations in the I-band filter. For the stellar density comparison, our star count is restricted to stars with I < 21 to match Mróz et al. (2019). For the efficiency-corrected event-rate comparisons, we restrict the events in the simulation to have Einstein crossing times 0.5 < t E < 300 days, source magnitude I src < 21, and impact parameter u 0 < 1 in order to match the completenesscorrected sample of Mróz et al. (2019; see also Table 4, column "Mock Mróz19," in Lam et al. 2020). In Table 2, the simulated stellar densities and event rates produced by each IFMR are compared to the observed event rates and stellar densities of Mróz et al. (2019).
Across different IFMRs the simulation event rates are comparable, as stellar lensing events dominate over compact object lensing. Overall the simulated event rates compare well with the completeness-corrected event rates observed by Mróz et al. (2019).
While the event rates from the simulation are in reasonable agreement with the observed rates, the stellar density of each field as reported in Mróz et al. (2019) is typically a factor of 2 higher than PopSyCLE predicts. This is a known issue found in other Galactic models, likely due to uncertainties in the length, angle, and overall structure of the Galactic bar, as well as variable extinction over small scales toward the bulge. The factor of 2 difference in star counts cannot be explained by accounting for stellar binarity or confusion (Abrams et al. in preparation). Additional details are presented in Appendix B.
We next consider whether differences in the IFMR are detectable from the observed (i.e., noncompleteness-corrected) distributions of microlensing events. To do this, we exclude simulated PopSyCLE microlensing events that are faint (baseline magnitude I base > 21), not substantially lensed (impact parameter u 0 > 2), or have a low observed amplification Δm < 0.1 mag. These cuts were based on the events reported by OGLE's Early Warning System (EWS; Udalski et al. 1994) from 2016 to 2018. See also
Table 4, column "Mock EWS," in Lam et al. (2020). Note. Observed versus simulated stellar density and efficiency-corrected event rates for nine select fields in the OGLE-IV survey. To calculate the efficiency-corrected event rates for the Spera15, Raithel18, and SukhboldN20 IFMR simulations, we apply the following cuts to match the completeness-corrected sample in M19: source magnitude I < 21 mag, maximum impact parameter u 0 < 1, and Einstein crossing time range 0.5 < t E < 300 days. The stellar densities for each field are also based on OGLE observability cuts (I < 21 mag), neglecting the effects of crowding. Stellar rather than compact object events dominate the microlensing event rate; thus, changing the IFMR only has a small impact on the overall event rates, generally within the uncertainties.
The mass distribution of BH lenses for the different IFMRs as detectable by a 10 yr OGLE Galactic Bulge survey is shown in Figure 5. The Spera15 IFMR has significantly more highmass lenses compared to either Raithel18 or SukhboldN20 IFMRs.
To explore whether these mass distribution differences are detectable, we investigate their impact on observable quantities measured by photometric microlensing surveys. All microlensing parameters depend on the Einstein radius θ E :
with d L the distance to the lens, d S the distance to the source, and M L the mass of the lens. However, θ E is generally not measurable with a photometric microlensing lightcurve; only quantities normalized by θ E can be measured. This includes the Einstein crossing time, t E , and the microlensing parallax, π E . The Einstein crossing time is
where μ rel is the magnitude of the relative source-lens proper motion. The Einstein crossing time characterizes the length of the photometric microlensing event.
The microlensing parallax is
where π rel is the relative parallax:
Microlensing parallax characterizes changes to the shape of the otherwise symmetric photometric lightcurve due to the Earth's motion around the Sun. It encodes information about the relative distance between the source and the lens. The Einstein crossing time and microlensing parallax scale with the lens mass as and ( ) p µ M 1 . 14 E L As shown by Lam et al. (2020) the most massive lenses (i.e., BHs) are characterized by a long Einstein crossing time and a small microlensing parallax.
The Spera15 IFMR allows a wider range of BH masses, and the more massive BHs produce longer Einstein crossing times and smaller microlensing parallaxes as compared to the SukhboldN20 and Raithel18 IFMRs, as shown in Figure 6. The Einstein crossing time is not as sensitive to changes in the lens mass as the microlensing parallax, since the relative proper motion also affects the Einstein crossing time. For these reasons, the Spera15 BH lens population in π E versus t E space does not change dramatically in t E , but includes events with much lower microlensing parallax as compared to the Raithel18 BH lensing population. Because the SukhboldN20 BH lens population does not have a significantly different mass distribution as compared to the Raithel18 BH lens population, it is difficult to distinguish them based on their resultant distributions of π E and t E . This reflects the fact that there are far fewer low-metallicity massive stars in the Milky Way then there are solar or supersolar metallicity ones. In the solar-and supersolar-metallicity regimes, the SukhboldN20 and Raithel18 IFMRs are similar (Figure 2).
An analysis of the distributions of π E and t E can be used to statistically constrain the Milky Way BH population (Golovich et al. 2022). The difference in the compact object mass distribution between the IFMRs is reflected in the difference between the distributions for Einstein crossing time and microlensing parallax, as shown in Figures 7 and 8. The differences between the Spera15, Raithel18, and SukhboldN20 IFMRs are largest at low microlensing parallax. The effect of increasing lens mass on Einstein crossing time is weaker. For this reason the best way to find events that are likely to be caused by higher-mass BHs is to select events not only with long Einstein crossing times but also with very low microlensing parallax (Sections 5.1 and 5.3.2).
We now consider whether these different IFMRs are statistically distinguishable with 10 yr of observations from an OGLE-like microlensing survey. In this case, we are not interested in the efficiency-corrected number, but rather the observed number of events. Taking all the observable events as outlined in Section 4.3 and then rescaling the number to have the simulated area and duration match those of the OGLE survey would result in an overestimation of the observed number of events, as those observational cuts do not capture sources of detection inefficiency such as observational gaps or sparse lightcurve coverage. We thus empirically rescale the number of events we observe in our simulated subset of the OGLE survey to the expected number of total events as follows. Mróz et al. (2017) published 2617 point-source point lens (PSPL) events in OGLE-IV's nine high-cadence bulge fields from 2010 to 2015. Mróz et al. (2019) published 5790 PSPL events in the remaining 112 OGLE-IV low-cadence bulge fields from 2010 to 2017. This corresponds to roughly 523 events yr -1 and 827 events yr -1 in the low-and high-cadence fields, respectively, for a total of 1350 events yr -1 total. This implies that over a 10 yr survey, an OGLE-like survey should observe around 13,500 events. We thus scale the total number of observable simulated events produced by the Raithel18 IFMR to 13,500. Applying the same scaling factor to the SukhboldN20 and Spera15 IFMR simulations result in 12,800 and 13,300 events, respectively.
As discussed in Section 4.3, we expect the Spera15 IFMR to have an excess of high Einstein crossing time events and low microlensing parallax as compared to the Raithel18 and SukhboldN20 IFMRs, since the Spera15 IFMR produces more high-mass BHs. Figure 9 shows a slight excess of long Einstein crossing time events for the Spera15 IFMR as compared to the Raithel18 IFMR and the SukhboldN20 IFMR. Figure 10 shows an excess of low microlensing parallax events for the Spera15 IFMR as compared to both the Raithel18 IFMR and the SukhboldN20 IFMR.
Next, we focus on whether the Spera15 and SukhboldN20 IFMRs can be distinguished via t E and π E . Integrating over the total Einstein crossing time distributions produced by each IFMR for events with t E > 100 days, we find a insignificant (2σ) difference in the number of events with Einstein crossing times greater than 100 days. Integrating over the total microlensing parallax distributions to find the number of events with π E < 0.02 produced by each IFMR, we find a significant ∼5σ difference. This indicates that low microlensing parallax is a more sensitive indicator of high lens mass (i.e., a BH lens) than a long Einstein crossing time.
Unfortunately, these small microlensing parallaxes are not possible to measure with ground-based photometry. Even for the long t E ≈ 100-300 days events where microlensing parallaxes are easier to constrain, OGLE is not very sensitive to π E 0.03 (Wyrzykowski et al. 2016; Figure 3). For shorterduration events the precision with which microlensing parallax can measured is even worse. This means that most BH lensing events, with their very small microlensing parallaxes, will actually have undetectable parallaxes and are often not reported Figure 5. Mass distribution of black hole lenses from OGLE-detectable microlensing events. The number of events is scaled to the number of events an OGLE-like survey would observe over a span of 10 yr, as described in Section 4.4. The Spera15 IFMR produces much more massive BHs than the Raithel18 and SukhboldN20 IFMRs, and this is reflected in the detected lens masses. The SukhboldN20 IFMR produces the fewest BH lensing events overall. (Karolinski & Zhu 2020). Thus, with current survey capabilities, different IFMRs cannot be distinguished reliably using photometric microlensing.
Comparing the overall Einstein crossing time distribution reported in Mróz et al. (2019) to the distributions produced by the PopSyCLE simulations (see Figure 9), we note that there is a slight discrepancy in the peak of the distribution, with our models predicting a peak at around 20 days while the Mróz et al. (2019) distribution peaks at around 25 days. This difference is not unexpected as our simulation does not exactly reproduce the OGLE field of view, and also does not yet include the effect of binary lenses and sources, which likely shift the peak of the Einstein crossing time distributions to longer timescales (N. S. Abrams et al. 2023, in preparation). Overall, a 5 day shift of the Einstein crossing time distribution will not significantly affect the fraction of longduration events (t E > 100 days) used to constrain the IFMR. For a more extensive discussion of the comparisons of PopSyCLE simulations to OGLE observations, please see Lam et al. (2020).
The Nancy Grace Roman Space Telescope is an upcoming NASA flagship mission with a 2.4 m telescope planned for launch in 2026-2027. As one of its Core Community Surveys, Roman will conduct a microlensing survey toward the Galactic bulge in the infrared to find thousands of cold exoplanets (Penny et al. 2019). In addition to photometric microlensing, Roman will also simultaneously be able perform astrometric microlensing measurements, allowing lens masses to be directly constrained. Lam et al. (2020) estimate that Roman's microlensing survey will be able to detect hundreds of BHs.
We now consider whether Roman can distinguish different IFMRs solely with photometric measurements. As an infrared space mission, Roman has advantages over an optical groundbased survey like OGLE. Infrared observations allow for a larger number of sources to be detected as most sources will be lowmass stars, which emit at longer wavelengths. Observations in the H band will also not suffer as much from interstellar extinction, which is especially important toward the Galactic bulge. Roman will also be able to see sources with baseline magnitudes less than 24th mag in the H band and as deep as 26th mag if the data are stacked (Penny et al. 2019) as compared to OGLE, which is limited to sources with baseline magnitudes brighter then 22 mag in the I band (Udalski et al. 1994).
Of the simulation fields in Table 2, field OGLE-IV-BLG500 is roughly located near the center of the proposed Roman exoplanet microlensing survey (Penny et al. 2019). Fields OGLE-IV-BLG506 and OGLE-IV-BLG675 are also near the proposed Roman fields. We thus use these three fields to represent the Roman sample used for this analysis. To replicate the observing capabilities of Roman, in the PopSyCLE simulation we use a diffraction-limited blending radius of 0 09 and make observations in the H-band filter. We define events detectable by Roman if they have baseline magnitude I base < 24, impact parameter u 0 < 2, and observed amplification Δm > 0.05 mag (see also Table 4, column "Mock WFIRST," in Lam et al. 2020). 9 The number of microlensing events detected by the PopSyCLE simulation is then scaled to equal a 1.97 deg 2 survey area and 5 yr duration to roughly match the number of events expected to be observed by Roman's microlensing survey. However, this method of determining Roman observable events neglects the effects of observational gaps, which lower the detection efficiency.
The fiducial Roman microlensing survey as presented in Penny et al. (2019) consists of six 72 day Galactic bulge observing seasons. The seasons are centered on the vernal and autumnal equinoxes, with an observational gap of about 3.5 months between available observing windows; this is because of spacecraft pointing limitations due to the observatory design. A complete analysis simulating Roman's detection efficiency is beyond the scope of this paper. Instead, we make a simple estimate of the detection efficiency. Based on 6 × 72 days of observations over a 5 yr window, this corresponds to a duty cycle of 24%. We thus consider 24% of the events detected by the criterion of the previous paragraph to be a sufficiently reliable estimate of the true number of microlensing events detectable by Roman. This ultimately corresponds to 27,300, 27,300, and 27,500 microlensing events detected by Roman's 5 yr microlensing survey for the Raithel18, SukhboldN20, and Spera15 IFMRs, respectively. For comparison, Table 2 of Penny et al. (2019) estimates 27,000 and 54,000 microlensing events with |u 0 | < 1 and |u 0 | < 3, respectively, to be detectable by Roman. Their simulation, detection criteria, and method for estimating the observable number of events are quite different from ours; this demonstrates that our simple estimates are reasonable. Thus we proceed with the assumption of 27,300-27,500 Roman-Figure 7. A histogram of Einstein crossing times for OGLE-detectable microlensing events as a function of lens type for each IFMR. The number of events is scaled to the number of events an OGLE-like survey would observe over a span of 10 yr, as described in Section 4.4. The largest differences between the IFMRs are found at long Einstein crossing times (150 days) as that is where the overall distribution is dominated by BH lenses. detectable microlensing events in the following results. Figures 11 and 12 show the expected distribution of Einstein crossing times and microlensing parallaxes detected by Roman assuming a Spera15 IFMR, SukhboldN20 IFMR, or Raithel18 IFMR.
Next, we focus on whether the Spera15 and SukhboldN20 IFMRs can be distinguished via t E and π E . Integrating over the total Einstein crossing time distributions produced by each IFMR for events with t E > 100 days we find a significant (∼5σ) difference in the number of events with Einstein crossing times greater than 100 days. Integrating over the total microlensing parallax distributions to find the total number of events with π E < 0.02 produced by each IFMR, we also find a significant ∼5σ difference.
Although it will have exquisite photometric precision, it is uncertain if Roman will be able to measure such small microlensing parallaxes accurately due to the observational gaps. In future work, we will explore how additional observations during Roman's bulge windows during its nonmicrolensing survey seasons, as well as observations from ground-based facilities like the Rubin Observatory, United Kingdom Infra-Red Telescope, or PRime-focus Infrared Mirolensing Experiment, could improve constraints on the microlensing parallax.
We next consider the effect of Maxwellian versus singlevalued birth kicks for NSs/BHs on the Einstein crossing time distribution. For one of the Roman fields, OGLE-IV-BLG500, we run a PopSyCLE simulation identical to the Raithel18 IFMR simulation described in Section 4.5, except we use single-valued birth kicks. The single-valued birth kick velocities are equal to the average velocity of the Maxwellian distribution, 350 km s -1 for NS and 100 km s -1 for BHs. A comparison of the Einstein crossing time distributions between the two simulations is shown in Figure 13.
We find that the effect of Maxwellian as opposed to single-value birth kick velocities has a fairly weak impact on the distribution of Einstein crossing times. This is because the kick velocities are added on top of existing stellar velocities to calculate the final remnant velocities (Section 2.4). Because there is already a significant amount of dispersion in the stellar velocity distribution, the effect of the additional dispersion from the Maxwellian kick distribution is diluted. The t E distribution for NSs is slightly wider, with a longer tail toward short t E , when using a Maxwellian distribution as compared to a single-kick velocity. However, within the uncertainties, the difference is not significant. For BHs and the full lens population, the t E distributions are identical within the Figure 8. A histogram of microlensing parallax for OGLE-detectable microlensing events as a function of lens type for each IFMR. The number of events is scaled to the number of events an OGLE-like survey would observe over a span of 10 yr, as described in Section 4.4. The largest differences between the IFMRs are found at small microlensing parallaxes (0.02) as that is where the overall distribution is dominated by BH lenses.
uncertainties. This means that even with 20,000 photometric microlensing events, different kick velocity distributions the same mean will likely be indistinguishable using Einstein crossing time distributions alone.
From Figure 6 it is apparent that, regardless of IFMR, BH lenses have the longest Einstein crossing times and the lowest microlensing parallaxes.
Astrometric follow-up can be used to measure the lens masses of candidate BHs (Lu et al. 2016;Lam et al. 2022;Sahu et al. 2022). As more astrometric follow-up campaigns are done it will be useful to quantify what percentage of events we should expect to be BH lenses based on our candidate selection criteria (either a minimum t E or a maximum π E ). Figures 14 and 15 show the ratio of BH events to total events for a 10 yr OGLE-like survey as a function of t E and π E , respectively.
Different IFMRs produce different BH fractions depending on the exact selection criteria. For example, a dearth of BH lenses in a sample of candidates chosen for astrometric followup based on t E > 100 days might indicate that the Milky Way IFMR is closer to SukhboldN20, and that the BH to NS formation ratio is lower then predicted by Spera15.
Using metallicity distributions and a star formation history from Galaxia, we create a present-day BH mass function (PDBHMF) for the Milky Way for each of the IFMRs. We find the total number of BHs produced by the Raithel18, Spera15, and SukhboldN20 IFMRs are 2.18 × 10 8 , 2.67 × 10 8 , and 1.72 × 10 8 respectively. Current estimates of the total number of stellar-mass BHs in the Milky Way range between 10 8 and 10 9 (Agol & Kamionkowski 2002). Fraction of BH events compared to total events in the OGLE-like survey sample as a function of t E . The number of events is scaled to the number of events an OGLE-like survey would observe over a span of 10 yr, as described in Section 4.4. The SukhboldN20 IFMR produces the fewest BHs and has the lowest fraction of BH events at long t E (>100 days). The Spera15 IFMR produces the most BHs and has the highest fraction of BH events at long t E . Regardless of IFMR, roughly at least a third of events with t E > 150 days are BH lenses.
Figure 15. Fraction of BH events compared to total events in the OGLE-like survey sample as a function of π E . The number of events is scaled to the number of events an OGLE-like survey would observe over a span of 10 yr, as described in Section 4.4. Because the Spera15 IFMR produces more massive BHs than Raithel18 and SukhboldN20, it has the highest fraction of BH events at small π E (<0.02). Although the SukhboldN20 IFMR produces more massive BHs than the Raithel18 IFMR, it also produces fewer BHs overall, so the fraction of BH events at small π E for those two IFMRs are similar.
Figure 16 compares the Milky Way PDBHMF predicted by each IFMR against the fiducial "Power Law+Peak" binary black hole (BBH) primary mass distribution derived using the Gravitational Wave Transient Catalog 3 (GWTC-3; The LIGO-Virgo-KAGRA Collaboration et al. 2021). For comparison purposes, the "Power Law+Peak" model has been vertically scaled so that its maximum roughly matches the maxima of Galactic BH distributions produced by the IFMRs. The GWTC-3 BBH population is not a direct analog to the simulated isolated Milky Way BH population, as the former is an extragalactic binary population at low metallicity. However, a comparison of the two populations' similarities and differences is still worthwhile.
Note that the four large spikes in the Spera15 IFMR PDBHMF are artifacts of the coarse metallicity bins used in this simulation (see Section 2.3 for more details). Finer metallicity bins would smooth over these spikes at the expense of increased computational time and the need to store many more SPISEA isochrones.
None of the PDBHMFs predicted by the different IFMRs match the fiducial BBH Power Law+Peak distribution. All three IFMRs have approximately flat PDBHMFs for masses 20 M e , while the BBH primary mass distribution follows a power law with spectral index α ∼ -3.5 (The LIGO Scientific Collaboration et al. 2021). None of the IFMRs reproduce the peak at 35 M e in the BBH mass distribution. The slope of the Spera15 IFMR PDBHMF is in reasonable agreement with the Power Law+Peak model for masses 35 M e . This is likely because the high-mass end of the Spera15 IFMR PDBHMF is composed of high-mass BHs that came from a low-metallicity progenitor population, similar to the progenitors of the GWTC-3 BBH population. In contrast, the low-mass end of the PDBHMF does not match. The SukhboldN20 IFMR PDBHMF has a similar slope to the power-law portion of the Power Law +Peak model for masses 15 M e , modulo the peak at ∼35 M e . Similar to the explanation for the Spera15 IFMR, the progenitors of the BHs in this mass range have low-metallicity progenitors The Raithel18 IFMR PDBHMF does not have the massive BHs in the BBH population as all progenitors are assumed to be solar metallicity.
Both the Raithel18 IFMR and the Power Law+Peak model favor a minimum BH mass above 5 M e , while both Spera15 and SukhboldN20 IFMRs allow for BH masses down to 3 M e . There is observational evidence for mass-gap BHs between 2 and 5 M e including two gravitational-wave merger remnants, one at ∼3 M e (Abbott et al. 2017) and another at ∼3.4 M e (Abbott et al. 2020), and the ∼2.6 M e merger component of GW190814 (Abbott et al. 2020). There has also been a detection of a ∼3 M e BH in a noninteracting binary system (Thompson et al. 2019). These results indicate that the IFMR should allow for at least a few BHs below ∼5 M e .
In this work we updated the way birth kicks are applied to NSs and BHs in PopSyCLE. Previously, each population was assigned a single-kick velocity value. Now, the kick values for each population are drawn from a Maxwellian distribution. A more accurate implementation would also have the birth kick velocities be a function of mass. In future versions of PopSyCLE we will consider different kick formalisms, such as having kick velocity decreasing with increasing remnant mass, or implementing predictions for birth kick velocities coming from gravitational-wave observations (O'Shaughnessy et al. 2017). If more massive BHs were given smaller birth kicks, this would likely make the Einstein crossing times slightly longer on average for massive BHs, and slightly shorter for less massive BHs.
The best BH candidates obtained from photometric microlensing surveys have both long Einstein crossing times and low microlensing parallaxes. Up until this point BH candidates for astrometric follow-up have mainly been selected based on t E alone (Lu et al. 2016;Lam et al. 2022;Sahu et al. 2022). However, even for t E > 120 days, around half of lenses are not BHs (Figure 14). Selecting on microlensing parallax in addition to Einstein crossing time would improve our ability to pick BH candidates for astrometric follow-up. However, it is unclear how well π E can be predicted prior to the photometric peak.
It may be possible to use Fisher matrix analysis to constrain π E over the course of the microlensing event, although it will be difficult in the case of BH lenses which have small π E (Karolinski & Zhu 2020). One of the next steps motivated by this work will be to determine exactly how well π E can be predicted prior to peak photometric amplification, in order to determine the viability of π E for use as a BH candidate selection criteria.
In this work we have assumed that all stars and compact objects are single. However, a significant fraction of stars are known to be in binary systems, which can also act as sources or lenses. Large microlensing catalogs such as Mróz et al. (2019) select only for microlensing events with a single lens and source. In some cases it is very clear that a microlensing signal is coming from a binary lens or source due to caustic crossings, significant asymmetry, or multiple peaks in the lightcurve. However, some lensing geometries may result in lightcurves that resemble single lens and source events. In other cases, observational gaps or data quality are not sufficient to rule out the possibility of binary microlensing. Abrams et al. (in preparation) explores the effect of adding binaries to the PopSyCLE microlensing simulation.
In addition to changes in the microlensing source and lens population, binarity can affect the IFMR itself. All the IFMRs considered in this paper come from models of single-star evolution, and so are only valid for single stars or binaries sufficiently wide enough to be noninteracting. Close binaries can allow for mass transfer between a massive star and its companion during the late stages of stellar evolution (Ivanova et al. 2013). Some binary star IFMRs simply assume hydrogen envelope loss for close binaries (Woosley et al. 2020), while others actually calculate the amount of mass and angular momentum transfer (Yoon et al. 2010). Population synthesis codes such as StarTrack (Belczynski et al. 2002(Belczynski et al. , 2008) that include binary stellar evolution have been used to investigate binary BH mergers (Belczynski et al. 2016). However, these have not been incorporated into microlensing simulations. A future version of PopSyCLE will incorporate the results of binary stellar evolution.
In this work we explored the possibility of using photometric microlensing to constrain the IFMR for massive stars.
1. We have added two metallicity-dependent IFMRs, Spera15 and SukhboldN20, to the SPISEA simple stellar population synthesis code. 2. Different IFMRs yield different Galactic BH mass distributions. This in turn affects the observed distributions of microlensing event parameters. We focus on long t E and small π E in this work as these are the regimes most sensitive to BH lenses. 3. Considering the number of events with t E > 100 days observed with an OGLE-like survey, the difference between the Spera15 and SukhboldN20 IFMRs is not statistically significant. Considering the number of events with π E < 0.02 observed with an OGLE-like survey, the difference between the IFMRs is statistically significant. However, ground-based surveys like OGLE are not able to constrain such small π E , making it not possible to distinguish the IFMRs photometrically. 4. Considering the number of events with t E > 100 days observed with a Roman-like survey, the difference between the Spera15 and SukhboldN20 IFMRs is statistically significant. Considering the number of events with π E < 0.02 observed with an OGLE-like survey, the difference between IFMRs is also statistically significant. Roman's excellent photometric precision should allow it to constrain small π E , but this is tempered by the large observational gaps in the survey. Further work to determine how best these gaps can be filled, either by Roman itself or by other facilities, will maximize the ability of Roman to detect small π E , and in turn constrain the IFMR and detect BHs. 5. Microlensing parallax is more sensitive to changes in lens mass than the Einstein crossing time. Because BHs are high mass, this means that it will be important that future surveys be able to accurately measure small microlensing parallax signals (π E < 0.02) in order to place meaningful constraints on the BH lens mass distribution. 6. Comparing astrometric BH lens detection efficiency based on candidate selection criteria to our predictions for each IFMR will be another way to place constraints on the Milky Way IFMR. Note. Observed versus simulated stellar density and efficiency-corrected event rates for three select fields in the OGLE-IV survey representative of the Roman microlensing fields. The observed stellar density and event rates are calculated as described in Table 2. Here, we consider two different Galactic bar models. In the v2 Galactic model, α = 11.1°and x 0 = 1.59 kpc. In the v3 Galactic model, α = 28°and x 0 = 0.7 kpc. Both the v2 and v3 models here use the SukhboldN20 IFMR. We only consider a single IFMR, as the difference in event rate across different IFMRs (Table 2) is smaller than the difference of changing the Galactic model.
SPISEA was formerly called PyPopStar, and is referred to as such inLam et al. (2020).
Available in SPISEA v2.1 and later.
http://web.pd.astro.it/mapelli/group.html
The Astrophysical Journal, 941:116 (18pp), 2022 December 20 Rose et al.
Specifically, we ran the simulation with a sampling cadence of 10 days. Note that in PopSyCLE, the sampling cadence is not equivalent to a real survey's observational cadence. This is why a sampling cadence of 10 days can detect events with t E < 10 days. See Section 4.3 and Figure3of Lam et al. (2020) for further details.
Note thatLam et al. (2020) allowed I base < 26; we chose a brighter limit of I base < 24 as a more conservative estimate that does not rely on stacking multiple observations together. On the other hand, we allow for smalleramplitude microlensing events by requiring Δm > 0.05 mag, whileLam et al. (2020) required Δm > 0.1 mag; we assume Roman's photometric precision and subhour sampling rate will trim out low-amplitude variables better than OGLE.