A red giant orbiting a black hole

ABSTRACT

We report spectroscopic and photometric follow-up of a dormant black hole (BH) candidate from Gaia DR3. The system, which we call Gaia BH2, contains a ∼1 M_⊙ red giant and a dark companion with mass |$M_2 = 8.9\pm 0.3\, {\rm M}_{\odot }$| that is very likely a BH. The orbital period, P_orb = 1277 d, is much longer than that of any previously studied BH binary. Our radial velocity (RV) follow-up over a 7-month period spans >90 per cent of the orbit’s RV range and is in excellent agreement with the Gaia solution. UV imaging and high-resolution optical spectra rule out plausible luminous companions that could explain the orbit. The star is a bright (G = 12.3), slightly metal-poor (⁠|$\rm [Fe/H]=-0.22$|⁠) low-luminosity giant (⁠|$T_{\rm eff}=4600\, \rm K$|⁠; |$R = 7.8\, R_{\odot }$|⁠; |$\log \left[g/\left({\rm cm\, s^{-2}}\right)\right] = 2.6$|⁠). The binary’s orbit is moderately eccentric (e = 0.52). The giant is enhanced in α-elements, with |$\rm [\alpha /Fe] = +0.26$|⁠, but the system’s Galactocentric orbit is typical of the thin disc. We obtained X-ray and radio non-detections of the source near periastron, which support BH accretion models in which the net accretion rate at the horizon is much lower than the Bondi–Hoyle–Lyttleton rate. At a distance of 1.16 kpc, Gaia BH2 is the second-nearest known BH, after Gaia BH1. Its orbit – like that of Gaia BH1 – seems too wide to have formed through common envelope evolution. Gaia BH1 and BH2 have orbital periods at opposite edges of the Gaia DR3 sensitivity curve, perhaps hinting at a bimodal intrinsic period distribution for wide BH binaries. Dormant BH binaries like Gaia BH1 and Gaia BH2 significantly outnumber their close, X-ray bright cousins, but their formation pathways remain uncertain.

1 INTRODUCTION

The Milky Way is very likely teeming with stellar-mass black holes (BHs). We can see today that it contains a few |$\times \, 10^4$| O stars with masses |$M\gtrsim 20\, {\rm M}_{\odot }$| (Garmany, Conti & Chiosi 1982; Reed 2003), a majority of which will likely leave behind BHs when they die. The lifetime of these massive stars is of the order of 10 Myr, or 0.1 per cent the age of the Milky Way. Notwithstanding variations in the Galactic star formation rate with cosmic time, this implies that the Milky Way should contain ∼1000 BHs for every luminous BH progenitor alive today, i.e. a few |$\times \, 10^7$| BHs.

Only a tiny fraction of this BH population has been observed. Observations of X-ray binaries have dynamically confirmed the presence of a BH in ∼20 systems (Remillard & McClintock 2006), and another ∼50 X-ray transients are suspected to contain BHs (Corral-Santana et al. 2016). In all these systems, there is ongoing mass transfer from a luminous star on to a BH, giving rise to an accretion disc that shines brightly in X-rays. BH X-ray binaries are an intrinsically rare outcome of binary evolution (e.g. Portegies Zwart, Verbunt & Ergma 1997), and empirical population models suggest that only about 1000 exist in the Milky Way (Corral-Santana et al. 2016). They nevertheless make up the bulk of the observed BH population, because they are easier to find than BHs that are not accreting.

The Gaia mission has opened a new window on the Milky Way’s binary star population. By precisely measuring the positions of ∼2 billion sources over time, Gaia can detect subtle wobbles caused by the gravitational effects of binary companions – even companions that do not emit any light themselves. As of the mission’s third data release (‘DR3’; Gaia Collaboration 2022c), nearly 500 million sources have been observed in more than 20 well-separated epochs. These data are ultimately expected to yield astrometric binary solutions for of the order of 10⁷ binaries (e.g. Söderhjelm 2004), but the quality cuts applied to solutions included in DR3 were stringent, so that only 1.7 × 10⁵ full astrometric orbital solutions have been published so far. The DR3 binary catalogue nevertheless represents more than an order of magnitude increase in sample size over previous samples of binary orbits and is therefore a promising data set to search for rare objects.

Gaia DR3 has thus far resulted in the identification of one unambiguous BH binary (Chakrabarti et al. 2022; El-Badry et al. 2023). That system, Gaia BH1, consists of a Sun-like star in a wide (P_orb = 186 d; a = 1.4 au) orbit with a |$\sim 10\, {\rm M}_{\odot }$| BH companion. The binary’s wide orbit is not easily explained by standard binary evolution models. At a distance of only d = 480 pc, Gaia BH1 is a factor of ∼3 closer than the next-nearest known BH, and so the Copernican principle suggests that widely separated BH binaries similar to Gaia BH1 are likely common compared to BHs in X-ray binaries. Most theoretical models also predict that wide binaries with non-accreting BHs should significantly outnumber their accreting cousins (e.g. Breivik, Chatterjee & Larson 2017; Langer et al. 2020).

Gaia data are well suited to identify such systems. While Gaia BH1 was the most obvious BH binary candidate published in DR3, El-Badry et al. (2023) also identified a second promising candidate, whose longer orbital period demanded longer term spectroscopic follow-up than was necessary for Gaia BH1. This paper is focused on that second object, which we refer to as Gaia BH2. Similar to Gaia BH1, the system contains a solar-mass star in a wide orbit around an object we suspect to be a BH. Unlike in Gaia BH1, the luminous star has left the main sequence and is ascending the giant branch. At a distance of 1.16 kpc, the object is likely the second-nearest known BH, after Gaia BH1.

The remainder of this paper is organized as follows. Section 2 describes the system’s discovery and Gaia orbital solution. Most of the observational data and modelling are described in Section 3, including light curves (Section 3.1), extinction (Section 3.2), UV photometry (Section 3.3), spectral energy distribution (SED) modelling (Section 3.4), mass inference of the giant (Section 3.5), spectroscopy (Section 3.6), orbit modelling (Section 3.7), spectral analysis (Section 3.8), Galactic orbit (Section 3.9), X-ray observations (Section 3.10), and radio observations (Section 3.11). We discuss the nature of the unseen companion in Section 4 and its possible formation history in Section 5. Prospects for detecting and characterizing additional BHs in binaries with Gaia are discussed in Section 6, and we conclude in Section 7. The appendices provide additional details about several aspects of the data and modelling.

2 DISCOVERY

Being a bright optical source (G = 12.3), Gaia BH2 appears in many archival photometric catalogues. To our knowledge, no spectroscopic observations of the source occurred prior to the Gaia observations. An orbital solution with a best-fitting period of 1352 d – somewhat longer than the ∼1000 d time baseline for observations included in Gaia DR3 – was published in Gaia DR3 (Gaia Collaboration 2022a) with source id 5870569352746779008. We discuss the cadence and phase coverage of the Gaia observations in Appendix  A, and the astrometric goodness of fit in Appendix  B. Both are unproblematic.

Unlike Gaia BH1, which has a purely astrometric solution (nss_solution_type = Orbital), Gaia BH2 has a joint astrometric + single-lined spectroscopic solution (nss_solution_type = AstroSpectroSB1). The reason for this difference is that individual-epoch radial velocities (RVs) are only derived for sources brighter than G_RVS = 12, where G_RVS is a source’s magnitude in the Gaia RVS bandpass. Gaia BH2 has G_RVS = 11.11, and RVs were measured (but not published) in 34 unique transits. Gaia BH1 has G_RVS = 12.83, and thus no individual-epoch RVs.

We identified Gaia BH2 as a promising BH candidate in the course of a broader program to spectroscopically follow-up Gaia DR3 binaries suspected to contain compact objects. The astrometric selection that first brought the system to our attention is described in Appendix E of El-Badry et al. (2023). The Gaia solution has a period of P_orb = 1352 ± 46 d. The astrometric mass function implied by the astrometric data alone is |$f\left(M_{2}\right)_{{\rm ast}}=\left(a_{0}/\varpi \right)^{3}\left(P_{{\rm orb}}/{\rm yr}\right)^{-2}=6.72\pm 0.5\, {\rm M}_{\odot }$|⁠; here a₀ and ϖ represent the semimajor axis of the photocentre and the parallax. Assuming a mass of |$1\, {\rm M}_{\odot }$| for the luminous star (Section 3.5), this implied an unseen companion mass of |$M_2 \approx 8.4\pm 0.5\, {\rm M}_{\odot }$|⁠. Since this far exceeds the apparent photometric mass of the luminous primary, we initiated a spectroscopic follow-up campaign, obtaining our first spectrum of the source in 2022 August.

In El-Badry et al. (2023), we reported that early follow-up data were consistent with the Gaia solution. However, since the data in-hand then covered only a small fraction of the orbit’s predicted RV range, we cautioned that continued RV monitoring was required to better assess the reliability of the Gaia solution. We have now carried out the necessary observations.

Gaia data for the source were also analysed by Tanikawa et al. (2022), who were the first to draw attention to the source as a BH candidate in a 2022 September preprint. They concluded that the Gaia orbital solution makes the source the most promising BH candidate among the ∼64 000 binary solutions in Gaia DR3 for which both time-resolved spectroscopic and astrometric data are available. Their work showed convincingly that the Gaia data suggest the presence of an unseen massive companion, which would be most simply explained by a BH. These authors did not obtain any follow-up data, so their conclusions depend entirely on the reliability of the Gaia solution.

3 PROPERTIES OF THE SOURCE

Gaia BH2 is a source in Centaurus (α = 13:50:16.7; δ = −59:14:20), near the Galactic plane (l = 310.4, b = 2.8), with G = 12.3. Its basic properties are summarized in Fig. 1. The field is relatively crowded, with 17 Gaia-detected neighbours within 10 arcsec, though all are at least 3.5 mag fainter than the source itself. The nearest Gaia-detected neighbour is a G = 17.95 source at a distance of 1.77 arcsec. This source is unresolved in the postage stamp shown in Fig. 1, which is from the DECaPS survey (Schlafly et al. 2018; Saydjari et al. 2023). Its astrometry is inconsistent with being bound to Gaia BH2.

The Gaia astrometric solution places Gaia BH2 at a distance d = 1.16 ± 0.02 kpc. We do not attempt to correct this distance for the Gaia parallax zero-point (Lindegren et al. 2021b), since there is little reason to expect that this zero-point – inferred from single-star astrometric solutions – applies to the binary solutions. On the colour–magnitude diagram, Gaia BH2 appears on the lower giant branch, just below the red clump. Our constraints on the star’s temperature and radius (Section 3.4) imply a mass near |$1\, {\rm M}_{\odot }$|⁠, as shown in the right panel of Fig. 1.

3.1 Light curves

