Constraining the binarity of black hole candidates: a proof-of-concept study of Gaia BH1 and Gaia BH2

Nearly a hundred of binary black holes (BBHs) have been discovered with gravitational-wave signals emitted at their merging events. Thus, it is quite natural to expect that significantly more abundant BBHs with wider separations remain undetected in the universe, or even in our Galaxy. We consider a possibility that star-BH binary candidates may indeed host an inner BBH, instead of a single BH. We present a detailed feasibility study of constraining the binarity of the currently available two targets, Gaia BH1 and Gaia BH2. Specifically, we examine three types of radial velocity (RV) modulations of a tertiary star in star-BBH triple systems; short-term RV modulations induced by the inner BBH, long-term RV modulations induced by the nodal precession, and long-term RV modulations induced by the von Zeipel-Kozai-Lidov oscillations. Direct three-body simulations combined with approximate analytic models reveal that Gaia BH1 system may exhibit observable signatures of the hidden inner BBH if it exists at all. The methodology that we examine here is quite generic, and is expected to be readily applicable to future star-BH binary candidates in a straightforward manner.


INTRODUCTION
Since the first discovery of GW150914 (Abbott et al. 2016), more than 90 candidates for binary black holes (BBHs) have been reported so far (The LIGO Scientific Collaboration et al. 2021).The formation and evolution of such BBHs are one of the important unsolved questions in astrophysics, and there are a variety of proposed scenarios; (1) final stage of isolated massive binary stars (e.g.Belczynski et al. 2002Belczynski et al. , 2007Belczynski et al. , 2012Belczynski et al. , 2016a,b;,b;Dominik et al. 2012Dominik et al. , 2013;;Kinugawa et al. 2014Kinugawa et al. , 2016;;Spera et al. 2019), (2) dynamical capture in dense clusters (e.g.Portegies Zwart & McMillan 2000;O'Leary et al. 2009;Rodriguez et al. 2016;Tagawa et al. 2016;Fragione et al. 2020;Trani et al. 2022), (3) binary formation channel of primordial black holes (e.g.Ioka et al. 1999;Bird et al. 2016;Sasaki et al. 2016Sasaki et al. , 2018;;Kocsis et al. 2018), and (4) evolution of wide binaries under the effect of galactic tide or cumulative flybys (e.g.Michaely & Perets 2016, 2019;Michaely & Naoz 2022).Regardless of those different formation scenarios, their progenitors are expected to have a longer orbital period.The subsequent dynamical evolution decreases their orbital energy and angular momentum, and eventually leads to the BBH merger events that are detectable using gravitational wave (GW) observations.Therefore, it is natural to expect that more abundant wide-separation BBHs remain undetected in the universe, or even in our Galaxy.
In order to search for BBHs with relatively long orbital periods that cannot be probed with GWs, Hayashi et al. (2020) and Hayashi & Suto (2020) pointed out that a BBH orbited by a tertiary star would be detectable in optical spectroscopic surveys from the radial velocity (RV) modulations of the tertiary star; the inner BBH produces the observable RV modulations of the star in short (a half of the orbital period of the BBH) and long (a nodal precession timescale and/or the von Zeipel-Kozai-Lidov timescale) terms.They examined the feasibility of the strategy from three-body simulations for hypothetical triple systems of an inner BBH and an outer tertiary star, and proposed that the methodology can distinguish between the single black hole (BH) and BBH when applied to the future star-BH binary candidates from on-going Gaia (Gaia Collaboration et al. 2016) and TESS (Ricker et al. 2014) surveys (e.g.Kawanaka et al. 2016;Breivik et al. 2017;Mashian & Loeb 2017;Yamaguchi et al. 2018;Masuda & Hotokezaka 2019;Shikauchi et al. 2020;Chawla et al. 2022).
Indeed, star-BH binary candidates, Gaia BH1 and Gaia BH2, recently discovered from Gaia DR3 astrometric data (Gaia Collaboration et al. 2022;El-Badry et al. 2023a,b;Chakrabarti et al. 2023;Tanikawa et al. 2023) would provide a good opportunity to directly check the methodology.Gaia BH1 is a binary of a ∼ 1M ⊙ main sequence star and a ∼ 10M ⊙ dark companion, with an orbital period P obs ∼ 190 days (El-Badry et al. 2023a;Chakrabarti et al. 2023; see also Rastello et al. 2023).Gaia BH2 was first discovered by Tanikawa et al. (2023) using Gaia astrometry, and later more robustly identified combining the follow-up RV observations by El-Badry et al. (2023b).Gaia BH2 is a binary of a ∼ 1M ⊙ red giant and a ∼ 9M ⊙ dark companion, with P obs ∼ 1300 days.The best-fit values of their system parameters are listed in Table 1.
Due to the limited precision and duration of the current spectroscopic monitoring observations, it is not possible to prove the presence of an inner BBH, instead of a single dark companion, in either system.Nevertheless, those systems are useful as a proof-of-concept in constraining the binarity of the dark companion for future star-BH binary candidates.
We first consider the short-term RV modulations on the timescale of half the inner orbital period.We next move on to the long-term RV modulations, which become important for inclined triples.We put constraints, as one application, on the binarity of dark companions in Gaia BH1 and Gaia BH2.For reference, Figure 1 shows the configuration of a triple that we consider in the present paper.
The rest of paper is organized as follows.Section 2 examines the short-term RV modulations.We first discuss the short-term semi-amplitude predicted from an analytic approximation for coplanar and circular triples.Then, we show that the outer eccentricity, such as those for Gaia BH1 and Gaia BH2, significantly increases the simple prediction by direct three-body simulations.Next, section 3, focuses on the long-term RV modulations induced by the nodal precession for moderately inclined triples.We also discuss analytic predictions first, and then examine their validity using three-body simulations.Section 4 considers more significantly inclined triples in which the von Zeipel- Kozai-Lidov (ZKL) oscillations (von Zeipel 1910;Kozai 1962;Lidov 1962) play an important role.Finally, we summarize the constraints on Gaia BH1 and Gaia BH2, and discuss a future prospect in section 5.

SHORT-TERM RV MODULATIONS
For a coplanar triple system, the inner binary efficiently induces short-term wobbles of the tertiary, with about a half the inner orbital period P in .For inclined triples, however, the additional long-term RV modulations are generated due to the misalignment between the inner and outer orbital angular momenta.This section focuses on coplanar triples, and discusses the amplitude of the short-term RV modulations using an analytic approximation and numerical simulations.The long-term RV modulations for inclined triples will be discussed in later sections.

Analytic estimates
. Schematic illustration of a triple system configuration that we consider in the present paper.The orbital angles are defined with respect to the reference Cartesian frame whose origin is set to be the barycenter of the inner orbit.The left panel shows the inner and outer orbits, and the right panel defines various angles specifying the orientations of the total (G tot ), inner (G in ) and outer orbital angular momenta of the system (G out ); the mutual inclination between the inner and outer binaries (i mut ), the orbital inclinations of the inner (i in ) and outer (i out ) binaries with respect to G tot , the inclination of the line-of-sight with respect to G tot (I los ), and the inclination of the outer orbit with respect to the observer's line-of-sight (I obs ).The time-dependent azimuthal angle of G out relative to the line-of-sight is ϕ(t).
The short-term RV modulations for coplanar and circular triples are, to the leading order of a in /a out , analytically approximated as (e.g.Hayashi et al. 2020) where ν in , ν out , θ in and θ out are the orbital frequencies and the true longitudes (the angle between a body's location and the ascending node) for inner and outer orbits, respectively.The true longitudes θ in and θ out in equation ( 1) are expressed as and for eccentric orbits, where f in , f out , ω in and ω out are the true anomalies and the pericenter arguments of inner and outer orbits, respectively, evaluated at the initial epoch.Incidentally, we define ν in and ν out as orbital frequencies of inner and outer orbits throughout the present paper, which should not be confused with the true anomalies (denoted by f here).
In equation ( 1), K short corresponds to a characteristic semi-amplitude of the short-term RV modulations defined as with G being Newton's gravitational constant, and a in , a out , P in , P out and I obs are the semi-major axes, orbital periods of inner and outer orbits, and the observed inclination, respectively.Equations ( 1), ( 4) and ( 5) can be derived from equations ( 21), ( 25)-( 28) in Morais & Correia (2008), while they use somewhat a different notation.
Since we assume a triple with an inner binary companion throughout the present analysis, the orbital parameters with the subscript "out" are interpreted to be those estimated for the star-BH binary (with the subscript "obs" in Table 1 for Gaia BH1 and Gaia BH2).Similarly, we assume that m 12 = m 1 + m 2 is equal to m c in Table 1.
Figure 2 plots the contours of K short in the q 21 ≡ m 2 /m 1 -P in plane; for Gaia BH1 (left) and Gaia BH2 (right).The shaded regions indicate those corresponding to dynamical instability condition for coplanar triples by Mardling & Aarseth (1999); Aarseth & Mardling (2001) (hereafter, MA01): The condition (6) turned out to be a good approximation for coplanar triples (i mut = 0 • ).Hayashi et al. (2022Hayashi et al. ( , 2023) ) examined the Lagrange stability timescales of triples in general, and found that the condition (6) needs to be improved especially for inclined triples that exhibit the ZKL oscillations.
Figure 2 indicates that the expected values of K short (dotted contours) are fairly small; at most O(10) m/s for Gaia BH1, and O(1) m/s for Gaia BH2.In reality, however, the observed semiamplitude should be sensitive to the mutual phases of the three bodies, in particular for eccentric outer orbits as in the cases of both Gaia BH1 and Gaia BH2.While K short in equation ( 4) is derived for circular orbits, we find that the effect of the outer eccentricity can be empirically taken into account by replacing a out by a out (1 − e out ) in equation ( 4), i.e., K short (1 − e obs ) −7/2 as plotted in red solid contours in Figure 2 (see Figure 3 for detail).
In the next subsection, we perform three-body simulations and show that the phase-dependent RV modulation amplitudes become even larger for Gaia BH1 and BH2 around the pericenter passages, due to their relatively large e obs .
at the initial epoch.Then, the initial true anomalies f j (j = in, out) in the simulations are computed from the mean anomalies M j using the standard relations among the true anomaly f , eccentric anomaly E, mean anomaly M , and orbital eccentricity e: and where variables with subscript j denote their values evaluated at the initial epoch.
In order to remove possible transient behavior due to the choice of initial phase angles, we first evolve the system over 100 outer orbital periods P out (= P obs ).The top panels of Figure 3 is the resulting RV curve for t = 100 P obs to 101 P obs .Then, we fit the simulated RV data with the public code, RadVel (see Fulton et al. 2018) so as to remove the overall Kepler motion of the tertiary star.The resulting residuals (middle and bottom panels) represent the short-term RV modulations.We perform the RadVel fitting only for a single outer orbital period.If we do so over longer periods, the residuals become significantly larger in most of the periods because a long-term trend due to the three-body effect cannot be removed as long as purely two-body dynamics is assumed (e.g.Hayashi & Suto 2020).In other words, the short-term modulations in the middle and bottom panels of Figure 3 are the robust residuals reflecting P in /2.
Left and right panels of Figure 3 correspond to Gaia BH1 with P in = 10 days, and Gaia BH2 with P in = 50 days, both of which are close to the dynamical stability limit (Figure 2).Red and blue curves show the results for the initial inner eccentricities of e in = 0 and e in = 0.2, respectively.The difference of e in produces a small phase shift of the total and residual RV, but does not affect their amplitudes in practice.
For reference, we plot the analytic short-term modulation semi-amplitude ±K short , equation ( 4), and also ±K short (1 − e obs ) −7/2 ; see magenta and cyan regions in middle and bottom panels of Figure 3. Clearly, K short significantly underestimates the simulated amplitudes.Indeed, the simulated RV modulations become even larger around the pericenter passage of the tertiary; the short-term RV modulations for Gaia BH1 and Gaia BH2 amount to ∼ 300 m/s and ∼ 100 m/s around the epoch.Those values are about 10-100 times larger than the analytic approximation K short , equation (4), and may be detectable for Gaia BH1 from the observed RV residuals according to Figure 4 of El-Badry et al. (2023a).

Possible degeneracy between an S-type planet and an inner binary black hole companion
An S-type planet (i.e., a planet around the stellar member of the binary) would induce short-term RV modulations that are similar to those produced by an inner BBH as we presented in the above.The possible degeneracy of the short-term RV modulations between an inner binary companion in a triple system and a planet orbiting a star in a binary system has been pointed by Schneider & Cabrera (2006), and discussed by Morais & Correia (2008).Although precise spectroscopic measurements can break the degeneracy of the short-term RV signals between an inner binary and a planet (Morais & Correia 2008) in principle, Hayashi et al. (2020) argue that the degeneracy is difficult to be broken under realistic observational noises and limited cadences.Therefore, it is interesting to consider the correspondence of the two interpretations, especially for the future Gaia BH1 signals.
According to Hayashi et al. (2020), the RV semi-amplitude K s induced by an S-type planet of mass M pl in a circular and coplanar orbit is .Mass ratio of inner BBH q 21 (dashed) and the corresponding planet mass M pl (solid) in RV signal semi-amplitude K s -period P s plane for Gaia BH1.We note that K s and P s are interpreted as RV semi-amplitude and period induced by a planet for the planet interpretation, while they are interpreted as K s = K short (1 − e obs ) −7/2 and P s = P in /2 for the BBH interpretation.
where P s is the orbital period of planet.Figure 4 shows contours of M pl (solid lines) and the mass ratio q 21 (dashed lines) for the planet and BBH induced RV modulations of the semi-amplitude K s and period P s , respectively.For the planet interpretation, we assume a planet with a circular and coplanar orbit, while we assume K s = K short (1 − e obs ) −7/2 and P s = P in /2 for the BBH interpretation.Figure 4 clearly indicates that a close-in hot jupiter induces a short-term RV modulation similar to that we expect from an inner BBH.While it is a possible false-positive for the inner BBH search, it may reveal a first-ever planetary system orbiting a BH.Given the observed frequency of the close-in planet, it is likely to find such a system in future even if Gaia BH1 is not the one.

LONG-TERM RV MODULATIONS FOR MODERATELY INCLINED SYSTEMS: NODAL PRECESSION
Consider next non-coplanar triples, i.e., the inner and outer orbits are mutually inclined.Hayashi & Suto (2020) pointed out that the long-term RV modulations of the tertiary body due to the nodal precession and the ZKL oscillations may carry interesting signatures of the hidden inner binary.The details of the inclined three-body dynamics are described in previous literature (e.g.Morais & Correia 2012;Naoz 2016).
In this section, we focus on the nodal precession in moderately inclined systems (i mut ≲ 50 • ).First, we consider analytic approximations for the timescale and the RV modulation amplitude of the nodal precession.Then, we perform three-body simulations to present more quantitative prediction, and discuss the observational feasibility.

Nodal precession timescale
If e in is initially small and i mut is moderate (i mut ≲ 50 • ), the outer ascending node Ω out regularly precess with the following timescale P Ω (e.g., Hayashi & Suto 2020): where C quad is the quadrupole strength coefficient: and G in , G out , and G tot are the amplitudes of inner, outer, and total angular momenta: In equations ( 12) and ( 13), µ in and µ out denote the reduced masses: where q 21 ≡ m 2 /m 1 is the mass ratio of the inner binary.
It is convenient to introduce the ratio of the amplitudes of inner to outer angular momenta: which is a key parameter that characterizes the long-term modulation due to the nodal precession.Figure 5 plots ξ against P in for Gaia BH1 and Gaia BH2 in solid and dashed lines; different colors correspond to (q 21 , e in ) = (1, 0), (1, 0.3), (0.1, 0), and (0.1, 0.3).As equation ( 17) indicates, ξ is sensitive to q 21 , but not to e in as long as e 2 in ≪ 1.The realistic range of ξ values is shown in this figure for a given set of q 21 and P in .Due to dynamical stability, ξ cannot exceed about 1.2 and 0.9 for Gaia BH1 and Gaia BH2, respectively.
By rewriting equation ( 10) in terms of ξ as equation ( 18) reduces to Gaia BH1 Gaia BH2 Figure 5.The ratio of the inner and outer angular momenta, ξ, plotted against P in for different q 21 ≡ m 2 /m 1 and e in .The solid and dashed lines correspond to Gaia BH1 and Gaia BH2, respectively. = 0.3  = 0.6  = 0.9 Figure 6.The normalized timescale of the nodal precession, P Ω /P out , plotted against cos i mut for different ξ.The solid and dashed lines correspond to Gaia BH1 and Gaia BH2, respectively.We fix q 21 = 1 and e in = 0 for simplicity.
Equation ( 19) implies that the nodal precession timescale is very sensitive to ξ. Figure 6 plots P Ω /P out as a function of cos i mut for Gaia BH1 (solid) and BH2 (dashed) with e in = 0 and q 21 = 1.0.The plot also shows that P Ω /P out is not sensitive to i mut as long as moderately inclined triples are considered.
. The normalized timescale of the nodal precession, P Ω /P out , plotted against ξ for i mut = 10 • (red) and 50 • (blue).The solid and dashed curves correspond to Gaia BH1 and Gaia BH2, respectively.The q 21 and e in are fixed as 1.0 and 0.0, respectively.
Figure 7 shows P Ω /P out as a function of ξ for Gaia BH1 (solid) and BH2 (dashed) with q 21 = 1 and e in = 0. Since P Ω is a strongly decreasing function of ξ, triples with a larger value of ξ are preferable for a successful detection of long-term RV modulations.

Relation between inclination angles
The long-term RV modulations due to the nodal precession are computed as a function of inclination angles illustrated in Figure 1.First, we derive a relation between i out and i mut , which proves to be useful in the later discussion.
The inner and outer inclinations i in and i out simply satisfy and If we further neglect the ZKL oscillations and simply consider the nodal precession alone, i mut is nearly constant.Therefore, the value of i out is determined from equations ( 20) and ( 21), and also stays constant in practice.For moderately inclined triples of i mut ≲ 50 • , we obtain Figure 8 plots the outer inclination i out against i mut for various values of ξ.The plot indicates that i out becomes asymptotically close to i mut as ξ increases.We intend to show Figure 8 25) for different ξ = 0.3, 0.6, 0.9, 1.2 and 10.For reference, we plot the dashed line that corresponds to the asymptotic limit of ξ → ∞.Figures 2 and 5 imply that ξ = 1.2 and 0.9 roughly correspond to the maximum possible values for Gaia BH1 and Gaia BH2, respectively.
relation between two angles.In the cases of Gaia BH1 and Gaia BH2, however, large values of ξ are nonphysical due to dynamical stability requirements; ξ ≲ 1.2 and ξ ≲ 0.9 for Gaia BH1 and Gaia BH2, respectively.Thus, rather small values of i out are allowed.For moderately inclined triples in which the nodal precession (i.e.Ω out precession) dominates the dynamics, Ω out and ϕ(t) (see Figure 1) change gradually from 0 • to 360 • with timescale P Ω .Thus, I obs (t) varies within the following range: We can insert in equation ( 24) the expression for i out in terms of i mut (< 90 • ) using the relation: derived from equations ( 22) and ( 23).
Figure 9 shows the constraints on the inclination angles of Gaia BH1 (left) and Gaia BH2 (right) from the observed value of I obs .If future observations detect any change of the RV semi-amplitude, or equivalently that of I obs , this plot is useful in inferring the geometric configuration of the corresponding triple system.
The observed RV semi-amplitude K(t) is proportional to sin I obs (t).In the case of the nodal precession alone, we obtain from Figure 1 cos I obs (t) = sin I los sin i out cos ϕ(t) + cos I los cos i out .
(26) Thus,  1).The bottom panels plot I los against i mut for given values of ξ.The region between two lines is permitted for each ξ, and the shaded regions correspond to ξ > 1.2 for Gaia BH1 and ξ > 0.9 for Gaia BH2, which are excluded from dynamical stability viewpoint.
where Γ ≡ cot i out cot I los , and the precession angle ϕ(t) varies from 0 • to 360 • periodically with the timescale of P Ω as Ω out precesses.
The above argument is simply summarized as and where V 0 is the RV semi-amplitude for an edge-on system: If future long-term RV monitoring (over the duration exceeding P Ω ) identifies the RV modulation of Gaia BH1 and BH2, equations ( 28) and (29) determine the inclination angles of the line-of-sight and the outer orbits, I los and i out , separately.If the dark companion is a single BH, instead of BBH, i out = 0 • and I los = I obs always (see Figure 1).The inner binarity of the dark companion may be revealed by i out ̸ = 0 • .Note that there is a parameter degeneracy of I obs ↔ 180 • − I obs in the RV observation, but the astrometry indeed breaks this degeneracy.Figure 10 summarizes the expected fractional change of the RV semi-amplitude, ∆ K , in the P in -P Ω plane for Gaia BH1 (left) and Gaia BH2 (right).Specifically, we define ∆ K using equations ( 28) and (29): For simplicity, we here assume e in = 0 and q 21 = 1 for both Gaia BH1 and Gaia BH2.In addition, we fix I los = 120 • and 30 • for Gaia BH1 and Gaia BH2, respectively, corresponding to the values close to their I obs ; see Table 1.
The light-blue regions correspond to the dynamically unstable region from MA01 (see equation ( 6)).We note that for high mutual inclination (i mut ≳ 50 • ), the analytic discussion based on the nodal precession becomes invalid since the ZKL oscillations become important.For moderate inclination, however, we can safely estimate ∆ K , and corresponding P in and P Ω .Figure 10 implies that ∆ K = 0.2-0.4variations are expected within 100 yrs for Gaia BH1 if P in = 5-10 days, while unrealistically long observational duration is required to detect the similar level of variations for Gaia BH2.

