Disentangling the Black Hole Mass Spectrum with Photometric Microlensing Surveys

From the formation mechanisms of stars and compact objects to nuclear physics, modern astronomy frequently leverages surveys to understand populations of objects to answer fundamental questions. The population of dark and isolated compact objects in the Galaxy contains critical information related to many of these topics, but is only practically accessible via gravitational microlensing. However, photometric microlensing observables are degenerate for different types of lenses, and one can seldom classify an event as involving either a compact object or stellar lens on its own. To address this difficulty, we apply a Bayesian framework that treats lens type probabilistically and jointly with a lens population model. This method allows lens population characteristics to be inferred despite intrinsic uncertainty in the lens class of any single event. We investigate this method’s effectiveness on a simulated ground-based photometric survey in the context of characterizing a hypothetical population of primordial black holes (PBHs) with an average mass of 30M ⊙. On simulated data, our method outperforms current black hole (BH) lens identification pipelines and characterizes different subpopulations of lenses while jointly constraining the PBH contribution to dark matter to ≈25%. Key to robust inference, our method can marginalize over population model uncertainty. We find the lower mass cutoff for stellar origin BHs, a key observable in understanding the BH mass gap, particularly difficult to infer in our simulations. This work lays the foundation for cutting-edge PBH abundance constraints to be extracted from current photometric microlensing surveys.

Unfortunately, BHs in isolation do not radiate measurable amounts of light, gravitational radiation or particles making them difficult to detect.Detectable emission from a BH is only produced through interaction with its environment.Massive, extra-galactic BHs can be detected through a strong gravitational interaction with another object causing gravitational radiation (e.g., Abbott et al. 2016Abbott et al. , 2021a) ) or through accretion, causing electromagnetic (EM) radiation (e.g., Akiyama et al. 2019;Event Horizon Telescope Collaboration et al. 2022;Fabbiano 2006).Studies using gravitational wave (GW) emission and EM observation can be effective for un-derstanding the extra-galactic BH population (e.g., Abbott et al. 2023Abbott et al. , 2021b,c;,c;Edelman et al. 2022;Roulet & Zaldarriaga 2019), provided that detection bias from observational selection effects is mitigated (Liotine et al. 2023).
Within the Milky Way, there are estimated to be ≈ 10 8 stellar origin BHs (SOBHs; e.g., Samland 1998).Despite this large expected abundance, only ∼50 SOBHs have been detected.The bulk of these BHs are found in X-ray binaries (e.g., Remillard & McClintock 2006;Corral-Santana et al. 2016), despite these systems being an intrinsically rare outcome of binary evolution (e.g., Kalogera 2001;El-Badry et al. 2023a).These systems are detectable due to bright X-ray emission from accretion of a luminous stellar companion onto the BH.Most recently, leveraging high-precision astrometry from Gaia (Gaia Collaboration et al. 2016, 2021), two nearby BHs that perturb the motion of their luminous binary companion have also been detected (El-Badry et al. 2023b,a;Chakrabarti et al. 2022).
Despite this diverse set of observational channels, they all require the BH, regardless of SOBH or PBH origin, to have a companion.None of these techniques are sensitive to detecting the population of isolated BHs within the the Galaxy.Gravitational microlensing is uniquely positioned to fill this detection blind spot, as it is the only practical method with which isolated BHs can be detected and characterized (Gould 2000;Bennett et al. 2002;Lam et al. 2022;Sahu et al. 2022;Chapline & Frampton 2016).Detecting a BH via microlensing only requires its close alignment with a distance background star.In addition to understanding BHs from single microlensing events, the characteristics of sets of microlensing events observed over the course of a survey can encode information about the underlying BH lens population (e.g., Lam et al. 2020;Rose et al. 2022;Mroz et al. 2021;Wyrzykowski & Mandel 2020;Wyrzykowski et al. 2016;Mróz & Wyrzykowski 2021).
However, robustly characterising the underlying lens population from microlensing surveys is challenging.This is because the photometric microlensing signal is degenerate in lens mass, distance, kinematics (Paczynski 1996), typical transient survey noise systematics (Golovich et al. 2022), and contains no direct lens mass or lens identity information.A microlensing event can also have an astrometric signal (Eddington 1919;Walker 1995;Hog et al. 1995;Miyamoto & Yoshii 1995;Rybicki et al. 2018), which can break the photometric degeneracies resulting in a direct measurement of the lens mass (e.g., Lu et al. 2016;Sahu et al. 2017;Kains et al. 2017;Zurlo et al. 2018;Lam et al. 2022;Sahu et al. 2022;McGill et al. 2023), and also differentiate lens subpopu-lations (e.g., Belokurov & Evans 2002;Lam et al. 2020;Pruett et al. 2022).However, currently there is no large sample of microlensing events with measured astrometry, which would be required to perform population inference, although this is set to change over the coming years (e.g., Gaia Collaboration et al. 2016;Spergel et al. 2015;Lam et al. 2020;Sajadian & Sahu 2023;Lam et al. 2023).
In the absence of astrometry, ∼10 4 photometric microlensing events have been detected over the past decades (e.g., Udalski et al. 2015;Kim et al. 2016;Jeong et al. 2015;Husseiniova et al. 2021) which can be used to constrain the underlying lens populations.Despite the degeneracies in the photometric signal, progress has been made in understanding the lens populations in the tails of the mass distribution i.e., Free Floating Planets (Mróz et al. 2017;Sumi et al. 2023) and SOBHs (e.g., Mroz et al. 2021) which effect the tails of the photometric microlensing event timescale distribution.However, current methods require manually pre-selecting or classifying events based on event characteristics, for example, assuming a set of candidate events with the longest timescales (e.g., Lu et al. 2016) and large parallax signals are caused by BHs (e.g., Wyrzykowski et al. 2016;Kaczmarek et al. 2022).
In the case of candidate BH lenses, auxiliary information can sometimes be used to constrain the identity of the lens.This information includes, baseline source astrometry (e.g., Wyrzykowski & Mandel 2020;Kaczmarek et al. 2022), testing if the event is consistent with the lens being dark, and assuming some model of the Galaxy which pins down the relative lens-source distance and kinematics (e.g., Wyrzykowski et al. 2016;Kaczmarek et al. 2022).However, conclusions about the lens identity are sensitive to unreliable source astrometry and distances, and assumptions about the location of the lens and source imposed by a given Galactic model (Mróz & Wyrzykowski 2021).Overall, definitively classifying the lens for a single microlensing event is difficult and can bias resulting inferences about the underlying lens population.
In this work we overcome the lens classification problem by extending the inference framework of Zevin et al. (2021) and Franciolini et al. (2022) and applying it to microlensing by treating the lens classification probabilistically.This method allows for all events to have some probability of belonging to each class (e.g, SOBH, PBH, or Star), effectively marginalizing over each possibility and bypassing the need to assume a single lens class.This approach allows the underlying lens population to be modelled jointly and in the absence of confident, individual event classifications.Our method gen-eralizes the work of Sumi et al. (2023) and Mroz et al. (2021) to comprehensively include survey selection effects, rate information and probabilistic lens classification jointly with an uncertain lens population model.
We apply this new framework with the goal of constraining and disentangling the mass spectra of the underlying lens population given a survey of photometric microlensing events.In this context, we focus on investigating if current photometric microlensing data can place constraints on the mass spectrum and abundance of PBHs in the Galaxy.Using simulated microlensing survey data, we evaluate the effectiveness of our method, including its ability to classify single lenses and marginalize over those classes to place constraints on the underlying lens populations.Through these exercises, we demonstrate the power of these methods on disentangling and constraining the mass spectra of PBHs and SOBHs and lay the foundation for these methods to be used in combination with Galactic simulations (e.g., Lam et al. 2020) to provide cutting edge constraints on the population of BHs in the Milky Way on current microlensing surveys like the Optical Gravitational Lensing Experiment (OGLE; Udalski et al. 2015).
This work also complements methods applied decades ago to study MAssive Compact Halo Objects (MA-CHOs) by the MACHO (e.g., Allsman et al. 2001), Expérience pour la Recherche d'Objets Sombres (EROS; Tisserand et al. 2007;Blaineau et al. 2022) and OGLE (e.g., Wyrzykowski et al. 2009Wyrzykowski et al. , 2010Wyrzykowski et al. , 2011a,b) ,b) collaborations.These projects all used photometric microlensing observations of the Magellanic Clouds to estimate the abundance and halo mass fraction due to MACHOs (of which, PBHs could be a specific realization).These experiments probed the galactic halo, conservatively attributing all microlensing detections to MACHOs and using optical depth calculations.The methods proposed here are designed for use with observations of the galactic bulge, which observe thousands of events and must be understood in the context of many lensing subpopulations, requiring different statistical methods.
We begin by describing modeling of single photometric microlensing lightcurves in Sec. 2. With the microlensing basics outlined, we describe our method in Sec. 3 including: accounting for observation bias (Sec.3.1), single event classification (Sec.3.2), and fully hierarchical inference (Sec.3.3).In Sec. 4 we describe our verification testbed which includes a set of simulated population models (Sec.4.1), a simulated microlensing survey (Sec.4.2), and performing inference at the single event (Sec.4.3) and population level (Sec. 4.4).With our model and simulation framework laid out, we then apply our method to a suite of simulated datasets in Secs. 5 and 6, focusing on single event and populationlevel inferences, respectively.We summarize our findings in Sec. 7.