We inspected the ASAS-SN V- and g-band light curves of the source (Kochanek et al. 2017), which contain more than 3000 photometric epochs over a 7-yr period, with a typical uncertainty of 0.02 mag. This did not reveal any significant periodic or long-term photometric variability: the optical source is constant at the ≳0.01 mag level.

Gaia BH2 was observed by TESS (Ricker et al. 2015) during sectors 11 and 38. The source is sufficiently bright (T ≈ 11.5) that it may be possible to measure asteroseismic parameters from the light curve (Stello et al. 2022). We extracted light curves from the TESS full frame images using eleanor (Feinstein et al. 2019) and analysed their power spectra following the methods described by Stello et al. (2022); see Appendix  C for details. We found a marginally significant detection of |$\nu _{\rm max}=61\, \mu \rm Hz$| in both sectors. If reliable, this value can be translated into a constraint on the giant’s surface gravity using asteroseismic scaling relations (e.g. Chaplin & Miglio 2013). Together with the spectroscopically measured effective temperature, it implies a surface gravity |$\log \left[g/\left({\rm cm\, s^{-2}}\right)\right]\approx 2.67$|⁠, consistent with the spectroscopic log g (Section 3.8). Given the marginal nature of the detection, we did not include it in our subsequent analysis.

3.2 Extinction

We estimate the extinction to Gaia BH2 using the 3D dust map from Lallement et al. (2022), which predicts a modest integrated extinction of |$A_{550\, \rm nm}=0.58\, \rm mag$| at d = 1.16 kpc, corresponding to E(B − V) ≈ 0.2 mag. The G-band extinction inferred by the GaiaGSP-Phot pipeline (Andrae et al. 2022) is A_G = 0.70, corresponding to E(B − V) ≈ 0.26, but we expect the value from the 3D dust map to be more reliable for stars in this evolutionary phase. We adopt E(B − V) = 0.2 ± 0.03 mag in our SED modelling.

3.3 Ultraviolet observation

Gaia BH2 does not have a published UV magnitude and is outside the published GALEX footprint. Measuring a UV magnitude or upper limit for the source is important for constraining possible flux contributions from a luminous secondary, so we obtained a 1625 s observation of the source with the UVOT photometer (Roming et al. 2005) on board the Neil Gehrels Swift Observatory in 2022 September (ToO ID 17935). We used the UVM2 band, with an effective wavelength of 2246 Å, and also obtained simultaneous X-ray observations using XRT (Burrows et al. 2005).

We analysed the UVM2 image using the uvotsource routine distributed with heasoft, using a 4 arcsec aperture centred on the Gaia source and an 8 arcsec background aperture in a nearby region with no obvious sources in the UVM2 image (RA = 207.558463, Dec. = −59.231955). Gaia BH2 is detected in the UVM2 image only with 2.7σ significance. A low-significance detection turns out to be sufficient for our purposes, since the main goal was to rule out a UV-bright companion. The AB magnitude in the UVM2 band is 22.20 ± 0.41 mag, corresponding to a flux density of |$(2.9 \pm 1.1)\times 10^{-17}\, \rm erg\, s^{-1}\, cm^{-2}$|⁠.

The source was not detected in the simultaneous XRT observation. We did not analyse the XRT data further, since our Chandra observations (Section 3.10) provide significantly deeper X-ray constraints.

3.4 SED model

We constructed the source’s broad-band SED by combining the Swift UVM2 magnitude, synthetic SDSS ugriz photometry constructed from Gaia BP/RP spectra (Gaia Collaboration 2022b), 2MASS JHK photometry (Skrutskie et al. 2006), and WISE |$W1\, W2\, W3$| photometry (Wright et al. 2010). We set an uncertainty floor of 0.02 mag in all bands to allow for photometric calibration issues and imperfect models. We then fit the SED with single-star models to infer the temperature and radius of the red giant. We predict bandpass-integrated magnitudes using empirically calibrated model spectral from the BaSeL library (v2.2; Lejeune, Cuisinier & Buser 1997, 1998). We assume a Cardelli, Clayton & Mathis (1989) extinction law with R_V = 3.1 and adopt a prior on the reddening E(B − V) = 0.20 ± 0.03. We use pystellibs.¹ to interpolate between model SEDs, and pyphot² 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.

The results are shown in the top panel of Fig. 2. A single-star model yields a reasonably good fit to the data, with χ²/N_data = 1.16, where N_data is the number of photometric points. The best-fitting SED model has |$T_{\rm eff}=4615\, \rm K$| and |$R = 7.77\, R_{\odot }$|⁠, with a bolometric luminosity of |$L\approx 25\, L_{\odot }$|⁠. Modelling of the source’s high-resolution spectrum (Section 3.8) leads to a very similar constraint on T_eff. The high precision of the Gaia parallax (ϖ/σ_ϖ = 47) allows us to constrain the star’s radius within |$\approx 3~{{\ \rm per\ cent}}$| and the luminosity within |$\approx 6~{{\ \rm per\ cent}}$|⁠.

The lack of significant UV excess leads to stringent limits on any possible main-sequence companions: if the companion were a normal star, it would be impossible to obtain a good fit to the SED with a single-star model (e.g. El-Badry & Rix 2022). Luminous stars that are massive enough to explain the observed astrometry and RVs are predicted to outshine the giant in the UV and contribute significantly to the optical spectrum. These scenarios are explored in the lower panels of Fig. 2 and discussed further in Section 4.

3.5 Luminous star mass and evolutionary state

We compare the measured radius and effective temperature of the luminous star to MIST single-star evolutionary models (Choi et al. 2016). We use solar-scaled models with Z = 0.012 to account approximately for the star’s measured sub-solar iron abundance and enhancement of α-elements (see Fu et al. 2018; Gaia BH2 has |$\rm [Fe/H]=-0.22$| and |$[\rm \alpha /Fe] = +0.26$|⁠, as described in Section 3.8). These models are shown in the right panel of Fig. 1.

Fitting the measured temperature and radius implies a luminous star mass of |$1.07\pm 0.19\, {\rm M}_{\odot }$|⁠. The uncertainty is dominated by uncertainty in the observed effective temperature, though uncertainties in the stellar models likely contribute at a similar level (e.g. Joyce et al. 2022). The corresponding age range is about 5–13 Gyr. In all plausible models, the star is on the lower red giant branch and has recently begun to expand following core hydrogen exhaustion; it is not a core helium burning red clump star that has already reached the tip of the giant branch and experienced a helium flash.

3.5.1 Could the star be a low-mass stripped giant?

We considered whether the giant might be a product of recent envelope stripping, in which case it might have significantly lower mass than implied by single-star evolutionary models (e.g. El-Badry et al. 2022; El-Badry & Rix 2022). We consider such a scenario unlikely given the binary’s wide and non-circular orbit, and the lack of evidence for ongoing mass transfer. In the current orbit, the star would overflow its Roche lobe at periastron only if it had a radius |$R\gtrsim 100\, R_{\odot }$|⁠. This essentially precludes a low-mass (⁠|$M\lesssim 0.5\, {\rm M}_{\odot }$|⁠) stripped giant, because a giant of that mass would not have reached a sufficiently large radius to overflow its Roche lobe in the Gaia BH2 system (e.g. Rappaport et al. 1995). A past period of mass transfer from the giant to the companion would also be expected to have circularized the orbit. Our analysis of the star’s spectrum reveals no excess absorption or emission in H α (Section 3.8), as has been commonly found in giants with ongoing or recently terminated mass transfer (Jayasinghe et al. 2021, 2022). Finally, the spectroscopic log g (Section 3.8) and the radius constraint from the SED fit imply a giant mass of |$M_\star = 1.17^{+0.88}_{-0.50}\, {\rm M}_{\odot }$|⁠, i.e. a 1σ lower limit of |$0.67\, {\rm M}_{\odot }$|⁠.

3.6 Spectroscopy

To measure atmospheric parameters, abundances, and RVs, we obtained multi-epoch spectra of Gaia BH2 using two instruments.

3.6.1 FEROS

We observed Gaia BH2 43 times with the Fiberfed Extended Range Optical Spectrograph (FEROS; Kaufer et al. 1999) on the 2.2m ESO/MPG telescope at La Silla Observatory (programs P109.A-9001 and P110.9014). The first several observations used 2 × 2 binning to reduce readout noise at the expense of spectral resolution; the remainder used 1 × 1 binning. The resulting spectra have resolution R ≈ 40 000 (2 × 2 binning) and R ≈ 50 000 (1 × 1 binning). Most of our observations used 1200 or 900 s exposures. The typical SNR at 5800 Å is 15–30.

We reduced the data using the CERES pipeline (Brahm, Jordán & Espinoza 2017), which performs bias-subtraction, flat fielding, wavelength calibration, and optimal extraction. The pipeline measures and corrects for small shifts in the wavelength solution during the course a night via simultaneous observations of a ThAr lamp with a second fibre. Several minor modifications to the pipeline were necessary to obtain good performance with 2 × 2 binning. We first calculate RVs by cross-correlating a template spectrum with each order individually and then report the mean RV across 15 orders with wavelengths between 4500 and 6700 Å. We calculate the uncertainty on this mean RV from the dispersion between orders, i.e. |$\sigma _{{\rm RV}}\approx {\rm std}\left({\rm RVs}\right)/\sqrt{15}$|⁠.

3.6.2 UVES

We observed Gaia BH2 five times with the Ultraviolet-Visible Echelle Spectrograph (UVES; Dekker et al. 2000) mounted on the 8.2m UT2 telescope at the VLT on Cerro Paranal (program 2110.D-5024). We used Dichroic 2 with central wavelengths of 437 and 760 nm on the blue and red arms, providing spectral coverage of the wavelength ranges from 3730–4990 and 5670–9450 Å. We took 600 s exposures and used 0.4 and 0.3 arcsec slits on the blue and red arms, respectively, which yielded resolutions of 80 000 and 110 000 with SNR ∼50 at 6500 Å.

We reduced the data using the ESO Reflex pipeline (Freudling et al. 2013) with standard calibrations. This performs bias-subtraction, flat fielding, wavelength calibration using afternoon ThAr arcs, optimal extraction, and order merging. We verified the stability of the wavelength solution using telluric absorption lines, measured RVs via cross-correlation of the merged red-arm spectra with a template, and adopted a per-epoch RV uncertainty of 0.1 |$\rm km\, s^{-1}$| to account for drifts in the wavelength solution between science exposures and afternoon arcs. Near-simultaneous FEROS and UVES observations allowed us to verify good agreement between the two RV scales at the |$\approx 0.1\, \rm km\, s^{-1}$| level.

3.6.3 Radial velocities

Our follow-up RVs are shown in Fig. 3 and listed in Table 1. In total, we obtained 48 RVs between 2022 August and 2023 March. All our observations have a precision better than 0.1 |$\rm km\, s^{-1}$|⁠. The precision of the FEROS RVs is somewhat higher than that of the UVES data because of the precise wavelength calibration enabled by simultaneous arcs.

