Gaia May Detect Hundreds of Well-characterized Stellar Black Holes

COSMIC generates binary populations from zero-age MS (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) ⁷ 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 Section 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 $\mathrm{log}({\text{}}Z/{\text{}}{Z}_{\odot })=-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) and then sampling semimajor axes that are uniform in log space up to 10⁵ 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.

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 that 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) that 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 MS, 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 that 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 supernova (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 proto−compact object that reinitiate 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 (NSs) and BHs between ∼3 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 σ = 265 km 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.

cosmic is designed to adaptively choose the size of a simulated population based on a convergence criterion that tracks how the shapes 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}_{\mathrm{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., Equation (6) of Chatziioannou et al. 2017) and defined as

$$\begin{eqnarray}&&\mathrm{match}=\displaystyle \frac{{\sum }_{k=1}^{N}{P}_{k,i}{P}_{k,i+1}}{\sqrt{{\sum }_{k=1}^{N}{P}_{k,i}{P}_{k,i}{\sum }_{k=1}^{N}{P}_{k,i+1}{P}_{k,i+1}}},\end{eqnarray} \tag{ 1 }$$

where P_k,i represents the height of bin k on the ith iteration (Breivik et al. 2020). The match criterion approaches unity as the shapes of each histogram converge. We use Knuth's rule to determine the bin widths 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⁻⁵ for the BH mass (M_BH), LC mass (M_LC), P_orb, and $Ecc$.

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

$$\begin{eqnarray}&&{N}_{\mathrm{BH}--\mathrm{LC},\mathrm{MW}}={N}_{\mathrm{BH}--\mathrm{LC},\mathrm{sim}}\displaystyle \frac{{M}_{{\rm{m}}12{\rm{i}}}}{{M}_{\mathrm{sim}}},\end{eqnarray} \tag{ 2 }$$

where M_m12i is the total amount of mass formed in galaxy m12i and ${M}_{\mathrm{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 realizations of the Milky Way population for each of the rapid and delayed models, producing 400 Milky Way realizations in total.

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 position-dependent dust map that 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). ⁸

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 constants:

$$\begin{eqnarray}&&{ \mathcal A }={a}_{\mathrm{LC}}(\cos \omega \ \cos {\rm{\Omega }}-\sin \omega \ \sin {\rm{\Omega }}\ \cos i)\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{ \mathcal B }={a}_{\mathrm{LC}}(\cos \omega \ \sin {\rm{\Omega }}\ +\sin \omega \ \cos {\rm{\Omega }}\ \cos i)\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{ \mathcal F }={a}_{\mathrm{LC}}(-\sin \omega \ \cos {\rm{\Omega }}\ -\cos \omega \ \sin {\rm{\Omega }}\ \cos i)\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{ \mathcal G }={a}_{\mathrm{LC}}(-\sin \omega \ \sin {\rm{\Omega }}\ +\cos \omega \ \cos {\rm{\Omega }}\ \cos i).\end{eqnarray} \tag{ 6 }$$

In the equations above, a_LC is the semimajor axis of the LC orbit defined as

$$\begin{eqnarray}&&{a}_{\mathrm{LC}}=a\displaystyle \frac{{M}_{\mathrm{BH}}}{{M}_{\mathrm{LC}}+{M}_{\mathrm{BH}}},\end{eqnarray} \tag{ 7 }$$

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

$$\begin{eqnarray}&&{x}_{\mathrm{project},\ \mathrm{LC}}={ \mathcal A }{ \mathcal X }+{ \mathcal F }{ \mathcal Y }\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{y}_{\mathrm{project},\ \mathrm{LC}}={ \mathcal B }{ \mathcal X }+{ \mathcal G }{ \mathcal Y },\end{eqnarray} \tag{ 9 }$$

where ${ \mathcal X }=\cos \,E-{Ecc}$, ${ \mathcal Y }=\sqrt{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

$$\begin{eqnarray}&&\alpha =\displaystyle \frac{{a}_{\mathrm{project},\ \mathrm{LC}}}{D},\end{eqnarray} \tag{ 10 }$$

where D is the distance to the binary.

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 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

$$\begin{eqnarray}&&\left|\displaystyle \frac{{\rm{\Delta }}({M}_{\mathrm{BH}}^{3}{M}_{\mathrm{tot}}^{-2})}{{M}_{\mathrm{BH}}^{3}{M}_{\mathrm{tot}}^{-2}}\right|=0.9\left(\displaystyle \frac{{\sigma }_{\xi }}{\alpha }\right){\left(\displaystyle \frac{N}{75}\right)}^{-1/2},\end{eqnarray} \tag{ 11 }$$

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 decl., each 05 apart between −90° and 90°. We divide R.A. for each of the above decl. into n_R.A. bins, where n_R.A. is the nearest integer $\geqslant 720\times \cos (\mathrm{dec})$. Our grid points ensure that the angular sky-projected area between any two adjacent grid points in R.A. and decl. is smaller than the field of view (072 × 069) for each of Gaia's telescopes (e.g., Lindegren et al. 2012). Using the Gaia Observation Forecast Tool (GOST), ⁹ we calculate and store the number of times Gaia would observe each grid point in 5 yr (N₅; between 2014 and 2019). We find the value of N₅ for any BH–LC binary in our Milky Way realizations simply by looking up this number corresponding to the nearest grid point from the R.A. and decl. of that particular model binary. Finally, we set N = 2N₅ 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: ¹⁰

$$\begin{eqnarray}\begin{array}{rcl}\left|\displaystyle \frac{{\rm{\Delta }}{M}_{\mathrm{BH}}}{{M}_{\mathrm{BH}}}\right| & = & 0.1\left|\displaystyle \frac{{\rm{\Delta }}{M}_{\mathrm{LC}}}{0.1{M}_{\mathrm{LC}}}\right|\left(\displaystyle \frac{2{M}_{\mathrm{LC}}}{{M}_{\mathrm{tot}}+2{M}_{\mathrm{LC}}}\right)\\ & & +\left|\displaystyle \frac{{\rm{\Delta }}({M}_{\mathrm{BH}}^{3}{M}_{\mathrm{tot}}^{-2})}{{M}_{\mathrm{BH}}^{3}{M}_{\mathrm{tot}}^{-2}}\right|\left(\displaystyle \frac{{M}_{\mathrm{tot}}}{{M}_{\mathrm{tot}}+2{M}_{\mathrm{LC}}}\right)\\ & = & 0.1\left|\displaystyle \frac{{\rm{\Delta }}{M}_{\mathrm{LC}}}{0.1{M}_{\mathrm{LC}}}\right|\left(\displaystyle \frac{2{M}_{\mathrm{LC}}}{{M}_{\mathrm{BH}}+3{M}_{\mathrm{LC}}}\right)\\ & & +0.9\left(\displaystyle \frac{{\sigma }_{\xi }}{\alpha }\right){\left(\displaystyle \frac{2{N}_{5}}{75}\right)}^{-1/2}\left(\displaystyle \frac{{M}_{\mathrm{BH}}+{M}_{\mathrm{LC}}}{{M}_{\mathrm{BH}}+3{M}_{\mathrm{LC}}}\right).\end{array}\end{eqnarray} \tag{ 12 }$$

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.

In order to gain insight for the underlying population of BH–LCs in the Milky Way at present, we first focus on 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 as part of Figure set 1 in the online journal). The P_orb distributions are clearly bimodal for both BH–MS and BH–post-MS (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 long-period 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 and 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (2 panels)

Tarball of figure set images

One of the most prominent differences between the two SN mechanisms manifests in the M_BH distribution; while the rapid prescription for core-collapse SNe does not allow production of BHs in the mass range 3–5 M_⊙ (often called the "mass gap" between 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 with MS (PMS) companions, and the rest come from core-collapse SNe. The M_BH distribution is 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.

Zoom In Zoom Out Reset image size

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 evolution 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–9 comes 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 that about 2% (15%) and 27% (76%) of BH–MS (BH–PMS) binaries with P_orb/yr ≤ 10 have near-circular ($Ecc$ ≤ 10⁻²) 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 circularization. In spite of this, the high fraction of eccentric BH–LC binaries primarily results 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 that are affected more strongly by tides after BH formation.

Zoom In Zoom Out Reset image size

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 as part of Figure set 4 in the online journal. 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 gigayears, 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 pileup at the last metallicity bins near $\mathrm{log}({\text{}}Z/{\text{}}{Z}_{\odot })=-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 Section 2.1). Nevertheless, the wide range in metallicities of our predicted present-day BH–LC binaries in the Milky Way indicates that if a sizable 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (2 panels)

Tarball of figure set images

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, BHs 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 grows. These stars live longer than normal owing 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.

We estimate that the Milky Way at present harbors ∼4 × 10⁵ (∼8 × 10⁴) BH–LC binaries according to the rapid (delayed) model. Of these, $17,{621}_{-144}^{+145}$ ($29,{616}_{-189}^{+181}$) have P_orb/yr ≤ 10. Of the BH–LC binaries with P_orb/yr ≤ 10, $9,{184}_{-133}^{+120}$ ($7,{091}_{-100}^{+113}$) BH–LC binaries are in detached configurations in the rapid (delayed) model (Table 1). ¹¹ We focus on the detached BH–LC binaries with P_orb/yr ≤ 10 in the Milky Way for Gaia's ability to resolve the LC's orbits throughout the rest of the paper.

Table 1. BH–LC Binaries in the Milky Way

Model	Type	Intrinsic			Resolvable
		All	Porb/yr < 10		Optimistic			Pessimistic
			All	Detached	Resolved	Av -corrected	BH Candidate	Resolved	Av -corrected	BH Candidates
	MS	369, 199	$14,{679}_{-145}^{+142}$	$8,{658}_{-123}^{+113}$	${383}_{-28}^{+33}$	${255}_{-26}^{+19}$	${76}_{-11}^{+9}$	${191}_{-19}^{+20}$	${139}_{-18}^{+12}$	${35}_{-7}^{+8}$
rapid	PMS	18, 990	$2,{947}_{-71}^{+52}$	${529}_{-27}^{+26}$	${88}_{-12}^{+11}$	${36}_{-7}^{+9}$	${31}_{-6}^{+8}$	${47}_{-9}^{+9}$	${14}_{-4}^{+5}$	${12}_{-4}^{+4}$
	Total	388, 189	$17,{621}_{-144}^{+145}$	$9,{184}_{-133}^{+120}$	${470}_{-30}^{+32}$	${292}_{-26}^{+21}$	${107}_{-12}^{+11}$	${239}_{-20}^{+21}$	${152}_{-17}^{+14}$	${47}_{-8}^{+9}$
	MS	70, 694	$24,{731}_{-175}^{+166}$	$6,{653}_{-104}^{+106}$	${125}_{-12}^{+15}$	${68}_{-9}^{+11}$	${13}_{-4}^{+5}$	${40}_{-8}^{+8}$	${26}_{-6}^{+7}$	${6}_{-3}^{+3}$
delayed	PMS	6, 355	$4,{886}_{-48}^{+47}$	${443}_{-28}^{+26}$	${57}_{-8}^{+10}$	${27}_{-7}^{+6}$	${11}_{-3}^{+5}$	${19}_{-5}^{+5}$	${9}_{-4}^{+3}$	${4}_{-2}^{+4}$
	Total	77, 049	$29,{616}_{-189}^{+181}$	$7,{091}_{-100}^{+113}$	${184}_{-14}^{+15}$	${95}_{-12}^{+12}$	${24}_{-5}^{+7}$	${60}_{-12}^{+8}$	${36}_{-8}^{+7}$	${10}_{-4}^{+4}$

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 realizations. "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 (Section 3.2). Any detached BH–LC binary with P_orb/yr ≤ 10 and resolved LC orbit is denoted as "Resolved." "A_v -corrected" denotes the expected number of resolved BH–LC binaries after correcting for extinction and reddening from interstellar dust (Section 3.1). A BH–LC is a "BH candidate" if Gaia's astrometry can constrain the BH's mass M_BH − ΔM_BH ≥ 3 M_⊙ and M_BH ≥ M_LC within the A_v -corrected resolved population (Section 3.3).

Download table as:  ASCIITypeset image

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 (Section 2.4). These realizations 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 realization. These differences provide an estimate of the level of statistical fluctuations in the expected BH–LC numbers resolvable by Gaia.

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 the 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.

Zoom In Zoom Out Reset image size

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$ (Figures 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.

Figures 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 the case of the rapid model, above 80% of resolvable BH–LC binaries are expected to contain an 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%; Section 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 HG during one of the CE episodes. Whether a binary can survive a phase of CE evolution with an 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 an 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 is 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 yr 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 realized, 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.

In this section we highlight some key properties of the resolvable BH–LC binaries. Comparison between the distribution of properties presented here and those presented in Section 4.1 sheds light on the selection biases of Gaia. Figures 8 and 9 show the corner plots for all resolvable binaries across all Milky Way realizations from our rapid and delayed models. Note that different Milky Way realizations can sample the same binaries multiple times (Section 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 an 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 SNe. 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 SNe 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 does 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 does 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 resolvable by Gaia is expected to exhibit an M_BH distribution that is very close to the intrinsic one in the Milky Way.

Zoom In Zoom Out Reset image size

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 cirularization 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 compared 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 stellar-mass BHs. Such constraints would be complimentary to those based on the locations of BH X-ray binaries away from the Galactic plane (e.g., Repetto et al. 2012; Repetto & Nelemans 2015) or those coming from detailed modeling of individual observed BH X-ray binaries (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 anticorrelation 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 apparent in Figures 8 and 9. For example, longer-period binaries are detected at further distances; all resolvable detached BH–LC binaries reside to the 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, and distance in the Milky Way is preserved (Section 2.1).

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 are described in Section 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 decrease the expected number of detections (Table 1). After correcting for reddening and extinction, we estimate that Gaia would resolve between ${292}_{-26}^{+21}$ (optimistic) and ${152}_{-17}^{+14}$ (pessimistic) BH–LCs in the rapid model, whereas the corresponding numbers for the delayed model are ${95}_{-12}^{+12}$ (optimistic) and ${36}_{-8}^{+7}$ (pessimistic). Our extinction- and reddening-corrected estimates for the BH–LC binaries resolvable by Gaia provide the most realistic estimates for detection yields. Thus, in the following sections we proceed with the extinction- and reddening-modulated resolvable systems.

Zoom In Zoom Out Reset image size

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 that 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 Section 3.3.

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 (Section 3.2) and effects of extinction and reddening (Section 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 a 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 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 number of resolvable BH–LC binaries with astrometric mass measurements strongly indicating a BH as the dark object is ${107}_{-12}^{+11}$ (${24}_{-5}^{+7}$) for the rapid (delayed) model. Based on the pessimistic threshold, the corresponding number is ${47}_{-8}^{+9}$ (${10}_{-4}^{+4}$) (Table 1).

Zoom In Zoom Out Reset image size

While the results shown in Figure 12 use the fiducial value ΔM_LC/M_LC = 0.1, we note that the precision in mass estimates of isolated single or binary stars can vary widely depending on the system, data quality, and measurement technique (e.g., Torres & Ribas 2002; Gallenne et al. 2016; Brogaard et al. 2018; Pavlovski 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 0 and 1, we find that the estimated number of BH candidates does 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 toward 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.

Zoom In Zoom Out Reset image size

The most readily available avenue to improve characterization of astrometrically identified candidates of interest is RV measurements by Gaia itself. Gaia contains an RV spectrometer that provides spectral resolution R ∼ 5000 (low) and R ∼ 11,500 (high) for stars with G ≤ 17 and G ≤ 11, respectively. We find that Gaia's high-R spectra may not be very useful in characterizing 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 the 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.

Zoom In Zoom Out Reset image size

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 Section 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 follow-up 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 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.

To date the discovered detached BH–LC binaries have used data from multiepoch spectroscopic surveys by the MUSE (e.g., Giesers et al. 2018, 2019) and APOGEE collaborations (e.g., Thompson et al. 2019; Jayasinghe et al. 2021). MUSE specifically focuses on 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 Gaia-resolvable BH–LC binaries in our models and apply two 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.

Zoom In Zoom Out Reset image size

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

$$\begin{eqnarray}&&{L}_{X}=\varepsilon \displaystyle \frac{{{GM}}_{\mathrm{BH}}{\dot{M}}_{\mathrm{acc}}}{{R}_{\mathrm{acc}}},\end{eqnarray} \tag{ 13 }$$

where ${\dot{M}}_{\mathrm{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 ${\dot{M}}_{\mathrm{acc}}\sim {10}^{-12}$–10⁻⁸ M_⊙yr⁻¹. Since ${\dot{M}}_{\mathrm{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. 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 and viscous heating parameter δ = 0.5, we find that ε ranges from 10⁻⁵ to 0.082 based on the ratio of ${\dot{M}}_{\mathrm{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–500 keV. Since X-ray observatories like Chandra and eRosita are sensitive in the energy band 0.1–10 keV, we evaluate the incident flux ${{ \mathcal F }}_{{\rm{X}}}$ in the energy band 0.1–10 keV assuming a power-law spectrum for X-rays,

$$\begin{eqnarray}&&N(E)\propto {E}^{-{\rm{\Gamma }}},\end{eqnarray} \tag{ 14 }$$

with photon index Γ = 2 (Yang et al. 2015). Furthermore, we correct ${{ \mathcal F }}_{{\rm{X}}}$ for interstellar extinction and reddening using

$$\begin{eqnarray}&&{{ \mathcal F }}_{{\rm{X}}}={\int }_{0.1}^{10}{e}^{{\sigma }_{\mathrm{ISM}}(E){N}_{H}}{EN}(E)\,{dE},\end{eqnarray} \tag{ 15 }$$

where σ_ISM(E) is the total photoionization cross section of the interstellar medium 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 position-dependent A_v obtained using mwdust:

$$\begin{eqnarray}&&\frac{{N}_{{\rm{H}}}}{{A}_{v}}=2.2\times {10}^{21}{\mathrm{cm}}^{-2}\end{eqnarray} \tag{ 16 }$$

(Güver & Özel 2050).

We show the reverse cumulative numbers of reddening-corrected Gaia-resolvable detached BH–LC binaries as a function of ${{ \mathcal F }}_{X}$ in Figure 16. We adopt the detection limits of Chandra and eROSITA as ${{ \mathcal F }}_{X}\geqslant 6\times {10}^{-15}\,\mathrm{erg}\ {{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-2}$ (based on CSC2.0 catalog Weisskopf et al. 2000) and ${{ \mathcal F }}_{X}\geqslant 2\times {10}^{-14}\,\mathrm{erg}\ {{\rm{s}}}^{-1}\,{\mathrm{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}_{-10}^{+9}$ (${24}_{-5}^{+6}$) and ${24}_{-6}^{+7}$ (${7}_{-3}^{+4}$) 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}_{-6}^{+8}$ (${8}_{-4}^{+4}$), respectively. Interestingly, although further characterization 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. 2022; Price-Whelan et al. 2020).

Zoom In Zoom Out Reset image size