PHOTOMETRIC MICROLENSING
Consider a point lens with mass M , and a more distant point source, at distances D L and D S from an observer, respectively.In the case of perfect lens-source alignment, the gravitational field of the lens deflects the light of the background source forming an Einstein ring of angular radius (Chwolson 1924;Einstein 1936), For imperfect lens-source alignments, two, usually unresolved, images of the source are formed (Liebes 1964;Refsdal 1964).As the lens passes between the source and observer the source images change brightness giving rise to an apparent amplification of the background source flux (Paczynski 1986), Here, u(t) is the magnitude of the lens-source angular separation vector in units of θ E .The relative lens-source trajectory, u(t), can be parameterized by (Gould 2004), Here, u 0 = |u 0 | is the magnitude of lens-source impact parameter in units of θ E , t 0 is the time of lens-source closest approach, t E = θ E /µ rel is the Einstein crossing time where µ rel is the relative lens-source proper motion vector.The first two terms of Eq. (3) make up the standard Paczynski (1986) rectilinear trajectory model.The third term in Eq. ( 3) accounts for the annual microlensing parallax signal which is caused by the acceleration of an Earth-based observer (Alcock et al. 1995).π E is the vector microlensing parallax which can be described by its magnitude is the relative lens-source parallax and ϕ is the angle between the ecliptic north and the direction of the lens-source relative proper motion in the heliocentric frame.
The expression for P (t; π E ) depends on the on-sky microlensing event location, and in all work that follows we use the results in Section 3.1 of Golovich et al. (2022) which are based on Gould (2013).Annual microlensing parallax can lead to typically subtle (e.g., Alcock et al. 1995;Golovich et al. 2022;Kaczmarek et al. 2022) but sometimes extreme (e.g., Wyrzykowski et al. 2016;Kruszyńska et al. 2022) asymmetrical deviations from the standard Paczynski (1986) lightcurve.
The flux of the unresolved blended images during the microlensing event can be written as, (4) F Base is the total baseline flux including both the unlensed source flux and all unresolved blended light.b sff is the fraction of unlensed source flux to the total base flux.Overall, we can describe the lightcurve with 7 parameters, θ = {F Base , b sff , t 0 , t E , u 0 , π E , ϕ}.
Examination of Eqs.(1-4) shows that the only parameters which can be inferred that contain any information on the lens mass, and therefore its identity, are t E and π E .However, both t E and π E are in units of θ E which cannot be inferred from the lightcurve in this simple scenario.Overall, this means that there is no direct lens mass information contained in the photometric signal for a single event -it is degenerate with the relative lens-source distance and velocity.
Prospects for understanding the nature and identity of lenses via photometric microlensing lensing does improve when a large sample of events can be detected over the course of a survey (e.g., Udalski et al. 2015;Kim et al. 2016;Husseiniova et al. 2021).In this case, different lens types (e.g., Stars, White Dwarfs, Neutron Stars, SOBHs, Free Floating Planets or PBHs) are expected to have differing population characteristics such as different mass distributions, kinematics, and spatial configurations in the Galaxy.These population level differences project down to populations of different lenses producing microlensing events with different characteristics (e.g., Mroz et al. 2021;Sumi et al. 2023).
Figure 1 shows a simulation of microlensing events in t E − π E space from the Population Synthesis for Compact object Lensing Events code (PopSyCLE; Lam et al. 2020) assuming an OGLE-IV-like microlensing survey.PopSyCLE combines galactic and evolutionary models (Kalirai et al. 2008;Sharma et al. 2011;Sukhbold et al. 2016;Raithel et al. 2018;Hosek et al. 2019) with microlensing survey characteristics to simulate detectable microlensing event catalogs for different lens populations.Figure 1 shows that the different lens types do occupy different, albeit overlapping, areas of t E − π E space.This separation is fundamentally caused by the scaling of these parameters with respect to the lens mass: These relationships result in the negative correlation between these two parameters with respect to a changing lens mass, all else being equal.In principle this suggests that given a survey of photometric microlensing events where we can measure  Lam et al. (2020).A strong correlation can be seen between the tE and πE parameters and mass (class) of the lens.This correlation will be be key to unraveling the subpopulation makeup of the total population of lensing objects in the galaxy.
t E and π E it is possible to make inferences about the different underlying subpopulations of lenses.

HIERARCHICAL INFERENCE WITH DETECTION BIAS
To robustly characterize subpopulations of lenses we must account for bias and uncertainty as rigorously as possible -from uncertainty in a single microlensing events' characteristics, to the uncertainty in the identity of the lens for a given event, to having an unknown lens population model.
We start with the concept of event detection probability in Sec.3.1 and its definition for a single event in the context of a population model.We then move on to assigning probabilistic lens classifications to single events in Sec.3.2.Finally, we put everything together in the context of a fully hierarchical population analysis in Sec.3.3.The rest of Sec. 3 first follows standard results in the literature (e.g., Loredo 2004;Vitale et al. 2020;Mandel et al. 2019;Taylor & Gerosa 2018) that are then extended to our specific class of models to improve computational tractability.
In what follows, θ are the parameters describing a single microlensing lightcurve, defined in Sec. 2. d is the lightcurve data for a single microlensing event -which is a collections of times, fluxes, and flux errors.Generally, {} denotes sets, for example, {d} and {θ} correspond to some set of lightcurves and events parame-ters, respectively.When considering population models, we will parametrize the full model as Λ, representing all parameters relevant to population modeling.When considering different subpopulations of lenses, we will denote each class by class a , where a ∈ [0, N pop ) (e.g., stellar, SOBH or PBH lenses).N is the total number of predicted microlensing events (detected or undetected) by the model.The parameters controlling the subpopulation distributions are {λ a }, and the parameters controlling the relative abundance of each subpopulation, {ψ a }, with a ∈ [0, N pop ) for N pop subpopulations.Examples of {λ a } are given in Sec. 4, including parameters of the mass spectrum of lenses.However, these parameters can represent any feature of a subpopulation of the population model, not just the mass spectrum.As {ψ a } are the relative abundances, a ψ a = 1.In summary, Λ = N ∪ {λ a } ∪ {ψ a }.

Detection Probabilities
Detection bias means that our observed set of events are not a fair sample of the true, underlying distribution.This is because some microlensing events are easier to observe than others.This effect can be accounted for by defining a "trigger", tr, and its probability.This means that once event data d has been recorded at the detector, it either produces a trigger, signifying it is a microlensing event, or not.This trigger is evaluated according to whether some deterministic criteria (ρ(d) > ρ threshold ) is met that is typically related to a signal-to-noise ratio (SNR) calculation; (5) Here, and in the work that follows we have assumed that the detection criterion is model-independent i.e., ρ only depends on d.
In the context of modeling a population of microlensing events, we can build on our definition of p(tr|d) to quantify event detection probabilities.First, we can calculate the probability of a trigger given an event exists with a certain set of parameters θ.As the concept of a trigger is inherently tied to detector noise and detector limitations, we will need to introduce a set of data d as a parameter and marginalize over all possible noise realizations consistent with the noise model p(d|θ), giving p(tr|θ) = p(tr, d|θ) dd , = p(tr|d)p(d|θ) dd . (6) This detection probability can be computed using Monte Carlo integration, averaging trigger probabilities over all possible data sets consistent with the single event model likelihood, p(d|θ), conditioned on the event properties θ.
Using Eq. ( 6), the detection efficiency for a population model given a set population parameters, commonly defined as α in the literature (e.g., Mandel et al. 2019;Vitale et al. 2020), is given by, p(tr|Λ) ≡ α = dd dθp(tr, θ, d|Λ) , = dd dθp(tr|d)p(d|θ)p(θ|Λ) .(7) This quantity can also be computed via Monte Carlo integration by simulating values of θ drawn from the population model and subsequently drawing a noise realization from the event likelihood, which is described in App.C. Eq. ( 7) shows that once the probability of a set of data is conditioned on θ, the probability of d is independent of the population model.Conceptually, α can be understood as the efficiency of a population model to produce detectable events.This can also be understood by noting that α can be related to the number of expected detectable events (N det ) and the total number of expected events (N ) by the relation α = N det /N (Loredo 2004;Vitale et al. 2020;Mandel et al. 2019;Taylor & Gerosa 2018), i.e., that α is the fraction of detectable events over the total number of events.