3.7 Orbit modelling

We jointly fit our follow-up RVs and the astrometry + RV constraints from Gaia with a model that has 14 free parameters: the five standard astrometric parameters (⁠|$\alpha , \delta ,\mu _{\alpha }^{*},\mu _{\delta },$| and ϖ); and the binary parameters period, eccentricity, inclination, angle of the ascending node Ω, argument of periastron ω, periastron time, centre-of-mass velocity, luminous star mass M_⋆, and companion mass M₂. For each call to the likelihood function, we predict the semimajor axis in physical units, |$a=\left[P_{\rm orb}^{2}G\left(M_{\star }+M_{2}\right)/\left(4\pi ^{2}\right)\right]^{1/3}$|⁠. We then calculate the angular semimajor axis of the photocentre (a₀), and that of the primary (a₁, i.e. the star whose RVs are being measured):

Here, q = M₂/M_⋆ is the mass ratio and ϵ is the G-band flux ratio of the companion to the primary. d is the distance, which we compute as d[kpc] = 1/ϖ[mas]. For a dark companion (ϵ = 0), a₁ = a₀. We then predict the six Thiele–Innes parameters:

Note that in the convention used here and in the Gaia archive, the astrometric parameters A, B, F, and G have angular units (mas), while the spectroscopic parameters C and H have physical units (au). In the Gaia data processing, C and H are constrained by the measured RVs, while A, B, F, and G are constrained by the astrometry (e.g. Pourbaix et al. 2022).

For each call to the likelihood function, we construct the predicted vector of Gaia-constrained parameters:

We then calculate a likelihood term that quantifies the difference between these quantities and the Gaia constraints:

Here, |$\boldsymbol{\mu }_{{\rm Gaia}}$| and |$\boldsymbol{\Sigma }_{\rm Gaia}$| represent the vector of the best-fitting parameters constrained by Gaia and their covariance matrix, which we construct from the corr_vec parameter reported in the Gaia archive.

We additionally predict the RVs of the luminous star at the array of times t_i at which we obtained spectra. This leads to an RV term in the likelihood,

where RV_i and RV_pred(t_i) are the measured and predicted RVs, with their uncertainties σ_{RV, i}. The predicted RVs are calculated by solving Kepler’s equation using standard methods. The full likelihood is then given by

We use flat priors on all parameters except M_⋆, for which we use a normal distribution informed by isochrones, |$M_{\star }/{\rm M}_{\odot }\sim \mathcal {N}(1.07,0.19)$|⁠. We sample from the posterior using emcee (Foreman-Mackey et al. 2013) with 64 walkers, taking 3000 steps per walker after a burn-in period of 3000 steps. The results of the fit are summarized in Table 2 (third block).

We initially considered including the flux ratio ϵ as a free parameter. But in this case we found the posterior for ϵ to pile up against 0, with a 1σ upper limit of ϵ < 0.11. That is, a luminous companion is disfavoured by the astrometry because it would produce too small of a photocentre semimajor axis, independent of any arguments about the expected mass–luminosity relation of plausible luminous secondaries. In subsequent analysis, we therefore fixed ϵ = 0 in order to avoid introducing a bias due to unrealistically large flux contributions from the secondary, which we suspect is truly dark.

3.7.1 Gaia-only constraints

In order to explore the relative importance of our follow-up RVs and the Gaia data for our final constraints, we also repeat the fit without the follow-up RVs, i.e. simply removing the ln L_RVs term from equation (12). We list the resulting constraints in Table 2 (fourth block) and also show them with cyan contours in Fig 4.

Fig. 3 (top panel) compares the measured RVs to the predictions of the Gaia-only solution. There is significant (±3 months) uncertainty in the predicted RV turnover time, presumably because the Gaia data did not cover a full period and were obtained two orbital cycles before our follow-up (see Appendix  A). The follow-up RVs are fully consistent with the Gaia-only solution.

Perhaps counterintuitively, the best-fitting period obtained from the Gaia-only fit is P_orb = 1300 ± 26 d, which is different at the ∼1σ level from the period constraint in the Gaia archive, P_orb = 1352 ± 45 d. Several other parameters are also slightly discrepant. The reason for this difference is that our modelling explicitly links the spectroscopic and astrometric Thiele–Innes parameters via the requirement that a₀ = a₁, and by the requirement that a single pair of (ω, i) can explain all six Thiele–Innes parameters. In contrast, the Gaia data processing only requires that the spectroscopic and astrometric orbits have the same period, eccentricity, and eccentric anomaly at fixed time; A, B, F, G, C, and H are fit to the astrometric and RV time series without enforcing constraints between them (see Pourbaix et al. 2022). Inspection of equations (3)–(8) reveals that the six Thiele–Innes parameters are functions of four Campbell elements (a₀, ω, Ω, i), or of five if a₁ is fit as a separate parameter. This makes the problem overconstrained. In the presence of noise, there is thus no guarantee that a single set of orbital parameters can actually match all the Thiele–Innes parameters, even for a perfectly well-behaved orbit. Our approach minimizes the difference between the predicted and measured parameters. As we show in Table 3, there does indeed exist a single set of orbital parameters that can reproduce all the Gaia constraints and our follow-up RVs.

3.7.2 Joint constraints

The middle and bottom panels of Fig. 3 compare the measured RVs to the predictions of the joint Gaia + RV fit. With the RVs included, the fit includes data taken over the course of approximately three orbital periods, resulting in a tightly constrained orbit. The follow-up RVs cover most of the best-fitting orbit’s dynamic range in RV and thus constrain the RV mass function more tightly than the pure-Gaia constraints (Section 3.7.3). We show the residuals of the measured RVs with respect to the best-fitting orbit in the top panel of Fig. 5. The fit is good, with χ²/N_RVs ≈ 0.7, suggesting that the typical RV uncertainty is slightly overestimated.

In Fig. 4, we compare the constraints from joint modelling to those based only on the Gaia constraints. Reassuringly, the two sets of constraints are consistent, but the constraints that include our RVs are tighter.

3.7.3 Companion mass and inclination

Our joint fit constrains the mass of the unseen companion to |$M_2 = 8.9\pm 0.3\, {\rm M}_{\odot }$|⁠. The dominant source of uncertainty in this estimate is the uncertainty in the red giant mass.

The astrometric constraint on the eccentricity plays a critical role in constraining the mass, because the inclination, i ≈ 35 deg, is relatively face-on. This means that the component of the star’s orbital velocity projected on to our line of sight is only |$\sin (35\, \rm deg)\approx 57~{{\ \rm per\ cent}}$| of the total orbital velocity. We show how the implied companion mass varies with inclination in Fig. 5, where we plot the implied companion mass as a function of assumed inclination, giving the RV mass function, |$f\left(M_{2}\right)=P_{{\rm orb}}K_{\star }^{3}\left(1-e^{2}\right)^{3/2}/\left(2\pi G\right)\approx 1.33\, {\rm M}_{\odot }$|⁠. If the orbit were edge-on, the companion mass implied by the period, RV semi-amplitude, and eccentricity would be |$M_2 \approx 2.6\, {\rm M}_{\odot }$|⁠. While this is much lower than the fiducial mass constraint when the inclination measurement is taken into account, it is still too massive for the unseen companion to be any plausible luminous star, whether a single star or a binary itself (Section 4).

Only 18 per cent of randomly oriented orbits are expected to have an inclination as low as 35 deg. We note that in contrast to RV surveys, astrometric binary surveys have a selection bias in favour of low inclinations (e.g. Gaia Collaboration 2022a). The good agreement between the RV prediction from the Gaia solution and our measured RVs, as well as the fact that RVs alone rule out a stellar companion even for an edge-on orbit (Fig. 2), suggests that the Gaia inclination constraint is reliable.

3.7.4 Pure RV constraints

Recognizing that the Gaia orbital solution could be wrong, we now explore constraints that can be obtained from our follow-up RVs alone.

We fit the FEROS and UVES RVs with the Joker (Price-Whelan et al. 2017), which uses rejection sampling to obtain unbiased posterior samples consistent with a set of RVs. We used a p(P_orb) ∝ 1/P_orb prior between 100 and 2000 d, and broad uniform priors on other orbital parameters. Even with 10¹⁰ prior samples, the Joker returned only a few samples that were marginally consistent with the RVs, and all of these samples were concentrated in a single posterior mode. This reflects the fact that the RVs strongly constrain the orbit, such that only a tiny fraction of orbital parameter space is consistent with the data. To better sample the posterior, we initialized the standard MCMC chains near the maximum-likelihood sample and drew more samples using emcee. The results are reported in the bottom block of Table 2 and are shown in Fig. 6. In the top panel, the grey lines show RV predictions from random posterior draws; the red line shows the prediction of the best-fitting joint Gaia + RV solution. In the bottom panels, the contours show constraints from the RVs alone, while the red square marks the best-fitting joint solution.

The RV-only constraints are fully consistent with those from the joint fit. Because our RVs cover only a small fraction of the period, the orbit is not fully constrained by RVs alone. In particular, the RVs allow for orbits with arbitrarily long periods combined with increasingly large eccentricities. For this reason, there are many posterior samples with periods near our adopted upper limit of 2000 d. However, with a prior of P_orb < 2000 d, the orbit becomes well-constrained. Most importantly, the RVs alone provide a well-defined minimum period and mass function, ruling out any P_orb ≲ 1150 d and any RV mass function |$f(M_2)_{\rm RVs}\lesssim 1.25\, {\rm M}_{\odot }$|⁠. This means that – independent of the Gaia data – the RV orbit and luminous star mass imply a companion mass of at least |$2.5\, {\rm M}_{\odot }$|⁠.

3.8 Spectral analysis

We analysed the UVES spectra to measure atmospheric parameters and abundances of the red giant and place further constraints on the light contributions of the secondary, using a combination of ab-initio and empirical modelling. Most of our analysis focused on the spectrum obtained on JD 2459935, but we verified that there are no significant differences between the rest-frame UVES spectra taken on different nights.

We first fit the UVES spectra using the Brussels Automatic Code for Characterizing High accUracy Spectra (BACCHUS; Masseron, Merle & Hawkins 2016) with the same set up as in El-Badry et al. (2023). BACCHUS enables us to derive the stellar atmospheric parameters, including the effective temperature (T_eff), surface gravity (logg), metallicity ([Fe/H]), and microturbulent velocity (v_micro) assuming Fe excitation/ionization balance, i.e. the requirement that lines with different excitation potentials all imply the same abundances. We use the 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, and 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₂ (T. Masseron, private communication). Finally, we also include the SiH molecular line list from the Kurucz linelists.³ 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). Once 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 lies in the range −0.6 dex < [X/Fe] < +0.6 dex, and performs a χ² 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 one 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) and Hawkins, Lucey & Curtis (2020b). We ran bacchus on the UVES spectrum after merging individual de-blazed orders and performing a preliminary continuum normalization using a polynomial spline. Further normalization is performed by bacchus. Spectral regions affected by telluric absorption were masked during fitting, and we only fit the portion of the spectrum with |$\lambda \lt 7600\,$| Å, where tellurics are less severe.

