We report discovery of a bright nearby G  13.8  d  480 rm pc Sunlike star orbiting a dark object. We identified the system as a black hole candidate via its astrometric orbital solution from the Gaia mission. Radial velocities validated and refined the Gaia solution and spectroscopy ruled out significant light contributions from another star. Joint modelling of radial velocities and astrometry constrains the companion mass of M2  9.62pm 0.18 mathrmModot . The spectroscopic orbit alone sets a minimum companion mass of M2gt 5 mathrmModot  if the companion were a 5 mathrmModot  star it would be 500 times more luminous than the entire system. These constraints are insensitive to the mass of the luminous star which appears as a slowly rotating G dwarf Trm eff5850 rm K log g  4.5 M0.93 mathrmModot  with nearsolar metallicity rm FeH  0.2 and an unremarkable abundance pattern. We find no plausible astrophysical scenario that can explain the orbit and does not involve a black hole. The orbital period PSUBorbSUB  185.6 d is longer than that of any known stellarmass black hole binary. The systems modest eccentricity e  0.45 high metallicity and thindisc Galactic orbit suggest that it was born in the Milky Way disc with at most a weak natal kick. How the system formed is uncertain. Common envelope evolution can only produce the systems wide orbit under extreme and likely unphysical assumptions. Formation models involving triples or dynamical assembly in an open cluster may be more promising. This is the nearest known black hole by a factor of 3 and its discovery suggests the existence of a sizable population of dormant black holes in binaries. Future Gaia releases will likely facilitate the discovery of dozens more.
['binaries: spectroscopic', 'stars: black holes', 'Astrophysics - Solar and Stellar Astrophysics', 'Astrophysics - Astrophysics of Galaxies']
A Sunlike star orbiting a black hole
{'A Sun-like star orbiting a black hole': 'Kareem El-Badry 1 GLYPH<150> 2 GLYPH<150> 3 ★ , Hans-Walter Rix 3 , Eliot Quataert 4 , Andrew W. Howard 5 , Howard Isaacson 6 GLYPH<150> 7 , Jim Fuller 5 , Keith Hawkins 8 , Katelyn Breivik 9 , Kaze W. K. Wong 9 , Antonio C. Rodriguez 5 , Charlie Conroy 1 , Sahar Shahaf 10 , Tsevi Mazeh 11 , Frédéric Arenou 12 , Kevin B. Burdge 13 , Dolev Bashi 11 , Simchon Faigler 11 , Daniel R. Weisz 6 , Rhys Seeburger 3 , Silvia Almada Monter 3 , Jennifer Wojno 3  \n- 1 Center for Astrophysics j Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA\n- 2 Harvard Society of Fellows, 78 Mount Auburn Street, Cambridge, MA 02138\n- 3 Max-Planck Institute for Astronomy, Königstuhl 17, D-69117 Heidelberg, Germany\n- 4 Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA\n- 5 Department of Astronomy, California Institute of Technology, Pasadena, CA 91125, USA\n- 6 Department of Astronomy, University of California, Berkeley, 501 Campbell Hall #3411, Berkeley, CA 94720, USA\n- 7 Centre for Astrophysics, University of Southern Queensland, Toowoomba, QLD, Australia\n- 8 Department of Astronomy, The University of Texas at Austin, 2515 Speedway Boulevard, Austin, TX 78712, USA\n- 9 Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Ave, New York, NY, 10010, USA\n- 10 Department of Particle Physics and Astrophysics, Weizmann Institute of Science, Rehovot 7610001, Israel\n- 11 School of Physics and Astronomy, Tel Aviv University, Tel Aviv, 6997801, Israel\n- 12 GEPI, Observatoire de Paris, Université PSL, CNRS, 5 Place Jules Janssen, 92190 Meudon, France\n- 13 MIT-Kavli Institute for Astrophysics and Space Research 77 Massachusetts Ave. Cambridge, MA 02139, USA', 'ABSTRACT': "We report discovery of a bright, nearby ( 𝐺 = 13 GLYPH<149> 8; 𝑑 = 480 pc) Sun-like star orbiting a dark object. We identified the system as a black hole candidate via its astrometric orbital solution from the Gaia mission. Radial velocities validated and refined the Gaia solution, and spectroscopy ruled out significant light contributions from another star. Joint modeling of radial velocities and astrometry constrains the companion mass to 𝑀 2 = 9 GLYPH<149> 62 GLYPH<6> 0 GLYPH<149> 18 𝑀 GLYPH<12> . The spectroscopic orbit alone sets a minimum companion mass of 𝑀 2 ¡ 5 𝑀 GLYPH<12> ; if the companion were a 5 𝑀 GLYPH<12> star, it would be 500 times more luminous than the entire system. These constraints are insensitive to the mass of the luminous star, which appears as a slowly-rotating G dwarf ( 𝑇 eff = 5850 K, log 𝑔 = 4 GLYPH<149> 5, 𝑀 = 0 GLYPH<149> 93 𝑀 GLYPH<12> ), with near-solar metallicity ( » Fe GLYPH<157> H … = GLYPH<0> 0 GLYPH<149> 2) and an unremarkable abundance pattern. We find no plausible astrophysical scenario that can explain the orbit and does not involve a black hole. The orbital period, 𝑃 orb = 185 GLYPH<149> 6 days, is longer than that of any known stellar-mass black hole binary. The system's modest eccentricity ( 𝑒 = 0 GLYPH<149> 45), high metallicity, and thin-disk Galactic orbit suggest that it was born in the Milky Way disk with at most a weak natal kick. How the system formed is uncertain. Common envelope evolution can only produce the system's wide orbit under extreme and likely unphysical assumptions. Formation models involving triples or dynamical assembly in an open cluster may be more promising. This is the nearest known black hole by a factor of 3, and its discovery suggests the existence of a sizable population of dormant black holes in binaries. Future Gaia releases will likely facilitate the discovery of dozens more.  \nKey words: binaries: spectroscopic - stars: black holes", '1 INTRODUCTION': "The Milky Way is expected to contain of order 10 8 stellar mass black holes (BHs), an unknown fraction of which are in binaries. The inventory of known and suspected BHs consists of about 20 dynamically confirmed BHs in X-ray binaries, an additional GLYPH<24> 50 Xray sources suspected to contain a BH based on their X-ray properties (e.g. McClintock & Remillard 2006; Remillard & McClintock 2006; Corral-Santana et al. 2016), a few X-ray quiet binaries in which a BHissuspected on dynamical grounds, and an isolated BH candidate discovered via microlensing (Sahu et al. 2022; Lam et al. 2022; Mroz et al. 2022). In X-ray bright systems, a BH accretes material from a close companion through stable Roche lobe overflow or stellar winds. Based on the distance distribution and outburst properties of  \nknown BH X-ray binaries, it has been estimated that of order 10 3 such systems exist in the Milky Way (Corral-Santana et al. 2016) only a tiny fraction of the expected total Galactic BH population.  \nBHX-raybinariesall have relatively short orbital periods. In Roche lobe overflowing systems with low-mass ( . 1 𝑀 GLYPH<12> ) main-sequence donors ('LMXBs'), mass transfer only occurs at periods 𝑃 orb . 1 day. Systems with massive or evolved donors can overflow their Roche lobes at somewhat longer periods. The longest orbital period for a known BH X-ray binary is 33 days, in GRS 1915+105 (Greiner et al. 2001), where the donor is a red giant. Binary population synthesis models predict that a large fraction of the total BH + normal star binary population may be found in wider binaries (e.g. Breivik et al. 2017; Chawla et al. 2022), where there is no significant mass transfer. These longer-period BH binaries are difficult to find because they are not X-ray bright, but they may represent the vast majority of the BH binary population. Searches for dormant BH binaries began  \neven before the identification of the first BHs in X-ray binaries (Guseinov & Zel'dovich 1966; Trimble & Thorne 1969). The intervening decades have witnessed the proliferation of vast spectroscopic and photometric surveys well-suited for finding dormant BHs, but only a few solid candidates have been identified (e.g. Giesers et al. 2018, 2019; Shenar et al. 2022a; Mahy et al. 2022).  \nGravitational wave observations have in the last decade also begun to detect large number of binaries containing stellar-mass BHs (e.g. The LIGO Scientific Collaboration et al. 2021). The BHs in these systems likely represent a similarly rare outcome of the binary evolution process to X-ray binaries, but their enormous gravitational wave luminosities during merger allow them to be detected all throughout the Universe. The formation channels of these merging BHs and their evolutionary relation to local BHs in X-ray binaries are still uncertain.  \nThe Gaia mission opens a new window on the Galactic binary population - including, potentially, BHs in binaries - through largescale astrometric orbit measurements. Gaia has been predicted to discover large numbers of BHs in binaries, with the predicted number varying by more than 4 orders of magnitude between different studies (e.g. Mashian & Loeb 2017; Breivik et al. 2017; Shao & Li 2019; Andrews et al. 2019; Wiktorowicz et al. 2020; Chawla et al. 2022; Shikauchi et al. 2022; Janssens et al. 2022). The large dispersion in these predictions reflects both inherent uncertainties in binary evolution modeling and different assumptions about the Gaia selection function.  \nThe first binary orbital solutions from Gaia were recently published in the mission's 3rd data release (Gaia Collaboration et al. 2022b,a), including about 170,000 astrometric solutions and 190,000 spectroscopic solutions. Early assessments of the BH candidate population in these datasets have been carried out for both spectroscopic solutions (e.g. El-Badry & Rix 2022; Jayasinghe et al. 2022) and astrometric solutions (e.g. Andrews et al. 2022; Shahaf et al. 2022). The DR3 binary sample represents a factor of GLYPH<24> 100 increase in sample size over all previously published samples of binary orbital solutions, and is thus a promising dataset in which to search for rare objects. At the same time, stringent SNR cuts were applied to the sample of orbital solutions that was actually published in Gaia DR3. The DR3 binary sample thus represents only a few percent of what is expected to be achievable in the mission's data releases DR4 and DR5.  \nThis paper presents detailed follow-up of one astrometric BH binary candidate, which we found to be the most compelling candidate published in DR3. Our follow-up confirms beyond reasonable doubt the object's nature as binary containing a normal star and at least one dormant BH. The remainder of this paper is organized as follows. Section 2 describes how we identified the source as a promising BH candidate. Section 3 presents the radial velocities (RVs) from archival surveys and our follow-up campaign. In Section 4, we constrain the mass of the unseen companion using the RVs and Gaia astrometric solution. Section 5 describes analysis of the luminous star's spectrum, including estimates of the atmospheric parameters and abundance pattern, and the non-detection of a luminous companion. Section 6 describes the system's Galactic orbit, and Section 7 discusses X-ray and radio upper limits. We compare the object to other BHs and BH imposters in Section 8, where we also discuss its evolutionary history. Section 9 discusses constraints on the occurrence rate of wide BH companions to normal stars. Finally, we summarize our results and conclude in Section 10. The Appendices provide further details on several aspects of our data and modeling.", '2 DISCOVERY': "In a search for compact object companions to normal stars with astrometric binary solutions from Gaia , we considered all 168,065 sources in the gaiadr3.nss\\_two\\_body\\_orbit catalog with purely astrometric ( nss\\_solution\\_type = Orbital ) or joint astrometric/spectroscopic ( nss\\_solution\\_type = AstroSpectroSB1 ) solutions. These solutions describe the ellipse traced on the sky by each source's 𝐺 GLYPH<0> band light centroid due to binary motion. In the AstroSpectroSB1 solutions, which are only available for bright ( 𝐺 . 13) sources, RVs and astrometry are fit simultaneously.  \nOur selection strategy and the candidates we considered are described in Appendix E. In brief, we searched for astrometric solutions with unusually large photocenter ellipses at fixed orbital period, exploiting the facts that (a) massive companions have larger orbits at fixed period due to Kepler's 3rd law, and (b) dark companions produce larger photocenter wobbles than would luminous companions of the same mass (e.g. van de Kamp 1975). This yielded 6 initially promising sources. Individual vetting and spectroscopic follow-up showed that in four of the six sources, the astrometric solutions are spurious, making a BH companion unlikely. In one case, the astrometric solution may be correct, but the luminous star is a giant and the orbital period is longer than the Gaia DR3 baseline, making the reliability of the solution and the nature of the companion difficult to assess without long-term spectroscopic monitoring. One candidate emerged as very promising, Gaia DR3 4373465352415301632. We colloquially refer to the source as Gaia BH1.", '2.1 Properties of the luminous source': "Basic observables of the luminous source are summarized in Figure 1. It is a bright ( 𝐺 = 13 GLYPH<149> 77) solar-type star in Ophiuchus ( 𝛼 = 17:28:41.1; 𝛿 = GLYPH<0> 00:34:52). The field (upper left) is moderately but not highly crowded. The nearest companion detected by Gaia is 5.1 magnitudes fainter at a distance of 3.9 arcseconds. None of the Gaia -detected neighbors have parallaxes and proper motions consistent with being bound to the source. In the color-magnitude diagram, the source appears as a solar-type main sequence star. The DR3 orbital astrometric solution puts it at a distance of 477 GLYPH<6> 4 pc.  \nThe Gaia astrometric solution has a photocenter ellipse with semimajor axis 𝑎 0 = 3 GLYPH<149> 00 GLYPH<6> 0 GLYPH<149> 22 mas. 1 Given the constraints on the source's parallax, this corresponds to a projected photocenter semimajor axis of 𝑎 0 GLYPH<157> 𝜛 = 1 GLYPH<149> 44 GLYPH<6> 0 GLYPH<149> 11 au. If the photocenter traced the semimajor axis (which it does for a dark companion in the limit of 𝑀 ★ GLYPH<157> 𝑀 2 ! 0), this would imply a total dynamical mass of 𝑀 tot = 11 GLYPH<149> 5 GLYPH<6> 2 GLYPH<149> 7 𝑀 GLYPH<12> . For larger 𝑀 ★ GLYPH<157> 𝑀 2 or a luminous secondary, the implied total mass would be larger. If we take the mass of the luminous star to be 𝑀 ★ = 0 GLYPH<149> 93 GLYPH<6> 0 GLYPH<149> 05 𝑀 GLYPH<12> (as implied by its temperature and radius; see below), the astrometric mass ratio function (see Shahaf et al. 2019) is A = 2 GLYPH<149> 32 GLYPH<6> 0 GLYPH<149> 17. This is much larger than the maximum value of A GLYPH<25> 0 GLYPH<149> 6 that can be achieved for systems with main-sequence components, including hierarchical triples. For a dark companion, this value of A implies a mass ratio 𝑞 = 𝑀 2 GLYPH<157> 𝑀 ★ = 14 GLYPH<149> 2 GLYPH<6> 2 GLYPH<149> 8.  \nWe retrieved photometry of the source in the GALEX NUV band (Martin et al. 2005), SDSS 𝑢 band (Padmanabhan et al. 2008), PanSTARRS 𝑔𝑟𝑖𝑧𝑦 bands (Chambers et al. 2016), 2MASS 𝐽𝐻𝐾 bands (Skrutskie et al. 2006), and WISE 𝑊 1 𝑊 2 𝑊 3 bands (Wright  \net al. 2010), and fit the resulting SED with a single-star model. We set an uncertainty floor of 0.02 mag to allow for photometric calibration issues and imperfect models. We predict bandpassintegrated magnitudes using empirically-calibrated model spectral from the BaSeL library (v2.2; Lejeune et al. 1997, 1998). We assume a Cardelli et al. (1989) extinction law with total-to-selective extinction ratio 𝑅 𝑉 = 3 GLYPH<149> 1, and we adopt a prior on the reddening 𝐸 ' 𝐵 GLYPH<0> 𝑉 ' = 0 GLYPH<149> 30 GLYPH<6> 0 GLYPH<149> 03 based on the Green et al. (2019) 3D dust map. We use pystellibs 2 to interpolate between model SEDs, and pyphot 3 to calculate synthetic photometry. We then fit the SED using emcee (Foreman-Mackey et al. 2013) to sample from the posterior, with the temperature, radius, metallicity, and reddening sampled as free parameters.  \nTheresults are shown in the upper right panel of Figure 1. A singlestar model yields an excellent fit to the data, with 𝜒 2 GLYPH<157> 𝑁 data = 0 GLYPH<149> 99, where 𝑁 data is the number of photometric points. The inferred temperature and radius correspond to a solar-type star near the main sequence, and evolutionary models then predict a mass of 𝑀 ★ GLYPH<25> 0 GLYPH<149> 93 GLYPH<6> 0 GLYPH<149> 05 𝑀 GLYPH<12> . The inferred mass is consistent between MIST (Choi et al. 2016) and PARSEC (Marigo et al. 2017) models within 0 GLYPH<149> 02 𝑀 GLYPH<12> . The fact that the source falls near the main sequence provides independent confirmation that the distance inferred from the Gaia astrometry is not catastrophically wrong, and suggest that there is no bright companion.  \nWe inspected the ASAS-SN 𝑉 GLYPH<0> and 𝑔 GLYPH<0> band light curves of the source (Kochanek et al. 2017), which contain 3300 photometric epochs over a 10-year period, with a typical uncertainty of 0.03 mag. This did not reveal any significant periodic or long-term photometric variability. The photometry from ZTF (Bellm et al. 2019), which is more precise but has a shorter baseline, also did not reveal any significant variability.", '3.1 Archival data': 'Gaia BH1 was observed twice by the LAMOST survey (Cui et al. 2012), which obtained low-resolution ( 𝑅 GLYPH<25> 1800) spectra in 2017 and 2019. The LAMOST spectra revealed a main-sequence G star with reported 𝑇 eff = 5863, log 𝑔 = 4 GLYPH<149> 36, and » Fe GLYPH<157> H … = GLYPH<0> 0 GLYPH<149> 29. The two RVs measured by LAMOST were 20 GLYPH<149> 0 GLYPH<6> 4 GLYPH<149> 0 and 8 GLYPH<149> 9 GLYPH<6> 5 GLYPH<149> 6 km s GLYPH<0> 1 , both not far from the mean RV of 23 GLYPH<149> 0 GLYPH<6> 2 GLYPH<149> 6 reported in Gaia DR3, 4 and consistent with no RV variation at all. However, the LAMOST observations both happened to occur at times when the astrometric orbit predicted the luminous star would be near apastron (Figure 2), and thus did not rule out the Gaia astrometric solution. Analysis of the Gaia scanning law (Appendix B) showed that most of the Gaia observations of the source also occurred near apastron. We thus initiated a spectroscopic follow-up campaign.  \n- 2 https://mfouesneau.github.io/pystellibs/\n- 3 https://mfouesneau.github.io/pyphot/', '3.2 Spectroscopic follow-up': "We obtained follow-up spectra using several instruments. Details about the observing setup, data reduction, and analysis for each instrument are listed in Appendix A. The first follow-up observation yielded an RV of GLYPH<25> 64 km s GLYPH<0> 1 , which was clearly different from the LAMOST and Gaia RVs and suggested that the Gaia astrometric solution might be correct. We obtained 39 spectra over the course of 4 months, as described in Appendix A and summarized in Figure 2. The follow-up RVs span most of the orbit's predicted dynamic range and broadly validate the Gaia solution.  \nWe measured RVs via cross-correlation with a synthetic template spectrum, which we took from the BOSZ library of Kurucz model spectra (Kurucz 1970, 1979, 1993) computed by Bohlin et al. (2017). The instruments we used have spectral resolutions ranging from 𝑅 GLYPH<25> 4300 to 𝑅 GLYPH<25> 55000, and in most cases the data have high SNR ( & 20 per pixel; see Table D1). Because we used several different instruments, with different line spread functions and wavelength coverage, the RV uncertainty is dominated by uncertainties in wavelength calibration and zeropoint offsets between different spectrographs. When possible, we performed flexure corrections using sky lines to minimize such offsets, and we used telluric and interstellar absorption lines to verify stability and consistency of the wavelength solutions.  \nIncluding calibration uncertainties, the per-epoch RV uncertainties range from 0.1-4 km s GLYPH<0> 1 . The RV zeropoint is set by telluric wavelength calibration of the Keck/HIRES spectra as described by Chubak et al. (2012). This brings the RV zeropoint to the Nidever et al. (2002) scale within GLYPH<24> 0 GLYPH<149> 1 km s GLYPH<0> 1 . Together with the two archival LAMOST RVs, our follow-up RVs are sufficient to fully constrain the orbit, even without including information from the astrometric solution (Section 4.3).", '4 COMPANION MASS': 'We explore three different approaches to constraining the orbit and companion mass: (a) simultaneously fitting the RVs and Gaia astrometric constraints, (b) using only the Gaia astrometric constraints, and (c) fitting only the RVs.', '4.1 Joint astrometric + RV orbit fitting': "The Gaia orbital solution is parameterized as joint constraints on 12 astrometric parameters: the 5 standard parameters for single-star solutions ( ra, dec, pmra, pmdec, parallax ) and 7 additional parameters describing the photocenter ellipse ( period, t\\_periastron, eccentricity, a\\_thiele\\_innes, b\\_thiele\\_innes, f\\_thiele\\_innes, g\\_thiele\\_innes ). The Thiele-Innes coefficients describe the ellipse orientation and are transformations of the standard Campbell orbital elements (e.g. Halbwachs et al. 2022). Individual-epoch astrometric measurements are not published in DR3.  \nWerefertothe vector of best-fit astrometric parameters constrained by Gaia as 𝝁 ast , and to the corresponding covariance matrix as 𝚺 ast. The latter is constructed from the corr\\_vec column in the gaiadr3.nss\\_two\\_body\\_orbit table. We include in our joint fit the 5 standard astrometric parameters as well as the period, eccentricity, inclination, angle of the ascending node Ω , argument of periastron 𝜔 , periastron time, center-of-mass velocity, luminous star mass 𝑀 ★ , and companion mass 𝑀 2 . For each call of the likelihood function,  \nDEC  \nFigure 3 shows the predicted astrometric orbit of the photocenter, which we calculate using pystrometry (Sahlmann 2019). The left panel shows the predicted motion of the source on the plane of the sky over a 6-year period. The black line (which is actually many narrow lines representing independent posterior samples) shows the  \n<!-- image -->  \n-  \nFigure 1. Properties of the luminous star. Upper left shows a 45 arcsecond wide 𝑧 / 𝑔 PanSTARRS postage stamp centered on the source. Upper right shows the broadband SED and a single-star fit. Red points show observed photometry; open black squares show the integrated fluxes predicted for the model spectrum (black line). Lower left compares the source to the rest of the Gaia astrometric binary sample on the color-magnitude diagram. Lower right compares the temperature and radius measured from the SED to single-star MIST evolutionary models with [Fe/H] = -0.2 (as measured from spectroscopy), which suggest a mass of 𝑀 ★ GLYPH<25> 0 GLYPH<149> 93 𝑀 GLYPH<12> .  \nwe then predict the corresponding vector of astrometric quantities, 𝜽 ast, and corresponding likelihood:  \nln 𝐿 ast = GLYPH<0> 1 2 GLYPH<0> 𝜽 ast GLYPH<0> 𝝁 ast GLYPH<1> | 𝚺 GLYPH<0> 1 ast GLYPH<0> 𝜽 ast GLYPH<0> 𝝁 ast GLYPH<1> GLYPH<149> (1)  \nWe neglect terms in the likelihood function that are independent of 𝜽 ast. The 5 single-star parameters for our purposes are nuisance parameters (they do not constrain the companion mass), but we include them in our fit as free parameters in order to properly account for covariances in the parameters constrained by Gaia . When predicting the parameters describing the photocenter ellipse, we assume the companion is dark. We predict the Thiele-Innes coefficients using the standard transformations (e.g. Binnendĳk 1960).  \nWe additionally predict the RVs of the luminous star at the array of times at which spectra were obtained, 𝑡 𝑖 . This leads to a radial velocity term in the likelihood,  \nln 𝐿 RVs = GLYPH<0> 1 2 ∑︁ 𝑖 GLYPH<16> RV pred ' 𝑡 𝑖 ' GLYPH<0> RV 𝑖 GLYPH<17> 2 𝜎 2 RV GLYPH<150> 𝑖 GLYPH<150> (2)  \nwhere RV 𝑖 and RV pred ' 𝑡 𝑖 ' are the measured and predicted RVs. The full likelihood is then  \nln 𝐿 = ln 𝐿 ast , ln 𝐿 RVs GLYPH<149> (3)  \nA more optimal approach would be to fit the epoch-level astrometric data and the RVs simultaneously, but this will only become  \npossible after epoch-level astrometric data are published in DR4. Since the Gaia astrometric fits are essentially pure likelihood constraints (i.e., they are not calculated with informative priors), it is possible to mix individual-epoch RV measurements and the Gaia astrometric constraints, without risk of multiplying priors.  \nWe use flat priors on all parameters except 𝑀 ★ , for which we use a normal distribution informed by isochrones, 𝑀 ★ GLYPH<157> 𝑀 GLYPH<12> GLYPH<24> N' 0 GLYPH<149> 93 GLYPH<150> 0 GLYPH<149> 05 ' . We sample from the posterior using emcee (ForemanMackey et al. 2013) with 64 walkers, taking 3000 steps per walker after a burn-in period of 3000 steps. The results of this fitting are visualized in Figures 2, 3, and C2. We infer an unseen companion mass 𝑀 2 = 9 GLYPH<149> 62 GLYPH<6> 0 GLYPH<149> 18 𝑀 GLYPH<12> .  \nFigure 2 compares the measured RVs to draws from the posterior. The scatter across posterior samples at fixed time is diagnostic of the uncertainty in the RV solution. Unsurprisingly, uncertainties are larger during phases where there are fewer measured RVs. Overall, the fit is good: we obtain a solution that both matches the observed RVs and predicts an astrometric orbit consistent with the Gaia constraints. This did not have to occur, and it speaks to the quality of the astrometric solution. Our constraints from the joint fit are listed in Table 1. Following the Gaia convention, the periastron time 𝑇 𝑝 is reported with respect to epoch 2016.0; i.e., JD 2457389.0.  \nFigure 2. Radial velocities. Points with error bars are measurements; gray lines are draws from the posterior when jointly fitting these RVs and the Gaia astrometric constraints. Top panel shows all available RVs, including observations by the LAMOST survey in 2017 and 2019; bottom panel highlights our follow-up in 2022. The best-fit solution has a period of 186 days, eccentricity 0.45, and RV semi-amplitude of 67 km s GLYPH<0> 1 . Together with the inclination constraint from astrometry, this implies a companion mass of 9 GLYPH<149> 62 GLYPH<6> 0 GLYPH<149> 18 𝑀 GLYPH<12> .  \n<!-- image -->  \nFigure 3. Left panel shows predicted motion of Gaia BH1's photocenter on the sky over a 6-year window. Middle and right panels show the predicted orbit with parallax and proper motion removed. Middle panel shows joint constraints from RVs and astrometry; right panel shows constraints from the Gaia astrometric solution alone.  \n<!-- image -->  \ntotal motion due to proper motion, parallax, and orbital motion. The dotted blue line isolates the contribution due to parallax and proper motion alone. The middle and right panels show the predicted orbital motion, with parallax and proper motion removed, based on the joint astrometry + RV fit (middle) or the pure astrometry fit (right; see below). The orbits predicted in the two cases are similar, but the fit with RVs included is better constrained, and has a slightly smaller photocenter semi-major axis.", '4.2 Purely astrometric constraints': 'We also explored fitting the orbit without the constraints from RVs; i.e., simply removing the ln 𝐿 RVs term from Equation 3. The results of this exercise are described in Appendix C, where we also compare the RVs predicted by the astrometry-only solution to the measured RVs. This comparison allows us to assess the reliability of the astrometric solution.  \nIn brief, we find that the measured RVs are consistent with the purely astrometric solution within its uncertainties. The joint RVs+astrometry fit yields a marginally lower companion mass and more face-on orbit than the best-fit pure astrometry solution, but the astrometric parameter constraints from our joint Gaia +RVs fit are consistent with the purely astrometric constraints at the 1.6 𝜎 level (Table C1). We also report the purely astrometric constraints in Table 1. The joint RVs+astrometry constraints are much tighter than those from astrometry alone.', '4.2.1 Astrometric uncertainties': 'Although we find no indication that the astrometric solution is unreliable, there is one way in which it is unusual: the uncertainties on the parameters describing the photocenter ellipse are significantly larger than is typical for a source with 𝐺 = 13 GLYPH<149> 77. This is reflected in the uncertainty on the photocenter ellipse semi-major axis, 𝜎 𝑎 0 = 0 GLYPH<149> 22 mas. The median uncertainty for sources with 13 GLYPH<149> 6 GLYPH<159> 𝐺 GLYPH<159> 13 GLYPH<149> 9 in the nss\\_two\\_body\\_orbit table is only 𝜎 𝑎 0 = 0 GLYPH<149> 026 mas, and only 1.4% of sources in that magnitude range have larger 𝜎 𝑎 0 than Gaia BH1.  \nThe unusually large 𝜎 𝑎 0 appears to be a result of two factors. First, astrometric solutions with larger 𝑎 0 have larger 𝜎 𝑎 0 (see Appendix G); the median 𝜎 𝑎 0 for sources in the same magnitude range with 2 GLYPH<149> 5 GLYPH<20> 𝑎 0 GLYPH<157> mas GLYPH<20> 3 GLYPH<149> 5 is 0.064 mas, and Gaia BH1 is only in the 85th percentile of 𝜎 𝑎 0 for such sources. Second, as we discuss in Appendix B, all Gaia observations of the source that contributed to the DR3 solution covered the same GLYPH<25> 50% of the orbit, even though they were spread over more than 5 orbital periods. Despite this, the combination of astrometry and RVs constrains the orbit well, and leads to a tighter constraint on 𝑎 0 than is achieved for typical Gaia sources at the same apparent magnitude. In any case, a catastrophic problem with the astrometric solution is firmly ruled out by the RVs, which are consistent with the astrometric solution and require a BH-mass companion even without the astrometric constraints.', '4.3 Solution based on radial velocities alone': "Although the Gaia astrometric solution provides strong constraints on the binary's orbit, it is useful to consider constraints based only on the measured RVs, which are insensitive to any possible problems with the Gaia solution. To this end, we fit the RVs with a model with the standard 6 orbital parameters for a single-lined binary. We first  \nTable 1. Physical parameters and 1 𝜎 uncertainties for both components of Gaia BH1. The inclination, 𝑖 GLYPH<25> 127 deg, is equivalent to 𝑖 GLYPH<25> 53 deg for a spectroscopic orbit. We compare constraints on the orbit based on both astrometry and RVs (3rd block), astrometry alone (4th block), and RVs alone (5th block).  \nsearch for the maximum-likelihood solution with simulated annealing (e.g. Iglesias-Marzoa et al. 2015) and then use MCMC to explore the posterior in the vicinity of the maximum likelihood solution. The dense sampling of RVs in the last 2 months of our observations allows us to robustly constrain the orbit. We first fit the RVs with no period constraint; this yielded a solution with 𝑃 orb = 184 GLYPH<149> 8 GLYPH<6> 0 GLYPH<149> 7 days. Since this is consistent with the Gaia solution, but the Gaia  \nsolution is based on a longer time baseline, we then fixed the period to 𝑃 orb = 185 GLYPH<149> 6 days, as inferred in Section 4.1. The constraints on the RV solution obtained in this way are reported in Table 1, and the best-fit RV curve is shown in Figure 4. The reduced 𝜒 2 is 0.54, suggesting that the fit is good and the RV uncertainties are overestimated on average.  \nFitting only the RVs yields an RV semi-amplitude and eccentricity consistent with the predictions of the purely astrometric and astrometric + RV fits. The main difference is that the inclination is unknown in the pure-RV solution, and the uncertainties on other parameters are somewhat larger. The periastron time inferred in this fit is 13 orbital cycles later than the one reported in the Gaia solution, as it corresponds to the periastron passage in August 2022 that is actually covered by our RVs.  \nThe constraint on the companion mass from RVs can be expressed in terms of the spectroscopic mass function,  \n𝑓 ' 𝑀 2 ' spec = 𝑀 2 GLYPH<18> 𝑀 2 𝑀 2 , 𝑀 ★ GLYPH<19> 2 sin 3 𝑖 = 𝑃 orb 𝐾 3 ★ 2 𝜋𝐺 GLYPH<16> 1 GLYPH<0> 𝑒 2 GLYPH<17> 3 GLYPH<157> 2 GLYPH<149> (4)  \n𝑃 orb , 𝐾 ★ , and 𝑒 are inferred directly from the RVs. Since both the terms GLYPH<16> 𝑀 2 𝑀 2 , 𝑀 ★ GLYPH<17> 2 and sin 3 𝑖 must be less than one, 𝑓 ' 𝑀 2 ' spec = 4 GLYPH<149> 08 GLYPH<6> 0 GLYPH<149> 08 sets an absolute lower limit on the mass of the unseen companion. This constraint is illustrated in the lower panel of Figure 4, where we plot the best-fit orbit and residuals, and the companion mass required to explain the mass function for different inclinations and 𝑀 ★ .  \nThe figure shows that the RVs cannot accommodate a companion mass below 4 𝑀 GLYPH<12> , independent of the astrometric solution and the mass of the G star. When a plausible mass for the G star is assumed, this limit increases to 5 GLYPH<149> 5 𝑀 GLYPH<12> . We also show the inclination implied by the Gaia orbital solution. Assuming this inclination and 𝑀 ★ = 0 GLYPH<149> 93 GLYPH<6> 0 GLYPH<149> 05 𝑀 GLYPH<12> leads to a companion mass of 9 GLYPH<149> 50 GLYPH<6> 0 GLYPH<149> 20 𝑀 GLYPH<12> , in good agreement with the value we inferred by simultaneously fitting the RVs and astrometry. The fact that this value is not identical to the 𝑀 2 = 9 GLYPH<149> 62 GLYPH<6> 0 GLYPH<149> 18 𝑀 GLYPH<12> inferred in the joint fit reflects the fact that the joint RV + astrometry fit implies a slightly larger 𝐾 ★ .", '4.3.1 Astrometric mass function': "To appreciate the additional information provided by the astrometric solution, we can consider the astrometric mass function,  \n𝑓 ' 𝑀 2 ' ast = GLYPH<16> 𝑎 0 𝜛 GLYPH<17> 3 GLYPH<18> 𝑃 orb yr GLYPH<19> GLYPH<0> 2 = 𝑀 2 GLYPH<18> 𝑀 2 𝑀 ★ , 𝑀 2 GLYPH<19> 2 (5)  \nwhere the 2 nd equality holds only for a dark companion. The quantity after the 2 nd equality is equivalent to 𝑓 ' 𝑀 2 ' spec except for a factor of sin 3 𝑖 . In the limit of small uncertainties, 𝑓 ' 𝑀 2 ' ast sets a dynamical limit on the mass of the companion that is more stringent than the spectroscopic mass function. In particular, 𝑓 ' 𝑀 2 ' ast is the mass that the orbit would imply if the G star were a massless test particle.  \nThe constraints from astrometry alone yield 𝑓 ' 𝑀 2 ' ast = 11 GLYPH<149> 2 GLYPH<6> 1 GLYPH<149> 9 𝑀 GLYPH<12> . Those from joint fitting of the RVs and astrometry yield 𝑓 ' 𝑀 2 ' ast = 8 GLYPH<149> 00 GLYPH<6> 0 GLYPH<149> 16 𝑀 GLYPH<12> . When we adopt 𝑀 ★ = 0 GLYPH<149> 93 GLYPH<6> 0 GLYPH<149> 05 𝑀 GLYPH<12> and solve for 𝑀 2 , the joint constraint yields 𝑀 2 = 9 GLYPH<149> 62 GLYPH<6> 0 GLYPH<149> 18 𝑀 GLYPH<12> .", '4.3.2 Possibility of underestimated uncertainties': "It is worth considering whether the astrometric uncertainties could be underestimated. Given that the Gaia single-star astrometric solutions have been shown to have uncertainties underestimated by GLYPH<24> 30%near  \nFigure 4. Top panel shows the best-fit orbital solution to the RVs, ignoring the astrometric solution. We show two periods; duplicated datapoints are transparent. Bottom panel shows the resulting constraints on the companion's mass. Black line assumes a mass of 0 GLYPH<149> 93 𝑀 GLYPH<12> for the luminous star as implied by the isochrones; red dashed line assumes (pathologically) that the star is a massless test particle. Shaded region shows the Gaia astrometric inclination constraint, which corresponds to a companion mass of 9 GLYPH<149> 50 GLYPH<6> 0 GLYPH<149> 20 𝑀 GLYPH<12> . Even if the Gaia astrometric solution were completely wrong and the luminous star had negligible mass, the implied mass of the companion would be ¡ 4 𝑀 GLYPH<12> , which is inconsistent with the observed SED and spectra for any luminous companion.  \n<!-- image -->  \n𝐺 = 13 GLYPH<149> 8 (e.g. El-Badry et al. 2021), and there is modest ( GLYPH<24> 1 GLYPH<149> 6 𝜎 ) tension between the joint fit and the pure-astrometry solution (Appendix C) this is not implausible. To explore the possible effects of underestimated astrometric uncertainties on our constraints, we multiplied all the astrometric uncertainties by 2, constructed the covariance matrix assuming the same correlation matrix (this is equivalent to multiplying the covariance matrix by 4), and repeated the joint fit of astrometry and RVs. This yielded a companion mass constraint of 𝑀 2 = 9 GLYPH<149> 54 GLYPH<6> 0 GLYPH<149> 26 𝑀 GLYPH<12> - simular to our fiducial result, and still a rather tight constraint.", '5 ANALYSIS OF THE G STAR SPECTRUM': "We analyzed the spectrum of the G star using several different methods. Most of our analysis used the Keck HIRES spectrum obtained on JD 2459791 using the standard California Planet Survey setup (CPS; Howard et al. 2010). We first fit the spectrum using the empirical SpecMatch-Emp code (Yee et al. 2017), which compares the continuum-normalized rest-frame spectrum to a library of HIRES spectra of FGKM dwarfs and subgiants observed with the standard CPS setup and analyzed with traditional spectroscopic methods. It estimates uncertainties based on the dispersion in parame-  \nTable 2. Parameters of the G star inferred by BACCHUS. 𝑁 lines is the number of independent absorption lines used to infer each element.  \nters of objects in the library with spectra most similar to an observed spectrum. SpecMatch-Emp yielded an effective temperature 𝑇 eff = 5801 GLYPH<6> 110 K and metallicity » Fe GLYPH<157> H … = GLYPH<0> 0 GLYPH<149> 34 GLYPH<6> 0 GLYPH<149> 09.  \nIn Figure 5, we compare the HIRES spectrum of Gaia BH1 to that of its closest match in the library, HD 33636. The parameters of that star, as measured by Valenti & Fischer (2005) from the HIRES spectrum and synthetic spectral models, are 𝑇 eff = 5904 K, log 𝑔 = 4 GLYPH<149> 43, 𝑣 sin 𝑖 = 3 GLYPH<149> 1 km s GLYPH<0> 1 , and » Fe GLYPH<157> H … = GLYPH<0> 0 GLYPH<149> 13. The region of the spectrum shown is centered on the Mg I b triplet and also contains strong lines of the elements Fe, Ni, Ti, Cr, Cu, Ca, and Y, many of which are labeled. The two spectra are strikingly similar. This simple comparison suggests that the luminous star in Gaia BH1 is a rather unremarkable solar-type star. Quantitative analysis of the spectrum yielded similar conclusions, as described below.  \nWealso fit the spectrum using the HIRES-adapted implementation of the Cannon (Ness et al. 2015) that was trained by Rice & Brewer (2020). TheCannon uses a library of spectra with known stellar labels (i.e., atmospheric parameters and abundances) to build a spectral model that smoothly interpolates between spectra, predicting the normalized flux at a given wavelength as a polynomial function of labels. This data-driven spectral model is then used to fit for the labels that can best reproduce an observed spectrum through a standard maximum-likelihood method. The Cannon fit yielded parameters 𝑇 eff = 5863 K, log 𝑔 = 4 GLYPH<149> 43, 𝑣 sin 𝑖 = 1 GLYPH<149> 34 km s GLYPH<0> 1 , and » Fe GLYPH<157> H … = GLYPH<0> 0 GLYPH<149> 29, similar to the values inferred by SpecMatch-Emp . It also returned constraints on the abundances of 14 metals beside iron. The full set of labels inferred by the Cannon can be found in Appendix F.  \nTo investigate constraints on the G star's projected rotation velocity 𝑣 sin 𝑖 , we convolved an 𝑅 = 300 GLYPH<150> 000 Kurucz spectrum from the BOSZ library with a range of rotational broadening kernels and with the 𝑅 GLYPH<25> 55 GLYPH<150> 000 HIRES instrumental broadening kernel, which we model as a Gaussian. We find that we can infer a robust upper limit  \nof 𝑣 sin 𝑖 GLYPH<159> 3 GLYPH<149> 5 km s GLYPH<0> 1 : a higher value would broaden the lines more than observed. However, the actual value of 𝑣 sin 𝑖 is not well constrained because at lower values of 𝑣 sin 𝑖 , instrumental broadening dominates, and the adopted value of 𝑣 sin 𝑖 has little effect on the predicted line profiles. For 𝑅 = 1 𝑅 GLYPH<12> , this translates to a rotation period 𝑃 rot & 14 GLYPH<2> ' sin 𝑖 ' days. Assuming sin 𝑖 GLYPH<24> 1 and using a canonical gyrochronological scaling relation for solar-type stars (e.g. Skumanich 1972; Mamajek & Hillenbrand 2008), this implies an age of & 1 Gyr. Application of such scaling relations is of course only appropriate if the G star was not previous tidally synchronized and spun down by interaction with the companion.", '5.1 Detailed abundances': "We fit the HIRES spectrum using the Brussels Automatic Code for Characterizing High accUracy Spectra (BACCHUS; Masseron et al. 2016) with the same set up as in Hawkins et al. (2020b). BACCHUS enables us to derive the stellar atmospheric parameters, including the effective temperature ( 𝑇 eff ), surface gravity (log 𝑔 ), metallicity ([Fe/H]) and microturblent velocity ( 𝑣 micro ) by assuming Fe excitation/ionization balance; i.e., the requirement that lines with different excitation potentials all imply the same abundances. For the set up of BACCHUS,wemakeuseofthe fifth version of the Gaia-ESO atomic linelist (Heiter et al. 2021). Hyperfine structure splitting is included for Sc I, V I Mn I, Co I, Cu I, Ba II, Eu II, La II, Pr II, Nd II, Sm II (see more details in Heiter et al. 2021). We also include molecular line lists for the following species: CH (Masseron et al. 2014), and CN, NH, OH, MgH and C 2 (T. Masseron, private communication). Finally, we also include the SiH molecular line list from the Kurucz linelists 5 . Spectral synthesis for BACCHUS is done using the TURBOSPECTRUM (Alvarez & Plez 1998; Plez 2012) code along with the line lists listed above and the MARCS model atmosphere grid (Gustafsson et al. 2008). The stellar atmospheric parameters derived from BACCHUS can be found in Table 2.  \nOnce the stellar parameters were determined, the model atmosphere is fixed and individual abundances were derived using BACCHUS' 'abund' module. For each spectral absorption feature, this module creates a set of synthetic spectra, which range between -0.6 GLYPH<159> [X/Fe] GLYPH<159> +0.6 dex, and performs a 𝜒 2 minimization between the observed and synthetic spectra. The reported atmospheric [X/Fe] abundances are the median of derived [X/Fe] across all lines for a given species. The uncertainty in the atmospheric [X/Fe] is defined as the dispersion of the [X/H] abundance measured across all lines for a species. If only 1 absorption line is used, we conservatively assume a [X/Fe] uncertainty of 0.10 dex. For a more detailed discussion of the BACCHUS code, we refer the reader to Section 3 of both Hawkins et al. (2020a,b). We ran BACCHUS on the full HIRES spectrum, after merging individual de-blazed orders and performing a preliminary continuum normalization using a polynomial spline. Further normalization is performed by BACCHUS. The resulting stellar parameters and abundances are listed in Table 2.  \nOverall, the abundance pattern of the G star is typical of the thin disk. When comparing the star's position in the [X/Fe] vs. [Fe/H] plane to the population of stars in the solar neighborhood with abundances measured by Bensby et al. (2014), it falls within the GLYPH<24> 1 𝜎 observed scatter for most elements, and within 2 𝜎 for all. Given that the G star is expected to have a thin convective envelope (making mixing of accreted material into the interior inefficient), this suggests that it did not suffer much pollution from its companion.  \nFigure 5. Comparison of the HIRES spectrum of Gaia BH1 (gold) and the standard star HD 33636 (blue). Black line shows the difference between the two spectra. Cutout is centered on the Mg I b triplet; other identifiable lines are labeled. The two spectra are very similar, indicating that the surface properties and abundances of the G star are normal for thin-disk stars in the solar neighborhood.  \n<!-- image -->  \nThe measured lithium abundance, 𝐴 ' Li ' = 2 GLYPH<149> 3 GLYPH<6> 0 GLYPH<149> 1, is also not unusual. For solar-type stars, lithium abundance can be used as an age indicator, because the surface lithium abundance is depleted over time (e.g. Skumanich 1972; Baumann et al. 2010). The age 𝐴 ' Li ' correlation is well-studied for Sun-like stars (e.g. Ramírez et al. 2012; Carlos et al. 2016), and for solar twins, 𝐴 ' Li ' = 2 GLYPH<149> 3 corresponds to an age of about 1 Gyr. However, 𝐴 ' Li ' varies strongly with both 𝑇 eff and [Fe/H] (higher 𝑇 eff and lower [Fe/H] both produce a thinner convective envelope and thus imply an older age at fixed A(Li)), so caution should be taken in applying the relation from solar twins to a star somewhat warmer and more metal-poor than the sun. More robustly, the measured abundance allows us to rule out youth: the star falls well below the Hyades (age GLYPH<25> 600 Myr) in the 𝐴 ' Li ' vs. 𝑇 eff plane (Takeda et al. 2013), so we can confidently rule out ages younger than this. Comparison of the G star's effective temperature and radius to isochrones (Figure 1) suggests an age of & 4 Gyr.  \nWe find no evidence for pollution of the G star's photosphere by 𝛼 -elements synthesized during the companion's death, as has been reported in some BH companions (e.g. Israelian et al. 1999; González Hernández et al. 2011). The barium abundance is slightly enhanced compared to the solar value, with » Ba GLYPH<157> Fe … = 0 GLYPH<149> 17 GLYPH<6> 0 GLYPH<149> 04. Strong barium enhancement is often interpreted as a result of accretion of the AGB wind of a binary companion ('barium stars'; McClure & Woodsworth 1990), usually in stars with white dwarf companions. However, the [Ba/Fe] we measure in Gaia BH1 is not high enough to qualify it as a barium star and is probably unrelated to the companion. We also do not find significant enrichment of other 𝑠 GLYPH<0> process elements, or of the 𝑟 GLYPH<0> process element europium.", '5.1.1 Limits on a luminous companion': 'As part of the CPS pipeline, the HIRES spectrum of Gaia BH1 was analyzed with the Reamatch tool (Kolbl et al. 2015), which searches for a second peak in the CCF after subtracting the best-fit template. For Gaia BH1, this yielded a null result, with a single narrow peak in the CCF and no evidence of another luminous star. Kolbl et al. (2015) found that for luminous binaries with typical solar-type primaries, reamatch has a ¡ 90%detection rate for companions that contribute ¡ 3% of the light in the optical in a single-epoch HIRES spectrum. These limits of course depend on the spectral type and rotation rate of the secondary, and a rapidly rotating secondary would be harder to detect. However, we find that even a secondary that contributes pure continuum cannot contribute more than 10% of the light in the HIRES spectrum, or it becomes impossible to achieve a good match to the observed line depths with SpecMatch-Emp . This can be recognized simply from the depth of the observed absorption lines (Figure 5): several lines have depths that reach 10% of the continuum flux. There cannot be another source contributing more than GLYPH<25> 10% of the light at these wavelengths, or the spectrum of the G star would have to reach negative fluxes in the absorption lines to produced the observed total spectrum.  \nWe also find no evidence of a luminous companion at longer or shorter wavelengths. Spectra from X-shooter provide good SNR over the full optical and NIR wavelength range, from 3100 to 24000 Å, and are consistent with a single G star. This and the GALEX + WISE photometry allow us to rule out pathological companions that are hot and small or cool and large.', '6 GALACTIC ORBIT': "To investigate the past trajectory of Gaia BH1, we used its parallax and proper motion from the Gaia astrometric binary solution, as well as the center-of-mass radial velocity inferred from the spectroscopic fit, as starting points to compute its orbit backward in time for 500 Myr using galpy (Bovy 2015). We used the Milky Way potential from McMillan (2017). The result is shown in Figure 6; for comparison, we also show the orbit of the Sun. The orbit appears typical of a few-Gyr old thin-disk star, with excursions above the disk midplane limited to GLYPH<6> 250 pc. The orbit also never comes very close to the Galactic center, and is not aligned with the orbits of any globular clusters (GCs). These conclusion are insensitive to the assumed Galactic potential and uncertainties in the source's proper motion or center-of-mass RV.", '7 X-RAY AND RADIO UPPER LIMITS': 'We checked archival all-sky X-ray and radio surveys to place upper limits on the flux from Gaia BH1. We used the second ROSAT all-sky survey (2RXS) catalog (Boller et al. 2016) and the Very Large Array sky survey (VLASS) Epoch 1 CIRADA catalog (Lacy et al. 2020; Gordon et al. 2020), respectively. At the position of Gaia BH1, the nearest ROSAT source (unassociated with any optical counterpart) is 2RXS J172552.5-003516, located 26\'4".1 away. The median separation between sources in ROSAT and SDSS sources is about 12" (Agüeros et al. 2009), which would place Gaia BH1 130 𝜎 away if a true association. Therefore, assuming no ROSAT detection, we place upper limits on the flux of Gaia BH1 to be the limit of the 2RXS catalog, GLYPH<24> 10 GLYPH<0> 13 erg s GLYPH<0> 1 cm GLYPH<0> 2 (Boller et al. 2016). This corresponds to a luminosity upper limit of 𝐿 X . 3 GLYPH<2> 10 30 erg s GLYPH<0> 1 in the 0.1-2.4 keV range.  \nIn the radio, the nearest VLASS source (also unassociated with an optical counterpart) is VLASS1QLCIR J172823.79-003035.0, located 6\'5".3 away. The beam size of VLASS is 2".5 (Lacy et al. 2020), which would place Gaia BH1 146 𝜎 away. Assuming no VLASS detection, we can place upper limits on the flux of Gaia BH1 to be the 3 𝜎 limit of a single epoch of the VLASS catalog, GLYPH<24> 0 GLYPH<149> 36 mJy, which in the 2-4 GHz range (S-band) corresponds to a flux limit of . 7 GLYPH<2> 10 GLYPH<0> 18 erg s GLYPH<0> 1 cm GLYPH<0> 2 (Lacy et al. 2020). This corresponds to a radio luminosity upper limit of 𝐿 R . 2 GLYPH<2> 10 26 erg s GLYPH<0> 1 in the S-band.', '7.1 Could Gaia BH1 be detected with X-ray or radio observations?': "We briefly consider whether an X-ray or radio detection is expected. A rough estimate of the expected accretion rate onto the BH can be obtained under the assumption that a spherically-symmetric wind from the G star is accreted at the Hoyle-Lyttleton rate:  \n/ 𝑀 acc = 𝐺 2 𝑀 2 BH / 𝑀 wind 𝑣 4 wind 𝑎 2 = 6 GLYPH<2> 10 GLYPH<0> 18 𝑀 GLYPH<12> yr GLYPH<0> 1 GLYPH<18> 𝑀 BH 10 𝑀 GLYPH<12> GLYPH<19> 2 GLYPH<18> / 𝑀 wind 10 GLYPH<0> 14 𝑀 GLYPH<12> yr GLYPH<0> 1 GLYPH<19> GLYPH<2> GLYPH<18> 𝑣 wind 600 km s GLYPH<0> 1 GLYPH<19> GLYPH<0> 4 GLYPH<16> 𝑎 1 au GLYPH<17> GLYPH<0> 2 GLYPH<149> (6)  \nHere / 𝑀 wind is the wind mass-loss rate of the G star and 𝑣 wind is the wind's velocity. The predicted / 𝑀 acc is very low: less than 10 GLYPH<0> 10 / 𝑀 edd , where / 𝑀 edd GLYPH<24> 2 GLYPH<2> 10 GLYPH<0> 7 𝑀 GLYPH<12> yr GLYPH<0> 1 ' 𝑀 BH GLYPH<157> » 10 𝑀 GLYPH<12> …' is the Eddington  \nrate with 10% efficiency. Assuming accretion onto the BH results in an X-ray luminosity 𝐿 X = 𝜂 / 𝑀 acc 𝑐 2 , the expected flux at Earth is  \n𝐹 X = 𝜂 / 𝑀𝑐 2 4 𝜋𝑑 2 = 1 GLYPH<149> 1 GLYPH<2> 10 GLYPH<0> 15 erg s GLYPH<0> 1 cm GLYPH<0> 2 GLYPH<16> 𝜂 0 GLYPH<149> 1 GLYPH<17> GLYPH<18> 𝑀 BH 10 𝑀 GLYPH<12> GLYPH<19> 2 GLYPH<18> / 𝑀 wind 10 GLYPH<0> 14 𝑀 GLYPH<12> yr GLYPH<0> 1 GLYPH<19> GLYPH<2> GLYPH<18> 𝑣 wind 600 km s GLYPH<0> 1 GLYPH<19> GLYPH<0> 4 GLYPH<16> 𝑎 1 au GLYPH<17> GLYPH<0> 2 GLYPH<18> 𝑑 500 pc GLYPH<19> GLYPH<0> 2 GLYPH<149> (7)  \nFor radiatively efficient accretion with 𝜂 GLYPH<24> 0 GLYPH<149> 1, this is within the flux limits of a deep Chandra ACIS observation, which are of order 10 GLYPH<0> 16 erg s GLYPH<0> 1 cm GLYPH<0> 2 . Unfortunately, models for advection-dominated accretion flows predict very low radiative efficiencies at the relevant accretion rates (e.g. Narayan & Yi 1995; Quataert & Narayan 1999). For example, extrapolating the models of Sharma et al. (2007) yields an expected 𝜂 GLYPH<24> 10 GLYPH<0> 6 , corresponding to 𝐹 X GLYPH<24> 10 GLYPH<0> 20 erg s GLYPH<0> 1 cm GLYPH<0> 2 . We conclude that an X-ray detection would be surprising, though it cannot hurt to look.  \nThe expected radio luminosity can be estimated by extrapolating the radio-X-ray correlation for quiescent BH X-ray binaries (e.g. Merloni et al. 2003; Gallo et al. 2006). This predicts a radio flux of order 10 𝜇 Jy at 8 GHz for 𝜂 = 0 GLYPH<149> 1 (below the VLASS limit, but detectable with the VLA), or GLYPH<24> 1 nJ for 𝜂 = 10 GLYPH<0> 6 (not detectable in the near future). Uncertainty in the G star's mass loss rate (which could plausibly fall between 10 GLYPH<0> 14 and 10 GLYPH<0> 13 𝑀 GLYPH<12> yr GLYPH<0> 1 ) and wind velocity near the BH (plausibly between 300 and 600 km s GLYPH<0> 1 ) leads to more than an order of magnitude uncertainty in these estimates, in addition to the uncertainty in 𝜂 .", '8.1 Comparison to recent BH imposters': "ManyrecentdormantBHcandidatesinbinarieshaveturnedoutnotto be BHs, but rather (in most cases) mass-transfer binaries containing undermassive, overluminous stars in short-lived evolutionary phases. Gaia BH1 is different from these systems in several ways.  \nFirst, the evidence for a BH companion does not depend on the mass of the luminous star, as it did in LB-1, HR 6819, and NGC 1850 BH1, which all contained low-mass stripped stars masquerading as more massive main-sequence stars (Shenar et al. 2020; Bodensteiner et al. 2020; El-Badry & Quataert 2021; El-Badry & Burdge 2022). In Gaia BH1, the mass of the luminous star could be zero, and the mass of the companion would still vastly exceed the limit allowed by the SED for normal stellar companion, the Chandrasekhar mass, and the maximum neutron star mass.  \nSecond, the BH nature of the companion does not depend on a low assumed inclination, as was for example the case in NGC 2004 #115 (Lennon et al. 2021; El-Badry et al. 2022a), LB-1, and NGC 1850 BH1. The inclination of Gaia BH1 could be edge-on, and RVs alone would still require the companion to be a BH. In addition, reliable constraints on the inclination are available from the Gaia astrometric solution; such constraints have not been available for any previous candidates.  \nThird, there is no evidence of ongoing or recent mass transfer (i.e., ellipsoidal variability or emission lines from a disk), as there was in all the candidates discussed above except NGC 2004 #115, as well as other recent candidates with red giant primaries (e.g. Jayasinghe et al. 2021; El-Badry et al. 2022b). The orbital period is  \nFigure 6. Galactic orbit of Gaia BH1, calculated backward for 500 Myr from the measured proper motion and center-of-mass RV. For comparison, we show the Sun's orbit calculated over the same period. The orbit is typical of a thin-disk star, ruling out a large natal kick or dynamical formation in a globular cluster.  \n<!-- image -->  \nlong enough, and the luminous star small enough, that there is no plausible evolutionary scenario in which it was recently stripped and either component is in an unusual evolutionary state. A stripped-star scenario is also disfavored by the unremarkable optical spectrum and surface abundances of the G star, and by its measured surface gravity (log » 𝑔 GLYPH<157> ' cms GLYPH<0> 2 '… = 4 GLYPH<149> 55 GLYPH<6> 0 GLYPH<149> 16), which implies a mass 𝑀 ★ = 10 log 𝑔 𝑅 2 GLYPH<157> 𝐺 = 1 GLYPH<149> 27 GLYPH<6> 0 GLYPH<149> 54 𝑀 GLYPH<12> .  \nFinally, the luminous star in Gaia BH1 has a luminosity of only 1 𝐿 GLYPH<12> , which is hundreds to thousands of times fainter than the luminous stars in all the BH imposters discussed above. This makes it vastly harder to hide any luminous companion in the SED: any normal star with a mass exceeding GLYPH<25> 0 GLYPH<149> 4 𝑀 GLYPH<12> would easily be detected in the spectrum.", '8.2 Nature of the unseen companion': "Considering just the RVs and the properties of the G star, it seems incontrovertible that the companion is a ¡ 5 𝑀 GLYPH<12> dark object. Given the good agreement between the observed RVs and those predicted by the astrometric solution, we are also inclined to trust the inclination constraint from the astrometric solution, and thus the most likely companion mass is GLYPH<24> 10 𝑀 GLYPH<12> .  \nThe companion has not been detected in X-rays or at radio wavelengths, and it lacks other observables associated with accreting BHs (e.g. outbursts or rapid flickering). An X-ray or radio detection would be unexpected given the weak stellar wind expected from the G star and wide separation. Constraints on the nature of the companion thus come down to probabilistic arguments conditioned on its inferred mass and low luminosity.  \nThe mass of the companion is consistent, for example, with a single 10 𝑀 GLYPH<12> BH, 2 lower-mass BHs, 5 neutron stars, 10 massive white dwarfs, or 200 brown dwarfs. Scenarios involving more than two objects inside the G star's orbit are difficult to assemble and unlikely to be dynamically stable. As we discuss in Section 8.4, a scenario in which the dark object is a close binary with total mass GLYPH<24> 10 𝑀 GLYPH<12> (containing 2 BHs or a BH and a neutron star) may be plausible. We thus proceed under the assumption that the unseen companion is either a single BH or a binary containing 2 compact objects, at least one of which is a BH.", '8.3 Comparison to other BHs': "Figure 7 compares Gaia BH1 to other stellar-mass BHs with known masses and orbits in the literature. Red and blue points show low- and high-mass X-ray binaries, whose parameters we take from Remillard &McClintock (2006) and the BlackCAT catalog of X-ray transients introduced by Corral-Santana et al. (2016). We take the most recent mass estimate for Cyg X-1 from Miller-Jones et al. (2021). We also show the binaries VFTS 243 (in the LMC; Shenar et al. 2022a) and HD 130298 (in the Milky Way; Mahy et al. 2022), which are both single-lined binaries containing GLYPH<24> 25 𝑀 GLYPH<12> O stars and unseen companions suspected to be BHs. Finally, we show in magenta two binaries in the GC NGC 3201 discovered with MUSE (Giesers et al. 2018, 2019) with suspected BH companions. For the detached systems besides Gaia BH1, the inclination is unknown, so only a lower limit on the BH mass can be inferred. Other detached BH candidates exist in the literature (e.g. Qian et al. 2008; Casares et al. 2014; Khokhlov et al. 2018; Thompson et al. 2019; Gomez & Grindlay 2021), but we do not include them because we consider their status as BHs more uncertain. Similarly, there are other X-ray bright binaries proposed to contain BHs that we do not include (e.g., Cyg X-3, SS 433) because a neutron star is not fully ruled out. We exclude IC 10 X-1 because the BH's mass is very uncertain (Laycock et al. 2015). We note that for many of the X-ray binaries included in the figure, mass estimates across different studies differ by significantly more than the reported uncertainties in individual studies.  \nUnsurprisingly, the X-ray bright systems are concentrated at short periods, as they involve accretion from the companion. All the systems with low-mass donors (red symbols) and 𝑃 orb & 0 GLYPH<149> 5 days have donors that are somewhat evolved and overflow their Roche lobes. The X-ray bright systems with OB-type companions are at slightly longer periods and are fed by wind accretion from donors that nearly fill their Roche lobes. The detached systems with OB primaries, VFTS 243 and HD 130298, are at only slight longer periods than the HMXBs, and thus might evolve to form systems similar to Cyg X-1.  \nThe orbital period and companion star in Gaia BH1 are most similar to in NGC 3201 #12560, the long-period BH candidate discovered in the GC NGC 3201 with MUSE (Giesers et al. 2018). That system has an orbital period of 167 days, eccentricity 0.61, and spectroscopic mass function of 3 GLYPH<149> 25 𝑀 GLYPH<12> , with an GLYPH<25> 0 GLYPH<149> 81 𝑀 GLYPH<12> luminous star. Since that system's inclination is unknown, the companion mass could be  \nFigure 7. Comparison of Gaia BH1 (black) to known BH X-ray binaries. Red and blue symbols correspond to accreting BHs with low- and high-mass companions. Magenta symbols show detached binaries in the globular cluster NGC 3201, and cyan points show detached binaries in which the luminous star is a high-mass ( & 20 𝑀 GLYPH<12> ) star. The properties of Gaia BH1 are most similar to the detached spectroscopic binary NGC 3201 #12560 (Giesers et al. 2018). However, that system is in a globular cluster and thus probably formed through a different channel. The closest system from an evolutionary standpoint is perhaps the LMXB GRS 1915+105, which has a red giant donor and was likely detached when it was a main-sequence star.  \n<!-- image -->  \nany value above 4 GLYPH<149> 5 𝑀 GLYPH<12> . The important difference, however, is that NGC 3201 #12560 is in a GC, where dynamical capture, exchange, and disruption of binaries are all efficient. Two apparent BHs were discovered in NGC 3201 among only 3553 stars with multi-epoch RVs. This implies a BH companion rate of order 1000 times higher than that in the local Galactic field (Section 9). Such a dramatic enhancement in the rate of BH companions in GCs can be most readily understood if the binaries are formed dynamically through capture or exchange (e.g. Fabian et al. 1975; Rasio et al. 2000; Kremer et al. 2018), and in any case, the long-term survival of a primordial binary with a 167 day orbit in such a dense stellar system is unlikely.  \nGiven the orbit of Gaia BH1 in the Galactic disk, formation from an isolated binary seems somewhat more likely than a dynamical formation channel. In this sense, the system is qualitatively different from any other binary shown in Figure 7. The most closely related system may be the BH + giant star 'microquasar' GRS 1915+105 ( 𝑃 orb = 33days; Castro-Tirado et al. 1992; Greiner et al. 2001), which may have formed through a similar channel to Gaia BH1 but with a significantly closer separation, or through a completely different channel.", '8.4 Evolutionary history': "The progenitor of the BH very likely had a mass of at least 20 𝑀 GLYPH<12> . For example, the models of Sukhbold et al. (2016) predict that at solar metallicity, an initial mass of 30 GLYPH<0> 50 𝑀 GLYPH<12> is required to produce a helium core mass of & 9 𝑀 GLYPH<12> (their Figure 19). Such a star would reach a radius of order 10 au as a supergiant if it were allowed to  \nevolve in isolation, which is significantly larger than the present-day separation of the G star and BH. This suggests that the two stars likely interacted prior to the formation of the BH. Given the extreme mass ratio, this interaction is expected to lead to a common envelope (CE) episode, in which the G star was engulfed by the envelope of the massive star, ultimately spiraling in and ejecting some or all of the envelope. Given the short lifetimes of massive stars, the G star would still be contracting toward the main sequence as the BH progenitor completed its evolution.  \nIt is not clear that a . 1 𝑀 GLYPH<12> star can survive a CE episode with a & 20 𝑀 GLYPH<12> companion. Indeed, many calculations in the literature find that ejection of the massive star's envelope by a low-mass star is impossible because there is not enough available orbital energy (e.g. Portegies Zwart et al. 1997; Podsiadlowski et al. 2003; Justham et al. 2006). The same challenge applies to modeling the formation of BH LMXBs, whose evolutionary channels are still poorly understood. We first review the 'standard' formation channel through a common envelope event, and then discuss possible alternatives.", '8.4.1 Common envelope channel': "In the CE scenario, the initial separation of the BH progenitor and the G star would have been in the range of 5-15 au (Portegies Zwart et al. 1997). The binary would have emerged from the CE episode as a GLYPH<20> 1 𝑀 GLYPH<12> star in a close orbit with a & 10 𝑀 GLYPH<12> helium core, which might also have retained some of its envelope. The standard 𝛼𝜆 CE prescription predicts that the ratio of the separations before and after  \nthe common envelope episode is (e.g. Webbink 1984)  \n𝑎 𝑓 𝑎 𝑖 = 𝑀 c 𝑀 1 GLYPH<18> 1 , 2 𝑀 e 𝛼𝜆𝑟 𝐿 𝑀 2 GLYPH<19> GLYPH<0> 1 GLYPH<150> (8)  \nwhere 𝑀 c and 𝑀 e are the mass of the BH progenitor's core and envelope, 𝑀 1 = 𝑀 c , 𝑀 e is the total mass of the progenitor, 𝑟 𝐿 is its Roche lobe radius in units of the semimajor axis, and 𝑀 2 is the mass of the G star. 𝛼 and 𝜆 are dimensionless parameters, respectively describing the fraction of the G star's orbital energy that goes into ejecting the BH progenitor's envelope, and the binding energy of the BH progenitor's envelope. If we take 𝑀 c = 10 𝑀 GLYPH<12> , 𝑀 e = 20 𝑀 GLYPH<12> , 𝑟 𝐿 = 0 GLYPH<149> 65, 𝑀 2 = 1 𝑀 GLYPH<12> , 𝛼 = 1, and 𝜆 = 0 GLYPH<149> 5 (a rather optimistic choice of 𝜆 , e.g. Klencki et al. 2021), we obtain 𝑎 𝑓 GLYPH<157> 𝑎 𝑖 = 0 GLYPH<149> 0025, such that for an initial separation of 𝑎 𝑖 = 2000 𝑅 GLYPH<12> , the final separation is 𝑎 𝑓 = 5 𝑅 GLYPH<12> . This is just wide enough that the final configuration might be stable and avoid a merger, but only for a narrow range of 𝑎 𝑖 . For lower values of 𝛼𝜆 , there are no initial separations that both lead to Roche lobe overflow of the BH progenitor and avoid a merger (e.g. Portegies Zwart et al. 1997). In any case, the predicted final separation is much closer than the observed 𝑎 GLYPH<25> 300 𝑅 GLYPH<12> .  \nStripped of its envelope, the helium core is expected to drive a prodigious Wolf-Rayet wind (e.g. Woosley et al. 1995), which will both diminish the final mass of the compact remnant and widen the orbit. Exactly how much mass is lost in this stage is uncertain (e.g. Woosley et al. 1995; Fryer & Kalogera 2001), and attempts to measure the relevant wind mass loss rates observationally are complicated by factors such as clumping (e.g. Smith 2014; Shenar 2022). In Gaia BH1, a strong wind cannot have persisted for very long, or the helium core mass would have fallen below the required 10 𝑀 GLYPH<12> . Eventually, the helium core will collapse to a BH, perhaps accompanied by some additional mass loss and/or an asymmetric natal kick.  \nThe wide orbit of Gaia BH1 and the unremarkable surface abundances of the G star strongly suggest that there has not been significant mass transfer from the G star to the BH. We can thus - in contrast with the situation for typical BH LMXBs (e.g. Podsiadlowski et al. 2002; Pfahl et al. 2003; Justham et al. 2006; Fragos & McClintock 2015) - rule out a scenario in which the luminous star was initially significantly more massive than today and has been stripped down to sub-solar mass by Roche lobe overflow. We can similarly rule out a scenario in which the initial mass of the BH was lower than today, and it grew significantly by accretion.  \nThe orbital eccentricity, 𝑒 GLYPH<25> 0 GLYPH<149> 45, places a constraint on mass loss and natal kicks during the BH's formation. If the natal kick was purely due to mass loss (and not an asymmetric stellar death) and the orbit was circular before the BH formed (e.g. Blaauw 1961), then the predicted final eccentricity would be  \n𝑒 = Δ 𝑀 BH 𝑀 BH , 𝑀 ★ GLYPH<150> (9)  \nsuch that mass loss of Δ 𝑀 BH GLYPH<25> 4 GLYPH<149> 8 𝑀 GLYPH<12> would be required to fully explain the observed eccentricity. Here Δ 𝑀 BH represents any mass lost during the BH progenitor's death that does not fall into the BH. The net velocity imparted to the center-of-mass in this case is 𝑣 sys = 𝑒 GLYPH<2> 𝑣 BHprog , where 𝑣 BHprog is the orbital velocity of the BH progenitor at the time of its death. Irrespective of the orbital period at the time, 𝑣 BHprog was likely rather low given the extreme mass ratio, and so (symmetric) mass loss alone is unlikely to have resulted in a system velocity larger than GLYPH<24> 10 km s GLYPH<0> 1 . Loss of & 4 𝑀 GLYPH<12> from the system at the time of the BH's formation seems somewhat unlikely in the CE channel because the core should collapse rapidly and most  \nof the envelope should have been removed by the CE episode, but the details of how the CE episode would proceed are uncertain.  \nOn the other hand, if a kick occurred mainly due to asymmetries during the BH progenitor's death (e.g., asymmetric ejecta or neutrino emission), the final eccentricity depends on the kick velocity and direction, and on the period at the end of the CE episode (e.g. Brandt & Podsiadlowski 1995; Hurley et al. 2002). In the limit of no mass loss, the systemic velocity induced by a kick is simply 𝑣 kick GLYPH<2> 𝑞 GLYPH<157> ' 1 , 𝑞 ' , which in the limit of 𝑞 GLYPH<29> 1 is close to 𝑣 kick . Here 𝑞 = 𝑀 2 GLYPH<157> 𝑀 ★ .  \nForming Gaia BH1 via a CE is quite challenging for two reasons: (1) the G star may not have had enough orbital energy to eject the envelope of its much more massive companion (this is generically a problem for BH LMXB formation models, but is particularly problematic here because the G star cannot have been born with a higher mass than observed today), and (2) if CE ejection was successful, the post-CE separation is expected to be significantly closer ( GLYPH<24> 5 𝑅 GLYPH<12> ) than the wide orbit observed today ( GLYPH<24> 300 𝑅 GLYPH<12> ) under standard assumptions for CE ejection efficiencies. A fine-tuned natal kick to the BH could in principle widen the orbit, but this would also make the final orbit highly eccentric and impart a large systemic velocity to the binary's center of mass, both of which are not observed. For example, assuming an optimistic post-CE separation of 20 𝑅 GLYPH<12> and random kick orientations, we find (e.g. Brandt & Podsiadlowski 1995) that kicks of order 200 km s GLYPH<0> 1 are required to widen the orbit to 𝑎 & 200 𝑅 GLYPH<12> . Such kicks result in a minimum eccentricity of 0.90 and a median of 0.95, both much larger than observed.  \nGiven the seemingly impossible nature of forming Gaia BH1 through the CE channel under standard assumptions, we performed a study of more extreme assumptions for CE and BH natal kicks to determine the required combination of assumptions to produce a binary with the observed properties of Gaia BH1. We follow the methods introduced in Wong et al. (2022) which uses COSMIC, a rapid binary population synthesis code, to explore all possible combinations of Zero Age Main Sequence (ZAMS) properties, CE assumptions, and natal kick assumptions; for an in-depth discussion of how COSMIC evolves binary-star populations, see Breivik et al. (2020).  \nInstead of selecting a single set of assumptions for CE and natal kicks, we use emcee to sample the posterior of the masses, orbital period, and eccentricity of the binary hosting Gaia BH1 with combinations of the ZAMS mass, orbital period, and eccentricity, as well as the CE ejection efficiency ( 𝛼 ) and the BH natal kick strength ( 𝑣 kick ) and direction ( 𝜃 , 𝜙 ) and the mean anomaly of the orbit at the time of core collapse ( M ) as model parameters. We place uniform priors on all ZAMS binary parameters as well as the CE and natal kick assumptions. We assign wide limits for the ZAMS primary mass prior between 1 𝑀 GLYPH<12> and 150 𝑀 GLYPH<12> and narrow limits for the ZAMS secondary mass prior between 0 GLYPH<149> 5 𝑀 GLYPH<12> and 2 𝑀 GLYPH<12> since the G star is not expected to have significantly gained or lost mass during the binary evolution. The orbital period and eccentricity priors are limited to 500 days and 6000 days and 0 and 1 respectively. We limit our prior on 𝛼 to be between 0 GLYPH<149> 1 and 20 and our prior on 𝑣 kick to be between 0 km s GLYPH<0> 1 and 300 km s GLYPH<0> 1 based on the relatively wide and moderately eccentric orbital configuration of Gaia BH1's binary host. We place uniform priors on the unit sphere for the natal kick angles, 𝜃 , 𝜙 , and all allowed values of the mean anomaly, M , at core collapse (0-360 deg ). Finally, we fix the metallicity of all simulated binaries to 𝑍 = 0 GLYPH<149> 63 𝑍 GLYPH<12> , where 𝑍 GLYPH<12> = 0 GLYPH<149> 02, consistent with » Fe GLYPH<157> H … = GLYPH<0> 0 GLYPH<149> 2 under the assumption that the G star follows the solar abundance pattern.  \nSince the age of the observed G star is not known to great precision, we require the properties of our simulated BH binaries to match the observed properties of the binary hosting Gaia BH1 just after the  \nformation of the BH. This choice does not affect our results because the G star's presence on the main sequence implies that tides are not expected to alter the binary's orbit between the formation of the BH and the present day.  \nWeinitialize 1024 walkers uniformly over the prior space in ZAMS mass, orbital period, eccentricity, natal kick strength and direction, and CE ejection efficiency. We evolve each walker for 100000 steps, thin the chains by a factor of 10, and retain the final 2000 steps of each chain to ensure a proper burn-in of our sampling.  \nWe find that the preferred evolutionary pathway is for binary to begin with an orbital period between 1600-4800 days depending on the initial eccentricity, which ranges between 0-0 GLYPH<149> 5. The BH progenitor's ZAMS mass is preferred to be between 65-75 𝑀 GLYPH<12> , while the G star's ZAMS mass is preferred to be between 0 GLYPH<149> 8-1 GLYPH<149> 1 𝑀 GLYPH<12> . From this initial configuration, the BH progenitor loses roughly 20 𝑀 GLYPH<12> due to strong stellar winds which widens the binary out to periods of GLYPH<24> 7000 days. As the BH progenitor evolves off the main sequence and begins core helium burning, its envelope expands and tides work to circularize and shrink the binary back to orbital periods between near GLYPH<24> 3000 days. Near the end of core helium burning, the BH progenitor fills it's Roche lobe and the mass transfer becomes dynamically unstable and forms a CE due to both the convective envelope of the BH progenitor and the orbital evolution which drives the binary components close together due to their highly asymmetric mass ratios. The BH progenitor's envelope is fully stripped in the CE, leaving behind a GLYPH<24> 20 𝑀 GLYPH<12> Wolf-Rayet star in a GLYPH<24> 20 day orbit with the G star. The Wolf-Rayet star's strong wind mass loss widens the orbit to roughly 80 days. Finally the BH progenitor reaches core collapse and explodes which both reduces the BH progenitor mass from GLYPH<24> 12 𝑀 GLYPH<12> to the mass of Gaia BH1 and imparts a natal kick that increases the eccentricity and orbital period to reach the present day observed values.  \nFigure 8 shows the combination of ZAMS parameters and binary evolution assumptions which produce Gaia BH1-like binaries whose masses, orbital period, and eccentricity match Gaia BH1's present day properties for the 50 th and 90 th percentiles of each distribution. Due to the moderately wide orbit of Gaia BH1, our models prefer low natal kick strengths with 𝑣 kick GLYPH<159> 45 km s GLYPH<0> 1 , which is consistent with Gaia BH1's orbit lying near the plane of the Galactic disk. There is a slight preference for kicks centered in the plane of the orbit ( 𝜃 GLYPH<24> 90 GLYPH<14> ), though there is support over nearly the entire allowed range. A correlation between 𝜙 and 𝑣 kick exists such that larger kicks prefer azimuthal angles near 90 GLYPH<14> and 270 GLYPH<14> . This is due to a higher likelihood of the binary to remain bound for larger natal kicks when the orbital and kick velocities are pointed in opposite directions. We do not show samples for the mean anomaly since there is a near equal preference for all values of M . This is not unexpected since the binary's orbit is circularized during the CE, thus fixing a constant orbital speed.  \nWhile the ZAMS masses are not correlated with the ZAMS orbital period, ZAMS eccentricity, or natal kick parameters, there is a correlation between the ZAMS primary mass and the CE ejection efficiency. This is expected since an increased primary mass will have a larger envelope which requires a larger, more efficient, 𝛼 to produce a successful CE ejection. The ZAMS orbital period is also correlated with 𝛼 ; longer initial periods require lower, less efficient 𝛼 's to shrink the orbit to the proper separation after the CE.  \nOur models prefer extremely efficient CE ejections with 𝛼 constrained to always be larger than 5 with a preference for 𝛼 = 14. We remind the reader that 𝛼 represents the fraction of the G star's orbital energy that goes into unbinding the BH progenitor's envelope. Values of 𝛼 significantly larger than 1 signify that the liberated orbital energy alone is insufficient to eject the envelope, and that some other  \nsource of energy is required. While there are large uncertainties in the binding energy of the BH progenitors in our models for which 𝜆 GLYPH<24> 0 GLYPH<149> 5, it is highly unlikely that such uncertainties will change 𝜆 by more than a factor of 10. Furthermore, while it is possible for 𝛼 values to exceed unity due to additional energy sources like recombination, it is unlikely that such sources will dominate the energy budget (Ivanova et al. 2015; Ivanova 2018). It is also possible that CE proceeds differently than the standard 𝛼𝜆 prescription (e.g. Hirai & Mandel 2022), though we leave this for a future study. We thus conclude that Gaia BH1 can be formed through standard binary evolution with a common envelope only under extreme (and likely unphysical) assumptions about how the common envelope evolution proceeds.  \nIn this sense, Gaia BH1 is reminiscent of the several known pulsars in wide, moderately eccentric binaries ( 𝑃 orb = 80 GLYPH<0> 400 days) with companions that appear to be low-mass main-sequence stars (e.g. Phinney & Verbunt 1991; Champion et al. 2008; Parent et al. 2022). These systems are also too wide to form via a common envelope, but too close to have avoided one. It is possible that the companions are white dwarfs, but their nonzero eccentricities are difficult to explain in this scenario. Perhaps they formed through a similar channel to Gaia BH1. We consider possible alternative formation scenarios below.", '8.4.2 Formation from a progenitor that never became a giant': "A CE event can be avoided if the BH progenitor never became very large. Models predict that at sufficiently high masses ( 𝑀 & 50 𝑀 GLYPH<12> , depending on wind prescriptions and rotation rates), stars do not expand to become red supergiants, but instead lose their hydrogen envelopes to winds and outbursts during and shortly after their mainsequence evolution, reaching maximum radii of order 50 𝑅 GLYPH<12> (e.g. Humphreys & Davidson 1994; Higgins & Vink 2019). A scenario in which the BH progenitor in Gaia BH1 never became a red supergiant has the attractive feature that it could seem to avoid a CE episode. It requires, however, that the G star formed uncomfortably close to the BH progenitor. The orbit would have expanded by a factor of GLYPH<24> 5 due to the BH progenitor's mass loss, as the quantity 𝑎 GLYPH<2> ' 𝑀 1 , 𝑀 2 ' is conserved under adiabatic orbit evolution (Jeans 1924). This in turn requires an initial separation of . 60 𝑅 GLYPH<12> and an initial period of less than 10 days. It seems somewhat unlikely that a premain sequence solar-mass star could form and survive so close to a massive star, but observational constraints on the existence of such companions to O stars - which would have luminosity ratios & 10 5 and angular separations of order 0.1 mas at typical Galactic distances - are extraordinarily difficult to obtain (e.g. Rizzuto et al. 2013; Sana et al. 2014).  \nA related set of models allow rapidly-rotating massive stars to avoid becoming red supergiants due to mixing, such that a chemical gradient never forms and a majority of the star is converted to helium ('quasi-homogeneous evolution'; e.g. Maeder 1987; Woosley & Heger 2006). Such a scenario could also avoid the problems with the CE channel, because the G star could have been born in an orbit close to the current separation, never interacting with the BH companion. The BH progenitor could have also had a somewhat lower initial mass, perhaps down to 15 𝑀 GLYPH<12> . However, most models predict that such evolution only occurs at low metallicity, because winds remove too much angular momentum at high metallicity (e.g. Yoon et al. 2006; Brott et al. 2011). The » Fe GLYPH<157> H … = GLYPH<0> 0 GLYPH<149> 2 we measure for the G star thus presents a challenge for this scenario.  \nFigure 8. The combination of sampled ZAMS component masses, initial orbital period, and initial eccentricity with uncertain binary evolution model parameters for CE ( 𝛼 ) and BH natal kicks ( 𝑣 kick , 𝜙 , 𝜃 ) which evolve to produce BHs in binaries with properties similar to Gaia BH1 and its binary host. We do not show distributions for the mean anomaly ( M ) because there is not a significant preferred angle nor any appreciable correlations with any other parameters. The contours show the 50 th and 90 th percentiles of each distribution. The diagram in the upper right corner illustrates the orbital configuration and natal kick direction for a Wolf-Rayet BH progenitor (blue) at core collapse in a circular orbit with a G star (green). Our models strongly prefer small natal kicks at BH formation and extremely efficient CE ejection efficiency.  \n<!-- image -->", '8.4.3 Formation through three-body dynamics': 'Mechanisms for forming BH LMXBs have been proposed involving tidal capture of inner binaries in hierarchical triples (e.g. Naoz et al. 2016)orverywideBH+normalstarbinariesperturbedbyencounters with field stars (e.g. Michaely & Perets 2016). In both cases, tidal effects can bring the binary into a close orbit only if a periastron  \npassage occurs within which the star comes within a few stellar radii of the BH. If this occurs, tides can efficiently dissipate orbital energy, locking the star into a close orbit and forming a LMXB. Such a scenario seems unlikely to have operated in Gaia BH1, because the observed periastron distance is too wide for tides to be important.', '8.4.4 Formation in a dense cluster': 'In dense environments like the cores of GCs, binaries with separations and eccentricities similar to Gaia BH1 can be formed via exchange encounters (e.g. Kremer et al. 2018). Formation in a GC is disfavored for Gaia BH1 due to its thin-disk like orbit and high metallicity.  \nAnother possibility is that the Gaia BH1 system was assembled dynamically in an open cluster. Recent work (Fujii & Portegies Zwart 2011; Sana et al. 2017; Ramírez-Tannus et al. 2021) suggests that many massive star binaries are initially formed in long-period orbits, but are subsequently hardened due to dynamical encounters over millions of years. In the case of Gaia BH1, the binary may have formed at longer periods, allowing the BH progenitor to evolve into a red supergiant and then collapse into a BH within GLYPH<24> 3 Myr. Dynamical interactions after the BH formation could then shrink the orbital separation to the value observed today. Alternatively, the BH progenitor may have evolved separately from the G star in Gaia BH1. After BH formation, a dynamical exchange within the birth cluster could have allowed it to capture the G star into the current orbit (Banerjee 2018a,b). The binary system could have simultaneously or subsequently been ejected from the cluster at low velocity, or simply remained bound after the cluster dissolved (Schoettler et al. 2019; Dinnbier et al. 2022).  \nShikauchi et al. (2020) recently investigated the dynamical formation of BH + main sequence binaries in open clusters and predicted that only a small fraction of such binaries observable by Gaia would have formed dynamically. However, their simulations did not include any primordial binaries, and thus did not account for binaries formed through exchange interactions. Further work investigating this formation channel is warranted.', '8.4.5 Formation in a hierarchical triple': "Given the problems with binary formation models noted above, triple or multi-body formation channels may have operated in Gaia BH1. One possibility is that the BH progenitor was originally in a compact binary ( 𝑃 orb . 20 d) composed of two massive stars. As the stars evolve, they transfer mass to each other, preventing either star from evolving into a red supergiant. A simple possibility is that the inner binary merged during unstable Case B mass transfer, and the merger product never swelled into a red supergiant (Justham et al. 2014), allowing the G star (originally an outer tertiary) to survive.  \nAlternatively, the inner binary may have undergone stable mass transfer. Consider a binary with 𝑀 1 GLYPH<24> 35 𝑀 GLYPH<12> and 𝑀 2 GLYPH<24> 20 𝑀 GLYPH<12> in an orbit with 𝑃 orb GLYPH<24> 4 days. As the primary evolves, its envelope is stripped during stable Case A mass transfer, producing a helium star which then collapses into a BH with 𝑀 BH GLYPH<24> 10 𝑀 GLYPH<12> . During the mass transfer, the orbital period and secondary mass only change slightly. A very similar evolution is thought to produce observed HMXBs such as Cyg X-1 (Qin et al. 2019). The secondary then evolves and has its envelope stripped during a second phase of stable Case A mass transfer, forming a helium star. Depending on the strength of its stellar wind and the core-collapse process, this helium star could collapse into a second BH, or explode and produce a neutron star, which may or may not remain bound. A similar model was proposed by Shenar et al. (2022a) for the binary VFTS 243. At the end of the process we are left with a BH and possibly a second compact object in a short-period orbit, with a still-bound outer tertiary.  \nTotest this idea, we constructed binary models with MESA (Paxton et al. 2011) using the problem setup described in Fuller & Ma (2019). Assuming non-conservative mass transfer (as expected from the high  \nmass transfer rates achieved during Case A mass transfer), we find very similar binary evolution histories to those discussed above. For example, if the second phase of mass transfer begins with a 20 𝑀 GLYPH<12> donor and a 10 𝑀 GLYPH<12> BH in an orbital period of 4 days, it ends with a 5 𝑀 GLYPH<12> helium star in a 9 GLYPH<149> 6 day orbit such that the orbit of the G star in Gaia BH1 could plausibly remain stable. As in the channels with single BH progenitors that never became giants, this scenario requires the G star to have formed quite close ( . 100 𝑅 GLYPH<12> ) from the inner binary. Given this and the expected orbital evolution of the inner binary due to mass transfer and winds, the parameter space of triples that can form a system like Gaia BH1, while not empty, is likely rather small. In the triple scenario it would also be somewhat surprising that the only good candidate we identified from Gaia DR3 has a (relatively) short period, because the effective search volume is larger for long periods (up to 1000 days; Section 9), and hierarchical triples with long outer periods would also be stable for a wider range of inner binary parameters.  \nWhile the triple model may seem unlikely given the rarity of BH HMXBs, it is important to note that the HMXB phase lasts at most . 10 6 yr, while the subsequent quiescent phase can easily last more than 10 10 yr. If a significant fraction of HMXBs host outer low-mass tertiary stars, this could account for the high space-density of objects like Gaia BH1 relative to HMXBs.  \nA prediction of this channel is that Gaia BH1 might still harbor a binary compact object, whose orbital motion would induce potentially observable perturbations to the RVs of the outer G star on half the orbital period of the inner binary. Hayashi & Suto (2020) show that the amplitude of these perturbations would scale as GLYPH<0> 𝑃 orb GLYPH<150> inner GLYPH<157> 𝑃 orb GLYPH<150> outer GLYPH<1> 7 GLYPH<157> 3 ; in Gaia BH1, the expected perturbation semi-amplitude ranges from GLYPH<24> 2 m s GLYPH<0> 1 for a 4 day inner period, to GLYPH<24> 20 m s GLYPH<0> 1 for a 10 day inner period, to GLYPH<24> 250 m s GLYPH<0> 1 for a 30-day inner period. Since the G star is bright and slowly-rotating, these amplitudes are within the realm of detectability. In addition, an inner binary would cause precession of the G star's orbit, which could be measurable via time-evolution of the orbit's orientation and RVs. The predicted minimum precession period ranges from about 35 years for a 30-day inner period, to 350 years for a 4-day inner period (and still longer for shorter inner periods). The RV semi-amplitude of this precession would be large (of order the observed semi-amplitude; 70 km s GLYPH<0> 1 ), but a long observing time baseline would be required to detect it.", '9 HOW MANY SIMILAR OBJECTS EXIST?': "Here we consider what the discovery of one dormant BH in a binary implies about the broader population. The precision of the conclusions we can draw is of course limited by the Poisson statistics of 𝑁 = 1. We will show, however, that Gaia has observed tens of millions of sources around which a BH companion could have been detected, if one existed. This makes possible robust inference about the occurrence rate of BH companions, as long as the selection function of the astrometric catalog can be understood.  \nAs described by Halbwachs et al. (2022), the cuts placed on astrometric solutions published in Gaia DR3 were quite conservative, designed to reduce false positives at the expense of completeness. Attempts to fit astrometric orbits were made only for sources satisfying ruwe > 1.4 (indicating a poor fit with a single-star solution) that were observed in at least 12 visibility periods. Sources with evidence of marginally resolved companions or contaminated photometry were likewise discarded. Sources for which a satisfactory fit could be achieved with an accelerating solution were not fit with an  \norbital solution. After orbital solutions were fit, those which did not satisfy all of the following three conditions were discarded:  \n𝜛 𝜎 𝜛 ¡ 20000 days 𝑃 orb (10)  \n𝑎 0 𝜎 𝑎 0 ¡ 158 √︁ 𝑃 orb GLYPH<157> day (11)  \n𝜎 𝑒 GLYPH<159> 0 GLYPH<149> 079 ln ' 𝑃 orb GLYPH<157> day ' GLYPH<0> 0 GLYPH<149> 244 (12)  \nHere 𝜎 𝑒 is the uncertainty on the orbital eccentricity. As shown by Halbwachs et al. (2022), Equation 10 removes a large fraction of otherwise acceptable solutions. After Equation 10 is applied, Equations 11 and 12 remove a relatively small number of additional solutions.  \nAll these cuts are more aggressive for binaries with shorter periods than those with longer periods. The astrometric solution for Gaia BH1 has 𝜛 GLYPH<157> 𝜎 𝜛 = 120 and 20000 GLYPH<157> 𝑃 orb = 107. That is, it just narrowly passes the cut in Equation 10, and would have been excluded if it were 10% more distant. Similarly, it has 𝑎 0 GLYPH<157> 𝜎 𝑎 0 = 13 GLYPH<149> 6; if this quantity were smaller than 158 GLYPH<157> p 𝑃 orb = 11 GLYPH<149> 6, it would have been excluded by Equation 11. The fact that Gaia BH1 narrowly escaped the quality cuts applied to the whole catalog has two implications.  \nFirst, there are likely other similar systems that could be discovered with only slightly looser cuts. We estimate the number of detectable systems explicitly in Section 9.1. Second, BH companions may be more common at shorter periods than at longer periods. Gaia DR3 was much more sensitive to binaries with long orbital periods (up to the observing baseline, GLYPH<25> 1000 days) than to those with short periods: Equation 10 will exclude systems with 𝜛 GLYPH<157> 𝜎 𝜛 GLYPH<159> 100 at 𝑃 orb = 200 days, but only those with 𝜛 GLYPH<157> 𝜎 𝜛 GLYPH<159> 20 at 𝑃 orb = 1000 days. At fixed absolute magnitude, this translates to a factor of GLYPH<25> 3 larger distance limit, and a factor of GLYPH<25> 10 larger survey volume (Section 9.1).  \nGaia BH1 has a shorter orbital period than 95% of all objects in DR3 with astrometric solutions. The fact that it appears to be the only credible BH with 𝑃 orb . 1000 days, while at longer-periods ( 𝑃 orb = 400 GLYPH<0> 1000 days), BH companions could have been detected around a much larger number of sources, thus points to a relative paucity of longer-period binaries consisting of a BH and a normal (low-mass) star.", '9.1 Effective survey volume of the astrometric binary sample': "Wenowexplore around how many sources Gaia could have detected a BH companion if one existed. Since 'detectability' depends on the astrometric uncertainties, this requires a noise model. We describe such a model in Appendix G.  \nThe actual uncertainties for a given source will depend on factors such as the orbital eccentricity, orientation, and Gaia scanning pattern, but we neglect these details. We assume that at 𝑃 orb GLYPH<159> 1000 days, the astrometric uncertainties 𝜎 𝜛 and 𝜎 𝑎 0 depend only on apparent magnitude and 𝑎 0 . We constrain the dependence on both parameters empirically using the 170,000 astrometric binary solutions published in DR3. The median astrometric uncertainty for sources with orbital solutions in DR3 ranges from 𝜎 𝜛 = 0 GLYPH<149> 020 mas at 𝐺 GLYPH<159> 14, to 0.045 mas at 𝐺 = 16, to 0.08 mas at 𝐺 = 17. The median uncertainty in 𝑎 0 is larger by a factor of GLYPH<25> 1 GLYPH<149> 5 on average. The astrometric uncertainties also depend on 𝑎 0 , increasing by about a factor of 2 from 𝑎 0 = 0 GLYPH<149> 3 mas to 𝑎 0 = 3 mas. This reduces the number of BHcompanions that could be detected, since BHs produce larger 𝑎 0 at fixed 𝑃 orb than luminous companions.  \nWe now assign a hypothetical 10 𝑀 GLYPH<12> dark companion to each  \nsource in the Gaia catalog and ask whether it could have been detected at a given period. When calculating the expected 𝑎 0 , we assign all sources on the main sequence a mass based on their absolute magnitude, and all giants a mass of 1 𝑀 GLYPH<12> . We apply the same cuts on apparent magnitude, photometric contamination, and image parameter metrics that were applied to sources in the actual Gaia binary catalog (Halbwachs et al. 2022). We calculated the expected uncertainties, 𝜎 𝜛 and 𝜎 𝑎 0 , based on the mean and observed scatter of these quantities among all astrometric solutions with the same 𝐺 and 𝑎 0 , as described in Appendix G. We assume that no binaries with periods within GLYPH<6> 15% of one year can be detected.  \nFigure 9 shows the number of sources in the gaia\\_source catalog for which a 10 𝑀 GLYPH<12> dark companion would likely pass different quality cuts. The black line corresponds to the detection thresholds actually employed in DR3. As discussed above, there is a strong selection bias against shorter orbital periods, which were only accepted at extremely high 𝜛 GLYPH<157> 𝜎 𝜛 .  \nIn DR4, Gaia is planning to publish epoch-level astrometry for all sources, irrespective of whether they satisfy cuts like Equation 10. This will make it possible to identify sources similar to Gaia BH1 at larger distances. The dotted red line in Figure 9 shows the effective search volume for the looser - but still relatively conservative threshold of a minimum astrometric SNR of 10. Using this threshold instead of Equations 10 and 11 would increase the search volume for 10 𝑀 GLYPH<12> companions at a period of 186 days by a factor of 20. Gains at shorter periods would be even more dramatic. Of course, removing these cuts is likely to result in a large number of spurious solutions - this is why they were applied in the first place. Inspection of the epoch-level astrometric measurements, coupled with spectroscopic follow-up of the most promising candidates, will be necessary and sufficient to eliminate spurious astrometric BH candidates.  \nIn addition, the factor of two improvement in observing baseline between DR3 and DR4 will make DR4 have more precise astrometry and be sensitive to orbits with periods as large as 2000 days, for which the search volume is also larger. While extrapolating future detection rates from the current search criteria and one detection is risky, we conclude that there is good reason to expect dozens more astrometric BHs from Gaia DR4. The number could also be higher if the occurrence rate of (wide) BH companions is significantly higher among massive stars than among solar-type stars, as is expected theoretically (e.g. Langer et al. 2020; Janssens et al. 2022). Few massive stars received orbital solutions in DR3, because most are too distant to pass the stringent parallax uncertainty cut of Equation 10.  \nOne can infer a rough lower limit on the number of objects similar to Gaia BH1 that exist in the Milky Way. Figure 9 shows that the effective search volume at 186 days is about 2 GLYPH<149> 5 GLYPH<2> 10 6 sources for which a 10 𝑀 GLYPH<12> companion would be detectable with DR3 cuts. Translating this into a constraint on the occurrence rate of BH companions in a given period range requires an assumed BH companion period distribution, which is obviously very uncertain. If we assume that a fraction 𝑓 of all low-mass stars have a 10 𝑀 GLYPH<12> BH companion with 100 GLYPH<159> 𝑃 orb GLYPH<157> days GLYPH<159> 300 and the period distribution is log-uniform in this range (d 𝑁 GLYPH<157> d log 𝑃 orb = const), then we find that 1 detection implies 𝑓 GLYPH<25> 4 GLYPH<2> 10 GLYPH<0> 7 ; i.e., an occurrence rate of order 1 in 2.5 million for BH companions to low-mass stars with 100 GLYPH<159> 𝑃 orb GLYPH<157> days GLYPH<159> 300. This is a constraint on the average incidence across a range of stellar masses, sampled in a way that is different from the mass function. A large majority of the sources in this search volume are mainsequence stars with masses between 0.5 and 2 𝑀 GLYPH<12> . Multiplying by GLYPH<24> 10 11 low-mass stars, this would translate to GLYPH<24> 40000 similar systems in the Milky Way - a factor of GLYPH<24> 40 larger than the expected  \nFigure 9. Number of sources in Gaia DR3in which a 10 𝑀 GLYPH<12> dark companion could be detected. Black line shows the cuts employed in DR3 (Equations 10 and 11), which remove most short-period binaries. Dashed red line corresponds to simply requiring 𝑎 0 GLYPH<157> 𝜎 𝑎 0 ¡ 10 and 𝜛 GLYPH<157> 𝜎 𝜛 ¡ 10; this results in a much larger searchable volume at short periods. Binaries with periods near 1 year are excluded, being difficult to detect with astrometry.  \n<!-- image -->  \nnumber of BH X-ray binaries (Corral-Santana et al. 2016), but still less than 0.1% of the expected total number of BHs.  \nWe can compare the occurrence rate of Gaia BH1-like systems to the number of stars formed that are massive enough to become BHs. For a Kroupa (2001) IMF, roughly 2 out of every 1000 stars formed has 𝑀 ¡ 20 𝑀 GLYPH<12> . Thus, if 1 out of every 5000 stars with 𝑀 ¡ 20 𝑀 GLYPH<12> has a solar-type companion in a Gaia BH1-like orbit when it becomes a BH, this could produce the population implied by Gaia BH1, with 1 in 2.5 million low-mass stars having a BH companion with 100 GLYPH<159> 𝑃 orb GLYPH<157> days GLYPH<159> 300. The '1 in 5000' factor could reflect a dearth of extreme mass ratio companions to massive stars, and/or the fact that a majority of such companions do not evolve to become intact BH + normal star binaries at these separations.  \nThe lack of a solid detection at 𝑃 orb = 400 GLYPH<0> 1000 days implies a lower BH companion incidence rate there. Assuming a log-uniform period distribution as above, the probability of detecting at least one BHcompanioninthisperiod range would be & 50%for 𝑓 ¡ 4 GLYPH<2> 10 GLYPH<0> 8 ; i.e., if the BH companion occurrence rate were more than 1 in 25 million in this period range. As above, this limit depends sensitively on the assumed period distribution.  \nLimits at shorter periods are weaker, given the implicit period bias of the selection cuts. For example, the number of sources for which a 10 𝑀 GLYPH<12> companionina60-dayorbitwouldbedetectedastrometrically in DR3 is only 0.3 million. At such short periods, the (lack of) detections in Gaia SB1(spectroscopic) solutions becomes a more stringent limit. DR3 contains 5.9 million sources satisfying rv\\_method\\_used = 1 (meaning that RVs were measured from individual epochs) and rv\\_nb\\_transits >10(enough epochs that an orbital solution could plausibly be inferred). About 2.5 million of these are near the main sequence. Orbital solutions with 𝐾 ★ ¡ 250 km s GLYPH<0> 1 were removed from DR3 (Gaia Collaboration et al. 2022a), but this will only remove BH companions to solar-type stars at 𝑃 orb . 3 days. This suggests that BHs could have been detected as companions to about 2.5 million sources at periods ranging from 3 to 1000 days, greatly exceeding the astrometric search volume at short periods. A non-detection of BH companions in the DR3 SB1 sample would seem to rule out a large increase in the BH companion incidence rate toward short peri-  \n, it is still unclear from the published information how SB1 solutions were processed and what cuts were ultimately applied to them. It is possible, for example, that close, tidally synchronized companions to BHs would have had too large 𝑣 sin 𝑖 for RVs to be reliably measured by the Gaia RVS pipeline.", '9.2 Prospects for detecting similar systems with spectroscopy and photometry': "Time-domain spectroscopic surveys provide another promising avenue for detecting more BH + normal star binaries like Gaia BH1. Indeed, Gaia BH1 was observed twice by the LAMOST survey; the source would have been easily recognizable as a promising BH candidate if it had been observed GLYPH<24> 10 times with suitable cadence. Gaia will eventually provide ¡ 20 epochs of spectroscopy for most bright, 'normal' stars with 𝐺 . 12 GLYPH<149> 5. The calculations above suggest that at most a few BHs similar to Gaia BH1 will satisfy 𝐺 . 12 GLYPH<149> 5. Manyepoch surveys from the ground are needed to reach targets that are fainter or harder to measure RVs for (e.g. hot and rapidly-rotating stars). Most current wide-field surveys do not obtain enough epochs per target to constrain the orbital period, but more targeted surveys focused on binarity do (e.g. Giesers et al. 2019; Wang et al. 2021; Shenar et al. 2022b). We note that, except in cases where orbital periods can be constrained with photometry or astrometry, spectroscopic surveys that obtain & 10 epochs for a modest number of targets are more useful for characterizing binaries than those obtaining 1-3 spectra for a larger number of targets.  \nBesides astrometry and spectroscopy, another avenue for detecting Gaia BH1-like binaries is through photometry. Ellipsoidal variation is not currently detectable (the expected variability amplitude is GLYPH<25> 5 GLYPH<2> 10 GLYPH<0> 7 ; Morris & Naftilan 1993) but will become detectable when the G star becomes a red giant. In the meantime, Gaia BH1like systems could be detected through self-lensing (e.g. Masuda & Hotokezaka 2019). If a binary like Gaia BH1 were viewed edge-on, the G star would brighten by GLYPH<25> 5% each time it moved behind the BH. The typical duration of this brightening would be 3-4 hours; i.e., a duty cycle of about 10 GLYPH<0> 3 . Self-lensing would only be detectable for the GLYPH<25> 0 GLYPH<149> 3% of binaries with sufficiently edge-on inclinations to be eclipsing, but given the large expected magnification amplitude, a self-lensing search could extend to significantly fainter targets than are accessible with spectroscopic surveys. Given that several selflensing binaries containing white dwarfs were detected by Kepler with longer periods and much smaller amplitudes than expected for a system like Gaia BH1 (Kruse & Agol 2014; Kawahara et al. 2018), future missions like Earth 2.0 (Ge et al. 2022) may detect rarer systems containing BHs.", '9.3 Other Gaia astrometric candidates': "Afew other recent works have searched Gaia DR3for dormant BH + normal star binaries. Andrews et al. (2022) analyzed the astrometric binary sample and identified 24 candidates for having BH or neutron star companions. They excluded Gaia BH1 from their sample as probably spurious because its orbital period is close to 3 GLYPH<2> the Gaia precession period. Our follow-up has shown that the Gaia solution was nevertheless robust.  \nThe candidates Andrews et al. (2022) do consider all have rather low masses for BHs. Only one, Gaia DR3 6328149636482597888, has an inferred companion mass above 2 𝑀 GLYPH<12> at the 2 𝜎 level. We have begun spectroscopic monitoring of this source, for which Andrews et al. (2022) reported 𝑀 2 = 2 GLYPH<149> 71 , 1 GLYPH<149> 50 GLYPH<0> 0 GLYPH<149> 36 𝑀 GLYPH<12> (2 𝜎 uncertainties), and  \nwe find it to be a metal-poor turn-off star ( » Fe GLYPH<157> H … GLYPH<25> GLYPH<0> 1 GLYPH<149> 5) on a halo-like orbit. Andrews et al. assumed a mass 𝑀 ★ = 1 GLYPH<149> 21 GLYPH<6> 0 GLYPH<149> 2 (2 𝜎 uncertainties) for the luminous star, but we find this to be an overestimate. Fitting the SED, we find 𝑇 eff = 6150 GLYPH<6> 100 K and 𝑅 = 1 GLYPH<149> 58 GLYPH<6> 0 GLYPH<149> 05 𝑅 GLYPH<12> , which for » Fe GLYPH<157> H … = GLYPH<0> 1 GLYPH<149> 5 MIST models, we find to correspond to 𝑀 ★ = 0 GLYPH<149> 78 GLYPH<6> 0 GLYPH<149> 05 𝑀 GLYPH<12> . Adopting this mass for the luminous star and fitting the Gaia astrometric solution, we find 𝑀 2 = 2 GLYPH<149> 25 , 0 GLYPH<149> 61 GLYPH<0> 0 GLYPH<149> 26 𝑀 GLYPH<12> (2 𝜎 uncertainties), consistent with a neutron star. The source's RVs thus far show a linear trend consistent with the astrometric solution. This is promising, but the Gaia orbital period is 736 days, so RV monitoring over an extended period is required to confirm the reliability of the Gaia solution, and more analysis is required to determine the nature of the companion. A BH seems unlikely.  \nAmong the lower-mass candidates identified by Andrews et al., which they propose the be neutron stars, a majority are significantly blueward of the main-sequence. This could indicate that (a) they have unresolved white dwarf companions, (b) they have low metallicities, or (c) the astrometric solutions are spurious. Our spectroscopic follow-up of these objects has shown that a majority do not have unusually low metallicities. Spurious astrometric solutions or white dwarf companions both remain plausible explanations. The candidates in their sample that are not blueward of the main-sequence are all sufficiently luminous ( 𝐿 & 1 𝐿 GLYPH<12> ) that white dwarf companions would likely escape detection. Among the neutron-star candidates, we find that when reliable mass estimates for the luminous stars are adopted from SED fitting, we in several cases infer companion masses significantly below the Chandrasekhar limit.  \nShahaf et al. (2022) also recently searched for compact object companions among the Gaia DR3 astrometric solutions. They identified 8 sources with inferred companion masses above 2 𝑀 GLYPH<12> , including Gaia BH1. Three of the remaining seven sources, including all the sources with inferred 𝑀 2 above 3 𝑀 GLYPH<12> , are ruled out by our spectroscopic follow-up, and likely have spurious astrometric solutions. The remaining candidates in their sample have sufficiently long periods that the astrometric solutions could not yet be tested robustly. Their sample also contains a few dozen neutron star candidates, most with periods ranging from 400 to 1000 days.  \nOne neutron star candidate, Gaia DR3 5136025521527939072, was presented by Gaia Collaboration et al. (2022a): it has an AstroSpectroSB1 solution which was validated by their RV followup, and also by our follow-up. However, we find that this source is also metal-poor and on a halo-like orbit. Fitting the SED and spectra, we find » Fe GLYPH<157> H … = GLYPH<0> 1 GLYPH<149> 3 GLYPH<150> 𝑇 eff = 6400 K GLYPH<150> 𝑅 = 1 GLYPH<149> 17 𝑅 GLYPH<12> , and 𝑀 ★ = 0 GLYPH<149> 77 GLYPH<6> 0 GLYPH<149> 05 𝑀 GLYPH<12> , lower than the 𝑀 ★ = 1 GLYPH<149> 2 𝑀 GLYPH<12> assumed by Gaia Collaboration et al. (2022a). Adopting this 𝑀 ★ and fitting the Gaia solution, we obtain 𝑀 2 = 1 GLYPH<149> 23 GLYPH<6> 0 GLYPH<149> 06 𝑀 GLYPH<12> , which is lower than the 𝑀 2 = 1 GLYPH<149> 5 𝑀 GLYPH<12> inferred by Gaia Collaboration et al. (2022a), and consistent with a white dwarf. This object and Gaia DR3 6328149636482597888 highlight another potential contaminant for astrometric compact object searches: moderately evolved, low-metallicity stars can be easily mistaken for higher-mass stars near the zero-age main sequence.  \nFurther analysis of the Andrews et al. (2022) and Shahaf et al. (2022) candidates is warranted, as some could indeed host neutron star or very low-mass BH companions. RV follow-up is essential for determining the reliability of the astrometric solutions. As our followup has already shown (Appendix E), the probability of encountering a spurious solution (or simply one that is in the tails of the noise distribution) is high when one selects objects in sparsely-populated regions of parameter space. For companions with masses consistent  \nwith neutron stars, other models (e.g., two white dwarfs) are also harder to rule out.  \nWeconclude that Gaia DR3 astrometric solutions contain one unambigous dormant BH, although we cannot exclude the possibility that up to a few low-mass ( . 3 𝑀 GLYPH<12> ) BHs also exist in the dataset. For BHs above 5 𝑀 GLYPH<12> , Gaia BH1 is the only solid candidate in the astrometric catalog among sources near the main-sequence (see Appendix E; there is one source with a red giant primary that may host a BH companion, but RV followup over a few-year period is required to vet it). Among Gaia DR3 spectroscopic binaries, there are almost certainly no solid BH companions to main-sequence stars (El-Badry & Rix 2022; Jayasinghe et al. 2022), but some neutron star candidates may turn out to be robust. BH companions to giants might still exist in the catalog, but are more difficult to vet.  \nWe note that if any of the neutron star or BH candidates discussed above turn out to be robust, their evolutionary histories will be a puzzle. Like Gaia BH1, they have orbits that are too wide to be understood as a result of a CE event, but too close to accommodate a red supergiant. In addition, their wide orbits could not survive the typical GLYPH<24> 250 km s GLYPH<0> 1 natal kicks expected for neutron stars. The existence of Gaia BH1 and these other candidates might thus point to a different dominant channel for forming moderately wide compact object + normal star binaries, such as the hierarchical triple scenario discussed in Section 8.4.5.  \nWhile this manuscript was under review, a preprint analyzing the same system was posted by Chakrabarti et al. (2022). Their conclusions are broadly consistent with ours. The best-fit BH mass they report, 𝑀 2 = 11 GLYPH<149> 9 , 2 GLYPH<149> 0 GLYPH<0> 1 GLYPH<149> 6 𝑀 GLYPH<12> , is marginally higher than the value we found, 𝑀 2 = 9 GLYPH<149> 62 GLYPH<6> 0 GLYPH<149> 18 𝑀 GLYPH<12> . We suspect that the difference in our results comes mainly from the fact that their RVs do not include any measurements near periastron, and thus do not constrain the RV amplitude well. This also leads to larger uncertainties in the BH mass. Their best-fit solution predicts an RV GLYPH<25> 20 km s GLYPH<0> 1 larger than our measurements near periastron, so we consider their constraints less reliable than those reported here. Nevertheless, it is encouraging that an independent analysis, based on independent RV data, reached similar conclusions about the nature of the Gaia BH1 system.", '10 SUMMARY AND CONCLUSIONS': "We have presented the discovery of Gaia BH1, a new binary system consisting of a GLYPH<24> 10 𝑀 GLYPH<12> dark object orbited by an otherwise unremarkable Sun-like star (Figure 1). The system was first identified as a promising BH candidate based on its Gaia astrometric solution. Our follow-up spectroscopic observations validated the astrometric solution and allowed us to place tight dynamical constraints on the mass of the dark object (Figure 2), which are summarized in Table 1. Even without consideration of the astrometric solution, the RVs imply a companion mass that can only be attributed to a BH (Figure 4; the 5 𝜎 lower limit on the spectroscopic mass function is 𝑓 ' 𝑀 2 ' spec ¡ 3 GLYPH<149> 68 𝑀 GLYPH<12> ). Scenarios not involving a BH are firmly ruled out by the object's large mass and stringent limits on the light contributions of a luminous companion.  \nGaia BH1 is different in several ways from other known BHs in binaries. It is nearby ( 𝑑 = 480 pc) and bright ( 𝐺 = 13 GLYPH<149> 8). Together with the lack of contamination from an accretion disk and the fact that the luminous star is cool and slowly rotating, this makes it possible to study the system in greater detail in the optical than is possible for BH X-ray binaries. For example, we measure chemical abundances for 22 elements in the G star's photosphere (Table 2) - and find them entirely unremarkable. The fact that the orbit is measured with both  \nastrometry and radial velocities (RVs) allows for tight constraints on the BH mass, with no inclination degeneracy. Several high-precision RV instruments are capable of measuring RVs for the G star with GLYPH<24> 1 m s GLYPH<0> 1 precision.  \nThe system's properties are difficult to explain with evolutionary models for isolated binaries. The maximum radius of the BH progenitor is expected to have been much larger than the G star's current orbit, leading to an episode of common envelope evolution in which the G star spiraled in through the BH progenitor's envelope and ultimately ejected it. However, a 1 𝑀 GLYPH<12> star is unlikely to have had enough orbital energy to eject the GLYPH<24> 20 𝑀 GLYPH<12> envelope of a red supergiant, and would be in danger of merging with the core under standard CE assumptions. This is a generic problem for formation models of BH low-mass X-ray binaries (LMXBs). For Gaia BH1, there is an addition challenge: the orbit today is much wider than expected after a common envelope episode, which would be predicted to leave the G star just a few solar radii from the BH progenitor. One cannot invoke strong natal kicks to the BH to widen the orbit again in this scenario, as this would produce a highly eccentric orbit and high space velocity.  \nThesystem's evolution may be better-explained in models in which the G star was initially a wide tertiary companion to a close binary containing two massive stars. In this case, interactions between the two stars could have prevented either one from expanding to become a red supergiant, such that the G star could have formed in an orbit similar to its current orbit (but somewhat tighter) and remained there ever since. In this case, the dark object could be a binary composed of two lower-mass BHs. High-precision RV follow-up offers the tantalizing possibility of testing this scenario. Alternatively, the binary may have formed through an exchange interaction. Its thin-disk-like Galactic orbit, modest eccentricity, and near-solar metallicity suggest that it did not form dynamically in a globular cluster, but dynamical formation in an open cluster that has since dissolved is more plausible.  \nGaia BH1 is about 3 times nearer to Earth than the next-closest BH known, the LMXB A0620-00 (McClintock & Remillard 1986), whose Gaia parallax implies a distance of 1 GLYPH<149> 44 GLYPH<6> 0 GLYPH<149> 25 kpc. This suggests that wide BH binaries like Gaia BH1, while being harder to detect, are significantly more common than BHs in close binaries with ongoing accretion. Given the conservative cuts applied to astrometric binary solutions published in Gaia DR3, it is likely that DR4 and DR5 will identify dozens of similar systems. Gaia BH1 will likely remain one of the nearest and brightest BHs at periods 𝑃 orb . 1000 days, but even nearer BH binaries may be discovered in longer-period orbits, which could not yet be characterized with DR3 data.  \nIt is both interesting and irksome that the first unambiguous BH discovered by Gaia seems to be somewhat of an oddball, and not easily accommodated by standard binary evolution models. One interpretation is that - among solar-type stars, which make up the bulk of the population to which Gaia is sensitive - common envelope evolution leads to the merger of most would-be BH + normal star binaries. In this case, only systems formed through alternative channels will survive. It is clear that BH X-ray binaries - which today are more compact versions of systems like Gaia BH1 - only narrowly escape a merger if they form through a common envelope, and represent a rare outcome of binary evolution. BH + normal star binaries with 𝑃 orb & 10 years can avoid a common envelope, but (a) Gaia is not yet sensitive to them, and (b) they may be disrupted by natal kicks. BH + massive star binaries are more likely to both survive common envelope and undergo stable mass transfer, but the number of massive stars accessible with Gaia astrometric orbits is thus far  \nmodest. In any case, population demographics of the detached BH + normal star binary population will become clearer in the next decade, as more systems are discovered.  \nBecause of its long orbital period and close distance, the angular size of Gaia BH1's orbit on the sky is more than 10 times larger than any other known BH binary. This makes the system an excellent target for interferometric follow-up. Observations at X-ray and radio wavelengths are also warranted. Although the accretion rate of the G star's wind by the BH is expected to be low and any detectable emission would be faint (Section 7), Gaia BH1 provides a unique opportunity to test imperfectly-understood models of accretion flows at very low accretion rates.  \nIn the future, the G star will become a red giant. After an X-ray bright period during which it transfers mass to the BH - perhaps appearing as a wider version of the red giant + BH binary GRS 1915+105 - the system will terminate its evolution as a wide BH + white dwarf binary. If the system turns out to host an inner BH binary, this evolution may be interrupted by its merger, with the resulting recoil potentially disrupting the outer binary.", 'ACKNOWLEDGEMENTS': "WethankMorganMacleod,E.SterlPhinney, Rohan Naidu, Lieke van Son, Mathieu Renzo, Tom Maccarone, Josh Simon, and Jeff Andrews for helpful discussions, and Vedant Chandra for help with data reduction. We are grateful to Yuri Beletsky, Sam Kim, Angela Hempel, and Régis Lachaume for observing help. HWR acknowledges the European Research Council for the ERC Advanced Grant [101054731]. This research made use of pystrometry, an open source Python package for astrometry timeseries analysis (Sahlmann 2019). This work made use of Astropy, 6 a community-developed core Python package and an ecosystem of tools and resources for astronomy (Astropy Collaboration et al. 2022).  \nThis project was developed in part at the Gaia Fête, held at the Flatiron Institute's Center for Computational Astrophysics in June 2022, and in part at the Gaia Hike, held at the University of British Columbia in June 2022. This work has made use of data from the European Space Agency (ESA) mission Gaia ( https: //www.cosmos.esa.int/gaia ), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos. esa.int/web/gaia/dpac/consortium ). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.  \nBased on observations obtained at the international Gemini Observatory, a program of NSF's NOIRLab, which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation on behalf of the Gemini Observatory partnership: the National Science Foundation (United States), National Research Council (Canada), Agencia Nacional de Investigación y Desarrollo (Chile), Ministerio de Ciencia, Tecnología e Innovación (Argentina), Ministério da Ciência, Tecnologia, Inovações e Comunicações (Brazil), and Korea Astronomy and Space Science Institute (Republic of Korea).", 'DATA AVAILABILITY': 'Data used in this study are available upon request from the corresponding author.', 'REFERENCES': '701454 (  \n), doi:10.1117/12.789972', 'A1 MagE': 'We observed Gaia BH1 4 times using the Magellan Echellette spectrograph (MagE; Marshall et al. 2008) on the 6.5m Magellan Clay telescope at Las Campanas Observatory. All observations were carried out with the 0.7 arcsec slit and exposure times ranging from 300 to 600 seconds. This yielded a typical SNR of 100 per pixel at 5,500 Å, spectral resolution 𝑅 GLYPH<25> 5 GLYPH<150> 400, and wavelength coverage of 3 GLYPH<150> 500 GLYPH<0> 11 GLYPH<150> 000 Å. To minimize drift in the wavelength solution due to flexure, we obtained a ThAr arc on-sky immediately after each science exposure.  \nWe reduced the spectra using PypeIt (Prochaska et al. 2020), which performs bias and flat field correction, cosmic ray removal, wavelength calibration, sky subtraction, extraction of 1d spectra, and heliocentric RV corrections. To assess the stability of the wavelength solution, we observed several standard stars, measured their RVs using the same procedure we used for the science targets, and compared the results to their measured Gaia DR3RVs. We found that the MagE-measured RVs were 1 GLYPH<149> 5 GLYPH<6> 2 km s GLYPH<0> 1 larger than those measured by Gaia , so we adopted a conservative uncertainty of 3 km s GLYPH<0> 1 for all the MagE RVs.', 'A2 GMOS': 'Weobserved Gaia BH1 4 times using the Gemini Multi-Object Spectrograph (GMOS; Hook et al. 2004) on the 8.1m Gemini-North telescope on Mauna Kea (program GN-2022B-DD-202). We used the R831\\_G5302grating with a 0.5 arcsecond slit and central wavelength of 5200 Å, leading to wavelength coverage from 4000 to 6400 Å and resolution 𝑅 GLYPH<25> 4 GLYPH<150> 400. We obtained a CuAr arc on-sky immediately  \nafter each science exposure. We used 342 second exposures, yielding a typical SNR of 90 per pixel at 6000 Å.  \nWe reduced the data using PypeIt , which required construction of a new template for the R831 grating at blue wavelengths. Crosscorrelation with a template yielded formal RV uncertainties of order 1 km s GLYPH<0> 1 , but we found variations of up to 10 km s GLYPH<0> 1 between RVs measured at the blue and red ends of the spectrum. We adopted a RV uncertainty of 4 km s GLYPH<0> 1 for the GMOS RVs, which was also validated by observations of Gaia BH1 with higher-resolution instruments that occurred nearly concurrently with some of the GMOS observations.', 'A3 X-shooter': 'We observed Gaia BH1 4 times using the X-shooter spectrograph (Vernet et al. 2011) on the 8.2m UT3 telescope at the VLT on Cerro Paranal (program 2109.D-5047). We used the 0.5, 0.4, and 0.4 arcsec slits on the UVB, VIS, and NIR arms. This setup yielded a resolution of 𝑅 GLYPH<25> 9 GLYPH<150> 700 in the UVB arm, 𝑅 GLYPH<25> 18 GLYPH<150> 400 in the VIS arm, and 𝑅 GLYPH<25> 11 GLYPH<150> 600 in the NIR arm, with near-continuous wavelength coverage from 3100-24000 Å. We reduced the spectra using pypeit .  \nWe generally measured RVs using the NIR spectra, whose wavelength solution is calculated from telluric OH lines in the science spectra, yielding a maximally stable wavelength solution. For the first epoch, we inadvertently used a 0.9 arcsec slit in the NIR arm, yielding a lower spectral resolution. We measured the RV from the VIS arm for this epoch. We adopt an uncertainty of 1 km s GLYPH<0> 1 on all the X-shooter RVs.', 'A4 HIRES': 'We observed Gaia BH1 9 times using the High Resolution Echelle Spectrometer (HIRES; Vogt et al. 1994) on the 10m Keck I telescope on Maunakea. The data were obtained and reduced using the standard California Planet Survey setup (CPS; Howard et al. 2010), including use of the C2 decker (0.86 arcseconds GLYPH<2> 14 arcseconds), which yields spectra with 𝑅 GLYPH<25> 55 GLYPH<150> 000 and wavelength coverage over most of 3700-8000 Å, with 2 gaps 100-200 Å wide. We used 600 second exposures, yielding SNR of order 40 per pixel at 6000 Å. The CPS reduction includes sky-subtraction using the long C2 decker.  \nAs part of the CPS pipeline, absolute RVs are measured using the method described by Chubak et al. (2012): to correct for instrumental shifts in the wavelength zero-point, the telluric A and B molecular oxygen absorption bands of each star are cross-correlated with those of RV standard stars, whose adopted RVs are matched to the scale of Nidever et al. (2002) and other RV standard catalogs. The resulting RVs are robust to 0 GLYPH<149> 1 km s GLYPH<0> 1 . We use these RVs throughout our analysis. Most of our observations were performed through an iodine cell, making it possible to measure higher-precision RVs from these data in the future, after a higher-SNR template spectrum has been obtained and the method of Butler et al. (1996) can be applied.', 'A5 FEROS': 'We observed Gaia BH1 17 times with the Fiberfed Extended Range Optical Spectrograph (FEROS; Kaufer et al. 1999) on the 2.2m ESO/MPG telescope at La Silla Observatory. We used 2 GLYPH<2> 2 binning in order to reduce readout noise (allowing us to measure RVs down to 𝐺 GLYPH<25> 16); this resulted in an effective resolution 𝑅 GLYPH<25> 37000 measured from sky lines. We took 2400 second exposures and achieved a typical SNR of 20 per pixel. We reduced the data with the standard ESOMIDASpipeline, which performs bias-subtraction, flat fielding,  \nwavelength calibration, and sky subtraction, and merges the spectra from different orders.  \nOn some nights, the wavelength solution drifted from the one established by the afternoon arcs by up to a few km s GLYPH<0> 1 . We corrected this by cross-correlating an arc spectrum (which is taken during each exposure with a second fiber) with a reference arc spectrum taken on the first night of our observations. Based on observations of an RV standard with the same setup, we adopted and uncertainty of 0 GLYPH<149> 3 km s GLYPH<0> 1 for all the FEROS observations, except in one case where poor conditions yielded a formal uncertainty of 0 GLYPH<149> 5 km s GLYPH<0> 1 . Nearconcurrent FEROS and HIRES observations of the source allowed us to verify that the two RV zeropoints are consistent within this tolerance.', 'A6 ESI': 'We observed Gaia BH1 one time with the Echellette Spectrograph and Imager (ESI; Sheinis et al. 2002) on the 10m Keck II telescope on Maunakea. We used a 300 second exposure with the 0.3 arcsec slit, yielding a resolution 𝑅 GLYPH<25> 12000 and SNR = 60, with useful wavelength coverage of 3900-10000 Å. We reduced the data using the MAuna Kea Echelle Extraction (MAKEE) pipeline, which performs bias-subtraction, flat fielding, wavelength calibration, and sky subtraction, and we refined the wavelength solution using sky lines. Based on RV measurements of RV standards on the same night, we adopted an uncertainty of 1 GLYPH<149> 5 km s GLYPH<0> 1 .', 'APPENDIX B: GAIA OBSERVATIONS': "As discussed in Section 4.2.1, the astrometric uncertainties for the parameters describing the orbit of Gaia BH1 are unusually large for the star's apparent magnitude. Here we show that this is likely a consequence of the Gaia scanning law.  \nEpoch-level astrometric data is not published in DR3, but one can determine when a source was observed using the Gaia observation scheduling tool (GOST). 7 Given a source position, GOST returns a list of observation times when the scanning law predicts that a source will transit across the Gaia focal plane. It is not guaranteed that Gaia actually obtains data on each source at the predicted times, as the tool does not account for gaps between CCDs and issues that cause temporary gaps in the datastream, such as micrometeor impacts. However, for Gaia BH1, GOST predicts 21 visibility periods (i.e., groups of observations separated from other groups by at least 4 days), and the gaia\\_source table reports that 21 visibility periods were indeed used in the astrometric solution. This suggests that the scan times predict by GOST are likely a good approximation of reality.  \nWe visualize these times in Figure B1. It is clear that although Gaia DR3 covered more than 5 orbital period of Gaia BH1, the phase coverage is far from uniform, and only about half of the orbit in phase was covered by astrometric data. This may be a consequence of the fact that the orbital period is a near-integer multiple of the Gaia precession period, and close to half a year. Our RV followup shows that despite this limitation, Gaia was able to measure an astrometric orbit that was basically correct. However, the incomplete phase coverage make it unsurprising that the uncertainties in the astrometric parameters are larger than usual for a source of this brightness.  \nGOST predicts that the north side of the orbit will remain poorly sampled even in Gaia DR4. We note that a few strategically timed observations with HST , JWST , GRAVITY+, or MICADO could fill in this part of the orbit and yield improved astrometric precision if analyzed together with the epoch-level Gaia astrometric data that will become available in DR4.", 'APPENDIX C: COMPARISON OF THE ASTROMETRIC AND RV SOLUTIONS': 'As discussed in Section 4, the tightest constraints on the Gaia BH1 orbit can be obtained through joint fitting of the RVs and astrometric constraints. However, it is also useful to compare the RVs predicted by the Gaia solution alone to the measured RVs, in order to assess whether or not the astrometric solution is consistent with the RVs. Figure C1 compares the measured RVs to predictions based on the Gaia astrometric solution, which are generated from posterior draws of the fit without RVs. Because the astrometric solution does not constrain the center-of-mass velocity, we set it to 45 km s GLYPH<0> 1 (consistent with the measured value) in all cases.  \nThe measured RVs are basically consistent with the Gaia -only solution, in the sense that some posterior draws go through all the RVs. Figure C2 and Table C1 also show that all inferred parameters are consistent within . 2 𝜎 . This suggests that the Gaia solution and its uncertainties are reliable, despite the incomplete orbit coverage discussed in Appendix B.', 'APPENDIX D: RADIAL VELOCITIES': 'Radial velocities are listed in Table D1. All times are reported as midexposure times in the heliocentric frame (HJD UTC). RVs are also in the heliocentric frame and their uncertainties include uncertainties in wavelength calibration. The reported SNR is at the wavelength range we used for RV measurement, which in most cases was near 5400 Å.', 'APPENDIX E: OTHER CANDIDATES': "Figure E1 illustrates our strategy for selecting BH candidates among the Gaia astrometric solutions. We show that binaries containing a BH and a low-mass main-sequence star are expected to have much larger photocenter wobble than binaries containing two luminous stars.  \nThe left panel shows a simulation. Here we begin with all sources in the gaia\\_source catalog with parallax\\_over\\_error > 50. We then attached a simulated binary population to these sources, drawing periods from the canonical lognormal distribution for solar-type stars (Raghavan et al. 2010) and mass ratios from a uniform distribution. For each simulated binary, we then calculate the expected photocenter semi-major axis and its uncertainty, and we retain only the binaries that satisfy the cuts imposed on the binary solutions published in DR3. We also generate a population of hierarchical triples, drawing both the inner and outer mass ratios from uniform distributions. We calculate the flux ratios using the 𝐺 GLYPH<0> band mass-magnitude relation tabulated by Janssens et al. (2022). Finally, we generate a population of BH + normal star and neutron star + normal binaries, by assigning 10 𝑀 GLYPH<12> and 1 GLYPH<149> 4 𝑀 GLYPH<12> dark companions to a random subset of the luminous sources, again drawing periods from the canonical lognormal. This is not meant to represent a realistic distribution of orbital periods, masses, or mass ratios, but only to show the expected photocenter orbit size at a given period.  \nIn the right panel, we show the same observables for all objects in DR3 with astrometric solutions ( Orbital and AstroSpectroSB1 solutions). As expected, most solutions fall in the region of parameter space that we expect to be populated by luminous binaries and triples. For each binary, we calculate ˜ 𝑀 , a quantity which is related to, but not equal to, the total dynamical mass:  \n˜ 𝑀 = 4 𝜋 2 𝑃 2 orb 𝐺 GLYPH<16> 𝑎 0 𝜛 GLYPH<17> 3 GLYPH<149> (E1)  \nFor a binary with a luminous star and a dark companion, 𝑎 0 𝜛 = 𝑎 𝑞 1 , 𝑞 , where 𝑞 = 𝑀 2 GLYPH<157> 𝑀 ★ and 𝑎 is the semi-major axis. Thus, ˜ 𝑀 =  \n𝑀 tot GLYPH<16> 𝑞 1 , 𝑞 GLYPH<17> 3 . For BH companions to low-mass stars, 𝑞 GLYPH<157> ' 1 , 𝑞 ' GLYPH<25> 1, and thus ˜ 𝑀 is only slightly less than 𝑀 tot. On the other hand, for a binary with two luminous stars, 𝑎 0 𝜛 = 𝑎𝛿 𝑞ℓ , where  \n𝛿 𝑞ℓ = GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> 𝑞 GLYPH<157> ' 1 , ℓ ' GLYPH<0> ℓ GLYPH<157> ' 1 , ℓ ' ' 1 , 𝑞 ' GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<150> (E2)  \nand ℓ = 𝐹 2 GLYPH<157> 𝐹 1 is the ratio of the two stars' 𝐺 GLYPH<0> band fluxes (e.g. Penoyre et al. 2022). For binaries with two main-sequence stars, 𝛿 𝑞ℓ . 0 GLYPH<149> 35 (e.g. Shahaf et al. 2019), and thus in general ˜ 𝑀 . 0 GLYPH<149> 05 𝑀 tot. ˜ 𝑀 is thus a good quantity with which to select massive dark companions to low-mass stars, though more sophisticated methods exist that are sensitive to lower-mass BHs and BH companions to massive stars (e.g. Shahaf et al. 2019, 2022).  \nWith this in mind, we selected objects from Gaia DR3 with ˜ 𝑀 ¡ 3 𝑀 GLYPH<12> . This yielded 6 objects, which are highlighted in the right panel of Figure E1 and enumerated in Table E. Sources brighter than 𝐺 RVs = 12 have multi-epoch RVs measured (but not yet published) by Gaia . The min-to-max range of these RVs is published as rv\\_amplitude\\_robust . Comparing this to the expected full RV amplitude implied by the Gaia solution provides a useful diagnostic of whether the solution is likely to be correct.  \nEach solution is accompanied by a goodness\\_offit value. As discussed by Halbwachs et al. (2022), this quantity is formally expected to follow a normal distribution, N' 0 GLYPH<150> 1 ' , with values significantly larger than 1 indicative of a problematic solution. However, this may not hold in practice, since the uncertainties in epoch-level measurements can be underestimated. We find that the distribution of goodness\\_offit values displays a sharp discontinuity at 𝐺 = 13. Fainter values approximately follow a N' 0 GLYPH<150> 1 GLYPH<149> 5 ' distribution, such that a value much larger than 1 is indeed unusual. However, for brighter solutions, the distribution peaks at GLYPH<25> 3 GLYPH<149> 5, and 13% of solutions of have goodness\\_offit > 10. We thus expect significantly larger goodness\\_offit values to be (potentially) consistent with a good solution for bright sources.  \nWe comment briefly on each source:  \n3640889032890567040 : Poor GoF and implausibly high companion mass. Ruled out by low rv\\_amplitude\\_robust .  \n4467000291193143808 : Period of 1647 days is significantly longer than Gaia DR3 baseline. Implausibly high companion mass; poor GoF for the source magnitude.  \n4373465352415301632 : The subject of this paper.  \n6281177228434199296 : Ruled out by our RV follow-up. Low rv\\_amplitude\\_robust also makes the astrometric solution implausible.  \n5870569352746779008 : Most plausible candidate beside Gaia BH1, but still uncertain. 𝑃 orb is longer than the Gaia DR3 observing baseline. The GoF is significantly larger than 1, but is not unusual for sources in the relevant range of 𝐺 magnitude. The CMD position is consistent with a lower giant, or possibly a red clump star.  \n2  \nFigure B1. Predicted observation times of Gaia BH1 from the Gaia observation scheduling tool (GOST). Black line shows the best-fit photocenter orbit. Red points show the predicted photocenter positions at the times when GOST predicts Gaia would have observed the source. Note that we do not have access to the actual measured Δ RA and Δ Dec values; only to the predicted scan times. GOST predicts that Gaia observed the source in 21 visibility periods which all fall on one side of the orbit. This may be a result of the fact that the orbital period is a near-integer multiple of Gaia's 63-day precession period, and almost exactly half a year.  \n<!-- image -->  \nTable C1. Comparison of the Gaia astrometric constraints on Gaia BH1's photocenter ellipse to constraints from the joint RVs+astrometry fit. All parameters are consistent at the 1.6 𝜎 level.  \nThe inclination is relatively low (35 GLYPH<6> 1 deg), so the nature of the companion hinges on the correctness of the astrometric solution. That is, if the Gaia solution is correct, RVs alone will only yield a spectroscopic mass function near 1 𝑀 GLYPH<12> , which is not high enough to unambiguously rule out a luminous companion to a red giant. The reliability of the astrometric solution will be easier to assess in DR4, when epoch-level astrometric data are published and a full orbit has been covered. The measured rv\\_amplitude\\_robust is not too far from the expected peak-to-peak amplitude. We have initiated a spectroscopic follow-up campaign, on which we will report in future work. Thus far, the RVs are consistent with the AstroSpectroSB1 solution. Further RV monitoring is warranted and necessary. This source was recently investigated in detail by Tanikawa et al. (2022), who found it to be the only credible BH candidate in Gaia DR3 with both astrometric and spectroscopic Gaia data, but did not carry out follow-up. If the solution turns out to be credible, the source will present a similar puzzle to Gaia BH1: the BH progenitor would likely not have fit inside the current orbit, but the orbit is much too wide to be a result of common envelope evolution. The current orbital separation is such that if the luminous star is a core helium burning red clump star, mass transfer via winds and/or stable Roche lobe overflow would have occurred prior to the helium flash.  \n3664684869697065984 : Ruled out by APOGEE RVs.  \nFor a 10 𝑀 GLYPH<12> BH companion, the cut of ˜ 𝑀 ¡ 3 𝑀 GLYPH<12> would be satisfied for luminous stars with 𝑀 ★ . 8 𝑀 GLYPH<12> . For highermass luminous stars, the AMRF, A statistic designed by Shahaf et al. (2019) is more sensitive. We verified that no objects in the  \ngaiadr3.nss\\_two\\_body\\_orbit catalog with photometric masses above 5 𝑀 GLYPH<12> have A ¡ 0 GLYPH<149> 6, as would be required to identify a dark companion with high confidence.", 'APPENDIX F: HIRES CANNON LABELS': 'The labels inferred by the Cannon fit to the HIRES spectrum of Gaia BH1 are listed in Table F1. The method is not expected to yield reliable formal uncertainties on labels, because at high SNR the errors are dominated by systematics in the spectral model rather than photon noise. We thus report the cross-validation error achieved by Rice & Brewer (2020) for each label when testing the model, which provides a rough lower limit on the uncertainties in each label.  \nThe atmospheric parameters and abundances reported by the Cannon are quite similar to those we find with BACCHUS (Section 5.1), which we adopt as our fiducial values. Note that the Cannon reports [X/H], while for the BACCHUS fit we reported [X/Fe] = [X/H] [Fe/H].', 'APPENDIX G: DEPENDENCE OF GAIA ASTROMETRIC UNCERTAINTIES ON OBSERVABLES': "In Figure G1, we show the mean uncertainty in 𝑎 0 and 𝜛 for DR3 astrometric binary solutions as a function of apparent magnitude and 𝑎 0 . We only include binaries with best-fit 𝑃 orb GLYPH<159> 1000 days, which  \nFigure C1. Comparison of measured RVs (same as Figure 2) to predictions from the Gaia astrometric solution alone. The center-of-mass RV (which is not constrained by astrometry) is set to the value inferred from RVs, but the cyan lines otherwise do not know about the measured RVs. The uncertainty in the predictions is larger than when RVs are included in the fit, but the measured RVs are generally consistent with the predictions, validating the Gaia solution.  \n<!-- image -->  \nis roughly the time baseline of observations used to calculate DR3 solutions. The figure shows that typical astrometric uncertainties especially the uncertainty in 𝑎 0 - depend both on apparent magnitude and on 𝑎 0 . At fixed apparent magnitude, binaries with 𝑎 0 GLYPH<25> 3 mas have uncertainties in 𝑎 0 and 𝜛 that are respectively about 2 GLYPH<149> 5 GLYPH<2> and 1 GLYPH<149> 7 GLYPH<2> larger than sources with 𝑎 0 GLYPH<25> 0 GLYPH<149> 3 mas.  \nIt is not immediately obvious whether this trend reflects an actual dependenceoftheastrometric uncertainties on 𝑎 0 amongallsolutions that were initially fit, or is mainly a result of the astrometric SNR cuts imposed on the published solutions (Equations 10-12). Figure G2 shows that removal of solutions with larger uncertainties and short periods likely has a significant effect on trends with 𝑎 0 : there is a 'wedge' of missing solutions with small 𝑎 0 and larger uncertainties.  \nTo determine whether an inherent dependence of the astrometric uncertainties on 𝑎 0 (or another parameter that scales with it) is required to explain the trends in Figure G1, we preformed the following experiment. Using the simulation setup described in Appendix E, we begin with all sources in the gaia\\_source catalog with parallax\\_over\\_error > 10 and create a simulated binary population that we then subject to the DR3 SNR cuts. To predict the uncertainties in 𝑎 0 and 𝜛 for each simulated binary, we first tabulate the mean and standard deviation of 𝜎 𝑎 0 and 𝜎 𝜛 as a function of apparent 𝐺 GLYPH<0> band magnitude, multiply this mean by an 𝑎 0 -dependent  \nscaling factor , and then draw uncertainties for each source from a Gaussian with the appropriate mean and standard deviation. Finally, we compare the simulated distribution of 𝜎 𝑎 0 and 𝜎 𝜛 for binaries surviving the quality cuts to the observed distribution shown in Figure G1.  \nWe parameterize the scaling factor as 𝑓 = 𝑎 log ' 𝑎 0 GLYPH<157> mas ' , 𝑏 . We find that when 𝑎 = 0 and 𝑏 = 1 (i.e., no dependence of uncertainties on 𝑎 0 ), the predicted correlation between 𝑎 0 and the astrometric uncertainties is weaker than observed. We can best reproduce the observed distribution for 𝜎 𝑎 0 when 𝑎 = 1 GLYPH<149> 3 and 𝑏 = 1 GLYPH<149> 4; for 𝜎 𝜛 , we find 𝑎 = 0 GLYPH<149> 64 and 𝑏 = 1 GLYPH<149> 2. We adopt these values when predicting the expected astrometric uncertainties of a simulated binary population in Section 9.  \nThis paper has been typeset from a T E X/L A T E X file prepared by the author.  \n/circledot  \n/circledot  \nFigure C2. Comparison between constraints obtained from joint fitting of RVs and astrometry (black) and constraints from astrometry alone (cyan). The two sets of constraints are generally consistent, but those with RVs are tighter. Numbers on the diagonals reflect the joint astrometry + RVs constraints. Both sets of constraints are listed in Table 1.  \n<!-- image -->  \nTable E1. BH candidates highlighted in Figure E1. GoF is the goodness\\_offit parameter, for which large values indicate a potentially problematic solution (see text for details). rv\\_amplitude\\_robust is the min-to-max range in the individual-epoch Gaia RV measurements (only available for bright sources); 'expected RV amplitude' is the min-to-max RV amplitude predicted by the astrometric solution. Sources are discussed individually in the text.", 'Table D1. Radial velocities.': '<!-- image -->  \nFigure E1. Parameter space within which we selected BH candidates. Left panel shows a simulation, in which we model a population of luminous (mainsequence) binaries, triples, and compact object + normal star binaries subject to the cuts imposed on binary solutions published in DR3. The relative number of objects in the three populations is chosen arbitrarily. Binaries containing dark companions fall at larger separations at fixed period than binaries containing two (or three) luminous stars. Right panel shows the same quantities for binaries published in DR3. Red star shows Gaia BH1; black points show other initially promising candidates, which we discuss in the text.  \n<!-- image -->  \nFigure G1. Mean uncertainty in photocenter semi-major axis ( 𝑎 0; left) and parallax (right) for all Gaia DR3 astrometric solutions with 𝑃 orb GLYPH<159> 1000 days. Fainter sources and sources with large 𝑎 0 have larger astrometric uncertainties. The correlation between 𝑎 0 and the astrometric uncertainties appears to be only partially a result of period-dependent quality cuts imposed on sources in DR3 (see text). We incorporate an explicit 𝑎 0-dependence when predicting the astrometric uncertainties of a hypothetical binary observed by Gaia .  \n<!-- image -->  \nFigure G2. Uncertainties in parallax (top) and photocenter semi-major axis ( 𝑎 0; bottom) for all sources with astrometric solutions published in DR3. Left column shows uncertainties as a function of 𝑎 0; right column shows them as a function of 𝑃 orb . Red star shows Gaia BH1. A correlation is apparent between the astrometric uncertainties and 𝑎 0. This is in part because sources with small 𝑎 0 and large uncertainties were removed from the catalog through period-dependent quality cuts (Equations 10-12). More structure is apparent as a function of 𝑃 orb : binaries with periods close to 0.5 years and 1 year, as well as harmonics of the 63 day precession period, have larger uncertainties. Uncertainties in 𝑎 0 also increase at 𝑃 orb ¡ 700 days, likely due to incomplete phase coverage.  \n<!-- image -->  \nTable F1. Parameters of the G star, from the HIRES cannon. The crossvalidation error is taken from Rice & Brewer (2020) and represents an estimate of the uncertainty in each label.'}