ESPRESSO Observations of Gaia BH1: High-precision Orbital Constraints and no Evidence for an Inner Binary

We present high-precision radial velocity observations of Gaia BH1, the nearest known black hole (BH). The system contains a solar-type G star orbiting a massive dark companion, which could be either a single BH or an inner BH + BH binary. A BH + BH binary is expected in some models where Gaia BH1 formed as a hierarchical triple, which is attractive because they avoid many of the difficulties associated with forming the system through isolated binary evolution. Our observations test the inner binary scenario. We have measured 115 precise RVs of the G star, including 40 from ESPRESSO with a precision of 3–5 m s−1, and 75 from other instruments with a typical precision of 30–100 m s−1. Our observations span 2.33 orbits of the G star and are concentrated near a periastron passage, when perturbations due to an inner binary would be largest. The RVs are well-fit by a Keplerian two-body orbit and show no convincing evidence of an inner binary. Using REBOUND simulations of hierarchical triples with a range of inner periods, mass ratios, eccentricities, and orientations, we show that plausible inner binaries with periods P inner ≳ 1.5 days would have produced larger deviations from a Keplerian orbit than observed. Binaries with P inner ≲ 1.5 days are consistent with the data, but these would merge within a Hubble time and would thus imply fine-tuning. We present updated parameters of Gaia BH1's orbit. The RVs yield a spectroscopic mass function fMBH=3.9358±0.0002M⊙ —about 7000σ above the ∼2.5 M ⊙ maximum neutron star mass. Including the inclination constraint from Gaia astrometry, this implies a BH mass of M BH = 9.27 ± 0.10 M ⊙.