Classification of a single event
With the definitions of detection probabilities in hand, we would like to know, given some value of the population parameters Λ, what is the probability that a single event belongs to a subcategory of the population?This leads us to the following posterior probability of an event belonging to a certain class This posterior probability is directly related to Bayesian model selection methods.If lens-classification is treated as a model selection problem, the ratio of these posterior probabilities for different classes is the posterior odds.
Taking the ratio when neglecting the prior probability of each class yields the Bayes' factor.All of these quantities (the normalized posterior probability, the posterior odds, and the Bayes' factor) are useful metrics at understanding the classification problem, but our focus will be on the normalized posterior probability.
The prior, p(class a |Λ), is the probability that a lens belongs to a given class without considering the data.The likelihood, p(d, tr|class a , Λ), can be simplified.We first re-write the likelihood using the product rule p(d, tr|class a , Λ) = p(tr|d, class a , Λ)p(d|class a , Λ) , = p(d|class a , Λ) , where we see that the selection effects completely drop out of the equation, i.e., p(d, tr|class a , Λ) = p(d|class a , Λ).This is because p(tr|d, class a , Λ) = p(tr|d) = 1 (Eq.5) for an event that has been detected.Going from the first line in Eq. ( 9) to the second, it is tempting to write p(tr|d, class a , Λ) as something related to the survey efficiency functions commonly published along with survey data.However, this is incorrect, as this probability is also conditioned on the data d, which takes precedence.While some microlensing surveys select events based on model parameters like u 0 and t E (e.g., Husseiniova et al. 2021), it should be noted that these are maximum likelihood estimates completely based on the characteristics of the data.These fitted parameters do not share the same meaning as the parameters θ, which are implicit when conditioning on class a in Eq. ( 9).Instead, we note that the correct interpretation is p(tr|d) = 1 for an event that is known to have a trigger, regardless of event model parameters or event classification.This means that the detection efficiency does not play a role in the likelihood of individual events1 .
Assuming that our set of considered lens subpopulations is complete, the evidence of a single lens (the denominator of Eq. ( 8)) is, where we are summing over the finite and complete set of lens classes.Here we have also used again the fact that tr only depends on d, and therefore p(d, tr|class a , Λ) = p(d|class a , Λ), from the arguments above.This leads to the following identification: p(d, tr|Λ) = p(d|Λ), or that the presence of a trigger does not carry any additional information not already contained in the stream of data itself.
After simplifying the terms in Eq. ( 8), the dependence on tr disappears, so that p(class a |d, tr, Λ) = p(class a |d, Λ).Fundamentally, this is because selection effects are the embodiment of factors that lead to some signals truly in the dataset being classified as an event (p(tr|d) = 1) and others to be missed (p(tr|d) = 0).When considering one event (already designated as a detection), selection effects play no role even when performing the analysis in the context of population models.However, as we will see in the population analysis in Sec.3.3, tr enters into the formalism when accounting for the fact that the full dataset being considered is incomplete.
We can now write Eq. ( 8) in a form that can be computed by introducing θ, Practically, we can compute the integral on the right hand side by importance sampling if we have S independent posterior samples θ c ∼ p(θ|d) drawn under some prior, π(θ), with wide support (Hogg et al. 2010), Here, the evidence for the single event analysis p(d) was absorbed into the updated evidence p(d|Λ) as an overall constant.With Eq. ( 12), we can leverage previously calculated posterior samples to address the question of lens classification for a single event.

Population Analysis
We now turn to inferring Λ using a set of N obs different microlensing events, {d i }, and detection information {tr}.The posterior probability density of Λ is well documented in the literature (e.g., Loredo 2004;Vitale et al. 2020;Mandel et al. 2019;Taylor & Gerosa 2018), so we state the result here and leave a detailed derivation to App. A. We have, Here, p(Λ) is the population parameter prior and p({d i }, {tr}, N obs ) is the evidence.The factor e −αN N N obs follows from assuming events are generated via an (inhomogeneous) Poisson process (e.g., Youdin 2011;Loredo 2004), which penalizes population models that do not predict the correct number of detected events (αN ).This factor accounts for selection bias by marginalizing over the unknown events in the data that fail to rise above the detection threshold.The quantity L obs i is the marginalized event likelihood, that is independent of selection effects for the reasons outlined in Sec.3.2.
The expression derived in Eq. ( 13) differs from those of past works, for example Mroz et al. (2021) (see their Eq.2) and Sumi et al. (2023) (see their Eqs. 4 and 10, where the total likelihood is the product of these two expressions).Both of these works neglect information about the overall rate of events (including the effects of Poisson statistics), which can be an important piece of information when disentangling population information.Eq. ( 13) generalizes these past methodologies such that differential rate information, which can differ dramatically between subpopulations of lenses, is formally included in the analysis.In the case that rate information is still deemed unnecessary, it can be marginalized out of Eq. ( 13) formally, with an appropriate choice of prior (Fishbach et al. 2018).Furthermore, we generalize the work of Sumi et al. (2023) by deriving a more nuanced treatment of selection affects, replacing the detection efficiency factor of their work with the integrated detection efficiency α, defined in Eq. ( 13).Mandel et al. (2019) showed that the treatment of Sumi et al. (2023) can lead to biased conclusions when considering systems with strong selection effects.While these past works were less susceptible to these differences because of a restricted focus on certain ranges of timescales, both of these extensions become increasingly important when considering a fully global analysis of the data, simultaneously considering all subpopulations with overlapping predictions of event parameters.
In addition to the full model in Eq. ( 13), we also explored the use of a restricted model.In this restricted model, we exploit the fact that we are using a mixture model for lens classes and fix the parameters of the lens subpopulations {λ a }, only allowing the subpopulation mixing fractions {ψ a } to vary.While this method is ≈ 10 2 times faster to compute than the full model, it can be susceptible to biased inferences for our problem.We detail the restricted model's exploration in App.B as it might be of use to other astrophysical problems or for future microlensing applications once these subpopulations are better understood.
With Eq. ( 11) and Eq. ( 13), we have the machinery to understand individual events in the context of their populations and extract hierarchical information from noisy, biased surveys robustly.

VERIFICATION DESIGN
In this section we describe the process of validating our proposed methods.To do this, we simulate microlensing events from our population models along with a microlensing survey.We then attempt to recover and disentangle the injected lens populations with our method.This self-consistent testbed enables the evaluation of our methods efficacy in an environment free of systematic bias.
As a specific test case, we consider five different population models denoted by {Λ 0 , Λ 1 , Λ 2 , Λ 3 , Λ 4 }, constructed to mimic features of the BH mass spectrum being reported through gravitational wave observation by the LIGO-Vigo Collaboration (LVC) collaboration (Abbott et al. 2021b,c).Namely, we utilize a BH mass spectrum with a power law component and a Gaussian component, as a model with these features yielded the best fit from LVC's analysis.With data generated from these population models, we use our methodology to determine if these features could be consistent with a subpopulation of PBHs and detected via microlensing.As a byproduct of these tests, we evaluate the ability of this method to constrain aspects of the stellar and SOBH subpopulations.Below, we outline the population models used in Sec.4.1, and our simulated microlensing survey in Sec.4.2.Finally, we describe the numerical methods used to construct and analyze a population realization in Sec.4.3 and Sec.4.4, respectively.

Population Models
We consider three intrinsic lens subpopulations (i.e., before observation) -Stars, SOBHs and PBHs, which only differ in their mass spectra.It is critical to note, however, that this framework can be extended to include inference on all hyperparameters of the lens subpopulations, including lens/source velocity distributions, spatial distributions, etc. Marginalizing over these uncertainties will be crucial to applying these methods to real data, but that extension is left to future work.Starting with the mass distribution adequately shows the effectiveness of this approach.
We structure every population model in the same way, only varying the number of PBH sources in the data.Common to all population models considered in the verification process, the stellar and SOBH subpopulations are represented by a Pareto type-II power law mass distribution implemented in SciPy (Virtanen et al. 2020) starting at 0.07M ⊙ and 5M ⊙ , respectively.These mass distributions have the form, Here, M , µ and σ are in solar masses.M min = µ+σ is the minimum mass.The parameter b is the tail parameter for a Pareto distribution, and is related to the spectral index of a power law through the relation b+1.The PBH subpopulation is described by a Gaussian centered at 30M ⊙ with mean (µ) and standard deviation (σ) parameters.The adopted parameter values and shapes of these mass spectra are shown in Table 1 and Fig. 2, respectively.For numerical stability we also require all lenses to have M < 1000M ⊙ which is well over the pair instability mass limit of SOBHs (Vink 2015;Heger & Woosley 2002;Heger et al. 2003).
In addition to lens masses, to generate microlensing events in all models we assume: D L , D S ∼ U(2000pc, 8000pc) where D L < D S , b sff ∼ U(0, 1), and ϕ ∼ U(0, 2π), which were chosen to be reasonably physical, simple distributions.We assume events have a baseline magnitude I Base ∼ N (µ = 21.066,σ = 1.780) with a reference point of F ref = 1 corresponding to I ref = 22 and log 10 µ rel ∼ N (µ = 0.81, σ = 0.21), which were chosen to be consistent with the PopSyCLE simulations (Lam et al. 2020).Finally, we assume events happen along a random line of sight, giving a variety of parallax orientations, and random peak magnification time t 0 uniformly distributed from 0 to 3650 days.For all models, we fix the number of detected stellar and SOBH lenses to N det Star = 3225 and N det SOBH = 27, respectively, to be reasonable approximations for what to expect from current microlensing surveys like OGLE while remaining computationally tractable.Averaged over many realizations, this corresponds to N Star = 15530 and N SOBH = 100 (detected or undetected) events given each subpopulation's population efficiency, α Star and α SOBH (also an average quantity), and noting the relationship between the two: N a = N det a /α a .The only difference between all the population models being described in this work is the relative contribution of the PBH Gaussian bump in the mass spectrum.For Λ 0 , PBHs contribute roughly 100% of the dark matter in our galaxy (f PBH = 1) with a PBH relative abundance of ψ PBH = 0.032 (Pruett et al. 2022).While this size of a PBH subpopulation is ruled out by observation and experiment for this range of masses, it provides a good point of reference for the analyses of the other population models.For Λ 1−4 , we progressively step down the   The true mass distributions used in our toy model universe.The distribution is comprised of three subpopulations, meant to reflect realistic distributions in nature.Note, the subpopulation mass distributions are normalized independently, so the relative amplitudes are not indicative of what was used to produce the data.contribution of the PBH subpopulation until f PBH = 0 (see Table 2).