The resulting stellar parameters and abundances are listed in Table 4. Fig. 7 compares portions of the observed spectrum to the best-fitting bacchus model. The fit is overall quite good: while there are some deviations between the model and data, the quality of the fit is typical for ab-initio (as opposed to data-driven) fits to high-resolution spectra. The remaining differences between data and model can be attributed to a combination of imperfect continuum normalization and imperfections in the line list.

The most important conclusion from this analysis is that the star’s iron abundance is subsolar (⁠|$[\rm Fe/H] = -0.22$|⁠), but the abundances of α-elements (i.e.; O, Mg, Si, Ca, and Ti) are enhanced relative to solar, with |$[\rm \alpha /Fe]=0.26\pm 0.05$|⁠. This abundance pattern is characteristic of old stars in the Galactic thick disc or metal-rich halo, which formed from material that was enriched by Type II SNe, with less enhancement of iron-peak elements by SNe Ia than typical stars in the thin disc. If Gaia BH2 formed in the thick disc, this would imply an age of 8–12 Gyr (e.g Nissen et al. 2020; Xiang & Rix 2022) and (assuming mass transfer was negligible) a stellar mass of |$M_\star \approx 1.0 \pm 0.07\, {\rm M}_{\odot }$| for the giant.

Beyond the strong α-enhancement, the abundances of the giant appear normal. We compared the star’s position in the [X/Fe] versus [Fe/H] plane to the population of stars in the solar neighbourhood with abundances measured by Adibekyan et al. (2012), Bensby, Feltzing & Oey (2014), and Battistini & Bensby (2016), finding that all elements fall within the |$\sim 2\, \sigma$| scatter of the observed local population. Lithium is not detected, as expected for a star on the giant branch that has undergone first dredge-up. There is no evidence of enhancement in neutron-capture elements relative to typical thick-disc stars.

3.8.1 Empirical comparison to other red giant spectra and constraints on luminous companions

To search for possible anomalous spectral features – which could arise, for example, due to a luminous companion, emission from an accretion disc, or absorption by an accretion stream – we compared the spectrum of Gaia BH2 to spectra of stars with similar stellar parameters and abundances observed by the GALAH survey (De Silva et al. 2015; Buder et al. 2021). To this end, we degraded the observed UVES spectrum to R = 28 000, shifted it to rest frame, and identified its nearest neighbour in pixel space among stars observed by GALAH with SNR > 100. The closest match was Gaia DR3 source 5909513401713377280, which in GALAH DR3 has |$T_{\rm eff} = 4591\pm 74\, \rm K$|⁠, log g = 2.45 ± 0.19, |$\rm [Fe/H]=-0.36\pm 0.05$|⁠, and |$[\rm \alpha /Fe]=0.28\pm 0.01$|⁠. These parameters are all similar to those we find for Gaia BH2, providing independent validation of our bacchus-inferred parameters. In Fig. 8, we compare the normalized and resolution-matched spectra of the two sources in the first and third GALAH wavelength windows, which, respectively, contain H β and H α. The spectra are nearly indistinguishable.

The similarity of the two observed spectra, and the good agreement between the data and the bacchus model spectrum in Fig. 7, speaks against the presence of a luminous secondary. The quantitative upper limit depends on the spectral type of the putative companion: a rapidly rotating secondary would be harder to detect. However, the depth of the observed absorption lines rules out even a companion that contributes only continuum. The total flux in the cores of the deepest lines is only |$\approx 2~{{\ \rm per\ cent}}$| of the continuum value (Fig. 7). A companion contributing more than |$\approx 2~{{\ \rm per\ cent}}$| of the total light at 4000–4300 Å would dilute the absorption lines beyond this value and is thus ruled out. We conservatively adopt a |$\lt 3~{{\ \rm per\ cent}}$| flux limit on the total flux; this is the origin of the blue upper limit in Fig. 2.

3.9 Galactic orbit

To investigate the Galactic orbit of Gaia BH2, we used its parallax and proper motion from the Gaia astrometric binary solution, as well as the centre-of-mass RV inferred from the joint 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 Fig. 9; for comparison, we also show the orbit of the Sun. The orbit is typical of a thin-disc star, with modest eccentricity and excursions above the disc mid-plane limited to ±80 pc.

3.9.1 Interpretation of the abundances and kinematics

The combination of thin-disc kinematics with α-rich chemistry is unusual. This is illustrated in Fig. 10, where we compare Gaia BH2’s chemistry and kinematics to the solar neighbourhood sample from Bensby et al. (2014). The left panel shows the Toomre diagram, while the right panel shows the abundance of titanium (a representative α-element) relative to iron. Most of the α-enhanced stars in the Bensby et al. (2014) sample have thick-disc kinematics. In contrast, Gaia BH2 is strongly enhanced in all α elements but has a Galactic orbit typical of the thin disc.

Given its thin-disc-like orbit, it is tempting to attribute the α-enhancement of Gaia BH2 to pollution from the BH progenitor. Indeed, similar α-enhancement in the donor stars of BH X-ray binaries has been interpreted as evidence for pollution during the supernova event. For example, 0.2–0.5 dex enhancement of some α-elements in the donor stars in GRO J1655-40 and V404 Cyg – the BH X-ray binaries whose detailed abundance patterns can be measured most readily – has been interpreted by Israelian et al. (1999) and González Hernández et al. (2011) as evidence that the BHs formed in a supernova. In those systems, the observed abundance patterns can plausibly be explained by pollution from a supernova with yields similar to those predicted for stars with initial masses of |$25\!-\!40\, {\rm M}_{\odot }$|⁠.

It is more challenging to apply a similar explanation to Gaia BH2. At the periastron distance of ≈2.4 au, a star with |$R_\star \approx 1\, R_{\odot }$| subtends only ≈10⁻⁶ of the sky as seen from the BH progenitor. The supernova ejecta are expected to consist mostly (⁠|$\approx 70~{{\ \rm per\ cent}}$|⁠) of oxygen (Thielemann, Nomoto & Hashimoto 1996) and must have had a total mass of |$\lesssim 10\, {\rm M}_{\odot }$| in order for the binary to remain bound. For a spherically symmetric explosion, we thus expect at most |$7\times 10^{-6}\, {\rm M}_{\odot }$| of oxygen to have been deposited on the companion during a supernova. During first dredge-up, this would have been mixed with the |$\approx 0.005\, {\rm M}_{\odot }$| of oxygen already present in the star’s convective envelope, and would not significantly increase the surface oxygen abundance observable today. Similar calculations apply to the predicted yields of other α-elements. It thus seems unlikely that the observed α-enhancement is a consequence of pollution by supernova ejecta, if these ejecta escaped the BH progenitor at high velocity and were deposited on the secondary directly.

It is possible, however, that some of the ejecta escaped the BH progenitor at low velocity, reached a separation of several au, and remained bound to the binary. In this case, a significant amount of the material may have eventually been accreted by the companion star, and it is possible that this supplied the |$\approx 0.005\, {\rm M}_{\odot }$| of α-elements required to explain the observed abundances. The expected total accretion on to the secondary in this scenario will depend on the energetics of the explosion and on the structure of the BH progenitor before its death; detailed calculations are required to explore this further.

3.10 X-ray observations

We observed Gaia BH2 for 20 ks using the Advanced CCD Imaging Spectrometer (ACIS; Garmire et al. 2003) on board the Chandra X-ray telescope on 2023 January 25 (proposal ID 23208881; PI: El-Badry). We used the ACIS-S configuration, with a spatial resolution of about 1 arcsec. The observations were timed to occur near the periastron passage, when the separation between the BH and the star was |$\approx 2.47\, \rm au$|⁠.

The Chandra data are shown in the upper left panel of Fig. 11. There is no obvious X-ray source coincident with the optical source, whose position is marked with a red circle. We used the chandra interactive analysis of observations (Fruscione et al. 2006) software to place upper limits on the X-ray flux. We first ran the chandra_repro tool to reprocess the observation; this creates a new bad pixel file and de-streaks the event file. We performed aperture photometry at the optical position of Gaia BH2 in the reprocessed event file using the srcflux tool. We detect no significant flux at the location of Gaia BH2 and obtain a background count rate of |$1.45 \times 10^{-4}\, \rm counts\, s^{-1}$|⁠. We assume a foreground hydrogen column density of N_H ≈ 1.36 × 10²¹ cm⁻² based on the optical extinction (e.g. Güver & Özel 2009) and use the Portable, Interactive Multi-Mission Simulator (PIMMS) tool to calculate the corresponding upper limit on the unabsorbed source flux. Based on the background count rate and N_H, and assuming a power-law spectrum with photon index 2, we obtain a 2σ flux limit in the 0.5–7 keV energy range of F_X < 5.1 × 10⁻¹⁵erg s⁻¹ cm⁻². This limit is shown as a horizontal dashed line in Fig. 11.

3.11 Radio observations

We observed Gaia BH2 for 4 h with the MeerKAT radio telescope in Lband (1.28 GHz) on 2023 January 13 (DDT-20230103-YC-01, PI: Cendes), when the separation between the BH and the star was |$\approx 2.54\, \rm au$|⁠. We used the flux calibrator J1939-6342 and the gain calibrator J1424-4913, and used the calibrated images obtained via the SARAO Science Data Processor (SDP)⁴ for our analysis.

The MeerKAT data are shown in the bottom left panel of Fig. 11. There is no detectable radio source coincident with Gaia BH2. We measured the flux density using the imtool package within pwkit (Williams et al. 2017) at the location of Gaia BH2. The RMS at the source’s position is 17 |$\mu$|Jy, so we report a 2σ (3σ) non-detection of <34 |$\mu$|Jy (<51 |$\mu$|Jy). The 2σ limit is shown with a horizontal dashed line in Fig. 11.

3.12 Expected X-ray and radio flux

3.12.1 Bondi–Hoyle–Lyttleton accretion

A rough estimate of the expected accretion rate on to the BH can be obtained under the assumption that a spherically symmetric wind from the giant is accreted at the Bondi–Hoyle–Lyttleton (BHL) rate:

The giant’s mass-loss rate, |$\dot{M}_{\rm wind}$|⁠, is uncertain. We consider the prescription from Reimers (1975):

where β_R is a dimensionless constant; we adopt β_R = 0.1 following Choi et al. (2016). For the wind velocity, a reasonable approximation is that it scales with the escape velocity at the stellar surface,

