Gaia May Detect Hundreds of Well-characterised Stellar Black Holes

Detection of black holes (BHs) with detached luminous companions (LCs) can be instrumental in connecting the BH properties with their progenitors' since the latter can be inferred from the observable properties of the LC. Past studies showed the promise of Gaia astrometry in detecting BH-LC binaries. We build upon these studies by: 1) initialising the zero-age binary properties based on realistic, metallicity-dependent star-formation history in the Milky Way (MW), 2) evolving these binaries to current epoch to generate realistic MW populations of BH-LC binaries, 3) distributing these binaries in the MW preserving the complex age-metallicity-Galactic position correlations, 4) accounting for extinction and reddening using three-dimensional dust maps, 5) examining the extended Gaia mission's ability to resolve BH-LC binaries. We restrict ourselves to detached BH-LC binaries with orbital period<10 yr such that Gaia can observe at least one full orbit. We find: 1) the extended Gaia mission can astrometrically resolve 30-300 detached BH-LC binaries depending on our assumptions of supernova physics and astrometric detection threshold; 2) Gaia's astrometry alone can indicate BH candidates for 10-100 BH-LC binaries by constraining the dark primary mass>3 Msun; 3) distributions of observables including orbital periods, eccentricities, and component masses are sensitive to the adopted binary evolution model, hence can directly inform binary evolution models. Finally, we comment on the potential to further characterise these BH binaries through radial velocity measurements and observation of X-ray counterparts.


INTRODUCTION
The recent discoveries of merging binary black holes (BHs) by the LIGO-Virgo and Kagra observatories have reignited the interest in understanding the astrophysical origins of BH binaries in short-period orbits (e.g., Abbott et al. 2016aAbbott et al. ,b, 2019Abbott et al. 2021a,b). A major hurdle in understanding the astrophysical implications of these detections, as well as modeling realistic populations of binary BHs can be attributed to uncertainties in how massive stars evolve and form compact objects (e.g., Woosley et al. 2002;Fryer et al. 2012;Sukhbold et al. 2016;Woosley 2017;Sukhbold et al. 2018;Pejcha 2020). chirag.chawla@tifr.res.in souravchatterjee.tifr@gmail.com It is expected that there are ∼ 10 7 − 10 9 stellarmass BHs in the Milky Way (e.g., Brown & Bethe 1994;Timmes et al. 1996;Samland 1998;Olejak, A. et al. 2020). However, discovering them using traditional methods such as X-ray or radio observations is notoriously difficult resulting in only ∼ 60 detections to date (e.g., Remillard & McClintock 2006). Only a small fraction of BH binaries are expected to be actively accreting at a detectable rate at any given time due to the stringent requirements on the orbital and stellar properties of accreting systems and their typically low duty cycles (Fragos et al. 2009;Fabbiano 2012;Gallo et al. 2014;Corral-Santana et al. 2016;Tetarenko et al. 2016). Furthermore, the observed population of BHs in low mass X-ray binaries may be biased toward lower masses due to observational selection effects (Jonker et al. 2021).
The population of BHs detected via gravitational waves (GWs) emitted from the inspiral and merger of binary BHs also have strong selection biases favoring distant and high-mass objects (e.g., Fishbach & Holz 2017). Thus, the historically standard methods for BH detections likely miss the bulk of the population of BH binaries in the Milky Way. Even when discovered, GW detections do not directly constrain age and metallicity of the progenitors of the merging BHs. Instead, such constraints come indirectly from population modeling of various astrophysical formation channels (e.g., Abbott et al. 2021a,b;Bavera et al. 2021;Boco et al. 2021). On the other hand, detailed modeling of complex accretion physics is needed to characterise the BH and the companion properties in X-ray binaries since the properties of the accretion disks are the observables in this case (e.g., Frank et al. 2002;Davis et al. 2005;Kreidberg et al. 2012).
The majority of BH binaries are expected to have stellar companions in orbits too wide for mass transfer via Roche-lobe overflow (RLOF) (e.g., Breivik et al. 2017). Detection of a detached BH binary with a luminous companion (LC) is interesting since such a source can effectively remove several limitations of the aforementioned methods. In a detached BH-LC binary, the properties, such as the age and metallicity of the BH progenitor can be assumed if the same properties can be measured for the LC. For example, knowledge of the brightness and distance of the LC can allow strong constraints for the luminosity and metallicity using well-understood stellar evolution models. If the BH-LC binary is primordial, as expected for most BH binaries in the Milky Way field, the metallicity and age of the BH progenitor is very likely the same as its companion. Furthermore, especially if the LC is a main-sequence (MS) star, the knowledge of luminosity and color can provide strong constraints on the mass of the star and as a result, the mass of the dark object can also be constrained (e.g., Andrews et al. 2019;Shikauchi et al. 2020). Interestingly, a fraction of detached BH-LC binaries in the field may also have been created dynamically inside star clusters which then get ejected from the host cluster (e.g., Kremer et al. 2018). Since star clusters consist of an effectively co-eval aggregate of stars, the same assumptions should be valid for dynamically produced systems as well to a large degree. 1 Detached BH-LC binaries are almost impossible to detect via traditional methods such as radio, GW, or X-ray emissions, with the only exception being wind-fed X-ray binaries. Nevertheless, several BH candidates in detached binary systems with a LC have already been detected by multi-epoch spectroscopic and photometric campaigns via measurement of the orbital motion of the LC (e.g., Giesers et al. 2018;Giesers et al. 2019; Thompson et al. 2019;Jayasinghe et al. 2021). Systematic surveys for radial velocities of all nearby stars by APOGEE or SDSS-V may provide more discoveries in future (e.g., Zasowski et al. 2017;Kollmeier et al. 2019).
Recently, several groups have proposed that large numbers of detached BH-LC binaries may be detected by resolving the orbital motion of LCs around dark objects using astrometry by Gaia (Gould & Salim 2002;Barstow et al. 2014;Mashian & Loeb 2017;Breivik et al. 2017;Yalinewich et al. 2018;Yamaguchi et al. 2018;Breivik et al. 2019;Wiktorowicz et al. 2020). Since Gaia will provide position and parallax measurements to ∼ µas precision, it will be relatively straightforward to estimate the LC's luminosity and temperature from magnitude and color without the need for additional followup observations. As a result, the age, metallicity, and mass of the LCs (especially if the LC is a main-sequence star) can be constrained using well-understood stellar modeling (e.g., Anders et al. 2019;Howes et al. 2019). In addition, a population of BHs detected this way, arguably will constitute the least biased detected to date, since all selection effects in this case depend primarily on the properties of the LC and not directly on the properties of the BH.
All of the aforementioned studies show that the astrometric motion of a large number of LCs in orbit around unseen dark objects should be resolved using Gaia. These studies also illustrate that while the basic idea is robust, the actual yield and the properties of the detected BH-LC binaries can vary widely depending on the assumptions for stellar evolution and binary interactions (e.g., Breivik et al. 2017), and how carefully one considers observability and selection effects. Hence, it is crucial to carefully consider the details related to the synthesis of the population of BH-LC binaries as well as their observation by Gaia. In this study we update our earlier works presented in Breivik et al. (2017), Breivik et al. (2019), and Andrews et al. (2019) by making several improvements in our binary population synthesis and observational considerations. In contrast to these previous works, we now include realistic stellar distributions in the Milky Way, and location-dependent stellar ages and metallicities based on the Ananke Framework of the Latte Suite of the FIRE-2 Simulations (Wetzel et al. 2016;Hopkins et al. 2018;Sanderson et al. 2020). Furthermore, we consider observational selection effects more accurately by taking into account the number of planned Gaia transits for each system and three dimensional extinction, both of which depend on the Galactic positions of the BH-LCs.
The paper is organised as follows. In section 2 we detail our numerical setup for BH-LC population synthesis and construction of synthetic Milky Way models. In section 3 we describe how we determine whether a BH-LC binary would be resolvable by Gaia. Here we take into account astrometric detectability, interstellar extinction and reddening, and possibility of mass constraints of the dark objects using astrometry alone. In section 4 we present our key results for the population of BH-LCs resolvable by Gaia's astrometry. In section 5 we explore promising avenues for followup studies. Finally, we summarize our results and conclude in section 6.