Survey Design and Selection Criteria
For our toy microlensing survey, we adopt OGLE-like (Udalski et al. 2015) characteristics; a ten year survey with a 3 day cadence and a magnitude measurement error of σ N = 0.1 mag for all magnitudes.We neglect gaps in the data from seasonal observations, leaving more realistic survey designs to future work.While this cadence is more suggestive of OGLE-III than OGLE-IV and might affect the short timescale end of the t E distribution, this work focuses on the long timescale end of the distribution.Small changes to the observation cadence in the range of hours to days should not drastically impact the long timescale end of the t E distribution, as   2. Each set of data Λi was created with the same models for the stellar and SOBH subpopulations, while the PBH subpopulation was logarithmically decreased from fPBH = 1 to fPBH = 0.
evidenced by a comparison of the OGLE-IV detection efficiency curve and the detection efficiency curve in these simulations shown in Fig. 18 above ∼10 days.Furthermore, the choice of using σ N = 0.1 mag is conservative when considering the estimated variance of the magnitude measurements from Mróz et al. (2019) (Fig. 9).While remaining in the white, Gaussian, and stationary noise limit, using a more accurate estimate of the magnitude measurement variance which varies with magnitude will only improve the conclusions of this work.For detection thresholds, we use a simplified version of the OGLE-IV (Mróz et al. 2019) criteria.Namely, for each lightcurve, we take the maximum flux, F max , and calculate the average baseline flux, F base , and variance of the baseline, σ 2 base , using the data greater than 360 days away from F max (cutting out a 720 day window).We calculate the significance through the χ 3+ parameter, defined as where i indexes the flux measurements and begins at F max , including all consecutive points above F base + 3σ base .If χ 3+ < 32, we classified the event as a nondetection.For an event to be detected we also require it to have baseline magnitude less than 21, and the corresponding change in magnitude between F max and F base to be > 0.1 mag.

Population Model Realization and Single Event Inference
To create a simulated catalog of microlensing events we draw samples of θ from the population model.For each of these microlensing events parameters we simulate a light curve according to Sec. 4.2, corrupting the data with white Gaussian noise with standard deviation of σ N = 0.1 mag.If the events meet the detection criteria in Sec.4.2, we add it to our microlensing event catalog.This process is continued until we have N det a events for each subpopulation, a ∈ {Star, SOBH, PBH} (as outlined in Table .2).This is regardless of what N a should be for each subpopulation as calculated by N det a /α a = N a , which is the average quantity over many realizations.By fixing N det a for each subpopulation instead of the intrinsic number N a , we are able to more directly compare the outputs of the various realizations.Fig. 4 shows an example of a simulated event that meets detection criteria with t E = 125 days and π E = 0.35.
To obtain posteriors samples of θ for the events in our simulated catalog we first transform certain parameters to log space ({log F Base , b sff , log t 0 , log t E , log u 0 , ϕ, log π E }) to increase sampling efficiency.We then use a custom Markov chain Monte Carlo (MCMC) sampler defined in (Perkins et al. 2021;Perkins & Yunes 2022), which has been validated in 15+ dimensional parameter spaces with jagged, multi-modal features in the posterior.This sampler is built on the concept of parallel tempering (Swendsen & Wang 1986;Earl & Deem 2005) to efficiently explore multi-modal posteriors, and utilizes Fisher information matrices to construct efficient proposal densities.In each run of the sampler, random draws from the prior are used at starting points to ensure we are not biased by starting at the known true  values and the sampler is run until ≈ 1000 independent samples are collected which is determined by the chain auto-correlation length.
For our prior distributions, we used priors we deemed to be appropriately uninformative for each parameter.This meant sampling uniformly in the flux F Base , the time of maximum magnification t 0 , the impact parameter u 0 2 , the blending fraction b sff , and angle ϕ.We sampled uniformly in the logarithm of t E and π E (i.e., uniform in scale for these parameters) to ensure proper exploration across the entire parameter space.The priors are summarized in Table 3.We assume a likelihood consistent with our simulated Gaussian white noise model 3 , ln L(d|θ) ∝ − 1 2 where the index t runs across the entire data L data and I(θ, t) is the prediction for the magnitude from our microlensing model, as a function of the event parameters and time.Fig. 4 shows the reconstruction of an example lightcurve in the synthetic catalog with the reconstruction from this analysis overlaid.Fig. 5 shows the inferred 2 We restrict our analysis to u 0 > 0 which neglects degeneracies leading to multi-modal posteriors (Gould 2004).We can do this in the current study as the simulated data was also restricted to u 0 > 0, ensuring no bias.Analysis of real data will have to be more careful to asses the effects of multi-modality (e.g., Kaczmarek et al. 2022).Although, its worth noting that work focused on compact objects is less sensitive to these issues as compact object preferentially fall into the low π E space, softening the degeneracy and leading to the joining of these different modes. 3This ensures no systematic bias for this preliminary study.Of course, future work could begin to relax this requirement to study systematics or more complicated noise models expected in actual data (Golovich et al. 2022) posterior in t E −π E space for the catalog of events drawn from Λ 0 and Λ 4 .

Population Inference
For the population-level parameters, we use the same MCMC software outlined in Sec.4.3 to obtain posterior samples.To perform the inference, we will need to forward model a population model.In this work, we will be using the same simple models outlined in Sec.4.1, but when using real data, population simulations like PopSyCLE will need to be implemented (although this pipeline is independent of the exact forward model method being employed).We perform inference with two separate population models: once using all three distributions (producing an unbiased estimate free of systematics) and once using only just Stars and SOBHs.The Λ 4 catalog is the only set of simulated data that can be perfectly modeled by only stellar and SOBH subpopulations, as it contains no PBHs.In the case of the other simulated datasets, Λ 0−3 , modeling with just stellar and SOBH subpopulations introduces systematic bias enabling us to understand how a PBH subpopulations signal could be detected (or evade detection).Overall, The priors used for Λ are detailed in Table 4 and were chosen to be uninformative, uniform, and with boundary values commensurate with the uncertainty in the prior understanding of the subpopulation of stars and SOBHs.Additionally, for the SOBH subpopulation, we stipulate that the minimum mass of the power law distribution must be in the range of [2, 6]M ⊙ , as there is observational evidence for this upper limit from GW observations, X-ray binaries, and radial/photometric observations (Abbott et al. 2023;Mapelli 2020;Özel et al. 2010;Thompson et al. 2019).

APPLICATION TO SINGLE EVENTS
We first explore the implications of treating the class of single events probabilistically and compare it to example, typical cuts in π E −t E (e.g., Golovich et al. 2022) to classify events.Specifically, for the purpose of identifying BHs, we compare our method with a linear cut in log 10 t E -log 10 π E space defined such that 50% of the events with posterior medians below the line are classified as BHs according to simulations from a population model, which is used by (Golovich et al. 2022).
Fig. 6 shows a comparison of these two lens classification methods in the π E − t E space for Λ 4 .Specifically,  .The data and reconstruction from our analysis for a specific event in the synthetic catalog.The top panel shows the raw data (in gray), the true signal in the data (dashed green) and the 90% confidence reconstruction of the signal from our posterior distribution (shaded orange).In the bottom panel, we show the residual, defined by the difference between the reconstruction and the data, divided by the average of the reconstruction and the data.