where β_wind is another dimensionless constant.

We then assume that the X-ray flux is emitted isotropically with a radiative efficiency η_X, such that |$F_{\rm X,BHL}=\frac{\eta _{X}\dot{M}_{{\rm BHL}}c^{2}}{4\pi d^{2}}$|⁠. This leads to a predicted X-ray flux for accretion at the BHL rate:

where we have scaled to parameters similar to Gaia BH2. Given the many approximations that enter this scaling, it is only expected to hold to order of magnitude. The most important prediction of equation (16) is the strong scaling with the radius and luminosity of the luminous star, |$F_{\rm X,BHL} \propto R_\star ^3 L_\star$|⁠. This means that at fixed mass, a giant leads to a much larger predicted X-ray flux than a dwarf. Physically, this occurs because (a) the giant has a significantly higher mass-loss rate, |$\dot{M}_{\rm wind}$|⁠, and (b) the wind from a giant is slower, allowing the BH to capture a larger fraction of it. As a result, the predicted X-ray flux from Gaia BH2 is ≈1000 times larger than that from Gaia BH1, even though Gaia BH2 is farther away and has a wider orbit.

We can also predict the radio flux assuming the source falls on the canonical X-ray–radio flux correlation. We take the version of this relation from Gallo et al. (2006):

where L_R and L_X represent the radio and X-ray luminosity. Following Gallo et al. (2006), we translate this into a radio flux density under the assumption of a flat radio spectrum (S_ν ∝ ν⁰) up to a maximum frequency |$\nu _{\rm max} = 8.4\, \rm GHz$|⁠. We show the expected X-ray and radio fluxes from Gaia BH2 in Fig. 11 with red-dashed lines. For each value of |$\dot{M}_{\rm BHL}$|⁠, we predict the radiative efficiency using the fitting function from Xie & Yuan (2012), who used simulations of hot accretion flows to calculate these efficiencies as a function of the accretion rate at the event horizon in Eddington units, for several values of the electron heating parameter, δ. We consider their predictions for δ = 0.5. For a plausible wind speed of |$v_{\rm wind}\approx 150\, \rm km\, s^{-1}$| (somewhat lower than the giant’s escape velocity, since the wind must escape the giant’s gravitational potential before reaching the BH), the predicted X-ray and radio fluxes are in tension with the observed non-detections.

3.12.2 Reduction in |$\dot{M}$| near the BH

The calculations above implicitly assume that the accretion rate near the event horizon is equal to the BHL rate. This may not be realized in practice if, for example, a significant fraction of the inflowing material is lost to winds (e.g. Blandford & Begelman 1999). Numerous simulations of radiatively inefficient accretion flows predict that |$\dot{M}$| decreases towards the event horizon roughly as |$\dot{M} \propto r^{0.5}$|⁠, where r is the distance from the BH (Pen, Matzner & Wong 2003; Pang et al. 2011; Yuan, Wu & Bu 2012; Ressler et al. 2021). There are several proposed physical explanations for this scaling. For concreteness, we consider the model from Yuan et al. (2012), in which |$\dot{M} \propto r^s$|⁠, where s ≈ 0.5 down to r ∼ 10R_s, and s ≈ 0 at smaller radii, where R_s is the Schwarzschild radius.

The BHL accretion rate corresponds to the accretion rate on scales of the BH’s accretion radius, |$R_{\rm acc} = GM_{\rm BH}/v_{\rm wind}^2 \approx 0.3\, \rm au$|⁠. The model described above thus predicts that the accretion rate at the event horizon (and all radii interior to |$\sim 10\, R_s$|⁠) is

i.e. a factor of ∼500 reduction in |$\dot{M}$| relative to the BHL rate!

A lower |$\dot{M}$| also leads to a lower predicted radiative efficiency. The expected accretion rate at the horizon is |$\dot{M} \approx 4\times 10^{-16} \approx 2\times 10^{-9}\dot{M}_{\rm edd}$|⁠, and the corresponding radiative efficiency from Xie & Yuan (2012) is predicted to be η_X ≈ 7.5 × 10⁻⁵, as opposed to η_X ≈ 0.004 for the BHL case. Substituting this value into equation (16) and finally multiplying by 0.002 (equation 19), we obtain a predicted |$F_X\approx 10^{-17}\, \rm erg\, s^{-1}\, cm^{-2}$| for Gaia BH2, which is more than two orders of magnitude below the Chandra limit. The corresponding radio flux is ≈0.2 |$\mu$|Jy, more than two orders of magnitude below the MeerKAT limit. The predicted X-ray and radio fluxes are shown with solid black lines in Fig. 11.

To summarize: if we take the accretion rate on to the BH to be the BHL rate predicted by equation (13), then the predicted radiative efficiency is η_X = 0.004, and the X-ray flux predicted by equation (16) is |$\approx 4\times 10^{-13}\, \rm erg\, s^{-1}\, cm^{-2}$|⁠, which would have been easily detected. On the other hand, if we reduce the predicted accretion rate to account for a reduction in |$\dot{M}$| between the accretion radius and the horizon (equation 18), the predicted accretion rate and radiative efficiency both fall significantly, and the expected X-ray and radio fluxes would not be detectable. The non-detection of Gaia BH2 thus supports models in which winds, convection, or other processes cause a significant reduction in the horizon accretion rate relative to the BHL rate. Such models have also been supported by polarization measurements of emission from Sgr A* and M87*, which show that the accretion rate near the event horizon is two to three orders of magnitude less than the Bondi value (Quataert & Gruzinov 2000; Kuo et al. 2014).

4 NATURE OF THE COMPANION

Since we have not detected the companion directly, we can infer its nature only from (a) dynamical constraints on its mass and (b) the fact that it emits little if any detectable light. The combination of Gaia astrometry and our RVs yields a robust and precise constraint on its mass: |$8.9\pm 0.3\, {\rm M}_{\odot }$|⁠. This relies on astrometric constraints on the orbital inclination, which are not accessible with RVs alone.

It is instructive to disregard the astrometric inclination constraint and consider only the RVs. The orbital solution from our joint fit has a spectroscopic mass function, |$f\left(M_{2}\right)=P_{{\rm orb}}K_{\star }^{3}\left(1-e^{2}\right)^{3/2}/\left(2\pi G\right)=1.33\pm 0.01\, {\rm M}_{\odot }$|⁠. This represents the absolute minimum companion mass compatible with the RV orbit, in the limit of an edge-on orbit and a luminous star with zero mass. Adopting |$M_{\star } = 1.07\pm 0.19\, {\rm M}_{\odot }$| and assuming an edge-on orbit (in conflict with the astrometric solution), this implies an minimum companion mass of |$2.64\pm 0.16\, {\rm M}_{\odot }$| (Fig. 5).

We now consider what kinds of companions could have the dynamically implied mass. As demonstrated in Fig. 2, a |$9\, {\rm M}_{\odot }$| main-sequence star would outshine the giant at all wavelengths. It would have a bolometric luminosity of |$L_2 \gtrsim 5,000\, L_{\odot }$|⁠, at least 200 times brighter than the giant. An evolved |$9\, {\rm M}_{\odot }$| star would be even more luminous, so any |$9\, {\rm M}_{\odot }$| star is ruled out by the non-detection of a second component in the UVES spectra and by the faintness of the source in the Swift UVM2 data. A |$9\, {\rm M}_{\odot }$| companion is also too massive to be a white dwarf or neutron star. While a small cluster of white dwarfs or neutron stars could escape detection, such a configuration would be difficult to assemble and would rapidly become unstable. We conclude that if the Gaia astrometric solution is correct, the companion must be a BH, or a close binary containing a BH.

What if the astrometric constraint on the inclination is spurious? In that case – if the period and eccentricity from the Gaia solution are still reliable – the companion mass could be as low as |$\approx 2.6\, {\rm M}_{\odot }$|⁠. This is still too massive to be a white dwarf or a neutron star, and a |$2.6\, {\rm M}_{\odot }$| main-sequence companion is inconsistent with the SED and spectra (third panel of Fig. 2). We also consider the possibility that the companion is an unresolved binary consisting of two stars with total mass |$2.6\, {\rm M}_{\odot }$|⁠. An equal-mass inner binary will yield the highest mass–luminosity ratio, and the bottom panel of Fig. 2 shows that this scenario is also ruled out by the faintness of the source in the UV and the spectral limit in the blue optical.

There is also an evolutionary problem for scenarios in which the companion contains one or more main-sequence stars: in all cases, the stars would be more massive than the giant itself, and should thus be more evolved than it is. But one or more red giant companions would outshine the giant in the optical.

4.1 Could the Gaia orbit be spurious?

Since our RV follow-up has covered only a limited fraction of Gaia BH2’s orbit, it is in principle possible that the Gaia orbit is simply wrong. In practice, our follow-up RVs alone constrain the orbit fairly precisely (Fig. 6), excluding any orbit with |$P_{\rm orb} \lesssim 1150$| d or |$f(M_2)_{\rm RVs} < 1.25\,M_{\odot}$|⁠. While we do not have access to the epoch-level astrometric data on which the Gaia solution is based, it it is very unlikely that the observed RVs would just happen to match the predictions of the Gaia solution if that solution were incorrect. The Gaia orbital solution appears unproblematic based on its goodness of fit metrics (Appendix  B), and the source was observed by Gaia almost 90 times over a 1000-d period, sampling most of the orbital phase (Appendix  A). Gaia BH2 has a AstroSpectroSB1 solution. In order to obtain an AstroSpectroSB1 solution in Gaia DR3, a source must first have independent spectroscopic and astrometric solutions, which can be combined in a manner that improves the goodness of fit relative to the purely astrometric solution (Pourbaix et al. 2022). This means that if the solution were spurious, the spectroscopic and astrometric orbits would independently have to be spurious in compatible ways – a scenario that is unlikely to occur for a source with such good astrometric phase overage. And finally, the minimum companion mass implied by the RVs alone – |$2.5\, {\rm M}_{\odot }$| – is still too high to be a white dwarf, or a main-sequence star that could escape spectroscopic and photometric detection.

More broadly, the fact that the measured RVs are in good agreement with the predictions of the Gaia solution, and the Gaia fit quality metrics are unremarkable, leaves us with little reason to doubt the reliability of the Gaia solution, including the inclination constraint. In this case, the inclination is robustly constrained to 34.9 ± 0.4 deg, and the companion mass to |$8.9\pm 0.3\, {\rm M}_{\odot }$|⁠. A companion of this mass must contain at least one BH. A single BH with |$M\sim 8.9\, {\rm M}_{\odot }$| thus seems to be the simplest explanation, although the data do not rule out an inner binary containing two BHs or a BH and another compact object or low-mass star.

4.2 Comparison to other known BHs