NUMERICAL SETUP
To study the properties of the BH-LC population Gaia may observe, we generate representative present day BH-LC binary populations using COSMIC (Breivik et al. 2020). COSMIC is a Python-based binary population synthesis suite that employs modified versions of the single and binary star evolution codes SSE/BSE (Hurley et al. 2000;Hurley et al. 2002). For a detailed discussion of the modifications to assumptions for binary interaction physics beyond the standard BSE release, see Breivik et al. (2020); Rodriguez et al. (2021). We detail the process to generate a representative Galactic population of BH-LC binaries in the following subsections.

Initialising the binary population
COSMIC generates binary populations from Zero Age Main Sequence (ZAMS) by assigning each system with an initial age, metallicity, primary mass, secondary mass, orbital period, and eccentricity. We assign ages and metallicities for each binary based on the distribution of star particles in the final snapshot of the simulated galaxy m12i in the Latte suite of the Feedback In Realistic Environments 2 (FIRE-2) 2 simulation suite (Wetzel et al. 2016;Hopkins et al. 2018). We use the ananke framework to assign three-dimensional Galactic positions for each binary by resampling from the smoothed density distribution of the star particles, as described in subsection 2.4 of Sanderson et al. (2020). These new star formation history assumptions are a major upgrade to our earlier work (Breivik et al. 2017) which assumed constant star formation rate and discrete metallicity values for the thin and thick disks. For a detailed discussion of the resemblance of the adopted star-formation history to that of the Milky Way see Sanderson et al. (2020). Since the stellar evolution fits employed in BSE and thus COSMIC are valid for metallicities in the range between log(Z/Z ) = −2.3 and 0.2, we limit the metallicity of our binaries to fall within this range and assign any metallicities outside of the range to the limiting values.
We assign ZAMS orbital parameters of the binary population by sampling primary masses, secondary masses, orbital periods, and eccentricities from observationally motivated probability distribution functions. In particular, the primary masses are drawn from the Kroupa (2001) initial stellar mass function (IMF) while the secondary companions are assigned masses following a uniform distribution of mass ratios with a range in mass from 0.08 M to the primary mass (Mazeh et al. 1992;Goldberg & Mazeh 1994). Initial eccentricities and orbital periods are drawn by first sampling eccentricities from a thermal distribution (Heggie 1975) then sampling semimajor axes that are uniform in log space up to 10 5 R (Han 1998). We reject any samples that produce binaries where one of the component stars fills more than half of its Roche radius. Finally, we assume that the initial binary fraction is 0.5 which places two of every three stars formed into a binary system.