Stars SOBHs
Figure 5.The median values of the posteriors for all the events in the synthetic catalog produced by population model Λ0 (left) and Λ4 (right) are scattered above, separated out by the subpopulation.For the SOBH and PBH subpopulations, we show the error bars (1σ), as calculated by our posteriors from the single event analysis.As a reminder, Λ4 does not include PBH lenses at all.Note the separation of heavier lenses (SOBH and PBH) down and to the right in tE-πE space, as expected.Roughly speaking, the distribution of events in this lower right quadrant of this space would give insight into the BH subpopulations of lensing objects in our galaxy.the 50% purity line is overlaid on the contours calculated by maximizing the class probability across the stellar and BH lenses (both SOBH and PBH).Fig. 6 reveals that our method captures high-order structure in the intrinsic uncertainty in the lens class predictions from the population models missed with the 50% purity line.
Firstly, Fig. 6 shows that there are regions of π E − t E space that do not trace the 50% purity line.Intrinsic to the population model itself, an event cannot have more than ∼ 50% lens class confidence of a BH vs a Star (dark regions), even if π E and t E are known perfectly.This has implications for allocating followup resources for events in progress, because if an event lies in one of the dark regions, then taking further high-cadence and high-precision photometric data as the event the event evolves to shrink the π E − t E posterior will not improve lens-class confidence.Conversely, events with diffuse constraints on π E − t E that are in areas of light contours could have their class confidence improved with followup observations.
To quantify the advantage of using our probabilistic lens classifications over purity cuts to identify BH candidates, we test both methods recovering the known BHs in the simulated datasets Λ 1−4 .Fig. 7 shows the purity of recovered BH candidates vs the fractions of PBHs in the simulated datasets.We test the purity cut method using three different priors on the individual event parameters -uniform in log π E and log t E (Table 3), a broad normal distribution for both parameters in Table 1 of Golovich et al. (2022), and uniform priors in t E and π E .For both methods, arbitrary thresholds need to be used to select lenses.In the case of the method proposed here, one must define the threshold probability, p threshold (class BH |d, Λ).For the linear method of Golovich et al. (2022), one must specify the target purity when calculating the line.To assess the impact of these two parameters, we consider two choices, namely a p threshold (class BH |d, Λ) of 0.5 and 0.9, and a target purity of 50% and 90%.However, defining a threshold probability, as we will see below, does not correspond to setting the final purity of the classification analysis.For all methods, Fig. 7 shows that more BHs are recovered and the purity of each sample increases as the number of PBHs in the dataset increases.This is due to PBHs increasing the abundance of BHs which populate the lower right corner of the π E − t E space (see Fig. 5 in Pruett et al. 2022), making it easy to separate from other lens classes.Fig. 7 also shows that the probabilistic lens class method outperforms the purity based methods across all simulated datasets as measured by the final purity of the sample.The largest performance gains are for simulated datasets with low numbers of PBHs.While not predicting quite as many correct BH candidates, the probabilistic method does lead to far 0.0 0.5 1.0 1.5 2.0 2.5 3.0 log 10 (t E ) (days)  2022).The contours are derived by taking the maximum probability between an event to include either a stellar or BH (both SOBH and PBH) lens.These are calculated assuming perfect measurement, i.e., that the posterior for the event parameters tE and πE are delta functions.The regions of dark shading illustrate regions in parameter space where intrinsic overlap in the predictions from different subpopulations fundamentally limits our ability to classify an event with these methods.With photometry alone, classifying events that fall in the dark green regions of parameter space will not be improved drastically even with infinite observational precision.The light regions reflect parts of parameter space where classification is highly certain.We note that the exact location of the contours fluctuate with numerical noise in our simulation, but the general structure is robust across different realizations.This is particularly true in regions with low expected rates, like high πE and high tE.
fewer false positives, as shown in the upper two panels of Fig. 7. Fig. 7 also shows the sensitivity of the purity cut methods to the prior distributions used when modeling each event.All the priors appear to be uninformative in different ways, however, the log-uniform prior selects at least an order of magnitude more BH candidates and a significantly less pure sample across all simulated datasets when compared to the prior used by Golovich et al. (2022) or the uniform prior.This is due to smaller objects, like stars, having poorly measured π E .In this case, the constraint on π E for stars is driven by the prior.Moreover, the purity cut methods relies on the posterior mean and not the full distribution.Overall, if the π E prior mean is in the region of π E − t E space which is dominated by BHs (as the log-uniform prior is), it will bias all the events to be classified as a BH when using the purity cut method.When using the uniform or normal distribution, we see a much better performance, where these two choices generally agree. .Top: the number of BH lenses correctly identified as BH candidates for each simulation according to each selection method (true positives): the method of this paper (green), the linear method using the priors from Golovich et al. (2022) (orange), the linear method using linear priors in tE and πE (purple), and the linear method using log-uniform priors in tE and πE (pink).The width of the band is calculated using two different tuning criteria.In the case of the method proposed in this work, we used a threshold probability of p(classBH|d, Λ) > 0.9 and p(classBH|d, Λ) > 0.5.In the case of the linear method, we used lines designed to have 50% and 90% purity.An important note is that this figure should not be interpreted as a method for measuring the abundance of BHs.The abundance is a population-level parameter and is more appropriately handled in the analysis of Sec. 6. Middle: the number of stellar lenses incorrectly identified as BH candidates for each simulation for each selection method (false positives).Bottom: the purity of these classifications, defined as the number of correctly identified BH candidates divided by the total number of BH candidate classifications.We see the probabilistic classification method of this paper outperforms the linear method with any of the single-event priors considered here, in terms of the highest purity.Furthermore, the probabilistic method is fundamentally independent of the priors used in the single event analysis and is generally robust to changes in the arbitrary threshold parameter (in terms of purity).
In contrast, Fig. 7 illustrates the insensitivity of the method in this paper to the arbitrary threshold probability.Changing the threshold from 0.9 to 0.5 increased the number of candidate events, but at an almost identical purity.The additional candidates gained by changing the threshold were equally likely to be a BH as a star.This comes from the distribution of p(class BH |d, Λ) for each survey, where the integral of this distribution ultimately determines the purity, not the lower boundary.On top of this insensitivity, we also note that this method is independent of the priors used when analyzing single events (see Eq. ( 12)), removing a possible source of systematic bias.
The above tests assume that we know the underlying true lens population model, which in reality is not true.To mitigate this, our probabilistic lens classification method can marginalize over a set of possible population models.In this case, instead of a point estimate for probability of the lens classification, we have a distribution of possibilities which captures the underlying lens population uncertainty.Fig. 8 shows the distribution of p(class BH |{λ a }, d) for a single event obtained when marginalizing over a set of restricted population models that only allow lens subpopulations mixing fractions to vary (see App. B).Fig. 8 shows that we were able to estimate the true p(class BH |{λ a }, d) accurately despite not perfectly knowing the underlying relative abundances of the different subpopulations.For illustration purposes, the bottom panel shows a graphical representation of the predicted distribution in t E -π E space for each subpopulation and the posterior of the specific event.
Two features in the results of this section provide compelling evidence for why the hierarchical analysis (results discussed in the next section) should incorporate the probabilistic nature of classification proposed in this work.First, our classification method never captures all the BHs in the data, as illustrated by the green band always falling below the markers representing the true number of BHs in each dataset in Fig. 7.The standard classification schemes and the formalism of this work either miss a large fraction of BHs in the data (with a high purity) or include many BHs but accompanied with many false positives (giving a low purity).This risks two types of bias when considering hierarchical inference: neglecting an important subset of BHs in the data or biasing results by including incorrectly classified events.Furthermore, the majority of the results in this section were achieved by conditioning on a specific population model.The impact of this choice is shown in Fig. 8, where the spread on classification confidence for this event can change an appreciable amount based on the uncertainty in the population model.This illus- Figure 8. Top: the probability distribution of a specific event to involve a BH lens (from Λ2, in this case), marginalized over the posterior distributions of mixing fractions for the underlying population model (assuming fixed shape parameters {λa}inj matching those used to create the data).The vertical line represents the probability for the event to involve a BH lens performing the analysis with the entire population model used to create the data (assuming both shape parameters and abundances fixed to those of the model which created the data, {λa}inj and {ψa}inj).Bottom: the posterior of the event in question (shown in black) as compared to the predicted distributions of each subpopulation in tE-πE space, for comparison.This particular event has a high probability of being a BH lens because of its large overlap with the SOBH and PBH distributions, despite uncertainty on the expected contribution of each subpopulation to the overall lens population.Critically, the bottom panel of this figure only represents a graphical representation of the likelihood.To determine the total probability a lens belongs to a certain class, one must also incorporate the prior probability.The contours are linearly spaced between 0 and 3.5.
trates the need to jointly infer the lens type along with the entire population model, simultaneously.In the next section, we demonstrate how our methods for hierarchical inference robustly address these issues and provide unbiased results.

APPLICATION TO POPULATIONS
To understand the lensing population model, the results of our analysis are broken down into several parts.We first consider the PBH subpopulation in the context of both its hyperparameter posteriors and through Bayes' factors.We then move on to consider the stellar and SOBH lens subpopulations.

