# ESPRESSO observations of Gaia BH1: high-precision orbital constraints and no
evidence for an inner binary

Pranav Nagarajan Department of Astronomy, California Institute of Technology,
1200 E. California Blvd., Pasadena, CA 91125, USA  Kareem El-Badry Department
of Astronomy, California Institute of Technology, 1200 E. California Blvd.,
Pasadena, CA 91125, USA  Amaury H.M.J. Triaud School of Physics and Astronomy,
University of Birmingham, Edgbaston, Birmingham B15 2TT, UK  Thomas A.
Baycroft School of Physics and Astronomy, University of Birmingham, Edgbaston,
Birmingham B15 2TT, UK  David Latham Center for Astrophysics $|$ Harvard &
Smithsonian, Harvard University, 60 Garden Street, Cambridge, MA 02138, USA
Allyson Bieryla Center for Astrophysics $|$ Harvard & Smithsonian, Harvard
University, 60 Garden Street, Cambridge, MA 02138, USA  Lars A. Buchhave DTU
Space, Technical University of Denmark, Elektrovej 328, DK-2800 Kgs. Lyngby,
Denmark  Hans-Walter Rix Max Planck Institute for Astronomy, Königstuhl 17,
D-69117, Heidelberg, Germany  Eliot Quataert Department of Astrophysical
Sciences, Princeton University, 4 Ivy Lane, Princeton, NJ 08544, USA  Andrew
Howard Department of Astronomy, California Institute of Technology, 1200 E.
California Blvd., Pasadena, CA 91125, USA  Howard Isaacson Department of
Astronomy, University of California, Berkeley, Berkeley, CA 94720, USA
Melissa J. Hobson Max Planck Institute for Astronomy, Königstuhl 17, D-69117,
Heidelberg, Germany

###### Abstract

We present high-precision radial velocity (RV) 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 are 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_{\text{inner}}\gtrsim 1.5$ days would have produced
larger deviations from a Keplerian orbit than observed. Binaries with
$P_{\text{inner}}\lesssim 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 $f\left(M_{\text{BH}}\right)=3.9358\pm 0.0002\,M_{\odot}$ — about
$7000\sigma$ above the $\sim 2.5\,M_{\odot}$ maximum neutron star mass.
Including the inclination constraint from Gaia astrometry, this implies a BH
mass of $M_{\text{BH}}=9.27\pm 0.10~{}M_{\odot}$.

Black Holes (162) — Multiple Stars (1081)

††facilities: ESPRESSO (VLT)††software: astropy (Astropy Collaboration et al.,
2013, 2018), REBOUND (Rein et al., 2023), RadVel (Fulton et al., 2020), emcee
(Foreman-Mackey et al., 2013), dynesty (Koposov et al., 2023)

##  1 Introduction

The Milky Way is expected to contain $\sim 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 $\sim$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 likely a rare outcome of binary
evolution, with population models estimating that only $\sim 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) binarity-induced wobble of $\sim
2$ billion stars. The mission’s 3rd data release (“DR3”) represents a factor
of $\sim 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 (Gaia
Collaboration et al., 2023a, b).

Recently, El-Badry et al. (2023a, hereafter E23) used data from Gaia DR3 to
identify Gaia BH1, a Sun-like star orbiting a dark companion with inferred
mass $M_{2}=9.62\pm 0.18~{}M_{\odot}$. At a distance of only 480 pc, the
system is the nearest known BH by a factor of $\sim$3\. While E23 identified
the companion as a likely dormant stellar-mass BH, they noted that the data
were also consistent with the unseen object being a close binary containing
two BHs with total mass $\approx 9.6\,M_{\odot}$. 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\approx 1.4$ AU, is significantly smaller than
the predicted maximum radius of a solar-metallicity progenitor of a $\sim
10~{}M_{\odot}$ 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).

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.

##  2 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.

###  2.1 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:

$K_{\text{short}}=\frac{m_{1}m_{2}}{m_{12}^{2}}\left(\frac{m_{12}}{m_{123}}\right)^{5/3}\left(\frac{P_{\text{inner}}}{P_{\text{outer}}}\right)^{7/3}\left(\frac{2\pi
Gm_{123}}{P_{\text{outer}}}\right)^{1/3}$ (1)