MA01 dynamically unstable
Excluded from geometry

MA01 dynamically unstable
Excluded from geometry Figure 10.Contours of the RV semi-amplitude changes ∆ K (red-solid) and the mutual inclination i mut (black-dotted) from nodal precession on P Ω -P in plane for Gaia BH1 (left) and Gaia BH2 (right).The lightblue regions correspond to the dynamically unstable region from MA01.The green regions are excluded from the observed value of I obs ; see equation ( 24).For relatively large mutual inclination, i mut ≳ 50 • , the predictions of these plots become less reliable because the ZKL oscillations become important.
The top panels of Figure 11 show the simulated RV semi-amplitudes against t/P Ω for Gaia BH1 (left) and Gaia BH2 (right), respectively.The red and blue curves indicate the envelope of the radial velocity K(t)/V 0 , i.e., neglecting the periodic changes over P out , which we define as RV max /V 0 and RV min /V 0 .The normalized RV semi-amplitude K/V 0 ≡ (RV max − RV min )/2V 0 is plotted in solid black curves, which should be compared with an analytic prediction, equation (31) with equations ( 28) and ( 29).In the plots, we show the analytically estimated ∆ K as magenta regions, and the expected semi-amplitude change from equation ( 27) as dotted green curves.We chose the initial phases ϕ(t = 0) so that they agree with the values from the simulations.
As expected, the ZKL oscillations are negligible for the present case (i mut = 20 • initially), and the mutual inclination is nearly constant over the period of P Ω .The simulated RV semi-amplitude changes almost sinusoidally with a period of ∼ P Ω (black curve), and its fractional change ∆ K is indeed in good agreement with the value predicted from the analytic approximation; ∆ K ≈ 0.2 for Gaia BH1, and ∆ K ≈ 0.3 for Gaia BH2, see Figure 10.
The example for Gaia BH1 indicates the RV semi-amplitude change of as large as 17 km/s, corresponding to ∆ K = 0.2, within P Ω /2 ≈ 26 yrs, depending on the phase.Furthermore, the zero-point of the RV curve also changes significantly.Thus, future long-term RV monitoring of Gaia BH1 should provide strong constraints on, or even detect, its inner BBH.On the contrary, the case of Gaia BH2 is very difficult because its P Ω is too long.