Binary stellar evolution
We use COSMIC to simulate binary evolution from ZAMS through to the present day with ages and metallicities assigned by the star particles in galaxy m12i. The present-day orbital characteristics of the BH-LC population depend strongly on the assumptions for the outcomes of RLOF mass transfer from the BH progenitor as well as natal kicks which can be imparted to the BH during its formation. We parameterize the stability of mass transfer using critical mass ratios defined in Belczynski et al. (2008) which delineate whether mass transfer remains dynamically stable or enters a common envelope (CE) which dramatically shrinks the orbit of the BH-LC progenitor. For stable mass transfer our treatment follows the treatment described in Hurley et al. (2002). We assume that the donor loses mass with a rate that steeply increases with the amount that the donor radius overfills its Roche radius. We assume that the accretor is able to accept mass at 10 times its thermally limited rate (a star's mass divided by its thermal time) during the main sequence, Hertzsprung gap, and core helium burning phase and at an unlimited rate during any giant-like phases. All mass that is not accreted is assumed to leave the binary with the specific angular momentum of the accretor.
COSMIC employs the αλ prescription for CE evolution where λ is a parameter which varies the binding energy of the donor's envelope based on its structure and α defines how efficiently orbital energy is used to eject the envelope (Paczynski 1976;Livio & Soker 1984;Tout et al. 1997). We set α = 1 which assumes that the initial orbital energy is converted into ejecting the envelope with 100% efficiency. We set the binding energy parameter, λ, as defined in the Appendix of Claeys et al. (2014). We assume the default choice in Claeys et al. (2014) such that no recombination energy from ionization in the donor envelope is considered in the envelope ejection (recombination contributes only a few percent to the overall energy budget for massive stars; Kruckow et al. 2016). All stable RLOF and CE interactions are assumed to perfectly circularize the binary.
The distribution of birth kicks BHs receive is quite uncertain because of the limited numbers of detected stellar BHs (e.g., Repetto et al. 2012;Repetto & Nelemans 2015;Atri et al. 2019) and the complexity of modeling SN explosions from first principles (e.g., Woosley et al. 2002;Sukhbold et al. 2018;Pejcha 2020). For this study, we consider the two most widely used mechanisms for the formation of the compact objects: the "rapid" and "delayed" mechanisms presented in Fryer et al. (2012). The two mechanisms differ in the growth time of instabilities in the convective region outside of the protocompact object which re-initiate the post core-bounce shock to complete the explosion. The primary difference between BH populations formed with the rapid and delayed mechanisms is the existence (rapid) or lack of (delayed) a mass gap separating neutron stars (NS) and BHs between ∼ 3 M and ∼ 5 M . Each prescription includes a recipe for the partial fallback (f fb ) of the stellar envelope. In addition to its effect on the BH masses, SN fallback reduces the strength of the natal kicks by a factor of 1 − f fb from the kick strength drawn randomly from a Maxwellian distribution with σ = 265km/s (Hobbs et al. 2005) and distributed isotropically. Throughout this work we test both models and call our simulated populations rapid and delayed after the name of the SN prescriptions used to create them.

Convergence of present-day binary properties
COSMIC is designed to adaptively choose the size of a simulated population based on a convergence criteria that tracks how the shape of parameter distributions for binaries, such as mass, orbital period (P orb ), and eccentricity (Ecc) change as the population grows. This is done by iteratively simulating populations of size N sim and adding each newly simulated population to the total population on each iteration. Histograms of the binary parameter data are also made at each iteration and compared bin by bin using a criterion inspired by matched filtering techniques (e.g. Eq. 6 of Chatziioannou et al. 2017) and defined as where P k,i represents the height of bin k on the ith iteration (Breivik et al. 2020). The match criteria approaches unity as the shapes of each histogram converge. We use Knuth's Rule to determine the binwidths of each binary parameter histogram for the i + 1 simulation iteration (Knuth 2019). At each iteration we record the total amount of ZAMS mass drawn in single and binary stars to scale the converged population up to a Galactic population. We continue to simulate binaries until match > 1 − 10 −5 for the BH mass (M BH ), LC mass (M LC ), P orb , and Ecc.

Creating synthetic Milky-Way population
We generate synthetic Milky-Way populations by sampling with replacement from the converged population of simulated binaries. To determine the number of BH-LC binaries in the Milky Way at present, we scale the number of BH-LC binaries in the converged population as where M m12i is the total amount of mass formed in galaxy m12i, and M sim is the mass of single and binary stars drawn to produce the simulated BH-LC population. We assign each of the N BH−LC,MW BH-LCs a complete set of stellar and orbital parameters including primary mass, secondary mass, age, metallicity, eccentricity, P orb , and luminosity from the converged population according to the unique binary id. The sampled binaries are then assigned a Galactocentric position by minimizing the difference in age and metallicity between the binary and star particles in the m12i galaxy. Since multiple BH-LCs can be assigned to a single star particle, the BH-LCs are distributed in a sphere centered around the star particle, where the radius of the sphere is determined using an Epanechnikov kernel with kernel size inversely proportional to the local density (Sanderson et al. 2020). Each BH-LC is also assigned an orientation with respect to the line-of-sight from Earth by sampling the inclination (i) uniformly in cos(i), and argument of periapsis (ω) and longitude of ascending node (Ω) uniformly between 0 and 2π.
The number of BH-LCs resolvable by Gaia (described in section 3) varies depending on the random assignment of the parameters described above. To incorporate the effects of this variance, we repeat the process described above to generate 200 realisations of the Milky-Way population for each of the rapid and delayed models, producing 400 Milky Way realizations in total.

DETECTABILITY OF BH-LC BINARIES USING Gaia
In this section we describe how we consider several important aspects that determine detectability of a BH-LC by Gaia, including the size of the sky-projected orbit of the LC with respect to the astrometric precision of Gaia, the number of Gaia observations at each sky position, and the effects of extinction and reddening on the target LC population.

Reddening
Extinction and reddening from interstellar dust have a significant effect on the observation of stars distributed in the Galactic disk. To account for the extinction due to interstellar dust we apply an extinction correction to all LCs in our population. We determine the extinction corrections using a complete three dimensional positiondependent dust map which combines three different dust maps (Drimmel et al. 2003;Marshall et al. 2006;Green et al. 2019) covering different regions of the Milky Way as implemented in mwdust (Bovy et al. 2016). 3

Resolving the sky-projected LC orbits by Gaia
For a BH-LC binary orbit to be resolved by Gaia we require that the extinction-corrected LC magnitude be brighter than Gaia's limit of G < 20. Gaia can observe stars fainter than G = 21, but in practice these are too faint to detect astrometric motion due to binarity. We further require that the projected angular size of the binary, α ≥ σ ξ ("optimistic") and α ≥ 3 × σ ξ ("pessimistic"), where σ ξ is Gaia's single position pointing error (Lindegren et al. 2018). In addition to the criteria for resolving the LC's orbital motion, we impose the condition P orb ≤ 10 yr to ensure that each binary completes a full orbit within the duration of Gaia's extended mission.
Since every BH-LC has a unique position and orientation in the Galaxy, we determine the projected angular size of the binary on the sky using the Thiele-Innes con-3 Green et al. (2019) includes data from Gaia DR2, Pan-STARRS 1, and 2MASS. Marshall et al. (2006) is based on only the 2MASS data. Drimmel et al. (2003) takes the data from COBE/DIRBE NIR observations.

stants:
A = a LC (cos ω cos Ω − sin ω sin Ω cos i) (3) B = a LC (cos ω sin Ω + sin ω cos Ω cos i) (4) F = a LC (− sin ω cos Ω − cos ω sin Ω cos i) (5) G = a LC (− sin ω sin Ω + cos ω cos Ω cos i). (6) In the equations above, a LC is the semimajor axis of the LC orbit defined as where a is the binary semimajor axis and M LC (M BH ) is the mass of the LC (BH). The projected cartesian position of the LC is then defined by where, X = cos E −Ecc, Y = √ 1 − Ecc 2 sin E, and E is the eccentric anomaly. Finally, we define the projected size of the LC orbit, a project,LC , on the sky in terms of the semimajor axis of the ellipse defined by Equation 8 and 9, as where D is the distance to the binary.

Astrometric mass measurement
The most promising way of discovering BH-LC binary candidates is through astrometric mass estimation of the BH, which does not suffer the inclination degeneracy present in radial velocity (RV) measurements (e.g., see discussion in Andrews et al. 2019). If the mass of the unseen companion to an LC is confirmed to be larger than the minimum mass of a BH then the unseen companion is likely a BH candidate. We showed in Andrews et al. (2019) that the astrometric mass function error can be approximated as where M tot = M BH +M LC and N is the number of times a binary is observed by Gaia.
We determine the position-dependent N for each BH-LC in the following way. We divide the sky into 164, 838 grids of equal angular area on the sky-plane. We take 360 equal intervals in declination, each 0.5 • apart between −90 • to 90 • . We divide RA for each of the above dec into n RA bins where n RA is the nearest integer ≥ 720 × cos(dec). Our grid points ensure that the angular sky-projected area between any two adjacent grid points in RA and dec is smaller than the field-of-view (0.72 • × 0.69 • ) for each of Gaia's telescopes (e.g., Lindegren et al. 2012). Using the Gaia Observation Forecast Tool (GOST) 4 we calculate and store the number of times Gaia would observe each grid point in five years (N 5 ; between 2014-19). We find the value of N 5 for any BH-LC binary in our Milky-Way realisations, simply by looking up this number corresponding to the nearest grid point from the RA and dec of that particular model binary. Finally, we set N = 2N 5 to account for a 10 yr observation duration.
If the LC mass is estimated using the astrometrically derived distance, then this mass measurement can be combined with the mass function measurement described above to determine the BH mass measurement accuracy ∆M BH . We use error propagation based on Equation 11 to estimate ∆M BH 5 : To determine the number of potential BH candidates, we assume that all LCs have mass measurements to 10% accuracy. We select BH candidates from our simulated population if M BH > M LC and the BH has a mass measurement such that M BH − ∆M BH ≥ 3 M . We impose the former condition to eliminate tricky scenarios where confirming a BH candidate is difficult because of the possibility that the LC may dominate the total light not because its companion is a dark remnant, but because it is a lower-mass and hence fainter star.

RESULTS
In this section we describe the BH-LC populations, both for the entire Milky Way and the subset resolvable by Gaia. We include both the overall detection rate as well as the stellar and orbital characteristics of the BH-LC binaries.

Present-day properties of simulated BH-LC binaries
In order to gain insight for the underlying population of BH-LCs in the Milky Way at present, we first focus on Orange, blue, and green histograms denote distributions for all BH-LC binaries, those that went through at least one CE evolution, and those that never went through a CE evolution. The short-period peak is dominated by BH-LC systems created via at least one CE episode, whereas, the long-period peak is created primarily by BH-LC binaries that never went through a CE episode. In our analysis we consider BH-LC binaries with P orb ≤ 10 yr (the vertical black dashed line to guide the eye), such that a 10-yr Gaia mission can resolve the full orbit. Interestingly, 10-yr Gaia mission potentially can detect BH-LC binaries created via the CE channel as well as those that never went through a CE evolution. (The equivalent figure for our delayed model is presented in the Appendix.) the observable properties of the BH-LC populations created in our models without imposing any observational selection effects. Figure 1 shows the overall distribution for the orbital periods (P orb ) of all BH-LC binaries in the Milky Way at present from the rapid model (section 2). The distribution of P orb for BH-LC binaries in the delayed model is very similar (shown in the appendix). The P orb distributions are clearly bimodal for both BH-MS and BH-PMS binaries for both of the adopted SN models. In both cases, the short-period peak is dominated by binaries created via at least one CE episode, whereas, the longperiod peak is created primarily with BH-LC binaries that never undergo a CE. The location of the trough for all cases is between P orb ∼ 1-10 yr. The majority (∼ 96%) of BH-LC binaries with P orb ≤ 10 yr have gone through at least one CE evolution independent of the adopted SN mechanism. Interestingly, the extended 10 yr mission of Gaia allows sampling of the separation between the CE and non-CE peaks. While outside the scope of this work, future studies could investigate how the observed shape and width of the trough can constrain uncertain CE physics.
One of the most prominent differences between the two SN mechanisms manifests in the M BH distribution; while the rapid prescription for core-collapse SN does not allow production of BHs in the mass range 3 to 5 M , (often called the 'mass-gap' between neutron stars (NSs) and BHs), the delayed prescription predicts no mass gap. Figure 2 shows the M BH distributions for BH-LC binaries at present day for the rapid and delayed models. The distributions are clearly different in the range between 3 and 5 M , as expected. All BHs with M BH between 3 and 5 M in the rapid model originate from accretion-induced collapse (AIC) of NSs. In contrast, in the delayed model, the AIC BHs contribute at ∼ 15 (55)% for BHs MS (PMS) companions and the rest come from core-collapse SN. The M BH distribution continuous all the way down to M BH /M = 3, our adopted mass that separates NSs from BHs in the delayed model, whereas, for the rapid model, there is a clear separation between the mass-gap BHs and the rest.
Interestingly, all BH-LC binaries with M BH /M 7 (10) in the rapid (delayed) model have P orb ≤ 10 yr leading to identical populations for the BH-LC binaries with lower-mass BHs regardless of P orb restrictions. This is because lower-mass BHs receive larger natal kicks in both SN prescriptions we have considered (Fryer et al. 2012) while for higher-mass BHs, natal kick magnitudes are reduced depending on the amount of mass that falls back onto the proto-compact object (e.g., Belczynski et al. 2008). The large natal kicks for low-mass BHs break their progenitor binaries with orbits P orb /yr 10, while, the more massive counterparts in wide orbits can remain intact. If Gaia is able to characterize BHs in binaries with M BH ∼ 10 M and P orb ≤ 10 yr, strong constraints can be placed on the strength of natal kicks for these systems.
A few other features, originating from a combination of the adopted initial stellar mass function, metallicity and age distributions, and the complex binary star evo- The blue and red curves represent MS and PMS companions, respectively. The faded and bright lines denote all BH-LC binaries without any constraints on P orb and those satisfying P orb ≤ 10 yr. The major difference between the two distributions come in the range MBH/M = 3-5. In case of the rapid model, all BHs within MBH/M = 3-5 come from AIC of NSs; we assume that the limiting mass for NSs is 3 M (section 2). Whereas, for delayed, BHs formed via core-collapse SN can have masses down to 3 M . As a consequence, the delayed model consists of a relatively higher fraction [25 (71) lution of the initial population of binaries, are noticeable and different in the two models. For example, in the rapid model a peak in the range M BH /M ≈ 7 to 9 come from the evolution of metal-rich (Z/Z ≥ 0.75) systems. Another less prominent and broader peak between M BH /M ≈ 13 and 19 consists of systems that are relatively metal poor (Z/Z ≤ 0.5). In contrast, in the delayed model the distribution peaks near M BH /M = 3 and again near M BH /M ≈ 15. Figure 3 shows the distributions of BH-LC orbital eccentricities for the rapid and delayed models. We find . In all cases there is a prominent peak for near-zero eccentricity, especially for BH-LC binaries with P orb /yr ≤ 10. This is because a large fraction of short-period BH-LC binaries go through CE evolution and are significantly affected by tides. The BH-MS (BH-PMS) binaries in the rapid and delayed models contain 0.8% (0.3%) and 7% (1%) P orb /yr ≤ 10 systems with Ecc ≥ 0.1, respectively. The delayed model contains higher fraction of eccentric BH-LC binaries compared to the rapid model. This can be attributed to typically larger natal kicks in case of delayed. The fraction of near-circular orbits is less in case of BH-MS compared to BH-PMS in both models since tidal circularisation is more effective in BH-PMS binaries than in BH-MS binaries.
that about 2% (15%) and 27% (76%) of BH-MS (BH-PMS) binaries with P orb /yr ≤ 10 have near-circular (Ecc ≤ 10 −2 ) orbits for the rapid and delayed models, respectively. The majority of the BH-LC systems with P orb ≤ 10 yr go through at least one CE evolution stage with the BH-progenitor's envelope which circularizes the binary before BH formation ( Figure 1). Moreover, even the small fraction of the short-period BH-LC binaries that never went through a CE phase can be significantly affected by tidal circularisation. In spite of this, the high fraction of eccentric BH-LC binaries primarily result from natal kicks during BH formation. The differences in the explosion mechanism between the rapid and delayed SN prescriptions lead to typically higher natal kicks in the delayed models. As a result, the fractions of short-period BH-MS (BH-PMS) binaries with significant Ecc ≥ 0.1 are 0.8% (0.3%) and 7% (1%) in our rapid and delayed models, respectively. The connection between larger eccentricities and stronger natal kicks is a robust prediction that is independent of our chosen SN prescription. For both SN mechanisms, BH-PMS binaries usually have a smaller fraction of eccentric systems compared to BH-MS binaries, since the former have larger LC envelopes which are affected more strongly by tides after BH formation. Figure 4 shows the distributions for age, metallicity, mass, and luminosity of the LCs for all BH-LC binaries and those satisfying P orb ≤ 10 yr for the rapid model. The corresponding figure for the delayed model is included in the Appendix. In contrast to the eccentricity and M BH distributions, the LC properties primarily depend on the stellar IMF, binary stellar evolution physics, and the star formation and metallicity evolution in the Milky Way; they do not strongly depend on the details of the assumptions of core-collapse SN physics. As a result, the distributions of LC properties are very similar in the rapid and delayed models. We find that the ages of the BH-LC progenitors spans several Gyr covering the full range of the star particle ages in the galaxy model m12i. As a consequence, the BH-LC progenitors also have a wide range in metallicities. The pile up at the last metallicity bins near log(Z/Z ) = −2.3 and 0.2 is simply because the star particles in the model galaxy m12i have a wider range in Z relative to the allowed Z range in COSMIC which is based on the single star fits of BSE (Hurley et al. 2000, see also the discussion in subsection 2.1). Nevertheless, the wide range in metallicities of our predicted present-day BH-LC binaries in the Milky Way indicates that if a sizeable population of BHs are detected, metallicity-dependent constraints on the properties of BHs may be obtained. Note that our earlier work made a simplifying assumption that used two fixed metallicity values corresponding to the thin and thick disks of the Milky Way (Breivik et al. 2017). Using more realistic metallicities that depend on the star-formation history in the Milky Way is one of the key improvements in the present study.
Almost all low-mass ( 0.1 M ) MS companions to the BH binaries show P orb ≤ 10 yr. This is likely related to the CE evolution. Orbits with lower LC masses (on an average) contain a smaller reservoir of orbital energy to eject the envelope. As a result, BH's with lower-mass LCs tend to have shorter post-CE orbits. Although, not as prominent, the same trend is observed for BH binaries with PMS companions; below M LC ∼ 1 M , all BH-PMS binaries have P orb ≤ 10 yr.
A small ( 0.1%) fraction of BH-MS binaries consist of MS stars in excess of 40 M . We find that in some cases, an initially high-mass (M ≥ 15 M ) MS star accretes from the BH progenitor through stable mass transfer via RLOF and grow. These stars live longer than normal due to rejuvenation (e.g., Hurley et al. 2002). Most of these systems are also very young ( 9 Myr). As expected, the mass range for the PMS companions is narrower compared that for the MS companions because the PMS phase is a shorter-lived stage of evolution. In addition, lower-mass LCs must be sufficiently old to advance off the MS leading to a higher relative number of BH-MS binaries with ages larger than 6 Gyr relative to BH-PMS binaries.

BH-LC binaries in the Milky Way resolvable by Gaia
In this section we introduce the realistic Milky-Way populations each consistent with the complex and interdependent age-metallicity-Galactic position observed in the Galaxy (subsection 2.4). These realisations are different from each other in the random orientations, the distance, and the Galactic locations of the BH-LC binaries. As a result there are differences in the number of detections by Gaia in each realisation. These differences provide an estimate of the level of statistical fluctuations in the expected BH-LC numbers resolvable by Gaia. Note-BH-LC numbers in the Milky Way predicted in our models. The numbers and errors denote the median and the spread between the 10th and 90th percentiles across the Milky-Way realisations. "Intrinsic" denotes the present-day population of BH-LC binaries in the Milky Way. "Resolvable" denotes the subset of BH-LC binaries resolvable by Gaia's astrometry over a 10 yr mission (subsection 3.2). Any detached BH-LC binary with P orb /yr ≤ 10 and resolved LC orbit is denoted as "Resolved." "Av-corrected" denotes the expected number of resolved BH-LC binaries after correcting for extinction and reddening from interstellar dust (subsection 3.1). A BH-LC is a "BH candidate" if Gaia's astrometry can constrain the BH's mass MBH − ∆MBH ≥ 3 M and MBH ≥ MLC within the Av-corrected resolved population (subsection 3.3). Figure 5 shows the cumulative number of BH-LC binaries for which Gaia can resolve the motion of the LC in orbit around the BH as a function of Gaia's G magnitude. The estimated total number of resolvable BH-LC binaries for the rapid model is between ≈ 240 (pessimistic) and 470 (optimistic). In case of the delayed model, the expected yield is lower, between about 60 (pessimistic) and 185 (optimistic) BH-LCs. The delayed SN prescription typically produces lower-mass BHs and higher natal kick magnitudes compared to the rapid SN prescription. As a result, a higher fraction of the progenitors of BH-LC binaries disrupt in our delayed model.
We now embark on an investigation of the relative importance of various formation channels of the resolvable BH-LC binaries. Identification of the relative importance of different channels can inform which binary evolution physics of any particular channel matters most in estimating the model populations. In addition, different formation channels can create BH-LC binaries with significantly different properties including P orb and Ecc (Figure 1, 3). Moreover, future BH-LC detections from Gaia can inform the relative abundances of various types of BH-LC binaries and shed light on the uncertain aspects of binary stellar evolution including the CE, and BH formation physics. Figure 6 and 7 show the detailed evolutionary pathways for the Gaia-resolvable BH-LC binaries and their relative importance. Below we mention the key findings. In case of the rapid model, above 80% of resolvable BH-LC binaries are expected to contain a MS star as a companion. Overall, ∼ 50% of all resolvable BH-LC systems come via CE evolution. Note that, the contribution from CE evolution in the intrinsic BH-LC population with P orb /yr ≤ 10 is much higher (∼ 96%; subsection 4.1). However, Gaia's astrometry can resolve the typically large orbits of the non-CE BH-LC binaries to larger distances. This selection bias increases the fraction of non-CE BH-LCs in the resolvable population. Looking at the resolvable BH-MS and BH-PMS populations separately, about 40% (95%) of the resolvable BH-MS (BH-PMS) binaries are expected to have undergone at least one CE phase.
Only about 2.4% of the resolvable BH-LCs contain BHs formed via AIC and all of them have a PMS companion. This is significant in two ways. First, if the mass-gap is real, it will be easily discovered from the resolvable BH-LC population without any significant contamination from BHs created via AIC of NSs. Second, the resolvable BHs in the mass gap are expected to have PMS companions only.
We find that in about 2% of the resolvable BH-LCs the BH progenitor is the donor while crossing the Hertzsprung Gap (HG) during one of the CE episodes. Whether a binary can survive a phase of CE evolution with a HG donor is uncertain(e.g., Ivanova & Taam 2004;Belczynski et al. 2008). In spite of the uncertainty of this channel, we do not exclude these systems from our analysis for the sake of completeness. Interestingly, if we do exclude systems that underwent a CE initiated by a HG donor, then the fraction of resolvable BH-LCs containing mass-gap BHs reduces to only 0.4%. In the delayed model, the relative importance of the formation channels are somewhat different. MS companions still dominate (69%) the population of resolvable BH-LCs. About 70% of the resolvable BH-LCs form via at least one CE episode and about 13% of the BH-LCs contain BHs created via AIC. The extension of Gaia's mission from 5 to 10 years allows Gaia to target BH-LCs up to P orb = 10 yr and still resolve the full orbit. According to our rapid (delayed) model, the total expected yield of resolved BH-LC binaries increases by 34% (16%) simply because of the mission extension and this increase in overall yield almost entirely comes from the increase in detectability of longer-period BH-LC binaries that do not go through CE evolution. We find that Gaia's 10 yr mission would identify significant fractions (50% and 30% for rapid and delayed models) of resolvable BH-LC binaries that do not come from the CE evolution channel. If such yields are realised, it may be possible to distinguish resolved BH-LC systems that form via CE evolution from those that do not simply using the P orb distribution.

Key Properties of resolvable BH-LC binaries
In this section we highlight some key properties of the resolvable BH-LC binaries. Comparisons between the distribution of properties presented here with those presented in subsection 4.1 sheds light on the selection biases of Gaia. Figure 8 and 9 show the corner plots for all resolvable binaries across all Milky Way realisations from our rapid and delayed models. Note that different Milky Way realisations can sample the same binaries multiple times (subsection 2.4). To avoid clutter, we plot resolvable binaries in the scatter plots only once and ignore the repeated draws. However, we create the histograms taking into account the repeated draws.
The difference between the rapid and delayed models manifests most prominently for BHs with 3 ≤ M BH /M ≤ 5. In the rapid (delayed) model, about 2.4% (65%) of the resolvable BH-LCs contain BHs in this mass range. If we discard the systems that have undergone CE involving a HG donor, the contrast is even starker. For resolved BH-LCs with BHs in this mass range, while in the rapid model we do not find any BH-MS binaries, in the delayed model BH-MS and BH-PMS contribute almost equally. Note that, even in the delayed model AIC BHs contribute only at about 13% within this M BH range and the rest come from core-collapse SN. We already highlighted similar trends in Figure 2 in the context of the intrinsic population of BH-LCs, however, it is interesting that this difference is prominent in the resolvable population as well. Detection of BHs with 3 ≤ M BH /M ≤ 5 with both MS and PMS companions in detached configuration would clearly indicate that the so-called mass-gap from core-collapse SN may not be real. Indeed, the recent discovery of BH-LC candidates with BH masses in this range (Thompson et al. 2019;Jayasinghe et al. 2021) and microlensing analyses (Wyrzykowski & Mandel 2020) provide observational evidence against the existence of a mass gap.
It is interesting that the shapes of the M BH distributions in our models do not change significantly between the intrinsic and resolvable populations of BH-LC binaries in the Milky Way independent of the adopted SN prescription. This indicates that the ability for Gaia to resolve a BH-LC binary do not strongly depend on M BH . Indeed, we find that the ratio between the total numbers of BH-LC binaries in the Milky Way and those that are resolvable by Gaia, do not strongly correlate with M BH (Figure 10). This is because Gaia's astrometric selection biases do not depend directly on the BH, rather they depend on the LC and orbital properties, such as G and P orb . As a result, the population of BHs resolv-  able by Gaia is expected to exhibit M BH -distribution that is very close to the intrinsic one in the Milky Way. The other potentially observable difference between the BH-LCs in the rapid and delayed models stems from the differences in the fallback-dependent natal kicks in the two respective SN mechanisms. The majority of the resolvable systems-about 90% (95%) in the rapid (delayed) model-go through dissipative binary evolution phases such as CE evolution and RLOF. The final eccentricities are almost entirely dependent on the natal kicks the BHs receive modulated later by tidal cirularisation or RLOF. Because of this, in both the rapid and delayed models, there are prominent peaks near circular orbits. We find that about 87% (33%) of the resolvable BH-LC binaries have Ecc ≤ 0.1 in the rapid (delayed) model. Similar to the intrinsic population, the resolvable BH-MS binaries in the delayed model show a higher fraction of eccentric systems com- pared to those in the rapid model. We find that the resolvable systems with Ecc ≥ 0.5 in both rapid and delayed models form predominantly from progenitor binaries that never underwent a CE evolution. Since, this group also consists of relatively wider orbits, they are also less affected by tides after BH formation. In both models, the fraction of eccentric BH-MS binaries is higher than that of eccentric BH-PMS binaries in the resolvable population since the BH-PMS binaries are more strongly affected by tides. Thus future observations of the eccentricity distribution of BH-LC binaries can put constraints on the natal kick strengths for stellarmass BHs. Such constraints would be complimentary to those based on the locations of BH XRBs away from the Galactic plane (e.g., Repetto et al. 2012;Repetto & Nelemans 2015) or those coming from detailed modeling of individual observed BH XRBs (e.g., Brandt et al. 1995;Gualandris et al. 2005;Fragos et al. 2009;Wong et al. 2014).
Similar to the intrinsic population in the Milky Way, the resolvable population also shows a wide range in metallicities.
We also find an expected anti-correlation between metallicity and M BH since metallicity-dependent stellar winds strongly reduce the maximum BH mass for near solar metallicities (e.g., Belczynski et al. 2010). Several other expected correlations are apprent in Figure 8 and 9. For example, longerperiod binaries are detected at further distances; all resolvable detached BH-LC binaries reside right of a line representing the limiting angular size of the orbit. Similarly, brighter systems are detected at higher distances. In contrast, no clear correlation is found between M BH and the distance to the resolvable binary. Some other relations, for example between distance and metallicity, simply come from the fact that we have used initial conditions such that the complex dependence between age-metallicity-distance in the Milky Way is preserved (subsection 2.1).

Effects of Reddening
In this section we evaluate the effects of reddening on our estimated yield of detached BH-LC binaries by Gaia. The details of our method to determine extinction and reddening ar described in subsection 3.1. Figure 11 shows the expected number of detections as a function of Gaia's G magnitude after taking into account the effects of extinction and reddening. We find that reddening and extinction significantly decreases the expected number of detections (Table 1). After correcting for reddening and extinction we estimate that Gaia would resolve between 292 +21 −26 (optimistic) and 152 +14 −17 (pessimistic) BH-LCs in the rapid model. Whereas, the corresponding numbers for the delayed model are 95 +12 −12 (optimistic) and 36 +7 −8 (pessimistic). Our extinction-and reddeningcorrected estimate for the BH-LC binaries resolvable by Gaia provide the most realistic estimates for detection yields. Thus, in the following sections we proceed with  Figure 10. The ratio between the Gaia-resolvable and the intrinsic numbers of detached BH-LC binaries with P orb /yr ≤ 10 in the Milky Way as a function of MBH. Orange and blue denote rapid and delayed models. The horizontal and vertical errorbars denote the MBH bin span for the calculation and range between the 10th and 90th percentiles, respectively. We do not find a strong correlation between Gaia's ability to resolve a BH-LC binary and MBH.
the extinction-and reddening-modulated resolvable systems.

Possible astrometric constraints on M BH
In this section we investigate whether Gaia's astrometry alone can clearly identify the nature of an unseen dark companion to an observed LC as a BH via mass constraints with sufficient accuracy. Note that, in our analysis, the LC exclusively dominates the emission observed by Gaia. Hence, a higher-mass unseen object is most likely to be a dark remnant. This essentially means that if Gaia's astrometry can provide mass constraints that clearly indicate that an unseen companion is more massive than any NSs or the LC mass, it must be a BH candidate.
By Gaia's final data release, the parallax to identified binaries will be supplied. As a result, the absolute luminosity and hence the mass of the LC can be constrained from the Gaia magnitudes using standard stellar evolution models (e.g., Anders et al. 2019;Howes et al. 2019). The constrained mass of the LC and the determined P orb from Gaia astrometry can then provide constraints for the unseen dark remnant in the binary. In Andrews et al. (2019) we showed that this is possible and the accuracy in the mass measurement for the dark component depends on the accuracy in the mass measurement of the LC and the accuracy of the along and across scans for the source. The analytic expression and the details of the calculation for mass errors of the BHs are given in subsection 3.3. Counts (< G) Figure 11. Cumulative number of Gaia-resolvable BH-LC binaries after incorporating the effects of extinction and reddening. The line-styles, colors, and the shades have the same meaning as in Figure 5. After correcting for extinction and reddening, the expected number of detections is between ≈ 36 (delayed model with pessimistic cut) 292 (rapid model with optimistic cut). We find a higher expected yield for the rapid model compared to the delayed model by a factor of about 3.0 (4.2) using the optimistic (pessimistic) cut. This difference can be attributed to typically lower-mass BHs and larger kicks in case of delayed.
In order to find the subset of BH-LC binaries where Gaia's astrometry alone can reveal that the dark remnant is a BH, in addition to all of Gaia's detectability considerations (subsection 3.2) and effects of extinction and reddening (subsection 3.1), we employ two additional filters, (i) M BH > M LC and (ii) M BH − ∆M BH ≥ 3 M . The former condition is to rule out confusion potentially created by very high-mass bright primary in orbit with another less bright secondary of mass above 3 M in observed systems. We denote this subset as BH candidates. Note that our constraint of M BH − ∆M BH ≥ 3 M is conservative and consistent with the assumed boundary between NSs and BHs in our simulations. We do not suggest that NSs can be as massive. The estimated masses of the heavy NSs to date  Figure 11, but after imposing the additional constraints, MBH − ∆MBH ≥ 3 M and MBH ≥ MLC on the extinction-and reddening-corrected detection numbers. We call this subset BH candidates. We estimate ∆MBH from Equation 12. We assume, on an average, 10% error in estimating MLC. These numbers indicate the number of cases where Gaia's astrometry of the LC alone can identify that the dark companion must be a BH.
are lower than 3 M (e.g., Cromartie et al. 2020;Biswas et al. 2021). Figure 12 shows the expected number of BH-LC binaries as a function of Gaia's G magnitude, where astrometry alone can determine that the mass of the unseen primary to a secondary LC is ≥ 3 M . We find that all of the above stringent criteria are satisfied for ≈ 37% (25%) of reddening-corrected resolvable systems in the rapid (delayed) model. According to our optimistic threshold for Gaia's astrometric precision, the total numbers of resolvable BH-LC binaries with astrometric mass measurements strongly indicating a BH as the dark object, is 107 +11 −12 (24 +7 −5 ) for the rapid (delayed) model. Based on the pessimistic threshold, the corresponding number is 47 +9 −8 (10 +4 −4 ) ( Table 1). While the results shown in Figure 12 uses the fiducial value ∆M LC /M LC = 0.1, we note that the preci- sion in mass estimates of isolated single or binary stars can vary widely depending on the system, data quality, as well as the measurement technique (e.g., Torres & Ribas 2002;Gallenne et al. 2016;Pavlovski et al. 2018;Brogaard et al. 2018;Rendle et al. 2019). Our Equation 12 provides a straightforward way to use any desired value of ∆M LC /M LC to ultimately update the estimated ∆M BH /M BH through Gaia's astrometry. Nevertheless, by varying ∆M LC /M LC between the full range of zero and one, we find that the estimated number of BH candidates do not depend strongly on the chosen value of ∆M LC /M LC (Figure 13). In fact, throughout the full range, the estimated numbers stay within the range due to statistical fluctuations for any one of the chosen values of ∆M LC /M LC including our fiducial value of 0.1. The lack of sensitivity towards the adopted value of ∆M LC /M LC can be easily understood from Equation 12. To identify BH candidates, we impose the condition M BH ≥ M LC . This makes the term involving ∆M BH /M BH less significant compared to the other term involving astrometric precision.

PROMISE FROM FOLLOWUP MEASUREMENTS BEYOND ASTROMETRY
While our primary focus of this work is on the ability of Gaia to detect and characterise detached BH-LC binaries using astrometry alone, of course, once identified by Gaia's astrometry, followup observations of these binaries can significantly improve the characterisation of the binary properties. In this section we highlight a couple of these promising avenues.

Radial velocities using Gaia's spectroscopy
The most readily available avenue to improve characterisation of astrometrically identified candidates of interest is RV measurements by Gaia itself. Gaia contains a RV spectrometer which provides spectral resolution R ∼ 5000 (low) and R ∼ 11500 (high) for stars with G ≤ 17 and G ≤ 11, respectively. We find that Gaia's high-R spectra may not be very useful in characterising BH-LC binaries. While the high-R allows measurement of lower RV semiamplitudes (K), the stronger constraint on G prohibitively limits the number of available BH-LC binaries in the Milky Way. Figure 14 shows the distribution RV semiamplitude (K) for the reddening-and extinction-corrected, resolvable BH-LC binaries with G < 17. We find that in the rapid model ∼ 35 BH-MS and 21 BH-PMS binaries could have high enough K to be resolved by Gaia's low-R spectroscopy ( Figure 14). In comparison, in the delayed model ∼ 14 BH-MS and ∼ 15 BH-PMS binaries are expected to have high enough K to be resolved by Gaia's low-R spectra. Using Gaia's high-R spectra, these numbers are at most a few in both rapid and delayed models. As can be seen in Figure 14, K for Gaia-resolvable BH-LC binaries is rather large. Hence, a compromise in the spectral resolution matters less than the loss of available systems in the Milky Way from a stronger constraint in G.
Overall, we find that a little over 35% of the detached Gaia-resolvable (reddening-corrected) BH-LC binaries will have some constraints on M BH either through astrometry (as described in subsection 4.5) or through RV measurements using Gaia's low-R spectra. Although, we limit our discussion to RV measurements using Gaia's spectrometer, of course, once the BH binary candidates are identified, a higher fraction of systems may be constrained via RV followup using other higher-R telescopes that can probe fainter stars. We also caution the readers that in our quick analysis we did not consider the details that affect the feasibility Counts (>log K) Figure 14. Reverse cumulative distribution of the RV semiamplitude (K) for detached BH-LC binaries resolvable by Gaia (using the optimistic threshold for astrometric resolution) with G < 17 after correcting for extinction and reddening. The vertical line denotes the level of K resolved by Gaia's low spectral resolution (R) data. Red and blue denote BH binaries with MS and PMS companions. Lines and shaded regions around the lines denote the median and the spread between the 10th and 90th percentiles.
of RV measurements of a particular LC which include availability of lines to observe within the band of the spectrometer, Galactic position, and stellar activity.
Till date the discovered detached BH-LC binaries have used data from multi-epoch spectroscopic surveys by the MUSE (e.g., Giesers et al. 2018;Giesers et al. 2019) and APOGEE collaborations (e.g., Thompson et al. 2019;Jayasinghe et al. 2021). MUSE specifically focuses in crowded fields such as globular clusters, whereas, primary APOGEE targets are in the field. Hence, it is interesting to compare Gaia's capabilities with that of APOGEE in the context of detecting BH-LC binaries using RV measurements. In order to find the number of Gaia-resolvable BH-LC binaries APOGEE may also detect, we take the subset of reddening-corrected Gaiaresolvable BH-LC binaries in our models and apply two  Figure 15. Distribution of the wavelengths (λ peak ) corresponding to the peak in the blackbody spectral-energy distribution for the resolvable rapid (red) and delayed (purple) binaries. Blue and orange vertical bands show the ranges of bandpass filters for Gaia (330-1050 nm) and APOGEE (1.51-1.70 µm). The wavelength band for Gaia is more favorable for the majority of the BH-LC binaries compared to that of APOGEE.
additional filters relevant for APOGEE. In particular, we use J − K s ≤ 0.5 and 12.2 ≤ H ≤ 14 (Zasowski et al. 2013;Price-Whelan et al. 2020). We find that APOGEE cannot detect any BH-MS binaries in our Gaia-resolvable populations and may detect ∼ 1 (∼ 3) BH-PMS binaries in the rapid (delayed) model. Note that, both APOGEE detections of BH candidates so far contain PMS companions (e.g., Thompson et al. 2019;Jayasinghe et al. 2021). The reason why Gaia is expected to be more prolific in discovering BH-LC candidates, especially, BH-MS binaries, can be understood by comparing the filters of Gaia and APOGEE ( Figure 15). APOGEE is sensitive to much cooler stars compared to Gaia. Most of the detached BH-LC binaries that Gaia can resolve are too faint in APOGEE.

X-ray counterparts
A fraction of the reddening-corrected, resolvable, detached BH-LC binaries may have X-ray counterparts if the stellar wind from the LC is accreted by the BH at a sufficiently high rate. We can estimate the bolometric X-ray luminosity for these wind-fed systems as- whereṀ acc is the accretion rate of the BH, R acc is the accretion radius, and ε is an efficiency parameter for the conversion of gravitational binding energy to radiation. For our calculations we assume, R acc ≈ 3 × R sc (R sc is the Schwarzschild radius) andṀ acc ∼ 10 −12 - The vertical pink lines denote FX = 6 × 10 −15 ergs s −1 cm −2 (dotted) and FX = 2 × 10 −14 ergs s −1 cm −2 (dash-dot), our adopted lower limits for detectability by Chandra and eROSITA, respectively. Red, blue, and purple denote BH-MS, BH-PMS, and all BH-LC binaries. Top and bottom panels are for the rapid and delayed models.
10 −8 M yr −1 . SinceṀ acc for the resolvable binaries is expected to be low, ≤ 10% of the Eddington rate, we adopt that the accretion process is in the regime of the advection dominated accretion flow (ADAF). We adopt the prescription given in Xie & Yuan (2012, their equation 11) to calculate the efficiency parameter. In particular, assuming the viscosity parameter α = 0.1, viscous heating parameter δ = 0.5, we find ε ranges from 10 −5 −0.082 based on the ratio ofṀ acc to the Eddington mass accretion rate.
We assume that the bolometric X-ray luminosity calculated using Equation 13 includes X-rays in the energy range 0.1 − 500keV. Since X-ray observatories like Chandra and eRosita are sensitive in the energy band 0.1 − 10keV, we evaluate the incident flux F X in the en-ergy band 0.1 − 10keV assuming a power law spectrum for X-rays, with photon index Γ = 2 (Yang et al. 2015). Furthermore, we correct F X for interstellar extinction and reddening using where σ ISM (E) is the total photoionization cross section of ISM taken from Wilms et al. (2000), and N H is the hydrogen column density between the source and the observer. We estimate N H simply from the positiondependent A v obtained using mwdust: (Güver &Özel 2009). We show the reverse cumulative numbers of reddeningcorrected Gaia-resolvable detached BH-LC binaries as a function of F X in Figure 16. We adopt, the detection limits of Chandra and eROSITA as F X ≥ 6×10 −15 ergs s −1 cm −2 (based on CSC2.0 catalog Weisskopf et al. 2000) and F X ≥ 2 × 10 −14 ergs s −1 cm −2 (Merloni et al. 2012), respectively. We find that the number of reddening-corrected, Gaia-resolvable, detached BH-LC binaries expected to have X-ray counterparts in eROSITA are 68 +9 −10 (24 +6 −5 ) and 24 +7 −6 (7 +4 −3 ) for the rapid and delayed models adopting the optimistic (pessimistic) astrometric threshold. The corresponding numbers for Chandra are 88 +11 −11 (33 +7 −7 ) and 25 +8 −6 (8 +4 −4 ), respectively. Interestingly, although further characterisation is needed to ascertain the nature of these sources several studies have reported possible BH-LC candidates by cross-matching Gaia astrometry and X-ray counterparts (e.g., Gandhi et al. 2020;Price-Whelan et al. 2020).

SUMMARY AND DISCUSSION
In this paper we have investigated the possibility for Gaia to detect a population of detached BH-LC binaries using astrometry. Using a state-of-the-art population synthesis code COSMIC (Breivik et al. 2020) and realistic stellar ages, metallicities, and Galactic positions adopted from the simulated Milky Way-mass galaxy m12i from the Latte suite of FIRE-2 simulations (Wetzel et al. 2016;Hopkins et al. 2018;Sanderson et al. 2020), which take into account the observed complex correlations between these parameters, we have created highly realistic present-day populations of BH-LC binaries expected to be found in the Milky Way adopting two widely used SN prescriptions, rapid and delayed (e.g., Fryer et al. 2012, ; section 2).
We have presented several relevant properties of our simulated present day populations of BH-LC binaries in the Milky Way, including the component masses, orbital properties, age, and metallicity (subsection 4.1). We summarize these findings below.
• The orbital period distribution, regardless of SN prescription, is bimodal with peaks corresponding to BH-LC populations with progenitors which have undergone at least one CE phase and those which have not experienced a CE. The extended Gaia mission lifetime of 10 yr allows for characterization of both populations.
• The BH mass distribution depends strongly on the SN prescription choice. The rapid prescription produces BH populations with a mass gap between 3-5M , while the delayed prescription does not.
• BHs with masses 10 M have orbital periods < 10 yr because of the correlation between natal kick strength and BH mass which unbinds BH-LC progenitors in wider orbits. If Gaia discovers BH-LCs with BH masses outside of this range, constraints on the correlation between natal kick strength and BH mass can be imposed.
• Eccentricity is a strong tracer of both natal kick strength and tidal circularisation. The majority of BH-PMS systems are circularized through tides, while ∼ 1% (7%) of BH-MS have Ecc > 0.1 in the rapid (delayed) models.
• The age and metallicity of BH-LC binaries in both the rapid and delayed models are broadly distributed, indicating the potential to observationally constrain a metallicity-dependent BH mass distribution.
In order to determine the subset of BH-LC population resolvable by Gaia, we create 200 realisations of the Milky Way sampling from our simulated intrinsic populations of BH-LC binaries for the rapid and delayed models. Each of these realisations preserve the observed correlations between the age, metallicity, and Galactic positions of stars in the Milky Way as well as the expected number and parameter distributions of the BH-LC population (subsection 2.4). For each of these realisations and each SN model, we investigate detectability using Gaia's astrometry (section 3). We summarise our findings below.
• Extinction and reddening from interstellar dust reduces the expected yield of Gaia-resolvable BH-LCs. The expected yield after reddening correction is 292 (95) for the rapid (delayed) model using the optimistic threshold; for the pessimistic threshold, these numbers are 152 (36).
• We find that Gaia's astrometry alone can identify BH candidates with some certainty for ≈ 107 (≈ 24) BH-LC binaries in our rapid (delayed) model. In addition to being resolvable (after reddening correction) by Gaia's astrometry, these systems satisfy M BH − ∆M BH ≥ 3 M and M BH ≥ M LC .
• Independent of the adopted SN prescription, ∼ 50% of the resolvable BH-MS progenitors undergo at least one CE phase while 90% of the resolvable BH-PMS experience at least one CE phase.
• The predicted intrinsic vs Gaia-resolvable Milky-Way populations of BH-LCs do not show a strong correlation with BH mass. This suggests that Gaia's observational selection does not strongly bias our view of the BHs in detached binaries in the Galaxy.
• Intrinsic correlations in BH mass and eccentricity with SN prescription and BH mass with metallicity persists in the Gaia-resolvable population, highlighting opportunities to provide constraints on uncertain SN physics with future population detections.
In this paper we have focused on the promise from Gaia's astrometry. Nevertheless, we point out that the orbital solutions and component characterisation of the BH-LC population potentially resolvable by Gaia astrometry can be improved significantly by analysing Gaia's RV data as well as RV followup using other telescopes. We find that the RV semi-amplitude K can be several to hundreds of km s −1 for the reddeningcorrected resolvable population of BH-LC binaries in our models (subsection 5.1, Figure 14). Gaia's own low-resolution spectroscopy can resolve the RV for ∼ 56 (∼ 29) BH-LC binaries within the reddening-corrected resolvable population in our rapid (delayed) model. In both models, about 35% of the extinction-and reddening-corrected Gaia-resolvable BH-LCs either have high enough K to be resolved by Gaia's low-R spectra or they can be identified as BH candidates through astrometric mass constraint independent of the adopted threshold for astrometric precision.
At present, APOGEE data provides the most relevant all-sky RV survey we can compare with. We find that the APOGEE RV survey may not be as prolific as what is expected of Gaia's RV data based on our models. This is because most resolvable BH-LC binaries, except for some BH-PMS binaries, are expected to be too faint in the APOGEE bands ( Figure 15).
Our models show that between ∼ 7 to 88 (depending on the adopted SN prescription and astrometric threshold) of the reddening-corrected Gaia-resolvable detached BH-LC binaries may also have X-ray counterparts as wind-fed systems potentially detectable by Chandra and eROSITA. In addition to these detached wind-fed systems, we estimate that the orbital motion of the PMS companion may be resolvable by Gaia astrometry for ≈ 76 BH-PMS binaries currently transferring mass via ROLF in both rapid and delayed models. We have ignored these binaries in this work to focus on the detached BH-LC binaries. We suspect that finding an orbital solution astrometrically for the non-detached systems will be significantly more challenging compared to doing so for the detached systems even if Gaia astrometry can resolve the orbital motion of the LC. However, these systems can be interesting candidates for multimessenger studies potentially creating a population of BH-LC binaries that connect the populations of BH Xray binaries and detached BH binaries potentially resolvable in large numbers by Gaia soon.
In conclusion, our models suggest that Gaia data could dramatically improve our understanding of the properties and demographics of BH binaries in the Milky Way. Such detections are expected to have few selections biases depending on the BH mass and would provide unprecedented constraints on BH-progenitor properties by constraining the stellar properties of the observed LC.
We thank the referee for constructive comments. CC acknowledges support from TIFR's graduate fellowship. SC acknowledges support from the Department of Atomic Energy, Government of India, under project no. 12-R&D-TFR-5.02-0200 and RTI 4002. JA acknowledges funding from CIERA and Northwestern University through a postdoctoral fellowship. The Flatiron Institute is supported by the Simons Foundation.

CODE AND DATA AVAILABILITY
The data and code used for this study is freely accessible on Zenodo Chawla et al. (2021).

APPENDIX
Here we present selected results from our delayed model. Figure 17 shows the P orb distribution for the delayed model for our converged population of BH-LC binaries. The distribution of P orb in the delayed model is very similar to the same for the rapid model (Figure 1). Figure 18 shows the relevant LC properties for the intrinsic population of BH-LC binaries in our delayed model. We have used propagation of error to estimate ∆M BH from the analytic expression for the astrometric mass function for the LC's motion (without the sin i degeneracy) given in Equation 11. Andrews et al. (2019) derived this expression using a suite of Markov-chain Monete Carlo parameter-estimation exercise using injected LC motion as would be seen by Gaia assuming a 5-yr mission duration. We have updated their expression for a 10-yr Gaia mission for our purpose. Below we show the steps used to derive Equation 12 from Equation 11. For simplicity, we write the astrometric mass function without the sin i term as-