In Figs 12, 13, and 14, we compare Gaia BH2 and BH1 to other known BHs. The red and blue points show low- and high-mass X-ray binaries, whose parameters we take from Remillard & McClintock (2006) and the BlackCAT catalogue of X-ray transients introduced by Corral-Santana et al. (2016). We also show the binaries VFTS 243 (in the LMC; Shenar et al. 2022) and HD 130 298 (in the Milky Way; Mahy et al. 2022), both of which are single-lined binaries containing |$\sim 25\, {\rm M}_{\odot }$| O stars and unseen companions suspected to be BHs. In magenta, we show two binaries with suspected BH companions in the globular cluster NGC 3201, discovered with MUSE (Giesers et al. 2018, 2019).

For the BHs in the Milky Way – which includes all the LMXBs but only two of the HMXBs – we additionally collected distance estimates and quiescent G-band magnitudes. Five systems – A0620-00, Cyg X-1, V4641 Sgr, GRO J1655-40, and HD 130298 – have reasonably precise Gaia DR3 parallaxes, with pallax_over_error > 5. For these sources, we use the geometric distance estimates from Bailer-Jones et al. (2021), which are informed by the parallax and a Galactic model prior. For the other sources, which are fainter and/or more distant, the parallaxes are not very constraining, and so we collect distance estimates from the literature that are based on the orbital period, donor spectral type, and apparent magnitude in quiescence. Where available, we take these estimates from Jonker & Nelemans (2004). For the objects not included in that work, we take the estimates from individual-object papers: Casares et al. (2009) for GS 1354–64, Heida et al. (2017) for GX 339–4, Mikołajewska et al. (2022) for MAXI J1820+070, and Mata Sánchez et al. (2021) for MAXI J1305−704. A few sources are too faint in the optical (in some cases due extinction) to be detected by Gaia; these are not shown in Fig. 14.

Gaia BH2 has the longest orbital period of all currently known and suspected BHs. Its orbital period is a factor of 7 longer than that of Gaia BH1, and ≈40 times longer than that of the longest period X-ray binary, GRS 1915+105. The masses of Gaia BH1 and BH2 are similar to the bulk of the BH X-ray binary population; only GRS 1009–45 and GRO J0422+32 host BHs of significantly lower mass.

Both Gaia BH1 and BH2 are closer, and have longer orbital periods, than any other BH binaries. They are brighter than any of the BH X-ray binaries with low-mass donors, but fainter than the HMXB Cyg X-1 (which hosts a |$\sim 400,000\, L_\odot$| luminous star). Given that their companion stars are cool and slowly rotating, their abundances and RVs can be measured with higher fidelity than in other BH companions. The two systems’ close distances imply that widely separated BH binaries significantly outnumber close and accreting BH X-ray binaries.

5 FORMATION HISTORY

As with Gaia BH1 (El-Badry et al. 2023), it is unclear whether Gaia BH2 formed from a primordial binary or via a more exotic channel. We first consider formation from a primordial binary and then discuss alternative channels in Section 5.5.

5.1 Natal kick constraints

The formation of a compact object after core collapse is expected to impart a natal kick to the compact object, modifying the orbital size and eccentricity. Given the wide current orbit of Gaia BH2, it is unlikely that tides have appreciably circularized the binary following the formation of the BH. Assuming that the system formed from a pre-existing binary, its semimajor axis and eccentricity may have been changed by an asymmetric supernova kick and by instantaneous mass-loss (e.g. Blaauw 1961). We used cosmic (Breivik et al. 2020) and emcee (Foreman-Mackey et al. 2013) to model the combined binary properties at core collapse and natal kick strength based on the present-day orbital period and eccentricity of Gaia BH2 using the 1σ measurement errors. We followed a similar procedure to the one described in Wong et al. (2022) and El-Badry et al. (2023). But instead of modelling the primordial binary properties, we initialized our model with the BH progenitor mass |$12\, {\rm M}_{\odot }$| just before core collapse, secondary mass |$1.07\, {\rm M}_{\odot }$|⁠, and zero eccentricity. This effectively marginalizes over all previous evolution up to the point of core collapse of Gaia BH2’s progenitor. Our choice of BH progenitor mass is defined by the Fryer et al. (2012) ‘delayed’ model and produces a BH with mass |$9\, {\rm M}_{\odot }$|⁠. The Blaauw kick associated with the instantaneous mass-loss imparts an eccentricity of ∼0.3, so a moderate natal kick is required to produce the observed eccentricity, assuming that the orbit was previously circularized by tides.

The strength of the natal kick is degenerate with the orbital separation at core collapse, a_cc, which we therefore include as a parameter in our model. We used uniform priors for the orbital separation at core-collapse between 1 and |$5\, \rm {au}$|⁠, the natal kick velocity between 0 and |$100\, \rm km\, s^{-1}$|⁠, and natal kick angles that are uniform on a sphere. We ran 1024 walkers for 50 000 steps, thinned the chains by a factor of 10, and retained only the last 2000 steps to ensure a proper burn in. We find that the present-day orbital period and eccentricity prefer a core-collapse semimajor axis |$a_{\rm {cc}} = 2.9\pm 0.93\, \rm {au}$| and natal kick speed |$v_{\rm {kick}} = 36^{+21}_{-11}\, \rm km\, s^{-1}$|⁠. Assuming Gaia BH2 formed as a primordial binary, this suggests that the semimajor axis was likely smaller in the past than its present-day value. The exact constraints depend on the assumed pre-supernova mass of the progenitor, which is uncertain. However, a generic result of this modelling is that the observed wide orbit and moderate eccentricity disfavour strong natal kicks to the BH.

5.2 BH progenitor mass

Given the BH’s mass of |$\sim 9\, {\rm M}_{\odot }$|⁠, its progenitor likely had an initial mass |$M_{\rm init}\gtrsim 25\, {\rm M}_{\odot }$| (e.g. Sukhbold et al. 2016; Raithel, Sukhbold & Özel 2018). Because more massive stars are predicted to have stronger winds, the relation between initial mass and BH mass for stars of near-solar metallicity is predicted to be non-monotonic, and – given current observational constraints – essentially any initial mass |$M_{\rm init}\gtrsim 25\, {\rm M}_{\odot }$| could plausibly have produced a |$9\, {\rm M}_{\odot }$| BH.

5.3 Problems with formation through common envelope evolution

If Gaia BH2 formed as a primordial binary with a |$\gtrsim 25\, {\rm M}_{\odot }$| primary and a ∼1M_⊙ secondary, mass transfer through Roche lobe overflow of the primary would almost certainly lead to an episode of common envelope evolution due to the extreme mass ratio. This is predicted to reduce the orbital separation by a factor of 100–1000, depending on the structure of the primary when mass transfer began (see the discussion in El-Badry et al. 2023). Since Roche lobe overflow could not have begun at a separation significantly larger than 10 au, this leads to a predicted post-common envelope separation at least a factor of 50 closer than observed. Given the system’s low space velocity and modest eccentricity, we also cannot invoke a fine-tuned natal kick to the BH to widen the orbit again after a common-envelope episode. The same considerations also disfavoured a common envelope for Gaia BH1, and Gaia BH2’s wider orbit exacerbates the tension.

5.4 Formation from a star too massive to become a red supergiant

We next consider whether Gaia BH2 could have formed from a massive star that never became a red supergiant and thus avoided mass transfer altogether. We modelled the evolution of the progenitor with single massive star evolutionary models of solar metallicity, calculated with mesa by Klencki et al. (2020). The orbital separation at the time of the progenitor’s death is unknown, since a natal kick of unknown magnitude and direction could have modified the orbit in a variety of ways. But as discussed above, the modest eccentricity suggests that the final pre-supernova separation was likely not too different from the current semimajor axis (∼5 au). We adopt a value of 3 au.

The massive star models undergo significant mass-loss through winds. Assuming this occurs isotropically on a time-scale that is long compared to the orbital period, the orbit will expand such that the product of the semimajor axis and total binary mass is conserved. This means that the orbit must have been tighter in the past in the absence of mass transfer. For each massive star model, we assume a final separation of 3 au, and then calculate the separation at each time-step based on the mass at that time. We also calculate the Roche lobe radius, using the fitting formula from Eggleton (1983) and assuming a |$1\, {\rm M}_{\odot }$| companion. Finally, we calculate the ratio of the massive star model’s radius to the Roche lobe. We plot the maximum value of this quantity across all time-steps in Fig. 15. If this ratio is larger than 1, it means that the BH progenitor would have overflowed its Roche lobe at some point if it were in the Gaia BH2 system, leading to a common envelope inspiral.

All models with |$M_{\rm init} \le 65\, {\rm M}_{\odot }$| would overflow their Roche lobes in the Gaia BH2 system by more than a factor of 2, meaning that they could not have avoided a common envelope. However, the most massive models behave quite differently: these models lose their hydrogen envelopes to strong winds while still on the main sequence and never become red giants (see e.g. Smith & Conti 2008). If the progenitor of Gaia BH2 followed such an evolutionary channel, the system could have formed as, e.g. a |$70+1\, {\rm M}_{\odot }$| primordial binary with separation |$\sim 0.5\, \rm au$|⁠, completely avoided mass transfer through Roche lobe overflow, and died as a Wolf–Rayet + |$1\, {\rm M}_{\odot }$| binary with separation ∼3 au. The mass threshold above which stars avoid a red supergiant phase varies significantly between models, and is as low as |$\sim 40\, {\rm M}_{\odot }$| in some models (e.g. Agrawal et al. 2020; Bavera et al. 2022).

An important caveat to these calculation is that the maximum radius of massive stars depends sensitively on parameters that are uncertain, particularly their wind mass-loss rates. Some recent studies (e.g. Beasor et al. 2020) suggest that winds of red supergiants with |$M\gtrsim 25\, {\rm M}_{\odot }$| are weaker than predicted by classical models. This may imply that the mass-loss rates of more massive stars are also lower, in which case envelope loss while on the main sequence might not occur. The mass threshold above which models avoid a red supergiant phase also depends significantly on assumptions about rotation (Maeder 1987). We conclude that there are significant uncertainties in the evolutionary models for massive stars that avoid expansion, but it is possible that Gaia BH2 formed through such a channel.

5.5 Other formation channels

It is also possible that Gaia BH2 did not form from a primordial binary. One possibility is that the system formed dynamically through an exchange interaction in a dense cluster. The binary’s thin disc-like Galactic orbit precludes formation in a classic globular cluster, but the system may have formed in an open cluster that has since dissolved. This scenario does not leave any smoking gun signatures and is thus difficult to test definitively for any single binary. N-body simulations of binaries in clusters (e.g. Kremer et al. 2018) that are tailored to systems like Gaia BH1 and BH2 are required to determine whether such dynamical exchange+ejection events in clusters are common enough to explain the observed rate of wide BH + normal star binaries. If the observed α-enhancement of the giant is a result of pollution from the BH progenitor (see Section 3.9.1), this would disfavour a dynamical formation channel.