LONG-TERM RV MODULATIONS FOR SIGNIFICANTLY INCLINED SYSTEMS: ZKL OSCILLATIONS
Finally, we consider triple systems whose inner binary orbit is significantly inclined, i mut > 50 • , relative to the outer orbit.In this case, analytic discussion is not easy due to the strong ZKL oscillations.Thus, we present examples of numerical simulations alone.
The middle and bottom panels of Figure 11 are the same as the top panels except that their initial mutual inclinations are i mut = 60 • and i mut = 90 • , respectively.Note that they are plotted against t/P out , and their long-term modulation period is roughly consistent with the quadrupole ZKL timescale (e.g.Antognini 2015): The middle panels, with i mut = 60 • initially, indicate that the amplitude of the ZKL oscillations are still modest in this example, and the resulting semi-amplitude change (black curves) is roughly sinusoidal as expected for nodal precession alone.Moreover, the analytic prediction of ∆ K , equation (31), agrees with the simulated value within ten percent.
In contrast, the bottom plots, with i mut = 90 • initially, show non-trivial RV curves, due to the strong ZKL oscillations.For most of the time, the systems stay at mutually orthogonal orbits, but suddenly move to i mut ≈ 40 • .While the change of the mutual inclination is very periodic roughly with the ZKL timescale T ZKL , equation (32), the corresponding RV semi-amplitude changes are no longer periodic.Therefore, long-term RV monitoring of such systems may detect the significant change of the RV semi-amplitude even for relatively short timescales, or barely no change for long duration, depending on the phase of the observation over the sporadic behavior represented in the bottom panels of Figure 11.