PBH posterior information
We begin by evaluating our ability to measure the relative abundances of the different subpopulations, focusing first on PBHs.Fig. 9 shows that our ability to detect PBHs varies with the number of PBHs in the data.When a significant number of PBHs exist in the data (Λ 0 ), we recover an accurate ψ PBH posterior inconsistent with zero (>4.5σ), providing evidence for a PBH subpopulation.As the number of PBHs in the simulated data are stepped down, we find strong (Λ 1 ), then mild (Λ 2 ), and finally no (Λ 3,4 ), constraint with the recovered ψ PBH posterior being inconsistent with zero.As the number of PBHs in the data decreases, our ability to measure a non-zero PBH abundance diminishes.However, even in the case of Λ 3,4 , we can still place an upper bound on the PBH subpopulation, which in this case, plateaus at ≈ 20-25% of the DM fraction (f PBH ) and at ≈ 1% of all lenses in the population (ψ PBH ).
The mixing fraction of SOBHs, ψ SOBH , also encodes information about the existence of a subpopulation of PBHs.Fig. 9 shows the ψ SOBH posterior for both the two (Star + SOBHs) and three (Star + SOBHs + PBHs) population configuration.When a substantial number of PBHs exist in the data (Λ 0 ) and only the Star + SOBHs population model is used, the SOBHs subpopulation absorbs the PBHs making ψ SOBH ≈ 5 times larger than than its true value.However, this signature of a PBH subpopulation is unlikely to be useful when applying this method to real data, due to it being completely dominated by the factor of ∼100 uncertainty on the expected number of SOBHs in the Milky Way (Samland 1998;Timmes et al. 1996;van den Heuvel 1992).
The Star + SOBHs model cannot, however, completely absorb and explain away the PBHs in the case of a large number of PBHs actually contained in the data (Λ 0 and Λ 1 ).If the SOBH and PBH subpopulations were perfectly degenerate, the two ψ SOBH posterior distributions would be wide but overlapping, extending from ∼0 to ∼0.05 for both classes of models.This is because the sum of the PBH and SOBH subpopulations would always account for the total contribution of both subpopulations.However, because the population model favors not having a large number of SOBH lenses but instead tends toward an SOBH subpopulation consistent with zero relative abundance when includ- .Shown above are the posterior distributions obtained from out simulated datasets for the mixing fraction or relative abundance of the PBH (ψPBH left) and SOBH (ψSOBH right) subpopulation.In the case of the PBH subpopulation, we also show its contribution to DM (fPBH).Vertical lines indicate the true value used to create the data.The prior on the absolute abundance of each subpopulation, N classa , was uniform between 1 and twice the number of lenses in the specific simulation, which translated to a prior on the relative abundance that was broad (between 0 and 1) but mildly peaked around 0.5.For the PBH subpopulation, we see an indication that the PBH subpopulation contributes meaningfully to the explanation of the data, due to the posteriors being inconsistent with zero in the case of datasets Λ0 and Λ1 and mildly contributes in the case of Λ2.In the case of the other datasets (Λ3 and Λ4), we can merely make statements about the maximum contribution of PBHs to the lensing population (∼ 1%) and their contribution to the DM fraction (∼ 20 − 25%).While datasets Λ3−4 peak slightly away from zero, the mean of the distribution is less than 1.5σ and is due to correlations with the SOBH abundance and the marginalization process.For the SOBH subpopulation, we see a separation of the posteriors when using two and three subpopulations.This indicates that including the PBH component leads to a better description of the data, showing preference for the more complex model.The (Star + SOBH) model is not flexible enough to explain the data when considering Λ0, Λ1, Λ2.
ing a PBH subpopulation in the modeling, a hierarchy emerges in the explaining power of each class of population model.In this case, using both the SOBH and PBH subpopulation to describe the entire BH subpopulation is more informative than the population model which neglects the PBH subpopulation.Fig. 10 shows that our method can recover the lens mass spectrum across all lens subpopulations.For the simulated datasets Λ 0−4 and for the two (Star + SOBHs) and three (Star + SOBHs + PBHs) population models, both the prior and posterior distributions on the lens mass spectrum are shown.In all cases, when the true model that generated the simulated data is used, we recover an unbiased lens mass spectrum in agreement with the true distribution to within 90% credibility.The disagreement between the two classes of populations models is greatest for the datasets with the most PBH lenses.Without the flexibility of the PBH subpopulation component, the SOBH power law subpopulation shifts to compensate for the missing category of lenses.This leads to fine-tuning of the SOBH mass spectrum in the population model that neglects the PBH subpopulation because only a narrow part of the SOBH subpopulation parameter space yields reasonable agreement with the data.This effect can been seen in Fig. 10 as a nar-rowing of the SOBH mass spectrum between 10−10 3 M ⊙ for Λ 0,1 .
Fine tuning is an aspect of model complexity that must be considered when evaluating competing models.For detecting a subpopulation of PBHs, we have a more flexible population model (Star + SOBHs + PBH) that we have to compare against less flexible population model (Star + SOBHs) that requires fine tuning to explain the data.Overall, Fig. 10 provides a diagnostic in the model selection problem of determining the evidence for the additional subpopulation of PBH lenses, where we see a systematic inability of the simpler model to accurately recover the true distribution.
As the number of PBHs in the data decreases (from Λ 1 through Λ 4 ), our ability to disentangle the subpopulation mass spectra drops.This is shown in Fig. 10 by the similar mass spectrum reconstructions between the two and three subpopulation models.This suggests that the extra flexibility of the higher dimension model is unwarranted or is fitting the noise.In the case of Λ 1,2 , there is still marginal evidence for the existence of the PBH subpopulation, although this is difficult to claim based purely on the mismatch between the posterior mass spectra between the two subpopulation models. .Shown above are the various reconstructions of the mass spectrum of lensing objects from our ten analyses performed on five sets of data.The range of reconstructions allowed by the prior is shown in the top panel.Each subsequent panel shows the 90% confidence reconstruction from the posterior probability of the inference analysis for each set of data, Λi.For each dataset, we conduct two analyses: one assuming two subpopulations (hatched) and one assuming three (unhatched).The true distribution used to create the synthetic data is shown as a black dashed line.From this figure, we can see that our analysis provides an unbiased reconstruction of the underlying true distribution, accurate to within 90% confidence when systematic bias is not present.The ability to disentangle the subpopulations is clear when a large subpopulation of PBHs is present in the data (Λ0), but these conclusions become increasingly more uncertain as the "strength" of the PBH subpopulation shrinks (Λ1 to Λ4).
Fig. 11 shows the posterior constraints on µ PBH and shows that for Λ 0,1 information can be inferred about the structure of the PBH mass spectrum.We see strong evidence for the PBH mass spectrum bump around the correct location of 30M ⊙ .The recovered PBH mass spectrum bump is always wider than true values, which The true values used to create the data are shown as solid vertical lines.For reference, the prior for this parameter was uniform between 1M⊙ and 80M⊙.For Λ0 and Λ1, the posteriors favor the true value.The other datasets are less informative, simply ruling out a high-mass component of the lensing population.The lower mass parts of the spectrum overlap with the stellar subpopulation (which dominate the catalog by a large margin), allowing the extra flexibility to be absorbed by this primary subpopulation.
indicates that our method is not sensitive to the width of the PBH mass spectrum bump compared to the location of its peak.For Λ 2−4 , Fig. 11 shows no constraint on µ PBH , therefore the data did not favor a high-mass PBH component in the lensing population, and only an upper bound can placed on the mass range of the PBH subpopulation.In these cases, the lower part of the PBH mass spectrum overlaps with the dominating stellar subpopulation, which can absorb the PBHs as noise.

Evidence for a simulated PBH subpopulation
In addition to examining the posteriors of the two subpopulation models, we can also compare their overall performance on the datasets directly.There are many statistics that can be used to compare competing models which all have their advantages and drawbacks.From χ 2 -based metrics (e.g., Wyrzykowski et al. 2016;Andrae et al. 2010), to information criteria (e.g., Kains et al. 2018) to Bayes' factors (e.g, Jenkins & Peacock 2011) and cross validation scores (e.g., Welbanks et al. 2023;McGill et al. 2023), these statistics can estimate and approximate different aspects of model performance.Here we compare models using the maximum likelihood and the Bayes' factor, where the Bayes' factor is estimated as a byproduct of the parallel tempering MCMC methods described in Section 4.3.The Bayes' factor is a widely used method for model selection and its advantages include its interpretation as the comparison of the posterior probability for each model, and that it penalizes model complexity not supported by the data.The Bayes' factor's main draw back is its sensitivity to prior distributions.Despite the Bayes's factor not being the perfect model comparison tool, we find it informative for our problem.Fig. 12 shows that for all simulated datasets the Bayes' factor always disfavours a PBH subpopoulation.This is driven by the wide, uninformative priors used in all models (Sec.4), and that the two subpopulation model can partly absorb the PBH subpopulation with some fine tuning, as discussed above.While sufficient evidence cannot be found for any of the data sets in isolation with the Bayes' factor, there is a strong trend in Fig. 12 showing that as the number of PBH lenses in the data increases, our ability to distinguish between the two subpopulation model classes improves significantly.The difference in the logarithm of the maximum likelihood increases by ≈ 4 − 5 (equivalent to the maximum likelihood value of the three subpopulation model increasing by a factor of ∼ 100) and the logarithm of the Bayes' factor increases by ∼ 7.While conclusive evidence for a subpopulation of PBH lenses cannot be claimed in this toy example with the Bayes' factor, the trend of its improvement from Λ 0−4 suggests that given sufficiently informative priors it could be used to determine the presence of a PBH subpopulation.

Stars and SOBHs
Beyond PBHs, we can asses the ability of the model to infer features of the SOBH and stellar mass spectra.Fig. 13 shows the posterior constraints on M min,Star and M min,SOBH .We find that for all population models and simulated datasets we are able to obtains a tight constraint on M min,Star of roughly 0.08M ⊙ ± 0.02M ⊙ .The large number of stellar lenses in the data suggests there is enough information to make robust claims about the minimum stellar mass.Conversely, we find M min,SOBH is never constrained due to its posterior distribution always approximately recovering the prior distribution.This suggests we are not able to probe the existence and properties of the mass gap between NSs and BHs (Farr et al. 2011;Özel et al. 2010;Fryer & Kalogera 2001).While this conclusion should be revisited for future surveys with tighter measurements and larger catalog sizes (e.g., Roman Space Telescope; Spergel et al. 2015) or when other data is taken into account (such as astrometry), our initial analysis using photometric only, OGLEtype data does not give confidence that this will be a measurable feature.The vertical lines indicate the true value used to create each dataset, while the hatched distributions refer to inference performed assuming two subpopulations and unhatched distributions were inferred assuming three subpopulations.The prior for Mmin,Star is a uniform distribution between 0.01M⊙ and 1M⊙.The prior for Mmin,SOBH is a uniform distribution between 2M⊙ and 6M⊙.We see that we can robustly measure the minimum mass of the stellar power law distribution, regardless of the data set or population model employed.However, we cannot infer the minimum mass of the SOBH power law distribution, regardless of the data set or population model employed.The data sets and type of data (photometry only) do not contain enough information to place robust bounds on the minimum mass of the SOBH subpopulation that would be useful in measuring things like the NS-BH mass gap (Farr et al. 2011;Özel et al. 2010;Fryer & Kalogera 2001).Posterior distribution for the shape parameter of the stellar (left) and SOBH (right) power law distribution (bStar and bSOBH).The vertical lines indicate the true value used to create the data.The two types of analyses (assuming two or three subpopulations) are shown as a hatched and non-hatched distribution, respectively.The prior for these parameters is a uniform distribution between 0.1 and 10.For the stellar subpopulation, the posterior is approximately the prior distribution because of correlations with the location parameter, µStar.This correlation is shown in Fig. 15.Considering the SOBH subpopulation, the posteriors are uninformative in the case of assuming three subpopulations.However, when only considering two subpopulations (thereby focusing on astrophysics as opposed to exotic physics), the posteriors contain significantly more information than the priors.The slight bias in the posteriors are connected to the same correlation shown in Fig. 15, but for the SOBH subpopulation model, and the process of marginalization.