It is also possible that Gaia BH2 formed as a triple, with the red giant an outer tertiary to an inner binary containing two massive stars. Interactions between the two stars could have prevented either one from becoming a red supergiant. This scenario is discussed further by El-Badry et al. (2023) for Gaia BH1; a similar scenario could also work for Gaia BH2. Indeed, a triple scenario likely requires less fine-tuning for Gaia BH2 than for BH1, because the wider orbit makes it possible to accommodate a wider range of inner binaries without the system becoming unstable. As with Gaia BH1, precision RV monitoring of the luminous star may eventually detect or rule out a dark inner binary, which would cause a short-period, low-amplitude modulation of the overall RV curve.

5.6 Future evolution

Over the next |$\approx 100\, \rm Myr$|⁠, the luminous star is expected to continue ascending the red giant branch, until it ignites core helium burning. The MIST models shown in Fig. 1 predict that if it evolved in isolation, the giant would reach a radius of about |$170\, R_{\odot }$| just before the helium flash. However, in the current binary orbit, the giant would overflow its Roche lobe at periastron when it reached |$R\approx 110\, \rm R_{\odot }$|⁠. The orbit will likely circularize as the giant expands, and there will likely be a brief period of mass transfer through stable Roche lobe overflow near the tip of the giant branch, when the system will be observable as a symbiotic X-ray binary. After helium ignition, the giant will contract and the binary will appear as a wide BH binary with a red clump or hot subdwarf star. Another episode of stable Roche lobe overflow is possible after core helium burning, depending on how much of the star’s envelope mass remains. Either way, the binary is expected to terminate its evolution as a wide BH + white dwarf binary, probably with a low-eccentricity orbit.

6 AN EARLY ASSESSMENT OF THE BH POPULATION REVEALED BY GAIA DR3

6.1 Completeness of follow-up

We identified both Gaia BH1 and BH2 as BH candidates using simple astrometric cuts designed to find massive and dark companions to low-mass stars (see appendix E of El-Badry et al. 2023). We have now followed-up all of the six candidates those cuts yielded. Two have been validated, and the remaining four turned out to have spurious Gaia orbital solutions.

A broader search for compact object companions among the Gaia DR3 astrometric solutions was carried out by Shahaf et al. (2023). Gaia BH2 was excluded from their sample both because the luminous source is not on the main sequence and because its orbital period is longer than 1000 d. Besides Gaia BH1, Shahaf et al. (2023) identified an additional seven BH candidates, whose Gaia orbital solutions appeared to imply a dark companion with |$M_2 \gt 2\, {\rm M}_{\odot }$|⁠. We are carrying out RV follow-up of these sources, which we will describe in detail elsewhere. Only for one of them, Gaia DR3 6328149636482597888, do the RVs inspire confidence in the Gaia solution. However, that source is a low-metallicity halo star, and so the masses of both the luminous star and the unseen companion are lower than inferred when the system is interpreted with solar-metallicity models (see El-Badry et al. 2023). The best-fitting companion mass is |$M_2 \approx 2.2\, {\rm M}_{\odot }$|⁠. This is consistent with a variety of astrophysical scenarios; we defer more detailed analysis of this source to future work.

In summary, we suspect that Gaia BH1 and BH2 are the only BHs with DR3 astrometric orbital solutions, at least among sources with implied companion masses |$M_2\gtrsim 2.5\, {\rm M}_{\odot }$|⁠. There are, however, a few dozen candidates with |$1 \lesssim M_2 \lesssim 2\, {\rm M}_{\odot }$| – most readily interpreted as neutron stars and white dwarfs – still to be followed up.

6.2 Is there a BH – neutron star mass gap?

It is conspicuous that there are now two vetted Gaia compact object binary solutions with 8 < M₂/M_⊙ < 10, a few dozen candidates with |$1 \lesssim M_2 \lesssim 2.5\, {\rm M}_{\odot }$|⁠, and no reliable solutions with 3 < M₂/M_⊙ < 8. The implied BH/NS companion mass distribution is reminiscent of the mass gap observed for BHs in X-ray binaries (e.g. Kreidberg et al. 2012), though the mean BH mass inferred for X-ray binaries is somewhat lower (⁠|$\sim 7\, {\rm M}_{\odot }$|⁠). While Gaia is somewhat more sensitive to high-mass companions, which produce a larger photocentre wobble at fixed orbital period, both Gaia BH1 and BH2 would plausibly have been detectable if the BH had a mass of, say, |$5\, {\rm M}_{\odot }$|⁠. The data are thus already in tension with a scenario in which low-mass BHs significantly outnumber |$\sim 10\, {\rm M}_{\odot }$| BHs, at least among BHs with stellar companions in wide orbits. The fact that Shahaf et al. (2023) identified several dozen wide neutron star candidates with dynamically implied masses of |$1-2\, {\rm M}_{\odot }$| also suggests that Gaia would be sensitive to lower mass BHs, but further follow-up is required to vet the neutron star candidates.

The next decade will likely see the discovery of both more wide BH + normal star astrometric binaries and isolated single BHs discovered via astrometric microlensing (e.g. Lam et al. 2022; Mroz, Udalski & Gould 2022; Sahu et al. 2022). Comparison of the mass distributions of BHs within these populations may shed light on the formation pathways of the BH binaries.

6.3 Orbital period distribution

In Fig. 16, we compare the orbital periods of Gaia BH2 and BH1 to the period distribution of all ≈186 000 astrometric binaries published in DR3. Most of these systems are solar-type binaries, whose intrinsic period distribution is well approximated by a lognormal that peaks at a period of 10⁵ d with a dispersion of 2.3 dex (Raghavan et al. 2010). To estimate the relative sensitivity of Gaia DR3 to binaries at different periods, we divide the observed period distribution by this lognormal. The latter decreases towards longer periods (for linear period intervals), so this rescaling gives more weight to binaries with long periods.

Fig.16 shows that Gaia DR3 was most sensitive to binaries with orbital periods between 1.25 and 2.5 yr, since shorter periods produce smaller astrometric wobbles at fixed mass, and longer periods could not be fully sampled during the ∼1000 d observing baseline for data published in DR3. The strongly reduced sensitivity to periods of 1 yr, due to confusion of orbital and parallax motion, is also evident. The fact that the two BHs discovered from Gaia DR3 have orbital periods on the long- and short-period edges of the observational sensitivity curve hints that the BH companion period distribution may be bimodal, with a valley at 1–3 yr periods. The strength of the conclusions we can draw is of course limited by the fact that we have detected only two BHs. We can nevertheless draw some conclusions. If the intrinsic period distribution of BH binaries were flat in linear orbital period, 80 per cent of the detected binaries would have periods between 1 and 3 yr. There is thus only a |$0.2^2 = 4~{{\ \rm per\ cent}}$| chance of detecting two BHs outside this range and none inside it if the underlying BH period distribution is flat.

As discussed by El-Badry et al. (2023), Gaia DR4 and DR5 will likely enable the discovery of dozens of BH + normal star wide binaries, both because they will be based on a longer observing baseline than DR3 and because epoch-level astrometric data will be published for all sources, without stringent cuts on signal-to-noise ratio. DR4 will be based on ∼2000 d of observations, and DR5 will likely be based on ∼4000 d of observations. Gaia will thus eventually enable the discovery and characterization of BH binaries with periods up to ∼15 yr (and perhaps longer if partial orbits can be fit uniquely; e.g. Andrews et al. 2021). The intrinsic period distribution of BH + normal star binaries is quite uncertain, but some models predict that it will rise steeply at P_orb ≳ 10 yr, where the BH binaries will have formed without the BH progenitor ever having overflowed its Roche lobe (e.g. Breivik et al. 2017; Chawla et al. 2022). This prediction depends sensitively on the strength of BH natal kicks, which can easily unbind the widest binaries. In any case, the empirical period and mass distributions of BH binaries will become much better sampled in future Gaia releases.

7 CONCLUSIONS

We have carried out follow-up observations of a newly identified BH + luminous star binary discovered via Gaia astrometry and spectroscopy. The system, which we call Gaia BH2, consists of an apparently unremarkable lower red giant branch star in a 3.5-yr orbit with an unseen massive companion. Basic properties of the system are summarized in Table 2. Our main conclusions are as follows:

• The luminous source: the luminous source is a bright (G = 12.3), nearby (⁠|$d=1.16\, \rm kpc$|⁠) star on the lower red giant branch, just below the red clump (Fig. 1). From evolutionary models, we infer a mass |$M_\star = 1.07\pm 0.19\, {\rm M}_{\odot }$|⁠. From a high-resolution spectrum (Fig. 7), we measured moderately subsolar metallicity (⁠|$\rm [Fe/H] = -0.22$|⁠) but strong enhancement of α-elements (⁠|$\rm [\alpha /Fe] = +0.26$|⁠). The detailed abundance pattern is otherwise unremarkable (Table 4). Because the system’s orbit is wide, the α-enhancement is most likely primordial, but it is also possible that it is a result of pollution from low-velocity ejecta during a failed supernova.

• Orbit and companion mass: our RV follow-up over a 7-month period validates the Gaia orbital solution (Fig. 3) and yields improved constraints on the companion mass (Fig. 4). The orbit has a 1277-d period – longer by a factor of 7 than any other known BH binary – and is moderately eccentric, with e = 0.52. The inclination is ≈35 deg, i.e. relatively face-on. The dynamically implied mass of the companion is |$M_2 = 8.9\pm 0.3\, {\rm M}_{\odot }$|⁠. This constraint depends only weakly on the assumed mass of the luminous star, but significantly on the inclination (Fig. 5).

• Limits on luminous companions: the SED from the UV to the near-IR (Fig. 2) and high-resolution spectrum (Figs 7 and 8) are both well fitted by single-star models. A luminous companion contributing more than |$\approx 3~{{\ \rm per\ cent}}$| of the light at 4000 Å is ruled out spectroscopically. The giant is cool (⁠|$T_{\rm eff} = 4600\, \rm K)$|⁠, and so the fact that the source is faint in the UV rules out any scenario in which the companion is one or more main-sequence stars (Fig. 2). These limits rule out any plausible luminous stellar companion if the Gaia solution is correct. Even if the inclination inferred from the Gaia orbit were spurious – in which case the RVs would allow a companion mass as low as |$2.6\, {\rm M}_{\odot }$| – a stellar companion would be detectable.

• Nature of the companion: the dynamically implied mass of the companion – |$8.9\pm 0.3\, {\rm M}_{\odot }$| – is too large to be a neutron star or white dwarf, or any plausible combination of white dwarfs, neutron stars, and faint stars. The simplest explanation is a single BH. As with Gaia BH1, it is possible that the companion is an unresolved binary containing at least one BH.