SUMMARY AND DISCUSSION
Triple systems are ubiquitous in the universe, and trigger a wide variety of interesting observable events in astronomy.While nearly a hundred of BBHs have been discovered from the GWs emitted at the final instance of their coalescence, there is no candidate for triples including two BHs yet.Needless to say, such triples are fascinating targets for observational astronomy.Furthermore, star-BBH or even triple BH systems may provide an important mechanism to accelerate the GW merger of the detected BBHs (e.g.Liu & Lai 2018;Trani et al. 2022).
Formation and evolution of stellar triples are fundamental, but theoretically challenging, problems in broad areas of astrophysics.Their proper understanding requires many complicated physical processes, including the evolution of common envelope phases, supernova explosions, and the subsequent dynamics of the resulting compact objects (e.g.Toonen et al. 2021).Thus, future discoveries of star-BH binaries and star-BBH triples that we consider in the present paper would shed complementary observational insights that are useful in constructing and testing theoretical models.Hayashi et al. (2020) and Hayashi & Suto (2020) have proposed a methodology to discover a hidden inner BBH in star-BH binary candidates from the radial velocity modulations of the orbiting For simplicity, we assume the equal-mass inner binary (q 12 = 1) in an initially circular orbit (e in = 0).For reference, the range of the analytic estimate of ∆ K , equation ( 31), is plotted in magenta regions, and the expected semi-amplitude change, from equation ( 27), is plotted as dotted green curves.Equation ( 10) diverges for i mut = 90 • , and we do not show the analytic predictions in those cases.
(tertiary) star.Recent discoveries of such systems, Gaia BH1 and BH2 (El-Badry et al. 2023a,b;Chakrabarti et al. 2023;Tanikawa et al. 2023), provide a great opportunity to examine the feasibility of their methodology in detail as a proof of concept.Even if the dark companion of Gaia BH1 and BH2 turn out to be a single BH instead of a BBH, the analysis presented here is readily applicable for future star-BH candidates that remain to be discovered.The results of our proof-of-concept study are summarized below.
(1) short-term RV modulations induced by the inner BBH: Inner BBH generates a smallamplitude modulation of period P in on the RV of the tertiary star.The semi-amplitudes based on an analytic approximation are O(10) m/s for Gaia BH1, and O(1) m/s for Gaia BH2, if the tertiary is on a coplanar and circular orbit.In reality, relatively large eccentricities of e obs ∼ 0.5 for both systems are expected to significantly increase the semi-amplitude.Our numerical simulations indicate that the semi-amplitude of the short-term RV modulation increases by more than a factor of (1 − e out ) −7/2 (≈ 11) near the pericenter passage.Thus, the resulting amplitudes amount to ∼ 300 m/s for Gaia BH1, and ∼ 100 m/s for Gaia BH2 at their pericenter passage phases.We conclude that high-cadence and precise RV followups near the pericenter passages of the star are promising to search for possible inner BBHs for star-BH candidates with large e obs .
(2) long-term RV modulations induced by the nodal precession: If the orbit of the inner BBH is moderately inclined relative to that of the tertiary, i mut ≲ 50 • , the nodal precession generates long-term modulations of the radial velocity, or equivalently of the inclination I obs of the tertiary relative to the observer's line of sight.Unlike the short-term RV modulation, the nodal precession changes the RV semi-amplitude of the tertiary by a factor of sin I obs .Thus, the change of the RV semi-amplitude, ∆ K V 0 , is significantly larger than that of the short-term modulation, but its modulation period P Ω may be unrealistically long.Our examples from three-body simulations (equal-mass circular BBH with P in = 10 days) predict the RV semiamplitude change of 17 km/s within ∼ 26 yrs for Gaia BH1, assuming that the line-of-sight inclination I los = 120 • is close to the observed inclination I obs = 126 • .6.For Gaia BH2, the nodal precession timescale is too long to be detectable within a reasonable observation duration.More importantly, we confirm that our simple analytic estimates of ∆ K and P Ω reproduce well the simulation results.
(3) long-term RV modulation induced by the ZKL oscillations: For highly inclined triples, the ZKL oscillations induce the drastic and non-periodic RV semi-amplitude change, and analytic approximation becomes less reliable than the case with the nodal precession alone.Thus, numerical simulations are required to make quantitative predictions.We confirm that the timescale of the corresponding RV modulations are consistent with the ZKL timescale T ZKL , which is roughly ∼ 100P out for our fiducial cases for Gaia BH1 (P out ∼ 190 days) and Gaia BH2 (P out ∼ 1300 days).Due to the rather sporadic and abrupt change of the RV semi-amplitude due to the ZKL oscillations, we may be able to detect the signatures of the long-term RV modulation depending on the observational phase.
We have demonstrated the feasibility of detecting an inner BBH from RV follow-ups of star-BH binary candidates, if some of them are indeed star-BBH triples.We studied the presently available best targets, Gaia BH1 and Gaia BH2, as a proof-of-concept, but found that future monitoring of Gaia BH1 may indeed detect an inner BBH within a reasonable timescale.Current gravitationalwave data seem to point to a possible mass gap of BHs between 3 ∼ 6M ⊙ , but its reality is still controversial.Thus, an inner BBH in Gaia BH1, if detected, has a huge impact on the formation scenarios of BHs in general.
The three observable signatures of the RV modulations of the tertiary discussed in the above summary are quite generic, and can be applied to more candidates from future Gaia data in a straightforward manner.We also mention that this method is applicable to tertiary pulsar -BBH triple systems using the pulsar timing analysis, instead of RV monitoring (Hayashi & Suto 2021).Furthermore, the long-term inclination modulation is sufficiently large to be identified from the orbital parameters determined by astrometry in ∼ 10 years from now.For instance, Liu et al. (2022) propose that the secular eccentricity variations induced by the apsidal precession resonances are potentially detectable with Gaia astrometry observations.
At this point, it is quite uncertain if such star-BBH and even pulsar-BBH triples exist within our reachable horizon.Nevertheless, we would like to conclude by referring to a universal principle that "everything not forbidden is compulsory" (White 1939;Thorne & Zytkow 1975;Sagan 1985).
Figure 3. Examples of simulated radial velocity curves, and the residuals after removing the best-fit Kepler motion.The fittings are performed with RadVel (Fulton et al. 2018).We assume an equal-mass inner binary (m 1 = m 2 = m c /2) for both systems.The left and right plots assume Gaia BH1 with P in = 10 days, and Gaia BH2 with P in = 50 days.The top panels shows the total RV curve with e in = 0 (red) and e in = 0.2 (blue), and the middle and bottom panels plot the corresponding RV residuals, respectively.
Figure4.Mass ratio of inner BBH q 21 (dashed) and the corresponding planet mass M pl (solid) in RV signal semi-amplitude K s -period P s plane for Gaia BH1.We note that K s and P s are interpreted as RV semi-amplitude and period induced by a planet for the planet interpretation, while they are interpreted as K s = K short (1 − e obs ) −7/2 and P s = P in /2 for the BBH interpretation.