CONCLUSIONS AND FUTURE WORK
In this work, we proposed and validated a methodology to conduct hierarchical inference simultaneously with the lens classification for individual events.The benefits of this framework over existing methods include properly accounting for Poisson statistics, measurement uncertainty, and selection bias while not assuming definite criteria for lens classification.This framework allows marginalization over population uncertainty, which is key to reliable inference given the current uncertainty about underlying lens population characteristics.
On a single event level, our method outperformed current purity-cut strategies used to classify lenses in search of BH candidates.We were able to recover purer samples of BH events while also quantifying population model uncertainty.Our probabilistic lens classification scheme also revealed and quantified intrinsic degenerate structure in the t E − π E -space.Further investigation of this structure could yield insights into reliably identifying BH microlensing events in real-time photometrically and efficiently allocating astrometric followup resources (e.g., Lu et al. 2016;Sahu et al. 2022;Lam et al. 2022).Although photometric microlensing parameters are difficult to constrain early in an event's evolution (Albrow 2004), our classification method could be used to classify an event at or after its photometric peak, where t E and π E are better constrained, to decide whether to use astrometric resources to measure the second astrometric peak.(e.g., Dominik & Sahu 2000).Further work on how well microlensing parameters can be constrained from a partial lightcurve would also benefit classification efforts (e.g., Dominik 2009).Real-time identification of black hole microlensing events is likely to become more important in the era of Rubin LSST survey planing (e.g., Street et al. 2023) and the integration of automated identification and followup planning methods into Target and Observation Management systems (e.g., Street et al. 2018;Coulter et al. 2023) our classification method could be used to classify an event at or after its photometric peak, where t E and π E are better constrained, to decide whether to use astrometric resources to measure the second astrometric peak.Further work on how well the microlensing parameters can be constrained from a partial lightcurve would be fruitful We find that our full hierarchical model leads to inference on the lens population mass spectrum and abundances that is accurate and effective while appropriately handling population uncertainty.In the context of characterizing a PBH subpopulation of lenses, our method produces posterior constraints on the PBH abundance inconsistent with zero at >4.5σ when considering f PBH = 1.For the more realistic case of f PBH ≲ 0.25 our ability to identify the subpopulation begins to deteriorate, and we are only capable of placing upper limits on PBH subpopulation.Moreover, a PBH subpopulation signature for any f PBH will likely be derived through a joint analysis of maximum likelihood measurements, Bayes' factors, and hyperparameter posteriors.
The results here can be compared to the constraints from past collaborations (e.g., Allsman et al. 2001;Tisserand et al. 2007;Wyrzykowski et al. 2009Wyrzykowski et al. , 2010Wyrzykowski et al. , 2011a,b;,b;Blaineau et al. 2022) which claim constraints on the DM fraction of MACHOs in the range of ≈ 2 − 20% for different mass ranges between 10 −7 −10 3 by studying events towards the Magellanic Clouds.However, direct comparison to the original microlensing MACHO constraints will have to be performed carefully due to two main reasons.Firstly, the effects of systematic noise arising for events detected towards Bulge (PBH confusion with SOBH and Stellar lenses) vs the Magellanic clouds (Galactic disk and self-lensing; Wyrzykowski et al. 2011b) are different.Secondly, the method presented in this work jointly infers a PBH mass spectrum and abundance which will have to be reconciled with the classic restrictive model of a delta mass function -a distinction in modeling assumptions shown to be important (Green 2016).Although modeling the PBH mass spectrum complicates comparison with the original Magellanic Cloud MACHO constraints, this key improvement is more general and will allow microlensing to contribute to other ongoing studies of the dark mass spectrum (e.g., Abbott et al. 2021b,c;Zevin et al. 2021;Franciolini et al. 2022).
For stellar lenses, the population constraints derived were largely independent of the number of BHs in the data, due to stellar lenses vastly outnumbering BH lenses.For SOBH lenses, the slope of the power law mass spectrum and their relative abundance can be accurately constrained if a PBH subpopulation is not present in the data.When PBHs were injected into the data, degeneracies between SOBH and PBH subpopulation models made it difficult to disentangle the characteristics of the two subpopulations.In all cases, we found the the minimum SOBH mass difficult to constrain.This suggests that extracting information about a possible SOBH mass gap via microlensing (e.g., Wyrzykowski & Mandel 2020) is likely difficult with current photometric surveys and will require more data (through larger and/or longer surveys), increased measurement precision or the incorporation of further information such as astrometric microlensing information.
There are multiple avenues of future research to be taken.The most immediate would be applying these methods to OGLE-IV data, beginning with the context of better understanding the BH population.To accomplish that, the simple population model presented here will need to be expanded to include more realistic information, such as the distributions of the flux blending fraction, the velocity distributions, the spatial distributions, other subpopulations like neutron stars and white dwarfs, etc. Implementing these extensions is purely a practical concern, as they can be formally integrated into the analysis through the framework presented here.These methods can also be extended to better understand other subpopulations of lenses.Variations have already been applied in the context of free floating planets (Sumi et al. 2023), but the framework outlined here can help to improve those constraints by marginalizing over other population uncertainties, accounting for Pois-son statistics and more carefully accounting for selection bias.
Finally, we also note that this methodology can easily be extended to heterogeneous data, i.e., incorporating simultaneous astrometric observations (e.g., for the Roman Space Telescope; Sajadian & Sahu 2023;Lam et al. 2023) or follow-up astrometric measurements from current space telescopes (e.g., Sahu et al. 2017;Kains et al. 2017;Zurlo et al. 2018;Lam et al. 2022;Sahu et al. 2022;McGill et al. 2023).The integration of both types of data will prove to be indispensable as they probe different event parameters and different distributions of events in the galaxy and can break photometric microlensing degeneracies.In the case of current astrometric measurements, the low number of events being followed up suggests a small impact to hierarchical inference.However, when considering future surveys like the Roman Space Telescope, the impact of this joint analysis remains an open question and should be investigated thoroughly.

APPENDIX
We will simply point out the simplifications that can be made if using the population model in Eq. (B15).Each event likelihood can be rewritten as where we note that the integral in L obs i,a can be pre-computed once for each event and each subpopulation, if the shape parameters λ a are fixed.
Similarly, the same argument can be used for α If α a can be pre-computed for each subpopulation, these two simplifications lead to Eq. (B18), allowing for drastic computational speedup at the cost of lost flexibility.We also analytically marginalize the overall rate, N , by assigning the prior p(N ) ∝ 1/N (Fishbach et al. 2018, Appendix A).This marginalization allows us to bypass modeling N , which is difficult due to observing effects like weather and observing schedules.This marginalizing process is possible and frequently useful when inferring the entire population model as well, independent of the mixture model we use here.However, we would like to compare our restricted model to the most flexible version of Eq. ( 13) for our initial validation.This leads us to the posterior of the restricted model, p({ψ a }|{d i }, {tr}, N obs , {λ a }) ∝ p({ψ a }|{λ a })α −N obs p({d i }, {tr}, N obs ) Here, the right hand side is now independent of N , and we have a likelihood that enables quick computation.For L obs i , we have leveraged Eq. ( 11) and can run our classification inference on all the events once, and then reuse those probabilities by combining them as weighted sums.Similarly, because we have fixed {λ a }, we can pre-compute α a values which can be added in weighted sums to calculate α quickly.
We can now test this restricted analysis by using it to perform inference on the relative abundances and compare the results to the full analysis.When sampling using the restricted analysis, we assign a prior consistent with the constraint a ψ a = 1.We will only sample ψ Star and ψ SOBH and use a Dirichlet prior informed by the original simulation (e.g., Golovich et al. 2022).The Dirichlet distribution is parametrized by a vector b, b = ( ψ − min( ψ)) + c , (B19) where c controls the width of the prior, and ψ is a vector of the relative abundances of the injected population model.We chose c = 1 to produce a broad, uninformative prior on the parameters.While more restricted in flexibility, this comes at the benefit of ≈ 10 2 computational speed up.The full population inference takes approximately 360 CPU hours on cluster-grade Intel Xeon E5-2695 nodes, and the restricted analysis takes approximately 1 CPU hour on a consumer-grade Intel Core i9-9980HK chip.Fig. 16 (top row) shows the posterior constraints abundances of the PBH and SOBH subpopulations compared with the full population model used in Sec. 6.We see good agreement between the two analyses, with constraints being tighter in the restricted Posterior distributions on the relative abundance for the PBH (left) and SOBH (right) subpopulations, for both the full analysis and the restricted analysis plotted as the shaded region and the hatched region, The prior for the mixing fractions in the case of the restricted analysis was a broad Dirichlet distribution.There is good agreement between the two methods in the absence of systematic effects, as the true distributions of each subpopulation are assumed to be those of the injection for the restricted analysis.Bottom Row: Posterior distributions on the relative abundance for the PBH (left) and SOBH (right) subpopulation, for both the unbiased restricted analysis using the correct distribution (unhatched) and a systematically biased restricted analysis (hatched).The assumed mass distributions for the biased analysis are shown in Fig. 17.The prior for the mixing fractions in the case of the restricted analysis was a broad Dirichlet distribution.While the conclusions are mildly consistent, a clear systematic bias can be seen in the posteriors.
model due its comparative lack of flexibility.This significantly more computationally efficient restricted analysis gives unbiased and similar constraints on the abundances to the full analysis when we condition on knowing the model that generated the data, however, it has to be validated in a more realistic setting before fully adopted.
To investigate the sensitivity of the restricted analysis to biased assumptions about the underlying population, we take 10-th percentile posterior sample from the full three subpopulation model for Λ 0−4 and run the restricted analysis with an alternate assumption for the underlying form of the mass spectra.This procedure injects bias into our analysis which simulates not knowing the subpopulation models exactly.While any number of alternative mass spectra could have been picked, we chose distributions that were different from the truth while still being consistent with the data in t E -π E space.Fig. 17 (left) shows the alternative mass spectra along the true spectrum used to generate the data.
Fig. 16 (bottom row) shows the posterior constraints on the PBH and SOBH abundances using the restricted model, both assuming an unbiased and the new, biased population model.The general disagreement between these posteriors highlights the fact that biased modeling assumptions can impact the conclusions drawn with the restricted model.
where ϑ is the set of parameters excluding t E such that ϑ = θ \ {t E }.In the first equation, we have also assumed that the distribution of the parameters ϑ are independent of t E , which is true in our toy problem.
We then construct a logarithmically spaced grid in t E that encapsulates the relevant parameter space (0.1 days < t E < 10 6 days).For each point in the grid (t E,l ), we calculate the partially-marginalized detection probability p det (t E,l ).To do this, we draw event parameters from the population distribution p(ϑ|Λ), from which we calculate a light curve.We then perform the second integral over data realization through another Monte Carlo integral, where the noise is drawn from the likelihood p(d|θ l , t E,l ).Given our assumption that the noise is white and stationary we have, (C22) This is now a one dimensional integral which can estimated using, for S independent samples.When considering the restricted hierarchical model, this can also be done for each subpopulation, independently.
Figure 2.The true mass distributions used in our toy model universe.The distribution is comprised of three subpopulations, meant to reflect realistic distributions in nature.Note, the subpopulation mass distributions are normalized independently, so the relative amplitudes are not indicative of what was used to produce the data.