where $m_{12}\equiv m_{1}+m_{2}$, $m_{123}\equiv m_{12}+m_{*}$, and
$P_{\text{inner}}$ and $P_{\text{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 $\left(P_{\text{inner}}/P_{\text{outer}}\right)^{7/3}$
and have a dominant period $\approx P_{\text{inner}}/2$ (Morais & Correia,
2008; Hayashi et al., 2020).

###  2.2 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 the
outer orbit is not circular (it has eccentricity $e_{\text{outer}}\approx
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 non-circular, non-coplanar
regime using REBOUND, a flexible N-body 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.

Table 1: Physical parameters and 1$\sigma$ 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).  Parameter | Description | Constraint (E23) | Median Constraint (this work) | MAP Solution (this work)
---|---|---|---|---
(1) | (2) | (3) | (4) | (5)
$P$ | Orbital Period |  $185.59\pm 0.05$ days |  $185.387\pm 0.003$ days |  $185.387$ days
$e$ | Eccentricity | $0.451\pm 0.005$ | $0.43230\pm 0.00002$ | $0.43230$
$i$ | Inclination | $(126.6\pm 0.4)^{\circ}$ | $(126.8\pm 0.2)^{\circ}$ | $126.832^{\circ}$
$\Omega$ | Longitude of Ascending Node | $(97.8\pm 1.0)^{\circ}$ | $(97.0\pm 0.7)^{\circ}$ | $96.9^{\circ}$
$\omega$ | Argument of Periastron | $(12.8\pm 1.1)^{\circ}$ | $(16.509\pm 0.003)^{\circ}$ | $16.510^{\circ}$
$T_{p}$ | Periastron Time (JD - 2457389) | $-1.1\pm 0.7$ | $2.07\pm 0.04$ | $2.07$
$M_{\text{BH}}$ | Black Hole Mass | $9.62\pm 0.18~{}M_{\odot}$ | $9.27\pm 0.10~{}M_{\odot}$ | $9.294~{}M_{\odot}$
$M_{\text{star}}$ | Luminous Star Mass | $0.93\pm 0.05~{}M_{\odot}$ | $0.93\pm 0.05~{}M_{\odot}$ | $0.933~{}M_{\odot}$
$\gamma$ | Center-of-Mass Radial Velocity |  $46.6\pm 0.6$ km s-1 |  $48.379\pm 0.001$ km s-1 |  $48.379$ km s-1
$\beta_{H}$ | HIRES Radial Velocity Offset | N/A |  $0.09\pm 0.03$ km s-1 |  $0.074$ km s-1
$\beta_{F}$ | FEROS Radial Velocity Offset | N/A |  $0.23\pm 0.01$ km s-1 |  $0.227$ km s-1
$\beta_{T}$ | TRES Radial Velocity Offset | N/A |  $0.51\pm 0.01$ km s-1 |  $0.516$ km s-1

####  2.2.1 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,\omega,K)$ basis. We use a generic Keplerian RV model and Gaussian
RV likelihood, and fit in the alternative
$(P,T_{c},\sqrt{e}\cos{\omega},\sqrt{e}\sin{\omega},K)$ basis, which imposes
flat priors on all orbital elements and helps speed MCMC convergence (Fulton
et al., 2018). Finally, we use the Powell minimization method with a tolerance
of $10^{-5}$ to derive the best-fit 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 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_{\text{inner}}/2$) variation in the RV
residual curve, which we attribute to precession in the orbital parameters of
the outer star.

Figure 1: Orbital configuration (left) and predicted RVs (right) of the
hierarchical triple scenario we seek to test. RVs are measured for a Sun-like
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 $P_{\text{inner}}=6$ days and eccentricity
$e_{\text{inner}}=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
$P_{\text{inner}}=6$ days, the amplitude of the RV residuals is $\sim$250 m
s-1 near periastron.

####  2.2.2 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, we sample the RV curve of the outer star at
$\sim$50 epochs spaced every few days, with increased cadence (i.e. every day)
near periastron.

We show the resulting RV residual curve for $P_{\text{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 $\chi^{2}_{\text{red}}\approx 666$. In the right panel, we show the
reduced $\chi^{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_{\text{inner}}\ll
1$ d, but rises steeply at $P_{\text{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_{\text{inner}}\gtrsim 1$ day.

####  2.2.3 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_{\text{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_{\text{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 $\alpha$)
in each panel.

Each row shows a different orbital configuration. In the top row, we simulate
the coplanar circular case, 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_{\text{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. This is primarily because the eccentric
outer orbit allows the inner BBH to get closer to the outer star at
periastron.

####  2.2.4 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_{\odot}$ 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
($\alpha=-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.

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 $\chi^{2}$
of best-fit Keplerian model as a function of period of inner binary. The
reduced $\chi^{2}$ is close to unity at $P_{\text{inner}}\lesssim 1$ day,
indicating a good fit. It rises rapidly at longer inner periods, indicating a
detectable RV perturbation. Figure 3: Left: RV residuals for a variety of
configurations of the hierarchical triple at $P_{\text{inner}}=6$ days. Right:
Power laws fitted to variation of amplitude of RV residuals with
$P_{\text{inner}}$. We recover the theoretical expectation of a power law
slope of $2.33$ in the coplanar circular case (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 $P_{\text{inner}}$. Bottom panel shows a
case where there is a single BH, but a Jupiter-mass planet orbits the star.

####  2.2.5 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_{\text{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 long-timescale RV residuals increases with the orbit of the
inner BBH, varying from about $50$ m s-1 at $P_{\text{inner}}=3$ days to $350$
m s-1 at $P_{\text{inner}}=9$ days. At small $P_{\text{inner}}$, the long-term
perturbation due to precession is larger than the short-timescale variations.
At longer $P_{\text{inner}}$, the two amplitudes are similar.

Figure 4: Variation of RV residual curve shapes with $P_{\text{inner}}$. 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 $P_{\text{inner}}$, 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 $P_{\text{inner}}$, but the short-timescale residuals grow
more quickly, such that short- and long-timescale residuals have comparable
amplitude at $P_{\text{inner}}=9$ days.

####  2.2.6 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.

Figure 5: Left: RV residuals at $P_{\text{inner}}=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 $P_{\text{inner}}$ 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.

####  2.2.7 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_{\text{inner}}\rightarrow\infty$, 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_{\text{inner}}$. In
almost all cases, the observed amplitude of short-term non-Keplerian RV
modulations is above our detectability threshold of $\approx 5$ m s-1 (see
Section 3) for $P_{\text{inner}}>1$ day. The only exceptions are the models
with $e_{\text{outer}}=0$ (which is ruled out since we know
$e_{\text{outer}}\sim 0.45$) or $q_{\text{inner}}=100$ (which is
astrophysically unlikely).

Figure 6: Evolution of amplitude of RV residuals as a function of
$P_{\text{inner}}$ with various orbital parameters of the inner BBH. We use
the dynamical stability criterion of 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 $P_{\text{inner}}$
significantly above 1 day, except in the (physically dubious) case of an inner
BBH mass ratio of $100$.

####  2.2.8 Dynamical stability

For sufficiently large inner orbits and high inner or outer eccentricities,
the triple is expected to become unstable on a short timescale. An approximate
dynamical stability criterion was derived by Aarseth & Mardling (2001):

$\frac{r_{p,\text{outer}}}{a_{\text{inner}}}>2.8\left[\left(1+\frac{m_{*}}{m_{1}+m_{2}}\right)\frac{(1+e_{\text{outer}})}{\sqrt{1-e_{\text{outer}}}}\right]^{2/5}$
(2)

where $a_{\text{inner}}$ is the semi-major axis of the inner binary and
$e_{\text{outer}}$ and $r_{p,\text{outer}}$ are the eccentricity and orbital
separation at periastron of the outer orbit. We use this criterion in Figure 6
to restrict the power law curves from entering regions where the hierarchical
triple system is unstable. Specifically, this is relevant for large
$P_{\text{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.

###  3.1 FEROS

We observed Gaia BH1 53 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, P110.A-9014, P111.A-9003, and
P112.2650). The first several observations used $2\times 2$ binning to reduce
readout noise at the expense of spectral resolution; the remainder used
$1\times 1$ binning. The resulting spectra have resolution $R\approx 40,000$
($2\times 2$ binning) and $R\approx 50,000$ ($1\times 1$ binning). Most of our
observations used 1800s exposures. The typical SNR per pixel at 5800 Å is 15.

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., $\sigma_{\text{RV}}\approx{\rm
std}\left({\text{RVs}}\right)/\sqrt{15}$. We used a Kurucz model spectrum
(Kurucz, 1979, 1993) with $T_{\text{eff}}=5750$ K, $\log g=4.5$, and
$\rm[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 $\approx 70$ m
s-1.

###  3.2 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\approx 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 cross-correlating 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., $\sigma_{\text{RV}}={\rm std}\left(\text{RVs}\right)/\sqrt{31}$. We used
the same Kurucz template from the BOSZ library that we used for the FEROS data
($T_{\text{eff}}=5750$ K, $\log g=4.5$, $\rm[Fe/H]=-0.25$). The median
uncertainty of the TRES RVs is $\approx 50$ m s-1.

###  3.3 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.

###  3.4 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\times 1$ binning, yielding a typical spectral resolution $R=140,000$,
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 all orders of the
extracted spectra with a G2 mask and then fits the resulting cross-correlation
function (CCF) with a Gaussian. 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 $\sim 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; uncertainty in the star’s micro/macroturbulent
velocities limits an even more precise measurement.

Figure 7: Typical ESPRESSO spectra of Gaia BH1 from the two epochs with
extremal RVs. Both spectra have SNR $\sim$30 per pixel, leading to $\sim 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\approx 140,000$) and the stability of the wavelength solution.

###  3.5 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 $\sigma_{\text{RV}}\gtrsim 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 principle 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_{\odot}$ and a standard deviation of
$0.05~{}M_{\odot}$. For all other parameters, we used uniform priors. The log-
likelihood assumes Gaussian uncertainties:

$\ln
L_{\text{RVs}}=-\frac{1}{2}\sum_{t}\left(\frac{\text{RV}_{\text{obs},t}-\text{RV}_{\text{pred},t}}{\sigma_{\text{RV},t}}\right)^{2}$
(3)

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 $\mu_{\text{ast}}$, their covariance matrix
$\Sigma_{\text{ast}}$, and the vector of predicted astrometric parameters for
a given likelihood call, $\theta_{\text{ast}}$, the astrometric log-likelihood
can be expressed as:

$\ln
L_{\text{ast}}=-\frac{1}{2}(\theta_{\text{ast}}-\mu_{\text{ast}})^{T}\Sigma_{\text{ast}}(\theta_{\text{ast}}-\mu_{\text{ast}}).$
(4)

We calculate $\theta_{\text{ast}}$ under the assumption that the companion is
dark. The total log-likelihood function includes both the astrometric and RV
terms:

$\ln L=\ln L_{\text{ast}}+\ln L_{\text{RVs}}$ (5)

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 $\theta_{\text{ast}}$ and $\mu_{\text{ast}}$, and remove the
corresponding rows and columns from $\Sigma_{\text{ast}}$. In predicting the
Thiele-Innes elements, we assume a parallax of $\varpi=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.

###  4.1 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$\sigma$ (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 of at most $\approx 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$,
$\Omega$, and $M_{\text{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_{\text{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_{\text{BH}}=9.27\pm
0.10~{}M_{\odot}$, corresponding to an RV semi-amplitude
$K_{\text{star}}=65.3785\pm 0.0009$ km s-1 and a spectroscopic mass function
$f\left(M_{\text{BH}}\right)=3.9358\pm 0.0002~{}M_{\odot}$. 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\sigma$, $3.7\sigma$, $3.4\sigma$, $4.5\sigma$ and $3.0\sigma$
respectively. The remaining orbital parameters are consistent with each other
to within 2$\sigma$.

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\sigma$ 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 $\chi^{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 $\gamma$ 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 $\sim 100$
m s-1 (see Lindegren & Dravins 2003 for a discussion).

Figure 8: 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\sigma$ level.
Figure 9: Corner plot for orbital parameters of Gaia BH1 from MCMC sampling
based on Gaia DR3 astrometry and our updated spectroscopic measurements,
assuming a two-body Keplerian orbit. The diagonal entries display the marginal
distribution of each parameter, while each of the other panels displays a
joint distribution. The new constraints (black) are plotted over the
constraints from E23 (blue). Our new data allow us to place much tighter
constraints on the orbital parameters than those obtained by E23. The best-fit
parameters are similar to those of E23 in an absolute sense, but several
parameters are in tension at the $3$–$4\sigma$ level (see Table 1).

###  4.2 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_{\text{inner}}\to 0$ or $q_{\text{inner}}\to 0$, we recover the two-body
Keplerian solution. For $0.5\lesssim P_{\text{inner}}/\text{day}\lesssim 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_{\text{inner}}\gtrsim 1.5$ days, the inner binary perturbs the outer orbit
too much, leading to a poorer fit than the two-body 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 5). 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_{\text{inner}}>0.5$ days because likelihood calls became prohibitively
expensive in the limit of $P_{\text{inner}}\to 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 $\lesssim 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 highest-probability
solution achieved across all such runs as the MAP 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_{\text{inner}}\lesssim 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 two-body 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.

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.  Parameter | Description | Median Derived Value (this work) | MAP Derived Value (this work)
---|---|---|---
(1) | (2) | (3) | (4)
$P_{\text{outer}}$ | Period of Outer Orbit |  $185.45_{-0.02}^{+0.01}$ days |  $185.4594$ days
$e_{\text{outer}}$ | Eccentricity of Outer Orbit | $0.43243_{-0.00005}^{+0.00005}$ | $0.43245$
$i_{\text{outer}}$ | Inclination of Outer Orbit | $(126.7_{-0.2}^{+0.2})^{\circ}$ | $126.653^{\circ}$
$\Omega_{\text{outer}}$ | Longitude of Ascending Node of Outer Orbit | $(98.2_{-0.9}^{+0.9})^{\circ}$ | $98.2^{\circ}$
$\omega_{\text{outer}}$ | Argument of Periastron of Outer Orbit | $(16.518_{-0.007}^{+0.008})^{\circ}$ | $16.521^{\circ}$
$T_{p,\text{outer}}$ | Periastron Time (JD - 2460000) of Outer Orbit | $-13.506_{-0.005}^{+0.005}$ | $-13.505$
$P_{\text{inner}}$ | Period of Inner Orbit |  $0.9_{-0.1}^{+0.1}$ days |  $0.98$ days
$e_{\text{inner}}$ | Eccentricity of Inner Orbit | $0.19_{-0.09}^{+0.08}$ | $0.205$
$i_{\text{inner}}$ | Inclination of Inner Orbit | $(108.6_{-5.0}^{+4.4})^{\circ}$ | $108.8^{\circ}$
$\Omega_{\text{inner}}$ | Longitude of Ascending Node of Inner Orbit | $(171.0_{-4.9}^{+3.4})^{\circ}$ | $172.4^{\circ}$
$\omega_{\text{inner}}$ | Argument of Periastron of Inner Orbit | $(8.6_{-6.2}^{+11.6})^{\circ}$ | $6.5^{\circ}$
$T_{p,\text{inner}}$ | Periastron Time (JD - 2460000) of Inner Orbit | $-0.7_{-0.2}^{+0.1}$ | $-0.70$
$q_{\text{inner}}$ | Mass Ratio of Black Hole Binary | $0.8_{-0.2}^{+0.1}$ | $0.80$
$M_{\text{BH, tot}}$ | Total Mass of Black Hole Binary | $9.24_{-0.08}^{+0.08}~{}M_{\odot}$ | $9.235~{}M_{\odot}$
$M_{\text{star}}$ | Mass of Luminous Star | $0.93_{-0.05}^{+0.05}~{}M_{\odot}$ | $0.9337~{}M_{\odot}$
$\gamma$ | Center-of-Mass Radial Velocity |  $48.375_{-0.002}^{+0.002}$ km s-1 |  $48.375$ km s-1
$\beta_{H}$ | HIRES Radial Velocity Offset |  $0.10_{-0.04}^{+0.03}$ km s-1 |  $0.101$ km s-1
$\beta_{F}$ | FEROS Radial Velocity Offset |  $0.23_{-0.01}^{+0.01}$ km s-1 |  $0.224$ km s-1
$\beta_{T}$ | TRES Radial Velocity Offset |  $0.52_{-0.01}^{+0.02}$ km s-1 |  $0.525$ km s-1

##  5 Discussion

###  5.1 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 $\chi^{2}$ value of a RadVel Keplerian fit to simulated
RVs for a hierarchical triple as a function of $P_{\text{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
$\chi^{2}_{\text{red}}=1.69$ from the updated Keplerian fit for Gaia BH1.

For most of the orientations and eccentricities shown in Figure 10, the
predicted reduced $\chi^{2}$ rises above the observed value for
$P_{\text{inner}}\gtrsim 1$ day, indicating that we can rule out an inner BBH
with $P_{\text{inner}}\gtrsim 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 $\chi^{2}_{\text{red}}=1.69$ for $P_{\text{inner}}\gtrsim 3$
days.

Figure 10: Reduced $\chi^{2}$ value of the best-fit Keplerian model as a
function of $P_{\text{inner}}$ 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 $\chi^{2}$ value generally increases
with increasing $e_{\text{inner}}$. At very high $e_{\text{inner}}$, the non-
uniform sampling of our RVs gives rise to oscillations in the reduced
$\chi^{2}$ value. Based on our best-fit model’s reduced $\chi^{2}$ value of
$1.69$, we can rule out most inner binaries with $P_{\text{inner}}\gtrsim 1.5$
days. Fine-tuned orientations exist that could hide inner binaries with
orbital periods up to $\sim 3$ days (e.g. lower right); we explore these in
more detail in Figure 11.

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 $\chi^{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_{\text{inner}}=1$ and $P_{\text{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_{\text{inner}}=0$ and $e_{\text{inner}}=0.6$. We also label
the orientations corresponding to the “best” and “worst” cases (in terms of
detectability, for $e_{\text{inner}}=0$) displayed in Figure 10.

The $\chi^{2}_{\text{red}}$ 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_{\text{inner}}\sim 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_{\text{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 $\chi^{2}$ value in each case,
assuming our observed cadence and RV uncertainties. We assume a uniform
$e_{\text{inner}}$ distribution and random orientations. Then, at each value
of $P_{\text{inner}}$, we compute the fraction of inner binaries that would
have been detected at our threshold of $\chi^{2}_{\text{red}}=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_{\text{inner}}\gtrsim 1.5$ days. This is consistent with the best-fit
derived value of $P_{\text{inner}}=0.9\pm 0.1$ days from our hierarchical
triple model.

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 ($\gtrsim 95\%$) inner
binaries with $P_{\text{inner}}>1.5$ days. Right: Dependence of the best-fit
Keplerian model’s reduced $\chi^{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 $P_{\text{inner}}=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
$\chi^{2}$ landscape is complex, but only a few fine-tuned inner binary
orientations allow $\chi^{2}_{\text{red}}<2$.

###  5.2 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$\sigma$.

We compare observed RV residuals for the best-fit two-body Keplerian and
hierarchical triple models in the bottom panel of Figure 12, once again
displaying only the high-precision 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 residuals is small, and the
improvement is not significant with respect to the RV uncertainties.

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$\sigma$. 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 $\sim$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.

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:

$\text{BIC}=k\ln{n}-2\ln{\hat{L}}$ (6)

where $k$ is the number of model parameters, $n$ is the number of data points,
and $\hat{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
models 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.

###  5.3 BH binary inspiral time and detectability with LISA

While an inner BH binary with $P_{\text{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_{\text{inner}}<1.5$ days, the BH binary merges in less than
the age of the Universe, even for cases with low orbital eccentricities.

Figure 13: Inspiral times of a equal mass ratio black hole binary
($M_{\text{BH},1}=M_{\text{BH},2}=4.65~{}M_{\odot}$) 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
$P_{\text{inner}}<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.

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_{\text{inner}}\lesssim 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_{\text{inner}}=0.5$ and $2.5$ days for $e_{\text{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

The 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_{\odot}$ white dwarf tertiary and a $0.2~{}M_{\odot}$ M dwarf
tertiary, both of which could have evaded spectroscopic detection. For the
tertiary star, we adopt an orbital 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 $\sim 5$ m
s-1, we find that we would detect a $1~{}M_{\odot}$ 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_{\odot}$ 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 $\gtrsim 0.2~{}M_{\odot}$ 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 $\sim 500$ AU, or periods
$<10^{6}$ days. A stronger limit of $\approx 0.1$ arcsec, corresponding to
periods $\lesssim 10^{4.5}$ days, could be set with speckle interferometry
(e.g. Howell & Furlan, 2022).

Figure 14: 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_{\odot}$ white dwarf and a $0.2~{}M_{\odot}$ 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 $\sim$ 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.

###  5.5 Implications for the formation history of Gaia BH1

We have shown that the $\sim 9.3~{}M_{\odot}$ dark object in the Gaia BH1
system is unlikely to be a binary with period longer than $\sim 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 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_{\odot}$
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.

Another formation channel that has recently been explored is dynamical
assembly in a star cluster. The G star’s near-solar 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 high-
mass 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 binary may not survive,
and (c) predict binaries with periods $P_{\text{orb}}\sim 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_{\text{inner}}=0.9\pm 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 highest-probability hierarchical triple solution. While improving
the upper limit of $P_{\text{inner}}\lesssim 1.5$ days to
$P_{\text{inner}}\lesssim 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).

##  6 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 $\sim
9.3\,M_{\odot}$. 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 $\gtrsim$ 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 high-precision RV measurements ($\sim
3$–$5$ m s-1 uncertainties) using ESPRESSO over one orbit of the Sun-like
star. We supplement these data with 75 (mostly new) medium-precision RVs
($\sim 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 $\chi^{2}_{\text{red}}=1.69$, we
set an upper limit of $P_{\text{inner}}\lesssim 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\pm 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 $\sim 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_{\text{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 fine-tuning. 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_{\odot}$ 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 emission 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.

We thank Toshinori Hayashi, Jim Fuller, and Dave Charbonneau for helpful
discussions. PN and KE were supported in part by NSF grant AST-2307232. AT and
TB are supported by a grant from the European Research Council (ERC) under the
European Union’s Horizon 2020 research and innovation program (grant agreement
number 803193/BEBOP). HWR acknowledges the European Research Council for
support from the ERC Advanced Grant ERC-2021-ADG-101054731. , LEGWORK (Wagg et
al., 2022a)

##  Appendix A Radial Velocity Data

We present all radial velocity measurements for Gaia BH1 used in this work in
chronological order in Table 3 below. The MagE, GMOS, XSHOOTER, and ESI data
are taken from E23. The remainder of the data (from the HIRES, FEROS, TRES,
and ESPRESSO instruments) are either updated in or new to this work.

Table 3: All RV measurements obtained via spectroscopic follow-up of Gaia BH1 using the MagE, GMOS, XSHOOTER, ESI, HIRES, FEROS, TRES, and ESPRESSO instruments.  HJD | Radial Velocity (km s-1) | Instrument
---|---|---
(1) | (2) | (3)
2459767.6226 | $63.8\pm 3$ | MagE
2459791.9186 | $131.90\pm 0.1$ | HIRES
2459795.6461 | $141.4\pm 3$ | MagE
2459796.4995 | $142.7\pm 3$ | MagE
2459798.8399 | $140.6\pm 4$ | GMOS
2459805.5101 | $127.7\pm 1.0$ | XSHOOTER
2459808.7388 | $118.0\pm 4$ | GMOS
2459813.6045 | $90.47\pm 0.15$ | FEROS
2459814.5874 | $86.23\pm 0.10$ | FEROS
2459815.5927 | $82.090\pm 0.049$ | FEROS
2459817.5278 | $74.475\pm 0.053$ | FEROS
2459818.5266 | $70.810\pm 0.054$ | FEROS
2459818.7870 | $67.8\pm 4$ | GMOS
2459819.5543 | $67.205\pm 0.058$ | FEROS
2459820.5465 | $63.950\pm 0.067$ | FEROS
2459821.5669 | $60.685\pm 0.059$ | FEROS
2459822.5745 | $57.665\pm 0.085$ | FEROS
2459823.5422 | $54.98\pm 0.11$ | FEROS
2459823.8525 | $53.76\pm 0.1$ | HIRES
2459824.5306 | $52.200\pm 0.033$ | FEROS
2459824.8516 | $51.18\pm 0.1$ | HIRES
2459825.5361 | $49.90\pm 0.16$ | FEROS
2459826.7920 | $46.59\pm 0.1$ | HIRES
2459828.5677 | $42.795\pm 0.073$ | FEROS
2459829.5373 | $42.1\pm 3$ | MagE
2459829.5768 | $40.520\pm 0.086$ | FEROS
2459830.6175 | $38.76\pm 0.15$ | FEROS
2459831.6223 | $37.34\pm 0.38$ | FEROS
2459833.7523 | $33.23\pm 0.1$ | HIRES
2459834.5509 | $31.81\pm 0.14$ | FEROS
2459834.7691 | $31.74\pm 0.1$ | HIRES
2459835.7678 | $30.14\pm 0.1$ | HIRES
2459838.7208 | $27.5\pm 4$ | GMOS
2459838.8082 | $26.35\pm 0.1$ | HIRES
2459840.7729 | $24.20\pm 0.1$ | HIRES
2459845.5069 | $19.4\pm 1.0$ | XSHOOTER
2459855.5012 | $14.2\pm 1.0$ | XSHOOTER
2459868.5128 | $9.3\pm 1.0$ | XSHOOTER
2459877.6978 | $10.5\pm 1.5$ | ESI
2459985.8801 | $139.05\pm 0.23$ | FEROS
2459989.8698 | $129.35\pm 0.08$ | FEROS
2459992.8690 | $117.79\pm 0.06$ | FEROS
2459996.8740 | $99.86\pm 0.09$ | FEROS
2459997.8709 | $95.49\pm 0.15$ | FEROS
2459999.8760 | $86.66\pm 0.09$ | FEROS
2460000.8793 | $82.43\pm 0.07$ | FEROS
2460001.8731 | $78.57\pm 0.11$ | FEROS
2460002.8546 | $74.67\pm 0.08$ | FEROS
2460003.8613 | $71.13\pm 0.06$ | FEROS
2460004.8639 | $67.51\pm 0.07$ | FEROS
2460005.8594 | $64.21\pm 0.04$ | FEROS
2460006.8705 | $60.87\pm 0.05$ | FEROS
2460012.8827 | $44.99\pm 0.04$ | FEROS
2460014.8826 | $40.84\pm 0.07$ | FEROS
2460016.8810 | $37.03\pm 0.05$ | FEROS
2460017.8634 | $35.43\pm 0.05$ | FEROS
2460018.8627 | $33.73\pm 0.05$ | FEROS
2460026.8575 | $23.50\pm 0.07$ | FEROS
2460037.9770 | $15.670\pm 0.059$ | TRES
2460040.8413 | $14.03\pm 0.05$ | FEROS
2460042.8033 | $13.0181\pm 0.0062$ | ESPRESSO
2460043.8640 | $12.87\pm 0.05$ | FEROS
2460046.7778 | $11.7593\pm 0.0049$ | ESPRESSO
2460049.7651 | $11.0815\pm 0.0044$ | ESPRESSO
2460050.8708 | $11.24\pm 0.08$ | FEROS
2460050.8927 | $11.275\pm 0.056$ | TRES
2460052.8000 | $10.5756\pm 0.0056$ | ESPRESSO
2460053.7976 | $10.60\pm 0.09$ | FEROS
2460054.8315 | $10.620\pm 0.089$ | FEROS
2460060.8260 | $10.1052\pm 0.0044$ | ESPRESSO
2460065.8700 | $10.3429\pm 0.0054$ | ESPRESSO
2460073.7926 | $11.75\pm 0.09$ | FEROS
2460075.8674 | $12.613\pm 0.047$ | TRES
2460076.7473 | $12.1620\pm 0.0047$ | ESPRESSO
2460076.7645 | $12.330\pm 0.066$ | FEROS
2460077.7304 | $12.725\pm 0.072$ | FEROS
2460080.8029 | $13.645\pm 0.083$ | FEROS
2460085.7878 | $15.200\pm 0.068$ | FEROS
2460087.8111 | $15.7071\pm 0.0032$ | ESPRESSO
2460089.9047 | $17.140\pm 0.049$ | TRES
2460093.8351 | $18.835\pm 0.052$ | TRES
2460094.6938 | $18.8160\pm 0.0051$ | ESPRESSO
2460101.7831 | $22.8400\pm 0.0046$ | ESPRESSO
2460101.8146 | $23.060\pm 0.037$ | FEROS
2460102.8134 | $24.065\pm 0.046$ | TRES
2460107.6790 | $26.8994\pm 0.0042$ | ESPRESSO
2460110.7907 | $29.47\pm 0.10$ | FEROS
2460111.7161 | $30.370\pm 0.054$ | FEROS
2460112.7846 | $31.0566\pm 0.0047$ | ESPRESSO
2460113.6682 | $32.065\pm 0.037$ | FEROS
2460113.8126 | $32.575\pm 0.059$ | TRES
2460117.7582 | $35.7466\pm 0.0045$ | ESPRESSO
2460121.7794 | $40.755\pm 0.072$ | TRES
2460133.5611 | $56.5368\pm 0.0087$ | ESPRESSO
2460133.7396 | $57.465\pm 0.065$ | TRES
2460142.6035 | $74.505\pm 0.054$ | FEROS
2460143.5587 | $76.516\pm 0.0044$ | ESPRESSO
2460153.4901 | $103.8949\pm 0.0083$ | ESPRESSO
2460162.6793 | $132.722\pm 0.095$ | TRES
2460163.5326 | $134.0932\pm 0.0037$ | ESPRESSO
2460164.4936 | $136.1711\pm 0.0049$ | ESPRESSO
2460165.5038 | $138.0140\pm 0.0043$ | ESPRESSO
2460167.5333 | $140.4032\pm 0.0039$ | ESPRESSO
2460168.5190 | $140.8277\pm 0.0037$ | ESPRESSO
2460169.5509 | $140.7025\pm 0.0041$ | ESPRESSO
2460170.5049 | $140.0269\pm 0.0065$ | ESPRESSO
2460172.4964 | $136.9670\pm 0.0034$ | ESPRESSO
2460173.5017 | $134.5784\pm 0.0041$ | ESPRESSO
2460175.6095 | $128.0143\pm 0.0042$ | ESPRESSO
2460176.5698 | $124.4401\pm 0.0047$ | ESPRESSO
2460177.5139 | $120.6559\pm 0.0040$ | ESPRESSO
2460178.4970 | $116.4871\pm 0.0031$ | ESPRESSO
2460179.4876 | $112.1362\pm 0.0038$ | ESPRESSO
2460181.4862 | $103.1152\pm 0.0036$ | ESPRESSO
2460182.5170 | $98.4518\pm 0.0058$ | ESPRESSO
2460183.4844 | $94.1418\pm 0.0034$ | ESPRESSO
2460187.6126 | $77.225\pm 0.063$ | FEROS
2460192.6551 | $59.970\pm 0.042$ | TRES
2460193.5093 | $56.9746\pm 0.0052$ | ESPRESSO
2460197.4998 | $46.5795\pm 0.0070$ | ESPRESSO
2460202.5055 | $36.3992\pm 0.0055$ | ESPRESSO
2460202.6309 | $36.623\pm 0.030$ | TRES
2460204.5154 | $33.0720\pm 0.0040$ | ESPRESSO
2460206.5202 | $30.1237\pm 0.0057$ | ESPRESSO
2460206.6351 | $30.410\pm 0.031$ | TRES
2460211.5028 | $24.0717\pm 0.0078$ | ESPRESSO
2460213.5108 | $22.0683\pm 0.0052$ | ESPRESSO
2460223.5072 | $15.295\pm 0.099$ | FEROS

## References

  * Aarseth & Mardling (2001) Aarseth, S. J., & Mardling, R. A. 2001, in Astronomical Society of the Pacific Conference Series, Vol. 229, Evolution of Binary and Multiple Star Systems, ed. P. Podsiadlowski, S. Rappaport, A. R. King, F. D’Antona, & L. Burderi, 77, doi: 10.48550/arXiv.astro-ph/0011514
  * Amaro-Seoane et al. (2017) Amaro-Seoane, P., Audley, H., Babak, S., et al. 2017, arXiv e-prints, arXiv:1702.00786, doi: 10.48550/arXiv.1702.00786
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Bekenstein (1973) Bekenstein, J. D. 1973, ApJ, 183, 657, doi: 10.1086/152255
  * Bohlin et al. (2017) Bohlin, R. C., Mészáros, S., Fleming, S. W., et al. 2017, AJ, 153, 234, doi: 10.3847/1538-3881/aa6ba9
  * Boisse et al. (2010) Boisse, I., Eggenberger, A., Santos, N. C., et al. 2010, A&A, 523, A88, doi: 10.1051/0004-6361/201014909
  * Brahm et al. (2017) Brahm, R., Jordán, A., & Espinoza, N. 2017, PASP, 129, 034002, doi: 10.1088/1538-3873/aa5455
  * Brown & Bethe (1994) Brown, G. E., & Bethe, H. A. 1994, ApJ, 423, 659, doi: 10.1086/173844
  * Buchhave et al. (2010) Buchhave, L. A., Bakos, G. Á., Hartman, J. D., et al. 2010, ApJ, 720, 1118, doi: 10.1088/0004-637X/720/2/1118
  * Chakrabarti et al. (2023) Chakrabarti, S., Simon, J. D., Craig, P. A., et al. 2023, AJ, 166, 6, doi: 10.3847/1538-3881/accf21
  * Chawla et al. (2022) Chawla, C., Chatterjee, S., Breivik, K., et al. 2022, ApJ, 931, 107, doi: 10.3847/1538-4357/ac60a5
  * Corral-Santana et al. (2016) Corral-Santana, J. M., Casares, J., Muñoz-Darias, T., et al. 2016, A&A, 587, A61, doi: 10.1051/0004-6361/201527130
  * Di Carlo et al. (2023) Di Carlo, U., Agrawal, P., Rodriguez, C. L., & Breivik, K. 2023, arXiv.  https://arxiv.org/abs/2306.13121
  * El-Badry & Rix (2018) El-Badry, K., & Rix, H.-W. 2018, MNRAS, 480, 4884, doi: 10.1093/mnras/sty2186
  * El-Badry et al. (2023a) El-Badry, K., Rix, H.-W., Quataert, E., et al. 2023a, MNRAS, 518, 1057, doi: 10.1093/mnras/stac3140
  * El-Badry et al. (2023b) El-Badry, K., Rix, H.-W., Cendes, Y., et al. 2023b, MNRAS, 521, 4323, doi: 10.1093/mnras/stad799
  * Fűrész (2008) Fűrész, G. 2008, PhD thesis, University of Szeged, Hungary
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306, doi: 10.1086/670067
  * Fulton et al. (2020) Fulton, B., Blunt, S., Hurt, S., et al. 2020, California-Planet-Search/radvel: Version 1.4.0, v1.4.0, Zenodo, doi: 10.5281/zenodo.3814717
  * Fulton et al. (2018) Fulton, B. J., Petigura, E. A., Blunt, S., & Sinukoff, E. 2018, PASP, 130, 044504, doi: 10.1088/1538-3873/aaaaa8
  * Gaia Collaboration et al. (2023a) Gaia Collaboration, Vallenari, A., Brown, A. G. A., et al. 2023a, A&A, 674, A1, doi: 10.1051/0004-6361/202243940
  * Gaia Collaboration et al. (2023b) Gaia Collaboration, Arenou, F., Babusiaux, C., et al. 2023b, A&A, 674, A34, doi: 10.1051/0004-6361/202243782
  * Gibson et al. (2016) Gibson, S. R., Howard, A. W., Marcy, G. W., et al. 2016, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9908, Ground-based and Airborne Instrumentation for Astronomy VI, ed. C. J. Evans, L. Simard, & H. Takami, 990870, doi: 10.1117/12.2233334
  * Giesers et al. (2018) Giesers, B., Dreizler, S., Husser, T.-O., et al. 2018, MNRAS, 475, L15, doi: 10.1093/mnrasl/slx203
  * Gray (2008) Gray, D. F. 2008, The Observation and Analysis of Stellar Photospheres
  * Guseinov & Zel’dovich (1966) Guseinov, O. K., & Zel’dovich, Y. B. 1966, Soviet Ast., 10, 251
  * Halbwachs et al. (2023) Halbwachs, J.-L., Pourbaix, D., Arenou, F., et al. 2023, A&A, 674, A9, doi: 10.1051/0004-6361/202243969
  * Hayashi & Suto (2020) Hayashi, T., & Suto, Y. 2020, ApJ, 897, 29, doi: 10.3847/1538-4357/ab97ad
  * Hayashi et al. (2023) Hayashi, T., Suto, Y., & Trani, A. A. 2023, ApJ, 958, 26, doi: 10.3847/1538-4357/acf4f6
  * Hayashi et al. (2020) Hayashi, T., Wang, S., & Suto, Y. 2020, ApJ, 890, 112, doi: 10.3847/1538-4357/ab6de6
  * Hirai & Mandel (2022) Hirai, R., & Mandel, I. 2022, ApJ, 937, L42, doi: 10.3847/2041-8213/ac9519
  * Howell & Furlan (2022) Howell, S. B., & Furlan, E. 2022, Frontiers in Astronomy and Space Sciences, 9, 871163, doi: 10.3389/fspas.2022.871163
  * Justham et al. (2014) Justham, S., Podsiadlowski, P., & Vink, J. S. 2014, ApJ, 796, 121, doi: 10.1088/0004-637X/796/2/121
  * Kaufer et al. (1999) Kaufer, A., Stahl, O., Tubbesing, S., et al. 1999, The Messenger, 95, 8
  * Koposov et al. (2023) Koposov, S., Speagle, J., Barbary, K., et al. 2023, joshspeagle/dynesty: v2.1.3, v2.1.3, Zenodo, doi: 10.5281/zenodo.8408702
  * Kozai (1962) Kozai, Y. 1962, AJ, 67, 591, doi: 10.1086/108790
  * Kurucz (1979) Kurucz, R. L. 1979, ApJS, 40, 1, doi: 10.1086/190589
  * Kurucz (1993) —. 1993, SYNTHE spectrum synthesis programs and line data
  * Lidov (1962) Lidov, M. L. 1962, Planet. Space Sci., 9, 719, doi: 10.1016/0032-0633(62)90129-0
  * Lindegren & Dravins (2003) Lindegren, L., & Dravins, D. 2003, A&A, 401, 1185, doi: 10.1051/0004-6361:20030181
  * Liu et al. (2022) Liu, B., D’Orazio, D. J., Vigna-Gómez, A., & Samsing, J. 2022, Phys. Rev. D, 106, 123010, doi: 10.1103/PhysRevD.106.123010
  * Luhn et al. (2020) Luhn, J. K., Wright, J. T., Howard, A. W., & Isaacson, H. 2020, AJ, 159, 235, doi: 10.3847/1538-3881/ab855a
  * Merritt et al. (2004) Merritt, D., Milosavljević, M., Favata, M., Hughes, S. A., & Holz, D. E. 2004, ApJ, 607, L9, doi: 10.1086/421551
  * Moe & Di Stefano (2017) Moe, M., & Di Stefano, R. 2017, ApJS, 230, 15, doi: 10.3847/1538-4365/aa6fb6
  * Morais & Correia (2008) Morais, M. H. M., & Correia, A. C. M. 2008, A&A, 491, 899, doi: 10.1051/0004-6361:200810741
  * Morais & Correia (2011) —. 2011, A&A, 525, A152, doi: 10.1051/0004-6361/201014812
  * Naoz (2016) Naoz, S. 2016, ARA&A, 54, 441, doi: 10.1146/annurev-astro-081915-023315
  * Pepe et al. (2021) Pepe, F., Cristiani, S., Rebolo, R., et al. 2021, A&A, 645, A96, doi: 10.1051/0004-6361/202038306
  * Peters (1964) Peters, P. C. 1964, Phys. Rev., 136, B1224, doi: 10.1103/PhysRev.136.B1224
  * Portegies Zwart et al. (1997) Portegies Zwart, S. F., Verbunt, F., & Ergma, E. 1997, A&A, 321, 207, doi: 10.48550/arXiv.astro-ph/9701037
  * Price-Whelan et al. (2017) Price-Whelan, A. M., Hogg, D. W., Foreman-Mackey, D., & Rix, H.-W. 2017, ApJ, 837, 20, doi: 10.3847/1538-4357/aa5e50
  * Rastello et al. (2023) Rastello, S., Iorio, G., Mapelli, M., et al. 2023, arXiv.  https://arxiv.org/abs/2306.14679
  * Rein & Liu (2012) Rein, H., & Liu, S. F. 2012, A&A, 537, A128, doi: 10.1051/0004-6361/201118085
  * Rein et al. (2023) Rein, H., Tamayo, D., Liu, S., et al. 2023, hannorein/rebound: 4.0.1, 4.0.1, Zenodo, doi: 10.5281/zenodo.10121021
  * Remillard & McClintock (2006) Remillard, R. A., & McClintock, J. E. 2006, ARA&A, 44, 49, doi: 10.1146/annurev.astro.44.051905.092532
  * Schwarz (1978) Schwarz, G. 1978, The Annals of Statistics, 6, 461 , doi: 10.1214/aos/1176344136
  * Shahaf et al. (2023) Shahaf, S., Hallakoun, N., Mazeh, T., et al. 2023, arXiv e-prints, arXiv:2309.15143, doi: 10.48550/arXiv.2309.15143
  * Shao & Li (2021) Shao, Y., & Li, X.-D. 2021, ApJ, 920, 81, doi: 10.3847/1538-4357/ac173e
  * Shenar et al. (2022) Shenar, T., Sana, H., Mahy, L., et al. 2022, Nature Astronomy, 6, 1085, doi: 10.1038/s41550-022-01730-y
  * Silsbee & Tremaine (2017) Silsbee, K., & Tremaine, S. 2017, ApJ, 836, 39, doi: 10.3847/1538-4357/aa5729
  * Speagle (2020) Speagle, J. S. 2020, MNRAS, 493, 3132, doi: 10.1093/mnras/staa278
  * Tanikawa et al. (2023) Tanikawa, A., Cary, S., Shikauchi, M., Wang, L., & Fujii, M. S. 2023, arXiv e-prints, arXiv:2303.05743, doi: 10.48550/arXiv.2303.05743
  * Toonen et al. (2020) Toonen, S., Portegies Zwart, S., Hamers, A. S., & Bandopadhyay, D. 2020, A&A, 640, A16, doi: 10.1051/0004-6361/201936835
  * Trimble & Thorne (1969) Trimble, V. L., & Thorne, K. S. 1969, ApJ, 156, 1013, doi: 10.1086/150032
  * Vogt et al. (1994) Vogt, S. S., Allen, S. L., Bigelow, B. C., et al. 1994, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 2198, Instrumentation in Astronomy VIII, ed. D. L. Crawford & E. R. Craine, 362, doi: 10.1117/12.176725
  * Wagg et al. (2022a) Wagg, T., Breivik, K., & de Mink, S. 2022a, The Journal of Open Source Software, 7, 3998, doi: 10.21105/joss.03998
  * Wagg et al. (2022b) Wagg, T., Breivik, K., & de Mink, S. E. 2022b, ApJS, 260, 52, doi: 10.3847/1538-4365/ac5c52
  * Yamaguchi et al. (2023) Yamaguchi, N., El-Badry, K., Fuller, J., et al. 2023, arXiv e-prints, arXiv:2309.15905, doi: 10.48550/arXiv.2309.15905