Figure 8 .
Figure8plots the outer inclination i out against i mut for various values of ξ.The plot indicates that i out becomes asymptotically close to i mut as ξ increases.We intend to show Figure8as a generic

Figure 9 .
Figure 9. Relations between i out and I los (top), and i mut and I los (bottom).The left and right panels correspond to Gaia BH1 and Gaia BH2, respectively.Shaded regions in the top panels are excluded due to the observed value of the inclination I obs (see Table1).The bottom panels plot I los against i mut for given values of ξ.The region between two lines is permitted for each ξ, and the shaded regions correspond to ξ > 1.2 for Gaia BH1 and ξ > 0.9 for Gaia BH2, which are excluded from dynamical stability viewpoint.

Figure 11 .
Figure 11.Example evolution of simulated RV semi-amplitudes and the mutual inclination angles for the initial values of i mut = 20 • (top), 60 • (middle), and 90 • (bottom).The left and right panels correspond to Gaia BH1 (P in = 10 days and I los = 120 • ) and Gaia BH2 (P in = 50 days and I los = 30 • ), respectively.For simplicity, we assume the equal-mass inner binary (q 12 = 1) in an initially circular orbit (e in = 0).For reference, the range of the analytic estimate of ∆ K , equation (31), is plotted in magenta regions, and the expected semi-amplitude change, from equation (27), is plotted as dotted green curves.Equation (10) diverges for i mut = 90 • , and we do not show the analytic predictions in those cases.

Table 1 .
Best-fit parameters for Gaia BH1 and BH2 systems