INTRODUCTION
The Milky Way is expected to contain ∼ 10 8 stellar mass black holes (BHs) (Brown & Bethe 1994;Chawla et al. 2022).A tiny fraction of this population has been observed, with the presence of a BH dynamically confirmed in only ∼20 X-ray bright systems to date (Remillard & McClintock 2006).While BH X-ray binaries are easier to find than wider binaries hosting non-accreting (i.e.dormant) BHs, they are intrinsically rare, with population models estimating that only ∼ 10 3 exist in the Milky Way (Portegies Zwart et al. 1997;Corral-Santana et al. 2016).
Searches for dormant BHs via spectroscopic techniques began even before the identification of the first BH X-ray binaries (Guseinov & Zel'dovich 1966;Trimble & Thorne 1969).In the intervening decades, there have been extensive spectroscopic searches for dormant BHs in Galactic binaries, but only in recent years have these searches begun to yield reliable BH detections (Giesers et al. 2018;Shenar et al. 2022).
Perhaps the most promising development in the search for dormant BHs thus far is the advent of precise, wide-field astrometry from the Gaia mission, which measures the positions, parallaxes, proper motions, and (potentially) binarityinduced wobble of ∼ 2 billion stars.The mission's 3rd data release ("DR3", Gaia Collaboration et al. 2023a,b) represents a factor of ∼ 100 increase in sample size over previous catalogs of binary orbits and is thus a promising dataset to search for various classes of rare binaries, including luminous stars orbiting dormant BH companions (Breivik et al. 2017;Janssens et al. 2022).
Recently, El-Badry et al. (2023a, hereafter E23) used data from Gaia DR3 to identify Gaia BH1, a Sun-like star in a 185-day orbit around a dark companion with inferred mass M 2 = 9.62 ± 0.18 M ⊙ .At a distance of only 480 pc, the system is the nearest known BH by a factor of ∼3.While E23 identified the companion as a likely dormant stellarmass BH, they noted that the data were also consistent with the unseen object being a close binary containing two BHs with total mass ≈ 9.6 M ⊙ .This raised the tantalizing possibility that Gaia BH1 could contain a wider cousin of the merging BH binaries now routinely detected at cosmological distances via gravitational waves.
Indeed, a BH binary in Gaia BH1 could solve some puzzles related to the system's formation.Isolated binary evolution models struggle to form "intermediate separation" BH + low-mass star binaries like Gaia BH1.The current orbital separation, a ≈ 1.4 AU, is significantly smaller than the predicted maximum radius of a solar-metallicity progenitor of a ∼ 10 M ⊙ BH, suggesting that the Sun-like star and BH progenitor would have interacted when the progenitor was a red supergiant.However, a Sun-like star is unlikely to successfully eject the envelope of a massive star, and is particularly unlikely to end up in an orbit as wide as 1.4 AU if it does survive.These hurdles could potentially be avoided if the system formed as a hierarchical triple, with the solar-type star orbiting two massive stars.In that case, the two massive stars could prevent one another from expanding to red supergiant dimensions, and after two episodes of mass transfer, the Sun-like star could find itself orbiting two BHs without ever having interacted with either one.This and related scenarios have been discussed as a possible formation channel for Gaia BH1-like binaries by several works (e.g.E23; Chakrabarti et al. 2023;El-Badry et al. 2023b;Hayashi et al. 2023;Di Carlo et al. 2023).Besides the triple scenario, other solutions for the current orbit have been proposed (see Section 5.5).
Binary population synthesis simulations predict that BH + BH binaries should be abundant, significantly outnumbering BH + star binaries (e.g.Shao & Li 2021).Whether many of these binary black holes (BBHs) form with distant tertiaries is uncertain, depending on modeling of complex processes in triple evolution (e.g.Silsbee & Tremaine 2017;Toonen et al. 2020).There is little doubt, however, that many massive stars form in triples (Moe & Di Stefano 2017), and so it is natural to search for the BH binary + normal star triples they may evolve to.
In this work, we explore and test the possibility that Gaia BH1 contains an inner BH binary.If an inner binary exists, its orbital motion would introduce non-Keplerian perturbations to the radial velocity (RV) curve of the outer star, which could be detectable with high-precision RV measurements (e.g.Hayashi et al. 2020;Hayashi & Suto 2020;Liu et al. 2022;Hayashi et al. 2023).While the expected perturbations are small, Gaia BH1 contains a relatively bright, Sun-like star, for which RVs can indeed be measured with very small uncertainties.
The remainder of this paper is organized as follows.In Section 2, we explore the orbital parameter space of possible inner BH binaries via simulations and make predictions for the observed amplitude of the RV residuals of the outer star with respect to the Keplerian case.In Section 3, we present our precision RVs for the outer star collected using ESPRESSO on the Very Large Telescope, along with a larger set of RVs for the outer star collected using various other lower resolution spectrographs.In Section 4, we provide an updated Keplerian fit to the RV curve of the outer star and compare our derived parameters to that of E23.While we find no convincing evidence of an inner BH binary, we also fit a hierarchical triple model to our spectroscopic data assuming the presence of a inner BBH, and derive the corresponding best-fit orbital parameters.In Section 5, we discuss the implications of our results for the potential formation channels of Gaia BH1.Finally, in Section 6, we summarize our main findings and consider avenues for follow-up study on the binarity of Gaia BH1.

EXPECTED RV SIGNATURES OF A BH BINARY
We begin by summarizing the expected RV perturbations due to a BH binary for a variety of inner binary periods, mass ratios, eccentricities, and orbital configurations.

Circular and coplanar orbits
In the case of a coplanar and circular hierarchical triple, the RV of an outer star (m * ) orbiting an inner binary (m 1 and m 2 ) observed by a distant observer has an approximate analytic solution given by perturbation theory (Morais & Correia 2008).While the mean motion of the outer star about the center of mass of the inner binary is approximately Keplerian, the short-term RV modulations are non-Keplerian, with characteristic semi-amplitude given by: (1) where m 12 ≡ m 1 + m 2 , m 123 ≡ m 12 + m * , and P inner and P outer are the orbital periods of the inner and outer orbits, respectively (Morais & Correia 2008;Hayashi et al. 2020).These short-term variations are smaller than the unperturbed Keplerian semi-amplitude by a factor (P inner /P outer ) 7/3 and have a dominant period ≈ P inner /2 (Morais & Correia 2008;Hayashi et al. 2020).For an equal-mass ratio inner binary with P inner = 6 days, the expected semi-amplitude for Gaia BH1 in the analytical case is ≈ 6 m s −1 .

Eccentric and inclined orbits
While the prescription above provides an analytic approximation to the RV perturbations of the outer star in the simplest possible scenario, physical binaries are rarely circular or coplanar.In the case of Gaia BH1, we know that the orbit of the G star is not circular (it has eccentricity e outer ≈ 0.45; see E23).Moreover, we expect the inner orbit (if it exists) to be eccentric and misaligned with the outer orbit because both BHs' natal kicks would have perturbed it.In general, we expect the orbit of the outer star to precess, increasing the complexity of the dynamics of the system.While Morais & Correia (2011) provide analytic approximations for the cases in which the orbits of the inner binary and outer star are circular and non-coplanar or eccentric and coplanar, detailed study of the general case requires numerical integrations (e.g.Hayashi & Suto 2020).
We calculate the expected RV perturbations in the noncircular, non-coplanar regime using REBOUND, a flexible Nbody integrator (Rein & Liu 2012).To begin, we use the default IAS15 integrator and uniformly sample 1000 epochs over one orbital period of the outer star (later, we will use the observing cadence of our measured RVs).In general, for each simulation, we fix the orbital parameters and mass of the outer star to the values determined by E23 (see Table 1), and leave the orbital parameters and mass ratio of the inner black hole binary as free parameters.We adopt default values for any REBOUND parameters that are not explicitly specified here.

Computing the RV residual curve
In most cases, the RV curve of the outer star is indistinguishable by eye from a simple Keplerian orbit, because the perturbations are small compared to the total RV amplitude of the outer orbit.To better visualize the non-Keplerian perturbations induced by the inner binary, we fit the predicted RV curve with RadVel (Fulton et al. 2018), a code for fitting Keplerian orbits, and plot the residuals of this fit.We initialize the parameters in the (P, T p , e, ω, K) basis (see Table 1 for descriptions of orbital parameters).We use a generic Keplerian RV model and Gaussian RV likelihood, and fit in the alternative (P, T c , √ e cos ω, √ e sin ω, K) basis, which imposes flat priors on all orbital elements and helps speed MCMC convergence (Fulton et al. 2018).Here, T c refers to the epoch of conjunction.Finally, we use the Powell minimization method with a tolerance of 10 −5 to derive the bestfit orbital parameters.
A visualization of the RV residuals that result from this process for a typical case, where the orbital parameters of the outer star are fixed to those determined for Gaia BH1 by E23 (see Table 1) and the inner equal-mass BH binary is eccentric and non-coplanar, is displayed in Figure 1.We set the period, eccentricity, and inclination (relative to the line of sight) of the inner orbit to be 6 days, 0.1, and 30 • respectively, with the remainder of the orbital elements matching those of the outer orbit.The qualitative behavior of the residuals is insensitive to the phase an orientation of the inner orbit.The maximum amplitude of the RV residuals, which occurs at periastron, is about 250 m s −1 .The characteristic short-timescale variability is more prominent in some portions of the orbit than others because the orbital plane of the outer star is not aligned with the orbital plane of the inner binary; consequently, the gravitational influence of the inner BBH on the outer star is stronger when the three bodies are aligned than when they are not.We also observe long-term (i.e. on timescales longer than P inner /2) variation in the RV residual curve, which we attribute to precession in the orbital parameters of the outer star.

Effects of RV uncertainties and finite observing cadence
Measuring the RV curve with the high cadence and precision shown in Figure 1 is currently infeasible.Here, we consider an observing strategy that is roughly representative of the observations we actually carry out for Gaia BH1 (see Section 3).To initialize the REBOUND simulation, we fix the orbital parameters of the outer orbit to those determined by E23 and assume the same orbital parameters for the inner binary as in Figure 1.Then, adopting a per-epoch RV uncertainty of 2.5 m s −1 (as might be expected for ESPRESSO observations), we sample the RV curve of the outer star at ∼50 epochs spaced every few days, with increased cadence (i.e.every day) near periastron.
We show the resulting RV residual curve for P inner = 6 days in Figure 2. At this inner binary period, the residuals are highly significant, with a maximum semi-amplitude of about 150 m s −1 .This is almost two orders of magnitude larger than the assumed uncertainties, leading to a very poor fit with reduced χ 2 ≈ 672.In the right panel, we show the reduced χ 2 value of the RadVel fit for inner binary periods ranging from 0.5 to 10 days.As expected, it is close to 1.0 for P inner ≪ 1 d, but rises steeply at P inner > 1 d, when the perturbations due to the inner binary begin to exceed the RV measurement uncertainties.This suggests that we can reasonably expect to detect any inner binary with P inner ≳ 1 day.

Dependence on inner binary mass ratio, eccentricity, and inclination
In Figure 3, we explore how the predicted deviations from a Keplerian orbit depend on the parameters of the inner binary.In the left panels, we show the residuals after subtracting the best-fit Keplerian orbit, always assuming P inner = 6 days.In the right panels, we show the RV residual amplitude (defined as the difference between the maximum and minimum RV residuals over the observing baseline) as a function of P inner , with other parameters held fixed.We fit a power law to each residual amplitude plot and report the best-fit power law index (denoted by α) in each panel.
Each row shows a different orbital configuration.In the top row, we simulate the case where the orbits are coplanar and both are circular, and find that the power law index is close to the theoretical expectation of 7/3 = 2.33 (see Section 2.1).In the next three rows, we fix the orbital parameters of the outer star to those determined by E23, and simulate cases in which the inner orbit is circular and coplanar with the outer orbit (2nd row), circular and perpendicular to it (3rd row), and eccentric and coplanar (4th row).In these cases, the trend is not a perfect power law, and the best-fit power law index is smaller than 7/3.This reflects the fact that the short-timescale and long-timescale perturbations have different scalings with P inner (Hayashi et al. 2023).In these cases, numerical integrations are critical, since the behavior of the RV residuals is not easy to explain analytically.Note that the RV residuals are always larger for the realistic case than for the circular + coplanar case on which most previous work has focused (Morais & Correia 2008;Hayashi et al. 2020).This is primarily because the eccentric outer orbit allows the inner BBH to get closer to the outer star at periastron.Orbital configuration (left) and predicted RVs (right) of the hierarchical triple scenario we seek to test.RVs are measured for a Sunlike star orbiting an inner BH + BH binary.In the general case, the inner and outer orbits are both eccentric and non-coplanar.The parameters of the outer orbit are fixed to those inferred for Gaia BH1 by E23.Here we assume an equal-mass inner binary with period Pinner = 6 days and eccentricity einner = 0.1.While the RVs of the Sun-like star are nearly consistent with a Keplerian orbit (upper right), subtraction of the best-fit Keplerian orbit reveals significant residuals (lower right).At Pinner = 6 days, the amplitude of the RV residuals is ∼250 m s −1 near periastron.

Distinguishing an inner binary from a planet
Finally, we simulate a case in which there is no inner BH binary, but there is a 10 −3 M ⊙ planet orbiting the outer star (5th row of Figure 3).Jettisoning considerations about how such a planet may have formed and survived until now, we assume it to have a circular orbit about the outer star and match the inclination and longitude of the ascending node of its orbit to that of the star's orbit.We find that the resulting RV residual curve is sinusoidal, and that the amplitude decreases with increasing period (α = −1/3), as expected.
This case is relevant because a planet could mimic a signal from an inner binary -indeed, perturbations due to an inner binary were considered as a false-positive for exoplanet searches before it became clear that exoplanets are more common than triple star systems with the relevant configurations (e.g.Morais & Correia 2008).One way to distinguish between the two cases is that the amplitude of the RV variation in the exoplanetary case does not significantly increase at the periastron of the outer orbit.

Short-and long-timescale RV perturbations
As is evident from the left panels of Figure 3, the predicted residuals due to an inner binary show features on both short and long timescales.The short-timescale features with dominant period P inner /2 reflect the dipole-like oscillations in the gravitational field felt by the star due to the inner binary's motion.The long-timescale residuals are a result of precession.Figure 4 shows example residuals for a range of inner binary periods.We fix the parameters of the outer orbit to those determined by E23 and constrain the inner orbit to be coplanar and circular.The amplitude of both short-and longtimescale RV residuals increases with the orbit of the inner BBH, varying from about 50 m s −1 at P inner = 3 days to 350 m s −1 at P inner = 9 days.At small P inner , the long-term perturbation due to precession is larger than the short-timescale variations.At longer P inner , the two amplitudes are similar.

Dependence on the observing time baseline
We also investigate how the amplitude of the RV residuals varies with the length of the observational baseline.In Figure 5, we again fix the orbital parameters of the outer star to those determined by E23 and constrain the inner orbit to be coplanar and circular.We then simulate observations over one, two, and three orbits.While the power law slope depends only weakly on the length of the observational baseline, extending RV coverage to two or three orbital cycles leads to significantly larger-amplitude RV residuals than coverage of one orbit.This is because precession of the outer orbit has an increasing cumulative effect when the observations cover multiple orbital cycles.We thus expect that tighter constraints on an inner binary can be obtained with a longer observing baseline.

RV residual amplitudes for plausible inner binary parameters
Figure 6 shows the expected RV residual amplitudes as a function of orbital period for simulations in which we systematically vary the inclination, eccentricity, and mass ratio of the inner BBH.To facilitate comparison with previous work that assumed a circular outer orbit, we also vary the outer orbit's eccentricity (though it is known for Gaia BH1).
Increasing the eccentricity of either orbit increases the amplitude of the observed residuals.This is because the outer star reaches a smaller minimum separation relative to the inner BBH at periastron.On the other hand, increasing the mass ratio of the BHs decreases the amplitude of the observed residuals; after all, in the limit of q inner → ∞, we would recover the two-body Keplerian orbit.We find that the mutual inclination of the orbits has a small effect on the observed residuals, with the residual amplitude varying non-monotonically with i inner .In almost all cases, the observed amplitude of short-term non-Keplerian RV modulations is above our detectability threshold of ≈ 5 m s −1 (see Section 3) for P inner > 1 day.The only exceptions are the models with e outer = 0 (which is ruled out since we know e outer ∼ 0.45) or q inner = 100 (which is astrophysically unlikely).
For sufficiently large inner orbits and high inner or outer eccentricities, the triple is expected to become unstable on a short timescale.We use the approximate dynamical stability criterion derived by Aarseth & Mardling (2001) to restrict the power law curves in Figure 6 from entering regions where the hierarchical triple system is unstable.Specifically, this is relevant for large P inner and highly eccentric orbits; an example of this is the truncation of the e = 0.75 power law curve in the lower left panel of Figure 6.
3. DATA We now describe the observed RVs.In brief, we analyze 115 RVs obtained over 432 days, including 40 from ESPRESSO with a typical precision of 3-5 m s −1 , and 75 from FEROS, TRES, and HIRES with a typical precision of 30-100 m s −1 .All of these RVs are listed in Table 3.
We reduced the data using the CERES pipeline (Brahm et al. 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 of a night via simultaneous observations of a ThAr lamp with a second fiber.We first calculate RVs by cross-correlating a synthetic 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., σ RV ≈ std (RVs) / √ 15.We used a Kurucz model spectrum (Kurucz 1979(Kurucz , 1993) ) with T eff = 5750 K, log g = 4.5, and [Fe/H] = −0.25 from the BOSZ library (Bohlin et al. 2017) as a template.
RVs from the first 17 observations were already presented by E23.However, we re-reduced those data with the CERES pipeline and achieved a significantly more stable wavelength solution compared to the ESO MIDAS reduction described by E23.The FEROS RVs analyzed here and listed in Table 3 thus supersede those measured by E23.The median uncertainty of the FEROS RVs is ≈ 70 m s −1 .

TRES
We obtained 13 spectra using the Tillinghast Reflector Echelle Spectrograph (TRES; Fűrész 2008) mounted on the 1.5 m Tillinghast Reflector telescope at the Fred Lawrence Whipple Observatory (FLWO) atop Mount Hopkins, Arizona.TRES is a fibrefed echelle spectrograph with a wavelength range of 390-910 nm and spectral resolution R ≈ 44, 000.Exposure times ranged from 1800 to 3600 seconds.We extracted the spectra as described in Buchhave et al. (2010).
As with the FEROS data, we measured RVs by crosscorrelating the normalized spectra from each of 31 orders with a template, and we estimate RV uncertainties from the dispersion between RVs measured from different orders; i.e., σ RV = std (RVs) / √ 31.We used the same Kurucz template from the BOSZ library that we used for the FEROS data (T eff = 5750 K, log g = 4.5, [Fe/H] = −0.25).The median uncertainty of the TRES RVs is ≈ 50 m s −1 .

HIRES
We analyze 9 telluric-calibrated RVs measured with the High Resolution Echelle Spectrometer (HIRES; Vogt et al. 1994) mounted on the 10 m Keck telescope at W. M. Keck Observatory at Mauna Kea, Hawaii.These are the same HIRES RVs that were analyzed by E23.We adopt an uncertainty of 100 m s −1 for all of these RVs.

ESPRESSO
We observed Gaia BH1 40 times with the Echelle SPectrograph for Rocky Exoplanets and Stable Spectroscopic Observations (ESPRESSO; Pepe et al. 2021) at the VLT (program 111.24GP.001).We used singleHR mode, in which the instrument can observe with any of the four 8.1 m Unit Telescopes (UTs).We used 2 × 1 binning, yielding a typical spectral resolution R = 140, 000 over a wavelength range of 380-686 nm, and used 900s exposures for all observations.The typical SNR at 600 nm was 20 and ranged from 15-25, depending on seeing and lunar phase.The ESPRESSO observations are spread over a 170 day period, spanning most of one orbit.The typical cadence was one observation every 3-6 days away from periastron, and nearly daily observations near periastron.Gaps in coverage near periastron were due to bad weather or scheduling constraints.
We reduced the data using version 3.0.0 of the ESPRESSO DRS pipeline maintained by ESO and executed through the EsoRex software.We subsequently measured RVs using the espda compu radvel routine within the ESPRESSO DAS package (version 1.3.7).This routine cross-correlates  We recover the theoretical expectation of a power law slope of 2.33 in the case where the orbits are coplanar and both are circular (Equation 1).In the other cases, the power law slope is lower, but the amplitude of the RV residuals is higher at all detectable values of Pinner.Bottom panel shows a case where there is a single BH, but a Jupiter-mass planet orbits the star.Here we fix the orbital parameters of the outer star to those determined by E23 and make the inner orbit coplanar and circular.At short Pinner, the residuals are dominated by a long-term undulation due to precession of the outer orbit.A lower-amplitude "wobble" on half the period of the inner binary is also present.Both the long-and short-timescale residuals become larger with increasing Pinner, but the short-timescale residuals grow more quickly, such that short-and long-timescale residuals have comparable amplitude at Pinner = 9 days.Aarseth & Mardling (2001) to restrict the power law curves from entering unstable regions in parameter space.In general, increasing the eccentricity of either orbit increases the amplitude of the observed residuals, while increasing the BBH mass ratio decreases the amplitude of the observed residuals.In addition, the residual amplitude varies non-monotonically with the orbital inclination of the inner BBH.The dashed horizontal line marks a residual amplitude of 5 m s −1 , roughly the sensitivity of our observations.This suggests that our observations are sensitive to essentially all plausible inner binaries with Pinner significantly above 1 day, except in the (physically dubious) case of an inner BBH mass ratio of 100.
all orders of the extracted spectra with a G2 weighted binary mask and then fits the resulting cross-correlation function (CCF) with a Gaussian.The mask gives maximum weight to strong lines and masks regions of the spectra affected by telluric absorption, chromospheric activity, and interstellar absorption (e.g.Pepe et al. 2002;Lafarga et al. 2020).We evaluated the CCFs on a 0.5 km s −1 grid.RV uncertainties are calculated from the curvature of the CCF (see Boisse et al. 2010, their Appendix A).The median RV uncertainty is 4.5 m s −1 .Examples of two typical spectra from our program are shown in Figure 7.
The high resolution of the ESPRESSO data allows us to set a more stringent limit on the projected rotation velocity of the G star than could be set previously: comparing a SNR ∼ 200 co-added rest-frame spectrum of the G star to Kurucz model spectra broadened with the rotational kernel from Gray (2008), we estimate v sin i < 2 km s −1 .The actual value of v sin i could be significantly lower than this; uncer-tainty in the star's micro/macroturbulent velocities limits an even more precise measurement.

Summary of RVs
We simultaneously analyze lower-precision (30-100 m s −1 ) RV measurements from FEROS, TRES, and HIRES, and higher-precision (3-5 m s −1 ) RV measurements from ESPRESSO, with coverage over 432 days (2.33 orbits).The RVs are shown in Figure 8 and listed in Table 3.
In addition to the FEROS and HIRES RVs discussed above, E23 also included several low-precision (uncertainties σ RV ≳ 1 km s −1 ) RVs from LAMOST, GMOS, X-SHOOTER, MagE, and ESI.We include these data in Figure 8 and Table 3 but do not include them in the fit because they are much less precise than the other data and would require 5 additional free parameters (i.e.instrumental offsets).
We also obtained 16 spectra with the Keck Planet Finder (Gibson et al. 2016), a new high-precision RV spectrograph installed on the Keck telescope.These data should in princi- ple allow RV measurement with uncertainties comparable to ESPRESSO.However, we have not yet succeeded in obtaining long-term stability in the instrument's wavelength solution due to software issues, and thus defer analysis of these data to future work.
4. ANALYSIS We built a predictive model for the RV variation of the outer star in Gaia BH1 in both the two-body Keplerian and hierarchical triple cases.In the Keplerian case, we included the center-of-mass RV, component masses, orbital elements (period, eccentricity, inclination, argument of periastron, longitude of the ascending node, and periastron time) of the outer orbit, and instrumental offsets (with respect to ESPRESSO) as free parameters.In the hierarchical triple case, we used REBOUND simulations in place of an analytical RV model, and additionally included the orbital elements and mass ratio of the inner BBH as free parameters.
In both cases, we placed a Gaussian prior on the mass of the Sun-like star, with a mean of 0.93 M ⊙ and a standard deviation of 0.05 M ⊙ .For all other parameters, we used uniform priors.The log-likelihood assumes Gaussian uncertainties: (2) In order to constrain the inclination and companion mass, we jointly fit the constraints on the star's plane-of-the-sky motion from Gaia.Our approach is very similar to E23.The Gaia astrometric solution is given by joint constraints on 12 astrometric parameters: right ascension, declination, proper motion in right ascension, proper motion in declination, parallax, orbital period, periastron time, eccentricity, and the four Thiele-Innes elements, A, B, F , and G (Halbwachs et al. 2023).Given the best-fit Gaia astrometric parameters µ ast , their covariance matrix Σ ast , and the vector of predicted astrometric parameters for a given likelihood call, θ ast , the astrometric log-likelihood can be expressed as: We calculate θ ast under the assumption that the companion is dark.The total log-likelihood function includes both the astrometric and RV terms: Our joint fitting of the astrometry and RVs in the two-body Keplerian case is almost identical to E23, with the following minor modification.E23 left the source's right ascension, declination, proper motions, and parallax as free parameters during fitting, to be constrained only by the Gaia astrometric solution.To reduce the number of free parameters, we excise these parameters from θ ast and µ ast , and remove the corresponding rows and columns from Σ ast .In predicting the Thiele-Innes elements, we assume a parallax of ϖ = 2.09 mas.We verified that this simplification (which reduces the dimensionality of the fit) speeds up convergence while having no significant effect on the constraints on the parameters of interest.
We began by fitting a two-body Keplerian orbit.Then, we added an additional set of orbital parameters, and tried a hierarchical triple fit instead.We finally considered whether the improvement in the log posterior probability was sufficient to warrant the extra free parameters (i.e.added model complexity).We describe this procedure and our results below.

Updated fit with astrometry for Gaia BH1
Assuming a Keplerian outer orbit and a single inner black hole, we used ensemble MCMC sampling (emcee; Foreman-Mackey et al. 2013) with 64 walkers and 8000 total iterations to derive updated best-fit parameters for Gaia BH1.We show the best-fit Keplerian RV curve and corresponding residuals in Figure 8.Compared to E23, our data now cover three orbital cycles, and our typical RV uncertainties are a factor of 100-1000 smaller, leading to much tighter constraints.
We show the resulting RV residuals at increasing levels of precision in the remaining three panels of Figure 8.At each level of precision, we see that most of the RV residuals are consistent with zero to within 1-2σ (given the reported uncertainties in our measurements).From the third level of the residual plots, we confirm that this also holds true for the highest-precision ESPRESSO data.The small scatter in the ESPRESSO RV residuals implies that the luminous star in Gaia BH1 has low RV jitter (due to e.g.convection, pulsations, or other systematics) of at most ≈ 3 m s −1 (see Luhn et al. 2020 for a discussion).Given the small residuals in the ESPRESSO data, it is clear that the scatter in the less-precise FEROS and TRES data acquired over the same time period is dominated by noise.However, these medium-precision RVs are still quite important for our results, because they cover three orbits and thus tightly constrain the orbital period, while the ESPRESSO data cover only one.Based on the lack of obvious structure in the residuals from a two-body orbit, there is no immediate evidence for deviations from a Keplerian orbit.
In Figure 9, we compare our constraints from the fit described above (black contours) to those from E23.The period, eccentricity, argument of periastron, periastron time, and center-of-mass RV of the orbit of the outer star are much more tightly constrained than they were by E23.The uncertainties in i, Ω, and M BH have decreased only slightly.These parameters are constrained most directly by astrometric data, which has not changed, but they are covariant with other parameters.The constraint on M star comes from our SED-informed prior and is unchanged.We report our updated best-fit values in Table 1.Compared to the E23 constraints, we find that the orbital period, eccentricity, and inferred companion mass have all decreased, while the argument of periastron, periastron time, and center-of-mass RV have increased.Our new value for the companion mass is M BH = 9.27 ± 0.10 M ⊙ , corresponding to an RV semiamplitude K star = 65.3785± 0.0009 km s −1 and a spectroscopic mass function f (M BH ) = 3.9358 ± 0.0002 M ⊙ .Our derived values for the orbital period, eccentricity, argument of periastron, periastron time, and center-of-mass RV differ from the values found by E23 by about 4.1σ, 3.7σ, 3.4σ, 4.5σ and 3.0σ respectively.The remaining orbital parameters are consistent with each other to within 2σ.
The small but significant tension between our best-fit orbital parameters and those of E23 suggests that the uncertainties reported by E23 were somewhat underestimated.There are a few possible reasons this could have occurred.E23 did not fit for instrumental RV offsets, and offsets between the seven spectrographs they used could be significant.However, the RVs analyzed by E23 are nearly all consistent with our updated solution at the 1-2σ level (2nd panel of Figure 8), suggesting that RVs are unlikely to be the main source of the disagreement.Another possibility is that the Gaia astrometric uncertainties are underestimated somewhat (see also Chakrabarti et al. 2023).The RVs obtained by E23 covered less than one orbital cycle, so the period constraint from their joint fit came primarily from astrometry.In contrast, our RVs cover three orbital cycles with high precision; so our constraint on the period and periastron time come almost entirely from the RVs.
Because our measured RVs are generally consistent with the best-fit RV curve within their uncertainties (i.e., we achieve a reduced χ 2 = 1.69), we are confident that the parameters we infer from RVs directly (e.g., period, eccentricity, and mass function) have robust uncertainties.The uncertainties on the parameters that depend mainly on the astrometric orbit -inclination, longitude of the ascending node, and BH mass -are harder to assess, and may be underestimated.It will become easier to robustly estimate uncertainties on these quantities when epoch astrometry is published in Gaia DR4.For now, the parameters reported in Table 1 supersede those reported by E23, and we adopt them for the two-body solution in the rest of this work.
Note that the reported γ in Table 1 represents the system's center-of-mass velocity relative to our adopted ESPRESSO zeropoint.That uncertainty in that zeropoint -which must be accounted for in e.g.kinematic analysis of the system's Galactic orbit or comparison with precision RVs from other instruments or templates -is significantly larger, on the order of ∼ 100 m s −1 (see Lindegren & Dravins 2003 for a discussion).

Hierarchical triple fit for Gaia BH1
Next, we tried fitting the data with a hierarchical triple model implemented with REBOUND.The model is relatively high-dimensional, with 19 free parameters (including instrumental offsets).We find the posterior to be multimodal and challenging to sample from.In the limits of P inner → 0 or q inner → 0, we recover the two-body Keplerian solution.For 0.5 ≲ P inner /day ≲ 1.3, there are multiple distinct combinations of inner binary parameters that yield slightly better fits (i.e., higher posterior probabilities) than the best-fit two-body solution (which, after all, has 7 fewer free parameters).For P inner ≳ 1.5 days, the inner binary perturbs the outer orbit too much, leading to a poorer fit than the twobody solution.Complex posteriors are common when fitting RV data, and robust algorithms have been developed to sample from them (e.g.Price-Whelan et al. 2017).Compared to the standard problem of fitting a two-body orbit, the hierarchical triple fit has the added challenges that (a) the space is higher-dimensional, with coupled inner and outer orbits, making brute-force approaches like rejection sampling infeasible, and (b) each call to the likelihood function requires running a REBOUND simulation, which is more expensive than evaluating the quasi-analytic equations required to predict RVs in the two-body case.
We experimented with a variety of sampling approaches, including directly running ensemble sampling with emcee, running ensemble sampling with emcee after using an optimizer to determine a favorable initialization, and directly running dynamic nested sampling with dynesty.Ultimately, we settled on the approach described below.
We used dynamic nested sampling (dynesty; Speagle 2020) to sample from the 19-dimensional posterior distribution.The free parameters of our fit are listed in Table 2. Compared to the two-body fit, the extra parameters are the period, eccentricity, inclination, longitude of the ascending node, argument of periastron, periastron time, and mass ratio of the inner orbit.As with the two-body fits, we added an astrometric term to the likelihood that compared the parameters of the outer orbit to the Gaia constraints (Equation 4).Since the Gaia solution assumes a two-body orbit, we treated the inner binary as a point mass when calculating the Thiele-Innes elements of the outer orbit.We used the random walk sampling method with 500 live points and set a maximum of 10 6 likelihood calls.We found that the sampler would hit the default stopping condition (defined to be when the estimated remaining evidence falls below the default threshold, see Speagle 2020) slightly before achieving the maximum number of likelihood calls.Consequently, we do not expect better performance with a larger value for the maximum number of likelihood calls.We assumed truncated normal priors on the parameters of the outer orbit based on the results of the two-body Keplerian fit (see Table 1), and flat, broad priors on orbital parameters of the inner binary.We set a lower limit of P inner > 0.5 days because likelihood calls became prohibitively expensive in the limit of P inner → 0.
We ran several dynesty runs with different sampling methods, initializations, and slight variations in the priors, and then compared the maximum posterior probabilities and marginalized constraints on inner binary parameters across runs.While the maximum probabilities were relatively stable across runs (with ln P max varying by ≲ 1), the best-fit parameters varied between runs by more than their formal uncertainties.This suggested that the runs were not fully converged.To explore the posteriors more thoroughly in the vicinity of the maximum probability solutions, we initialized an emcee chain with 64 walkers and 3125 steps at the maximum-probability sample from each dynesty run.These chains achieved slightly higher posterior probabilities than the best dynesty samples.We report the highestprobability solution achieved across all such runs as the MAP (maximum a posteriori) solution in Table 2, where we also report the marginalized median and middle 68% constraints from the corresponding emcee run.The MAP solution falls within the marginalized middle 68% range for all parameters.The best-fit solution is consistent with an inner BH binary with P inner ≲ 1.5 days.
We emphasize that the posterior distribution of the hierarchical triple fit is complex, and that -because our sampling is likely not fully converged -there is no guarantee that our reported solution corresponds to the best absolute solution, or that the reported uncertainties encompass all possible solutions.However, given the small residuals of both the twobody and three-body fits (Figure 12), we consider it unlikely that a significantly better three-body solution exists.For the best-fit solution we report, the improvements over the best two-body solution are marginal.

Limits on period of inner BH binary
To explore limits on the period of the inner binary (if indeed it exists), we present the reduced χ 2 value of a RadVel Keplerian fit to simulated RVs for a hierarchical triple as a function of P inner for various orbital configurations of the inner binary in Figure 10.In contrast to our simulations in Section 2, these simulations assume the exact observing cadence and uncertainties of our measured RVs.
As already demonstrated in Figure 6, the predicted RV signal of an inner binary depends on its period, eccentricity, mass ratio, and orientation.In Figure 10, we provide predictions for four different inner binary orientations, with each panel representing a different combination of inclination and longitude of the ascending node.We show a range of periods and eccentricities for the inner binary in each panel, setting the argument of periastron and periastron time to zero and the mass ratio to unity.We also plot a black dashed line indicating our observed value of reduced χ 2 = 1.69 from the updated Keplerian fit for Gaia BH1.
For most of the orientations and eccentricities shown in Figure 10, the predicted reduced χ 2 rises above the observed value for P inner ≳ 1 day, indicating that we can rule out an inner BBH with P inner ≳ 1 day for such orientations.However, for some inner binary orientations -such as the one labeled "worst case" in the bottom right panel -the predicted RV signature of an inner binary is weaker, only rising above reduced χ 2 = 1.69 for P inner ≳ 3 days.
We explore the sensitivity of our constraints to the orientation of the inner binary more thoroughly in the right panels of Figure 11, where we plot the reduced χ 2 value of the best-fit Keplerian model as a function of inclination and longitude of the ascending node of the inner orbit.We assume q inner = 1 and P inner = 2 days in this exploration, as might be expected for an inner binary just above our detection threshold.As before, we adopt the exact observing cadence and uncertainties of our measured RVs, and set the argument of periastron and periastron time of the inner orbit to zero.We show predictions for two values of the inner binary's eccentricity, e inner = 0 and e inner = 0.6.We also label the orientations corresponding to the "best" and "worst" cases (in terms of detectability, for e inner = 0) displayed in Figure 10.
The reduced χ 2 landscape shows the expected rotational and reflection symmetries.Beyond this, it is complicated, bearing imprints of both the observational properties (i.e.cadence and uncertainties) of the RVs, and the orientation of the outer orbit relative to our line of sight.For both circular and eccentric inner orbits, the most easily detectable inner binary is close to (but not exactly) coplanar with the outer orbit, while the hardest-to-detect orientations are closer to perpendicular.
Given the complicated dependence on inner binary eccentricity, mass ratio, and orientation, it is difficult to determine a single upper limit on the inner binary's period.For the most fine-tuned orientations, periods as long as P inner ∼ 3 days could escape detection in our data with an equal-mass inner binary, and even longer inner periods are possible if the inner binary mass ratio is allowed to be arbitrarily large.To probe further, we generate 1000 random eccentricities, inclinations, and longitudes of the ascending node at each value of P inner on a grid between 0.5 and 3 days.We set the argument of periastron and periastron time of the inner orbit to zero, and calculate the reduced χ 2 value in each case, assuming our observed cadence and RV uncertainties.We assume a uniform e inner distribution and random orientations.Then, at each value of P inner , we compute the fraction of inner binaries that would have been detected at our threshold of reduced χ 2 = 1.69.We plot the resulting curve in the left panel of Figure 11.We conclude that, for typical inner binary orientations and eccentric inner orbits -as are expected in the presence of natal kicks and Kozai-Lidov oscillations (Lidov 1962, Kozai 1962, Naoz 2016) -our observed RVs rule out most inner binaries with P inner ≳ 1.5 days.This is consistent with the best-fit derived value of P inner = 0.9 ± 0.1 days from our hierarchical triple model.

Comparison of the best-fit binary and triple models
We provide a comparison of the observed RV residuals relative to the best-fit two-body Keplerian model with the difference between the best-fit hierarchical triple and two-body Keplerian models in the top panel of Figure 12, focusing on the high-precision ESPRESSO data.We find that the residuals are consistent with the continuous curve representing the difference between our best-fit models to within 1-2σ.
We compare observed RV residuals for the best-fit twobody Keplerian and hierarchical triple models in the bottom panel of Figure 12, once again displaying only the highprecision ESPRESSO data.The log posterior probability (which includes the RVs from all instruments, the Gaia astrometry, and the priors) for the best-fit three body model is −82, which is slightly favorable compared to the log posterior probability for the best-fit two-body model of −94.While the scatter of the ESPRESSO RV residuals around zero for the hierarchical triple model appears to be slightly smaller than in the Keplerian case, the observed amplitude of The reduced χ 2 value generally increases with increasing einner.At very high einner, the non-uniform sampling of our RVs gives rise to oscillations in the reduced χ 2 value.Based on our best-fit model's reduced χ 2 value of 1.69, we can rule out most inner binaries with Pinner ≳ 1.5 days.Fine-tuned orientations exist that could hide inner binaries with orbital periods up to ∼ 3 days (e.g.lower right); we explore these in more detail in Figure 11.
the residuals is small, and the improvement is not significant with respect to the RV uncertainties.
From Figure 12, we conclude that the hierarchical triple model provides a marginally better fit to the observed RVs compared to the two-body Keplerian model.This is because the extra model parameters allow for a long-term variation (due to precession) that can absorb some of the scatter in the residuals.However, this improvement is small compared to the individual uncertainties; at this level of precision, the difference in the log posterior probability is most likely due to overfitting the noise in the observed data.
To emphasize this, we turn to the Bayesian information criterion (BIC), a model selection metric that penalizes increasing model complexity (Schwarz 1978).The BIC is defined to be: where k is the number of model parameters, n is the number of data points, and L is the maximum value of the likelihood function for the model in question.We find that the two-body Keplerian model and hierarchical triple mod-els have BIC values of 246 and 258 respectively.Since lower BIC values are generally preferred, the evidence appears to favor the two-body model over the three-body one.

BH binary inspiral time and detectability with LISA
While an inner BH binary with P inner < 1.5 days could possibly exist in Gaia BH1, this would be a fine-tuned scenario, as such a close binary would have merged in less than a Hubble time.To demonstrate this, we use the results derived by Peters (1964) to calculate the inspiral time for an equal-mass inner BH binary as a function of orbital period at various eccentricities in Figure 13.These calculations account for evolution of the eccentricity (i.e.circularization) during the inspiral.We make no attempt to account for perturbations due to the outer star, which would accelerate the merger via Kozai-Lidov oscillations for some orbital configurations.Based on our best-fit models, we choose a total binary BH mass of 9.3 solar masses.We confirm that the BBH inspiral time is shorter at higher orbital eccentricities.Furthermore, at P inner < 1.5 days, the BH binary merges in less Figure 11.Left: Fraction of simulated inner binaries that are detectable with our data as a function of inner period.We assume random inner binary orientations and a uniform inner binary eccentricity distribution, and we define "detectable" binaries as those which produce a worse Keplerian two-body fit than we observe.The observed RVs rule out most (≳ 95%) inner binaries with Pinner > 1.5 days.Right: Dependence of the best-fit Keplerian model's reduced χ 2 on the orientation of the inner BBH, for two values of the inner binary eccentricity.For each choice of inner binary inclination and longitude of the ascending node, we simulate observations of a triple with the exact observing cadence and uncertainties of our measured RVs and then fit a Keplerian two-body orbit.We assume Pinner = 2 days, as might be expected for an inner BBH just above our detection threshold.The "best" and "worst" cases shown in Figure 10 are labeled.The reduced χ 2 landscape is complex, but only a few fine-tuned inner binary orientations allow reduced χ 2 < 2.
than the age of the Universe, even for cases with low orbital eccentricities.
At sufficiently short periods or high eccentricities, an inner BH binary would be detectable via gravitational waves (GWs) with the LISA space observatory (Amaro-Seoane et al. 2017).We use LEGWORK (Wagg et al. 2022b) to estimate the SNR of an equal-mass BH binary in Gaia BH1 observed by LISA for a range of inner periods and eccentricities, assuming a 4-year mission duration.We find that binaries with P inner ≲ 0.15 days would be detectable with SNR > 5 for any eccentricity.For the same SNR threshold, the maximum detectable inner orbital period rises to 0.33 days for e inner = 0.5 and 2.5 days for e inner = 0.9.At higher eccentricities, the detection is primarily due to the GW signal from higher harmonics of the orbital period.Figure 13 demonstrates that catching any individual BH binary at a period where it would be detectable by LISA but had not already merged requires some fine tuning.Of course, LISA will have the advantage of being sensitive to all close BH binaries in the Milky Way, including those not orbited by a luminous tertiary.5.4.Limits on the presence of distant tertiary star fact that our RVs are well-fit by a Keplerian two-body orbit also places constraints on the presence of objects in wider orbits in the Gaia BH1 system.We investigate the perturbation that would result from a distant tertiary star by plotting the Keplerian RV residual amplitude of the Sun-like star as a function of orbital period of the tertiary in Figure 14.We consider the cases of 1 M ⊙ white dwarf tertiary and a 0.2 M ⊙ M dwarf tertiary, both of which could have evaded spectroscopic detection.For the tertiary star, we adopt an or-bital eccentricity of 0.5 and set both the inclination and longitude of the ascending node to zero.In addition, we assume the same observing time baseline as in Figure 1.For each tertiary, we randomly sample ten orbital phases and plot the resulting curves in Figure 14.
Adopting a detectability threshold of a RV residual amplitude of ∼ 5 m s −1 , we find that we would detect a 1 M ⊙ tertiary at orbital periods less than about 2700 days, or approximately 7 years, in the residuals of a two-body Keplerian fit.For a 0.2 M ⊙ tertiary, the lower limit on the orbital period is about 1500 days, or approximately 4 years, instead.We conclude that sensitivity to a low-mass tertiary star would require spectroscopic follow-up over an extremely long observational baseline.
Tertiaries that are bright enough to be detected by Gaia (corresponding to masses ≳ 0.2 M ⊙ for main-sequence stars at the distance of Gaia BH1) are already ruled out.Assuming a 1 arcsec Gaia resolution (which is somewhat optimistic for the faintest companions due to blending; El-Badry & Rix 2018), this rules out separations wider than ∼ 500 AU, or periods < 10 6 days.A stronger limit of ≈ 0.1 arcsec, corresponding to periods ≲ 10 4.5 days, could be set with speckle interferometry (e.g.Howell & Furlan 2022).

Implications for the formation history of Gaia BH1
We have shown that the ∼ 9.3 M ⊙ dark object in the Gaia BH1 system is unlikely to be a binary with period longer than ∼ 1.5 days.Shorter inner periods cannot be ruled out because they would lead to deviations from a Keplerian orbit that could not be detected with the current data.However, such short inner binaries would merge within a Hubble time The scatter of the residuals for the hierarchical triple model is somewhat smaller than in the Keplerian case, with an improvement in log posterior probability of ∼12.This reflects the fact that the more flexible model allows for a long-term trend in the residuals due to precession.However, the improvement is small compared to the RV uncertainties and appears consistent with overfitting.
and thus require significant fine tuning.Hence, the simplest -and arguably most plausible -explanation is that the dark object is a single, 9.3 M ⊙ BH.It is in principle possible that the dark object was initially a BH + BH or BH + NS binary and has already merged.However, we consider this scenario unlikely because anisotropic gravitational wave emission produces a kick during compact object mergers, with a typical magnitude above 100 km s −1 (e.g.Bekenstein 1973;Merritt et al. 2004).Such a kick would most likely have disrupted the outer binary, or imparted a larger eccentricity and space velocity on it than is observed.It is also possible that the system started as a triple but the inner binary merged before one or both stars died, as some evolutionary models predict merger products to explode as blue supergiants or Wolf-Rayet stars without ever expanding to red supergiant dimensions (e.g.Justham et al. 2014).The range of initial triple conditions for which such an evolutionary scenario is feasible in the Gaia BH1 system is, however, likely rather limited.A related possibility is that the BH could have formed from a progenitor with mass ≳ 40 M ⊙ that never became a red supergiant (e.g.Humphreys & Davidson 1979;Davies et al. 2018;Gilkis et al. 2021), but it is still uncertain whether very massive stars briefly expand to red supergiant dimensions in their postmain sequence evolution.
Another formation channel that has recently been explored is dynamical assembly in a star cluster.The G star's nearsolar metallicity and disk-like orbit disfavor a globular cluster, but the system may have formed in a cluster in the Galactic disk that has since dissolved.This scenario is difficult to test because it makes no predictions for present-day observables that could distinguish a dynamically-assembled BH binary from a primordial one.Several recent works have predicted that dynamical formation in intermediate-or highmass clusters could dominate the formation rate of Gaia BH1-like binaries (Rastello et al. 2023;Di Carlo et al. 2023;Tanikawa et al. 2023).However, the models considered so far (a) appear rather fine-tuned, involving multiple protagonists, stellar mergers, and a calibrated BH natal kick, (b) do not actually avoid a common envelope event, which the The BH binary merges faster for higher orbital eccentricities.The fact that equal-mass BH binaries with Pinner < 1.5 days merge within a Hubble time implies that hierarchical triple models for Gaia BH1 with inner BH binaries at these orbital periods (though consistent with the observed RV residuals from our two-body Keplerian fit) require some fine tuning.binary may not survive, and (c) predict binaries with periods P orb ∼ 100 days to form inefficiently compared to both shorter and longer periods (Rastello et al. 2023).Further work is required to explore the sensitivity of these predictions to initial cluster mass and primordial binary population.Another possibility is that Gaia BH1 formed through isolated binary evolution via a channel not captured in vanilla population synthesis models (e.g.Hirai & Mandel 2022).The large populations of wide white dwarf + main sequence (Yamaguchi et al. 2023;Shahaf et al. 2023) and neutron star + main sequence (El-Badry et al., in prep) binaries discovered by Gaia with orbits similar to that of Gaia BH1 lends some credibility to this possibility, but more work is required to explore it.5.6.Can further spectroscopic follow-up detect a inner BH binary in Gaia BH1?
As shown in Section 2.2.6, sensitivity to an inner binary improves significantly as the observing time baseline increases.This is primarily because inner binaries cause the outer orbit to precess, and the cumulative effects of precession grow over time.
Extending our analysis in Section 5.1, we find that acquiring 20 additional ESPRESSO RVs of Gaia BH1 near the August periastron passages in 2024 and 2025 would improve our constraint on the orbital period of an inner BBH to an upper limit of about 0.75 days.In this work, our best-fit hierarchical triple model for Gaia BH1 suggests the possibility of an inner BH binary with P inner = 0.9±0.1 days.Thus, future work on high-precision spectroscopic follow-up of Gaia BH1 would already be able to rule out (or confirm) our highestprobability hierarchical triple solution.While improving the upper limit of P inner ≲ 1.5 days to P inner ≲ 0.75 days may seem incremental, this would reduce the implied minimum inspiral time from the age of the Universe to only 2 Gyr (see Figure 13).

CONCLUSIONS
We have presented high-precision RV follow-up of Gaia BH1, a nearby binary system containing a Sun-like star orbiting a dark object with mass ∼ 9.3 M ⊙ .The system was recently discovered via Gaia astrometry and further characterized with low-precision RV follow-up.In their discovery paper, E23 interpreted the dark object as a likely dormant stellar-mass black hole (BH).They noted, however, that their data were also consistent with it being a close BH + BH or BH + NS binary, and that a hierarchical triple origin for the system could potentially alleviate some of the challenges associated with explaining its formation.
If Gaia BH1 hosts an inner black hole binary (BBH), the luminous star's RVs will display subtle deviations from a Keplerian orbit.High-precision RVs can thus put the hierarchical triple model to the test.In parallel with our work, the feasibility of such a test was recently explored by Hayashi et al. (2023) using idealized simulations.Here, we report new RV data and simulations matched to that data.Our main conclusions are as follows: • Sensitivity to inner BH binaries: We explore the parameter space of expected RV signatures of a hierarchical triple in the non-circular, non-coplanar regime using the N-body integrator REBOUND (Figure 1).We find that (a) increasing the eccentricity of either the inner or outer orbit increases the amplitude of observed RV residuals, (b) increasing the inner binary mass ratio decreases the amplitude of observed RV residuals, and (c) varying the mutual inclination non-monotonically changes the amplitude of observed RV residuals (see Figure 6).We determine that with optimal observing cadence we can reasonably detect an inner BH binary in the RV residuals of a Keplerian fit for inner periods ≳ 1 day (see Figure 2), and that residuals in realistic orbital configurations are always larger than in the circular and coplanar case (Figure 3).
We show that with observations over one period of the outer orbit, the short-timescale RV variations induced by an inner BH binary are of comparable magnitude to the long-timescale RV variations due to precession (see Figure 4).Extending the observational baseline to multiple orbital cycles helps detect the cumulative effect of precession in the outer star's orbit (Figure 5).
We find that it is possible to distinguish between perturbations due to an inner BH binary and an exoplanet orbiting the outer star with precise RV measurements.
A clear difference between the two cases is that the amplitude of RV variations does not increase significantly at the periastron of the outer star's orbit in the exoplanetary case (last panel of Figure 3).
• High-precision RV data: We obtained 40 highprecision RV measurements (∼ 3-5 m s −1 uncertainties) using ESPRESSO over one orbit of the Sunlike star.We supplement these data with 75 (mostly new) medium-precision RVs (∼ 30-100 m s −1 uncertainties) measured with FEROS, TRES, and HIRES over three orbits (Figure 8 and Table 3).We concentrated the high-precision RVs near a periastron passage, where the expected perturbations due to an inner binary are largest.We combine the spectroscopic and astrometric data to tightly constrain the parameters of a two-body Keplerian model and update the constraints from E23 (Table 1 and Figure 9).
• Two-body and three-body fits: The observed RVs are well-fit by a Keplerian two-body orbit.From our observed value of reduced χ 2 = 1.69, we set an upper limit of P inner ≲ 1.5 days on the period of any inner BBH (see Figures 10 and 11).
To assess whether a model with an inner BH binary can better explain the data than a two-body Keplerian model, we use REBOUND to fit a three-body model.We derive the best-fit orbital parameters for both the outer star and the inner BH binary in this case (see Table 2).We find that the putative inner BH binary would have an orbital period of 0.9±0.1 days, consistent with our limits from the Keplerian fit.While the hierarchical triple model can provide a marginally better fit to the observed RVs than the two-body Keplerian model (Figure 12), the improvement is not sufficient to justify the model's higher complexity, and we conclude that the improvement in log posterior probability is most likely due to overfitting the noise.
Acquiring ∼ 20 additional precise RV observations of Gaia BH1 around periastron in 2024 and 2025 would improve our constraint on the orbital period of the inner BBH to an upper limit of about 0.75 days.
• Merger timescale: While it is possible that an inner BH binary with P inner < 1.5 days exists, such a BH binary would be expected to merge in less than a Hubble time (see Figure 13).Hiding a BH binary in the period range not ruled out by RVs would thus require finetuning.Much of the short-period inner binary regime not ruled out by our RVs will eventually be probed by LISA (Section 5.3).
• Implications for the formation of Gaia BH1: Our results imply that Gaia BH1 is unlikely to currently host an inner BH binary: a single 9.3 M ⊙ BH is the most likely companion.It is also unlikely that it previously hosted such a binary that has now merged, because the kick due to anisotropic gravitational wave emis-sion would likely have unbound the outer orbit.We have not ruled out a triple scenario in which the inner binary merged while its components were still on the main sequence.More work is required to investigate alternative formation channels, such as dynamical assembly through exchange interactions, nonstandard treatments of common envelope evolution, and efficient wind mass loss to avoid a red supergiant phase in the BH progenitor.

Figure 1 .
Figure1.Orbital configuration (left) and predicted RVs (right) of the hierarchical triple scenario we seek to test.RVs are measured for a Sunlike star orbiting an inner BH + BH binary.In the general case, the inner and outer orbits are both eccentric and non-coplanar.The parameters of the outer orbit are fixed to those inferred for Gaia BH1 by E23.Here we assume an equal-mass inner binary with period Pinner = 6 days and eccentricity einner = 0.1.While the RVs of the Sun-like star are nearly consistent with a Keplerian orbit (upper right), subtraction of the best-fit Keplerian orbit reveals significant residuals (lower right).At Pinner = 6 days, the amplitude of the RV residuals is ∼250 m s −1 near periastron.

Figure 2 .
Figure 2. Left: RV residuals for outer star (bottom panel) after fitting a RadVel Keplerian orbit to 50 noisy RVs sampled from a REBOUND simulation using realistic cadence (top panel).The orbital parameters of the outer star are fixed to those derived for Gaia BH1 by E23, and the orbital parameters of the inner BBH are the same as those used for Figure 1.Right: Reduced χ 2 of best-fit Keplerian model as a function of period of inner binary.The reduced χ 2 is close to unity at Pinner ≲ 1 day, indicating a good fit.It rises rapidly at longer inner periods, indicating a detectable RV perturbation.
this case, Pinner refers to the period of the exoplanetary orbit.

Figure 3 .
Figure3.Left: RV residuals for a variety of configurations of the hierarchical triple at Pinner = 6 days.Right: Power laws fitted to variation of amplitude of RV residuals with Pinner.We recover the theoretical expectation of a power law slope of 2.33 in the case where the orbits are coplanar and both are circular (Equation1).In the other cases, the power law slope is lower, but the amplitude of the RV residuals is higher at all detectable values of Pinner.Bottom panel shows a case where there is a single BH, but a Jupiter-mass planet orbits the star.

Figure 4 .
Figure 4. Variation of RV residual curve shapes with Pinner.Here we fix the orbital parameters of the outer star to those determined by E23 and make the inner orbit coplanar and circular.At short Pinner, the residuals are dominated by a long-term undulation due to precession of the outer orbit.A lower-amplitude "wobble" on half the period of the inner binary is also present.Both the long-and short-timescale residuals become larger with increasing Pinner, but the short-timescale residuals grow more quickly, such that short-and long-timescale residuals have comparable amplitude at Pinner = 9 days.

Figure 5 .
Figure5.Left: RV residuals at Pinner = 6 days for observational baselines of increasing duration in the case where we fix the orbital parameters of the outer star to those determined by E23 and constrain the inner orbit to be coplanar and circular.Right: Power laws fitted to variation of amplitude of RV residuals with Pinner for each observational baseline.Covering more than one orbital cycle of the outer star can increase the observed RV residual amplitude significantly, because the cumulative effects of precession grow with the observing baseline.

Figure 6 .
Figure6.Evolution of amplitude of RV residuals as a function of Pinner with various orbital parameters of the inner BBH.We use the dynamical stability criterion ofAarseth & Mardling (2001) to restrict the power law curves from entering unstable regions in parameter space.In general, increasing the eccentricity of either orbit increases the amplitude of the observed residuals, while increasing the BBH mass ratio decreases the amplitude of the observed residuals.In addition, the residual amplitude varies non-monotonically with the orbital inclination of the inner BBH.The dashed horizontal line marks a residual amplitude of 5 m s −1 , roughly the sensitivity of our observations.This suggests that our observations are sensitive to essentially all plausible inner binaries with Pinner significantly above 1 day, except in the (physically dubious) case of an inner BBH mass ratio of 100.

Figure 7 .
Figure7.Typical ESPRESSO spectra of Gaia BH1 from the two epochs with extremal RVs.Both spectra have SNR ∼30 per pixel, leading to ∼ 4 m s −1 RV uncertainties.This cutout is centered on the Fraunhofer Na "D" lines.Both photospheric and interstellar lines are present.The photospheric lines trace the RV of the Sun-like star.The interstellar lines (labeled "ISM") are very sharp and stable between epochs, highlighting the high resolution of the data (R ≈ 140, 000) and the stability of the wavelength solution.

Figure 8 .
Figure8.Best-fit RV curve for Gaia BH1 based on Gaia DR3 astrometry and our updated spectroscopic measurements, assuming a two-body Keplerian orbit.The top panel shows the observed data points over 50 RV curves randomly sampled from the posterior.The bottom panels show the residuals relative to the MAP solution plotted at three levels of precision, with the last panel focusing on the latest ESPRESSO measurements.At all levels of precision, the RVs are generally consistent with the best-fit model at the 1-2σ level.

Figure 10 .
Figure10.Reduced χ 2 value of the best-fit Keplerian model as a function of Pinner for simulated data assuming the exact observing cadence and uncertainties of our measured RVs.Different panels show different orientations of the inner binary, and lines within each panel show different inner binary eccentricities.The reduced χ 2 value generally increases with increasing einner.At very high einner, the non-uniform sampling of our RVs gives rise to oscillations in the reduced χ 2 value.Based on our best-fit model's reduced χ 2 value of 1.69, we can rule out most inner binaries with Pinner ≳ 1.5 days.Fine-tuned orientations exist that could hide inner binaries with orbital periods up to ∼ 3 days (e.g.lower right); we explore these in more detail in Figure11.

Figure 12 .
Figure 12.Top: Points with error bars show the observed ESPRESSO RV residuals relative to the best-fit two-body Keplerian model.The continuous curve represents the difference between the best-fit hierarchical triple model and that Keplerian model.The residuals are consistent with the continuous curve to within 1-2σ.Bottom: Observed ESPRESSO RV residuals relative to the best-fit two-body Keplerian model (purple; same as top panel) and the best-fit hierarchical triple model (red).The scatter of the residuals for the hierarchical triple model is somewhat smaller than in the Keplerian case, with an improvement in log posterior probability of ∼12.This reflects the fact that the more flexible model allows for a long-term trend in the residuals due to precession.However, the improvement is small compared to the RV uncertainties and appears consistent with overfitting.

Figure 13 .
Figure13.Inspiral times of a equal mass ratio black hole binary (MBH,1 = MBH,2 = 4.65 M⊙) as a function of orbital period for various eccentricities.The BH binary merges faster for higher orbital eccentricities.The fact that equal-mass BH binaries with Pinner < 1.5 days merge within a Hubble time implies that hierarchical triple models for Gaia BH1 with inner BH binaries at these orbital periods (though consistent with the observed RV residuals from our two-body Keplerian fit) require some fine tuning.

Figure 14 .
Figure14.RV residual amplitudes predicted for the Sun-like star when the star + BH binary is orbited by a distant (outer) tertiary.Blue and orange curves correspond to a 1 M⊙ white dwarf and a 0.2 M⊙ M dwarf, respectively, with high-precision RVs obtained over one orbit.Different curves show different random orbital phases.Based on our detectability threshold of an RV residual amplitude of ∼ 5 m s −1 , we can rule out a tertiary white dwarf at all orbital periods less than about 2700 days, and a tertiary M dwarf at periods less than about 1500 days.

Table 1 .
Physical parameters and 1σ uncertainties derived for the orbit of the outer Sun-like star in Gaia BH1 assuming a two-body Keplerian model.We present both the values measured by E23 and the updated values inferred in this work (see Section 4).

Table 2 .
Best-fit parameters for the orbit of the outer Sun-like star in Gaia BH1 assuming a hierarchical triple model.Error bars on the median solution are derived from the 16th and 84th percentiles.