• Formation history: the current orbit of Gaia BH2, with a periastron distance of |$2.4\, \rm au$|⁠, is too close to accommodate BH progenitors with initial masses 20 ≲ M_init/M_⊙ ≲ 65, which are predicted to become red supergiants with radii larger than the current orbit (Fig. 15). One possibility is that the BH progenitor had a sufficiently high mass (⁠|$M_{\rm init}\gtrsim 65\, {\rm M}_{\odot }$|⁠, though the precise value varies between models) to lose its envelope while on the main sequence, avoiding a red supergiant phase. It is also possible that the system was formed dynamically through an exchange interaction in a dense cluster, or through triple evolution.

• X-ray and radio observations: we observed Gaia BH2 with both X-ray and radio facilities but did not detect it (Figs 11). Expressed in terms of the BHL accretion rate, the non-detection implies a radiative efficiency ≲10⁻⁴. The non-detection is consistent with model predictions for radiatively inefficient accretion flows in which most of the inflowing material is lost to winds or accretion is suppressed by convection, such that only a small fraction of the BHL accretion rate actually reaches the event horizon. For typical radiative efficiencies, the non-detection is inconsistent with models in which a mass flux equal to the BHL rate reaches the event horizon.

• The Gaia BHs in context: Gaia DR3 has thus far resulted in the identification of two widely separated BH + normal star binaries, whose properties are quite different from previously known BH binaries (Figs 12–14). Both have longer periods and are closer to Earth than any previously known BHs. This suggests that widely separated BH binaries significantly outnumber close systems with ongoing mass transfer. Both have masses between 8 and 10 M_⊙. There do not appear to be any BHs with astrometric solutions published in DR3 that have masses between 3 and 8 M_⊙. The two BHs have periods at the short- and long-period edges of the Gaia sensitivity curve (Fig. 16), perhaps hinting at a bimodal intrinsic period distribution for BH binaries.

ACKNOWLEDGEMENTS

We thank Frédéric Arenou, Jim Fuller, and Andy Gould for useful discussions, and Maren Hempel, Paul Eigenthaler, and Régis Lachaume for carrying out the FEROS observations. We are grateful to the ESO, MeerKAT, Chandra, and Swift directorial offices and support staff for prompt assistance with DDT observations.

HWR acknowledges the European Research Council for the ERC Advanced Grant [101054731]. This research used pystrometry, an open source Python package for astrometry timeseries analysis (Sahlmann 2019). This work used astropy,⁵ a community-developed core Python package and an ecosystem of tools and resources for astronomy (Astropy Collaboration 2022), and the ciao software (Fruscione et al. 2006) provided by the Chandra X-ray Center (CXC). We acknowledge the use of public data from the Swift data archive.

This project was developed in part at the Gaia Fête, held at the Flatiron Institute’s Center for Computational Astrophysics in 2022 June, and in part at the Gaia Hike, held at the University of British Columbia in 2022 June. This work has used 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.

The MeerKAT telescope is operated by the South African Radio Astronomy Observatory, which is a facility of the National Research Foundation, an agency of the Department of Science and Innovation. The scientific results reported in this article are based in part on observations made by the Chandra X-ray Observatory, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of the National Aeronautics Space Administration under contract NAS8-03060.

DATA AVAILABILITY

Data used in this study are available upon request from the corresponding author.

Footnotes

References

APPENDIX A: GAIA OBSERVATIONS

A1 Predicted astrometry

Fig. A1 shows the plane-of-the-sky motion of the photocentre of Gaia BH2 as it would be observed by Gaia. In the left panel, the dashed blue line shows the motion expected due to parallax and proper motion alone, while the black and grey lines show the combined effects of orbital motion, parallax, and proper motion. Lines are coloured darker in the time interval covered by Gaia observations included in DR3, between 2014 July and 2017 May. The interplay between parallax and proper motion leads to rather complex motion on the plane of the sky, which – given the limited time baseline of DR3 – does not yet appear periodic. None the less, the Gaia data alone constrain the orbit fairly precisely.

In the right panel, the individual cyan lines show draws from the posterior of the fit to the pure Gaia solution, with proper motion and parallax motion removed. Although there is noticeable dispersion in these predictions owing to observational uncertainties, the size and shape of the orbital ellipse are well-determined. The middle panel shows draws from the posterior of the joint fit of the Gaia data and RVs. As expected, the resulting orbital ellipse is consistent with the Gaia-only fit, but its parameters are more tightly constrained.

A2 Observation times

The AstroSpectroSB1 solution for Gaia BH2 is based on both astrometric and RV measurements. Single-epoch astrometric and RV data are not published in DR3, but some insight into when a source was observed can be provided by the Gaia observation scheduling tool (GOST).⁶ Given a source position, GOST provides a list of observation times when the scanning law predicts that a source will transit across the Gaia focal plane. For Gaia BH2, GOST predicts 97 CCDs transits during the observation window that contributed to DR3, distributed over 27 visibility periods (i.e. groups of observations separated by at least 4 d). Meanwhile, the gaia_source table reports matched_transits = 93 and visibility_periods_used = 26. This suggests that a majority of the planned transits indeed contributed to the DR3 solution, but a few were omitted, as is common. Fewer transits contributed RVs, with rv_nb_transits = 32 and rv_visibility_periods_used = 12. This is also expected, since only a fraction of focal plane transits cover the RVS spectrometer (Gaia Collaboration 2016).

In Fig. A2, we show the predicted astrometric positions and RVs at observation times predicted by GOST. It is important to note that (a) epoch-level data are not published in DR3, so there is no guarantee that the actual observed positions and RVs match the model predictions, (b) we have not simulated noise, (c) actual Gaia astrometric measurements are 1D, with the scan angle varying from visit to visit, and (d) while a majority (93/97) of the planned transits actually contributed data to the astrometric solution, only about one-third (32/93) contributed an RV.

We can nevertheless draw several conclusions from Fig. A2. The predicted astrometric measurement times span 1000 d–79 per cent of the binary’s orbital period–and cover most of the predicted plane-of-the sky motion. The typical actually realized along-scan uncertainty for each astrometric measurement for a source with G = 12.3 is 0.2 mas (Lindegren et al. 2021a), which is comparable to the size of the red points in Fig. A2. The expected per-epoch RV uncertainty for a source of Gaia BH2’s apparent magnitude and spectral type is also about 1.5 |$\rm km\, s^{-1}$| (see e.g. El-Badry et al. 2019, their fig. F1), also comparable to the size of the points in the lower right panel. The reported rv_amplitude_robust = |$37\, \rm km\, s^{-1}$| (representing the total range in measured epoch RVs, after clipping of outliers) suggests that RVs were measured over most of the orbit’s dynamic range in RV, though the two transits predicted near the periastron passage in 2016 likely did not result in an RV measurement.

We conclude that the Gaia measurements likely are adequate to constrain the orbit robustly, though it would obviously be preferable to have access to the actual epoch-level data.

APPENDIX B: GAIA GOODNESS OF FIT

One of the quality diagnostics calculated for Gaia orbital solutions is the goodness_of_fit estimator, F₂, which is defined (e.g. Halbwachs et al. 2022) as

Here, ν is the number of degrees of freedom and χ² is the sum of the uncertainty-normalized squared residuals between model and data. When there is good agreement between model and data and the uncertainties are reliably estimated, one expects F₂ to be normally distributed, i.e. |$F_2 \sim \mathcal {N}(0,1)$|⁠. Values of F₂ ≫ 1 thus indicate a formally poor fit. Gaia BH2 has F₂ = 3.07, raising potential concern about the quality of the Gaia solution.

In Fig. B1, we compare the F₂ values of Gaia BH2 and BH1 to those of all astrometric binary solutions published in DR3 (including both Orbital and AstroSpectroSB1 solutions). This comparison reveals that Gaia BH2 actually has a somewhat lower F₂ value (indicative of a better fit) than typical sources at its apparent magnitude: almost all sources with G < 13 have F₂ > 1. This most likely is a result of the change in Gaia window class at G = 13 (Rowell et al. 2021): epoch astrometry for sources with G < 13 is derived from a 2D window fitted with a point spread function, while epoch astrometry for fainter sources is calculated from 1D data. This leads to discontinuities in many astrometric quantities at G = 13 (e.g. El-Badry, Rix & Heintz 2021; Lindegren et al. 2021b).

The fact that Gaia BH2 has F₂ > 1 does not necessarily imply that the uncertainties on the AstroSpectroSB1 solution – which also enter into our joint fit – are underestimated. In calculating the astrometric solution, the uncertainties of the epoch-level data are inflated by multiplying by a constant, |$c=\sqrt{\chi ^{2}\nu ^{-1}\left(1-\frac{2}{9\nu }\right)^{-3}}$|⁠, after which F₂ = 0 (Halbwachs et al. 2022; Pourbaix et al. 2022).

APPENDIX C: TESS LIGHT CURVES

The eleanor light curves of Gaia BH2 (TIC 207885877) are shown in Fig. C1. We extract these using a 15 × 15 pixel target pixel file, with the background estimated in a 31 × 31 pixel region. We use the systematics-corrected flux in the default aperture chosen by eleanor, i.e. the aperture with the lowest combined differential photometric precision (CDPP) on 1-h time-scales. The expected flux contamination from background stars in the aperture is of the order of 10 per cent.

The source is clearly not strongly variable on short time-scales. Much of the apparent variability is noise; some variability on longer time-scales may be the result of imperfect background subtraction and removal of systematics. The light curve from sector 38 is noisier than the one from sector 11, mainly because sector 38 used three times shorter exposure times.

Power spectra for both sectors are shown in the right panels of Fig. C1. The red line shows the result of smoothing the power spectra with a Gaussian kernel with |$\sigma = 5\, \mu \rm Hz$|⁠. Several potentially significant peaks are evident in the smoothed power spectra. We tentatively identify the peak with |$\nu = 61\, \mu \rm Hz$| as the global asteroseismic parameter ν_max, since it is the strongest feature evident in data from both sectors. We emphasize that the identification is tentative since the time baseline and signal-to-noise ratio of the TESS data is still insufficient to unambiguously constrain ν_max in a source of this apparent magnitude (e.g. Stello et al. 2022). The implied surface gravity is |$\log g\approx \log g_{\odot }+\log \left(\nu _{{\rm max}}/\nu _{{\rm max,\odot }}\right)+\frac{1}{2}\log \left(T_{{\rm eff}}/T_{{\rm eff},\odot }\right)\approx 2.68$|⁠, where g is in cgs units and log g_⊙ ≈ 4.43, |$\nu _{\rm max,\odot }=3090\, \rm \mu Hz$|⁠, and |$T_{\rm eff,\odot }=5777\, \rm K$|⁠.