Figure 3 .
Figure3.The mass probability densities for each of the subpopulation models described in Table2.Each set of data Λi was created with the same models for the stellar and SOBH subpopulations, while the PBH subpopulation was logarithmically decreased from fPBH = 1 to fPBH = 0.
for the model containing just Stars and SOBHs and for the model containing all lens subpopulation, respectively.

Figure 4
Figure4.The data and reconstruction from our analysis for a specific event in the synthetic catalog.The top panel shows the raw data (in gray), the true signal in the data (dashed green) and the 90% confidence reconstruction of the signal from our posterior distribution (shaded orange).In the bottom panel, we show the residual, defined by the difference between the reconstruction and the data, divided by the average of the reconstruction and the data.
Figure7.Top: the number of BH lenses correctly identified as BH candidates for each simulation according to each selection method (true positives): the method of this paper (green), the linear method using the priors fromGolovich et al. (2022) (orange), the linear method using linear priors in tE and πE (purple), and the linear method using log-uniform priors in tE and πE (pink).The width of the band is calculated using two different tuning criteria.In the case of the method proposed in this work, we used a threshold probability of p(classBH|d, Λ) > 0.9 and p(classBH|d, Λ) > 0.5.In the case of the linear method, we used lines designed to have 50% and 90% purity.An important note is that this figure should not be interpreted as a method for measuring the abundance of BHs.The abundance is a population-level parameter and is more appropriately handled in the analysis of Sec. 6. Middle: the number of stellar lenses incorrectly identified as BH candidates for each simulation for each selection method (false positives).Bottom: the purity of these classifications, defined as the number of correctly identified BH candidates divided by the total number of BH candidate classifications.We see the probabilistic classification method of this paper outperforms the linear method with any of the single-event priors considered here, in terms of the highest purity.Furthermore, the probabilistic method is fundamentally independent of the priors used in the single event analysis and is generally robust to changes in the arbitrary threshold parameter (in terms of purity).
Figure10.Shown above are the various reconstructions of the mass spectrum of lensing objects from our ten analyses performed on five sets of data.The range of reconstructions allowed by the prior is shown in the top panel.Each subsequent panel shows the 90% confidence reconstruction from the posterior probability of the inference analysis for each set of data, Λi.For each dataset, we conduct two analyses: one assuming two subpopulations (hatched) and one assuming three (unhatched).The true distribution used to create the synthetic data is shown as a black dashed line.From this figure, we can see that our analysis provides an unbiased reconstruction of the underlying true distribution, accurate to within 90% confidence when systematic bias is not present.The ability to disentangle the subpopulations is clear when a large subpopulation of PBHs is present in the data (Λ0), but these conclusions become increasingly more uncertain as the "strength" of the PBH subpopulation shrinks (Λ1 to Λ4).

Figure 11 .
Figure11.Shown above are the various posterior distributions on the location parameter for the PBH Gaussian bump, µPBH.The true values used to create the data are shown as solid vertical lines.For reference, the prior for this parameter was uniform between 1M⊙ and 80M⊙.For Λ0 and Λ1, the posteriors favor the true value.The other datasets are less informative, simply ruling out a high-mass component of the lensing population.The lower mass parts of the spectrum overlap with the stellar subpopulation (which dominate the catalog by a large margin), allowing the extra flexibility to be absorbed by this primary subpopulation.

Figure 12 .
Figure12.Top: the logarithm of the Bayes' factor between the two and three subpopulation class of models for each set of data Λi.Bottom: the logarithm of the ratio of the maximum likelihood values between the two and three subpopulation class of model for each set of data Λi.While the Bayes' factor statistic does not suggest strong evidence for a PBH subpopulation in any dataset, the ability to differentiate between the two models (two or three subpopulations) clearly improves as the number of PBH lenses actually contained in the data increases (from Λ4 to Λ0).

Figure 13 .
Figure13.Shown above are the posterior distributions for the minimum mass of the stellar (left) and SOBH (right) power law distribution(Mmin,Star and Mmin,SOBH).The vertical lines indicate the true value used to create each dataset, while the hatched distributions refer to inference performed assuming two subpopulations and unhatched distributions were inferred assuming three subpopulations.The prior for Mmin,Star is a uniform distribution between 0.01M⊙ and 1M⊙.The prior for Mmin,SOBH is a uniform distribution between 2M⊙ and 6M⊙.We see that we can robustly measure the minimum mass of the stellar power law distribution, regardless of the data set or population model employed.However, we cannot infer the minimum mass of the SOBH power law distribution, regardless of the data set or population model employed.The data sets and type of data (photometry only) do not contain enough information to place robust bounds on the minimum mass of the SOBH subpopulation that would be useful in measuring things like the NS-BH mass gap(Farr et al. 2011;Özel et al. 2010;Fryer & Kalogera 2001).

Figure
Figure14.Posterior distribution for the shape parameter of the stellar (left) and SOBH (right) power law distribution (bStar and bSOBH).The vertical lines indicate the true value used to create the data.The two types of analyses (assuming two or three subpopulations) are shown as a hatched and non-hatched distribution, respectively.The prior for these parameters is a uniform distribution between 0.1 and 10.For the stellar subpopulation, the posterior is approximately the prior distribution because of correlations with the location parameter, µStar.This correlation is shown in Fig.15.Considering the SOBH subpopulation, the posteriors are uninformative in the case of assuming three subpopulations.However, when only considering two subpopulations (thereby focusing on astrophysics as opposed to exotic physics), the posteriors contain significantly more information than the priors.The slight bias in the posteriors are connected to the same correlation shown in Fig.15, but for the SOBH subpopulation model, and the process of marginalization.

Figure 15 .
Figure15.The two dimensional joint posterior on bStar and µStar, as inferred using a two subpopulation model and the Λ0 data set.The green (dotted) contours and histograms refer to the posterior, while the black (solid) contours and histogram refer to the prior distribution for these parameters.The solid black lines represent the true values of the parameters used to generate the data.The strong correlation causes the one dimensional, marginalized posteriors on bStar and µStar to be very broad, despite there being plenty of information about certain linear combinations of these parameters in the data.
work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344.The document number is LLNL-JRNL-852673-DRAFT.This work was supported by the LLNL-LDRD Program under Project 22-ERD-037.M.F.H. acknowledges the support of a National Aeronautics and Space Administration FINESST grant under No. ASTRO20-0022.C.Y.L. and J.R.L. acknowledge support from the National Science Foundation under grant No. 1909641 and the Heising-Simons Foundation under grant No. 2022-3542.C.Y.L. also acknowledges support from NASA FINESST grant No. 80NSSC21K2043, a research grant from the H2H8 Foundation, and a Carnegie Fellowship.The authors would like to thank Alex Geringer-Sameth, Michael Schneider, James Barbieri and James Buchanan for helpful discussions.
c,j ) , d c,j ∼ p(d|ϑ ,j , t E,l ) , ϑ ,j ∼ p(ϑ|Λ) .(C21)pdet (t E ) calculated at each grid point, we can interpolate across the detection probability values, giving the typical detection efficiency curves published by many surveys.The results of this calculation are shown in Fig.18, compared to a published detection probability curve for a specific field from OGLE for reference.With this interpolated function, we can now finish the rest of the detection efficiency calculation repeated for every population model we consider, α = dd dθp(tr|d)p(d|θ)p(θ|Λ) , = dt E p det (t E )p(t E |Λ) .
Shown above are the event parameters tE and πE from a microlensing simulation of the Milky Way bulge produced by PopSyCLE, as published by

Table 1 .
Above are shown the true model parameters used by each subpopulation (stellar, SOBH and PBH), where the stellar and SOBH subpopulations are modeled as power law distributions and the PBH subpopulation is modeled as a Gaussian distribution.The values were picked to reflect realistic scenarios, either inspired by PopSyCLE simulation or from the literature.

Table 2 .
Pruett et al. (2022)s the number of detected sources in the final catalog of each population model used in this study, along with the corresponding relative abundance.The final column shows the correpsonding DM fraction fPBH given the assumptions ofPruett et al. (2022).We denote each population model as Λi for the i-th model, where the only difference in the model used to generate the data was the number of PBH sources.Note that the relative abundances shown here are for the intrinsic population model, not the fraction of detected sources.

Table 3 .
The prior distributions utilized in this study for the individual event model parameters.The ranges are listed in the same units as the parameters in the first column (e.g., we sample in log tE and the prior is log-uniform for tE, but the range is written here in days).

Table 4 .
The prior ranges used in the inference of the simulated population data.All priors are uniform, with boundaries set by the values in this table.The quantity N data represents the total number of events (detected or not) predicted by the true population model.This number ranges between 15630 and 16160.