Uncovering a hidden black hole binary from secular eccentricity variations of a tertiary star

We study the dynamics of a solar-type star orbiting around a black hole binary (BHB) in a nearly coplanar system. We present a novel effect that can prompt a growth and significant oscillations of the eccentricity of the stellar orbit when the system encounters an “ apsidal precession resonance, ” where the apsidal precession rate of the outer stellar orbit matches that of the inner BHB. The eccentricity excitation requires the inner binary to have a nonzero eccentricity and unequal masses, and can be created even in noncoplanar triples. We show that the secular variability of the stellar orbit ’ s apocenter, induced by the changing eccentricity, could be potentially detectable by Gaia . Detection is favorable for BHBs emitting gravitational waves in the frequency band of the Laser Interferometer Space Antenna, hence providing a distinctive, multimessenger probe of the existence of stellar-mass BHBs in the Milky Way.

BHB progenitors are expected to be numerous but remain undetected as an abundant population in our Universe.Searching for these inspiraling BHBs is of great importance to understand the origin of gravitational-wave (GW) sources.If the BHBs are not accreting, these quiescent sources can be detected via GWs.Since BHBs in the inspiral phase are still far from merger, the associated GWs are in the low-frequency band that can be explored by future spaceborne GW observatories, such as LISA [62], TianQin [63], Taiji [64], B-DECIGO [65], Decihertz Observatories [66], and TianGO [67].Also, since a significant fraction of compact BHBs may be members of hierarchical systems [68][69][70], the motion of a nearby visible object (such as a star or a pulsar) can be used to search for BHBs.In this scenario, the inner BHB can perturb the outer orbit, either inducing short-term orbital oscillations (tertiary orbit becomes quasi-Keplerian), or causing long-term oscillations of the eccentricity and the orientation of the angular momentum when the tertiary orbit is highly inclined [71][72][73].If the tertiary object * Electronic address: bin.liu@nbi.ku.dk is bright enough, radial velocity measurements could be used to determine the short and long term deviations from a Keplerian orbit [74][75][76].
Currently, observations show that most triple-star systems are less inclined or nearly coplanar [77][78][79].Here, we consider a solar-type star orbiting around a BHB and study the secular evolution of the stellar orbit when the triple system is coplanar.We show that the outer binary may experience an eccentricity growth driven by the "apsidal precession resonance" [80,81].Compared to the complex evolution of the orientation of the angular momentum (i.e., precession or nutation), the secular change of the eccentricity of the outer stellar orbit could provide distinctive evidence to reveal the presence of a BHB.

II. APSIDAL PRECESSION RESONANCE
We consider an inner BHB with masses m 1 , m 2 , and a solar-type star (m = 1 M ) that moves around the center of mass of the inner bodies.The reduced mass for the inner binary is µ in ≡ m 1 m 2 /m 12 , with m 12 ≡ m 1 + m 2 .Similarly, the outer binary has µ ≡ (m 12 m )/(m 12 + m ).The semi-major axes and eccentricities are denoted by a in , a and e in , e , respectively.The orbital angular momenta of two orbits are thus given by When the triple system is less inclined or nearly coplanar, the ZLK oscillations are not allowed to occur, but a significant eccentricity excitation of the inner binary may still be induced [82,83].A secular, "apsidal precession resonance" plays a dominant role if the total apsidal precession of the inner binary matches the precession rate of the outer binary [80,81].Precession of both inner and outer binaries is driven by Newtonian and general relativistic (GR) effects .Such resonance allows efficient angular momentum exchange between the inner and outer binaries.
Here, we extend our previous studies to an "inverse" in , e > 0) and the maximum change of eccentricity ∆e max of the outer stellar orbit as a function of the semimajor axis ain.The cross-hatched region corresponds to dynamically unstable triple systems [85].The parameters are the same as in the left panel, except ain is relaxed to a range of values.The labeled GW frequency is the peak frequency at pericenter [47].The numerical results are obtained by integrating the SA secular equations over several timescales (as labeled), and the analytical result is given by the energy and angular momentum conservation laws (dashed blue line).Right panel: ∆e max induced by the apsidal precession resonance in the (m2/m1 − ain) plane.We use the same parameters as in the middle panel, taking into account the full range of mass ratios of the BHB.The three black contours (solid, dashed and dot-dashed) specify ∆e max = 0.1, 0.2 and 0.5, respectively.The gray crosses indicate the significant change of e that leads to instability.secular problem, and address the question of how the apsidal precession resonance modifies the eccentricity evolution of the tertiary for the first time.Since we are interested in the long-term orbital evolution, we adopt the single-averaged (SA; only averaging over the inner orbital period) secular equations of motion, taking into account the contributions from the Newtonian effect up to the octupole level of approximation and the leading order GR effect in both inner and outer orbits.The explicit SA equations are provided in [52,84].
In Figure 1, the left panels show the evolution of eccentricities of both inner and outer orbits.Starting with the eccentric inner binary, we see that e can be excited from the circular orbit and undergoes oscillations.In the middle panel, we find that a large fraction of systems can develop eccentricities (0.04 a in /AU 0.08), and an evident peak, ∆e max 0.8, can be resolved when the evolution is sufficiently long ( 10yrs).The right panel illustrates the level of e −excitation in the (m 2 /m 1 −a in ) plane.The eccentricity of the stellar orbits can be excited for a in 0.04 AU, and the systems with smaller mass ratio tend to have larger ∆e max .This is because the evolution of e is determined by the octupole-order secular interactions, which can be quantified by terms proportional to [45]: ( We see that the eccentricities of some outer orbits reach significantly large ∆e max close to unity, leading to un-bound orbits. For coplanar ( L = L ), non-dissipative (no gravitational radiation) systems, the secular dynamics can be understood analytically.When e in , e 1, the evolution of e in and e is governed by the linear Laplace-Lagrange equations [86,87].If we define the complex eccentricity variables as E in ≡ e in exp(i in ) and E ≡ e exp(i ), where in , are the longitude of pericenter of the inner and outer orbits, then the evolution equations are reduced to with where ) )] are the mean motion and the GR-induced pericenter-precession frequency of the inner (outer) binary for m 12 m , respectively.The maximum change of eccentricity as a function of the inner binary semimajor axis.We consider the inner BHB with the total mass m12 = 50 M and mass ratio m2/m1 = 0.2 (same as the middle panel of Figure 1), and the stellar orbits with different orbital period (as labeled).The results (solid lines) are obtained by the numerical integrations of SA equations for 25 yrs (left panel), 150 yrs (middle panel) and 300 yrs (right panel), and the results (dashed lines) are all from 10 yrs.To compare with the fiducial example (black lines), we fix the initial e 0 (e 0 in ) in the upper (lower) panels (as labeled).
Starting with e in = e 0 in , e = 0 at t = 0, Equation ( 2) can be solved to determine the time evolution of e (t) [see 80], which oscillates between 0 and e max , where Clearly, e max attains its peak value when ω in = ω occurs, and then we have Note that the linear theory is valid for the low−e systems.For the example in Figure 1, non-zero e 0 in may lead to unphysically large e peak ; however, Equations (7) and (8) are useful in the sense that we can expect: i) e −excitation appears when ω in ω (see the middle panel of Figure 1; the resonance occurs with finite e ); ii) e −excitation becomes stronger with increasing e in (see also Figure 2

below).
For finite eccentricities, Equation (2) breaks down.However, in the case of exact coplanarity, the maximum eccentricity of a triple can still be calculated algebraically, using energy and angular momentum conservations [88].This method works well for two orbits with arbitrary eccentricities, but it cannot show the time evolution and resonance features.A full derivation can be found in [81] and the solution is presented in Figure 1.
Figure 2 shows e −excitation with the arbitrary initial eccentricities for coplanar triples.We again consider the fiducial example, in which the inner BHB has the total mass m 12 = 50 M with mass ratio m 2 /m 1 = 0.2.We choose three values of the semimajor axis of the outer orbit (a ) and consider a range of a in that satisfies the stability criterion.To illustrate the role of the eccentricity, we start with the same initial configurations, namely, the argument of periapse, the longitude of ascending nodes and the true anomaly of the outer orbit are set to be the same at t = 0.Each system is evolved for a long timescale (to achieve the highest value of e ) and a short timescale (10 yrs), respectively.The maximum change of e is picked only for the system remains stable.
In the upper left panel of Figure 2, we see that all the initial circular outer orbits can become eccentric when the resonance occurs, and the maximum change ∆e max grows as e 0 in increases.In the lower left panel, e −excitation can still occur for the initially eccentric outer orbits.However, due to the stability, only a fraction of systems (with a narrow range of a in ; when e 0 0.6) may undergo significant eccentricity oscillations.Note that for the wider stellar orbits (P = 50, 100 days), as shown in the middle and right panels, the system has to be evolved for a sufficiently long time.This is because the timescale of e −excitation is of the order of [89] T e ein,e 1 When the inner and outer orbits are mutually inclined, no simple analytical result can be derived, and the longterm evolution of the outer orbits can only be studied numerically.Previous studies [71][72][73] showed that in inclined triple systems, the eccentricity of the outer orbit can oscillate moderately and the angular momentum (L ) undergoes nodal precession/nutation around the inner one (L in ).Moreover, L might experience a flip if the tertiary is a test particle.
Figure 3 presents the results of the triples with a series of initial inclinations.We see that in the upper panels, regardless of the values of the initial inclinations, e −excitation due to resonance always occurs, and the resonance location shifts when I 0 changes (I 0 is the initial inclination angle between L in and L ).In the lower panels, we find that the inclination varies for a wider range of a in compared to the change of e .In particular, ∆I always undergoes an additional excitation when ∆e max approaches the peak value for the inclined systems.

III. RESONANCE IN STELLAR ORBITS AND DETECTABILITY
We now focus on dependence of the change in eccentricity on parameters of the outer stellar orbit, considering the inner BHB is a LISA source.To explore the observability of this effect, we consider the detectability of a maximum eccentricity change, ∆e peak , within a certain timescale.
We initialize systems with e 0 in = 0.9 and e 0 = 0 in a coplanar configuration.For the inner BHB, we choose three values of the total mass (m 12 /M = 20, 50, 100), and allow the mass-ratio range to take on all values such that both masses are consistent with being a BH, m 1 , m 2 5M .To develop the joint detection with GW detectors in the coming future, we focus on the BHBs radiating GWs in the LISA frequency band.Thus, the semimajor axis a in is chosen from a uniform distribution that satisfies f GW ≥ 10 −4 Hz and with merger time (due to GW emission) larger than 10 3 yrs [90].Then, for the (outer) stellar orbits, we sample the semimajor axis a from a range of P ∼ [1, 180] days, considering only systems that are dynamically stable.Each system evolves to 10, 30 or 100 yrs, using the SA equations of motion.The maximum change of the e is recorded if the star remains gravitationally bound and stable during the evolution.Finally, since the orbital evolution relies on the initial geometry of the triples, to cover all possibilities, we randomly sample the longitude of pericenter of the inner orbit and the true anomaly of the outer orbit.Each triple system is evolved with 100 different initial geometries.
For a given set of parameters, the criterion of apsidal precession resonance (ω in = ω ) provides a good estimate for the resonance radius: resonance occurs at the location a = a Res for a given a in .In Figure 4, the upper panel clarifies the resonance locations when m 12 = 50M .The region of interest where the outer orbit potentially undergoing e −excitation due to the apsidal precession resonance is located within a wide range of a (or P ).
The lower panels of Figure 4 present the results of the averaged maximum changes of the outer eccentricity (∆e peak ), over all 100 runs, as a function of a for different masses of BHBs.A remarkable eccentricity excitation is achieved for all tested values of the BHB mass.For the m 12 = 20M case, the maximum eccentricity growth is found for the smallest a and decreases for larger a .For the BHBs (m 12 /M = 50, 100), the perturbation becomes stronger, and the induced peak eccentricity can be so high (because of e peak ∝ µ in ) that some of the systems with small a can become unbound.The peak of ∆e peak appears to shift to larger a for the larger BHB masses.The decrease of ∆e peak with increasing a is the result of the secular timescale of generating the eccentricity growth (See Figure 3).Note that, in principle, the resonance can occur as a ≤ a Res (the vertical blue line; P 160 days).As shown, more systems can have larger ∆e peak when the evolving time is longer (See Figure 2).We find that e −excitation is mainly contributed by the inner BHBs with relatively small mass ratios (m 2 /m 1 0.5).
The secular variability of the orbital eccentricity in- where the apsidal precession resonance occurs.We consider the BHB with m12 = 50 M , m2/m1 = 0.2 and the initial e 0 in = 0.9, e 0 = 0.The blue region is given by Equations 3-4 including the finite e 0 in and e (as labeled).The dashed lines characterize the frequency of GWs emitted by the inner BHB.Lower panels: Maximum of ∆e max (i.e., ∆e peak ) as a function of a for different m12.The solid and dashed lines are obtained by the numerical integrations for 10, 30 and 100 yrs, respectively.The vertical blue line to a = a Res at e = 0 via fGW = 10 −4 Hz.The cross-hatched region in the bottom panel indicates that the systems eventually become unbound when evolving for a longer timescale.duces a change in the projected orbit that may be probed via astrometric monitoring with surveys such as Gaia [91].To determine detectability with Gaia, we compute a signal-to-noise ratio (SNR) for astrometric detection, ρ = θ signal /θ Gaia .Assuming the signal can be wellapproximated by the change in apocenter, the maximum such signal for a system at distance D is approximated by To compute noise, we follow [92], which draws on [93][94][95][96][97] to evaluate the Gaia astrometric precision for a singlescan.Here, we assume that our source is a solar-type star with absolute V-band magnitude of −26.8, and V − I c = 0.688 [98].Following [96], the single-scan precision is computed from the end-of-mission, sky-averaged parallax  10)-( 11) and ∆e peak is given by the data from Figure 4 for 10 yrs.The subfigure in the top panel illustrates an example of detectability where the source is at D = 900 pc, and the shaded region corresponds to the systems with θ signal ≥ θGaia.The white dashed lines specify the region with ρ > 2.
uncertainty θ eom (D), as, The pre-factors account for an average 140 Gaia visits over 10 years, a geometrical averaging factor of 2.15, and a contingency margin of 1.1 [96].Note that the contingency margin is 1.1 instead of the value of 1.2 chosen in [96].This is based on newest information from Gaia EDR3 uncertainties as expressed in Section 1 of [99].To compute the end-of-mission astrometric precision, we use the most up-to-date fitting formula from the Gaia expected science performance document [99], Z ≡ max 10 0.
Here the 0.527 prefactor is for a 10 year mission (referred to as Gaia DR5 in [99]), the conversion between Gaia-Gand V-magnitudes, in the last line above, is given in Table A2 of [100], and m V (D) = −26.8+ 2.5 log 10 [(D/AU) 2 ] is the apparent magnitude of a sun-like star at distance D in Astronomical Units (AU). Figure 5 summarizes the optimal-characteristic SNR as a function of distance to the source and outer orbital period, for three different inner binary masses.Here, we use the data from Figure 4 that gives the maximum change of eccentricity for a range of a over 10 yrs (for m 12 = 100 M , we assume ∆e peak 0.3 due to the long term instability).We see that overall SNR improves D decreases, and the boundary of detectability is marked at ∼ 10 3 pc given by ρ = 2 [97].Two peaks of high SNR are the results of the pure e −excitation (Figure 4) and the enhancement of θ signal from wide binaries.Although the detectability here is evaluated based only on ∆e peak , the outer pericenter argument can vary in the actual detection [75].When the triple system is inclined, a combination of predicted changes in orbital inclination and line-of-sight orbital projection can also increase or decrease our estimate here.Our results are characteristic of the optimal astrometric signatures over the course of the Gaia mission.Note that the average time between visits is ∼ 26 days, hence, the shorter period systems may be sampled at suborbital frequencies.However, because this is a secular effect, even a longer than orbital sampling rate could probe the longer-timescale secular change of the projected orbit.Follow-up analysis must simulate mock systems to determine the required orbital sampling rate and SNR needed to reliably detect this effect.Further work is also required to determine in which cases the current Gaia pipeline would flag such secularly evolving systems as binaries with an unseen companion, or miss identify them due to difficulty in fitting to a Keplerian orbit.

IV. SUMMARY AND DISCUSSION
We have studied a novel secular dynamical effect of a solar-type star around a compact BHB in a nearly coplanar triple configuration.We point out that the stellar orbit may undergo significant oscillations of the orbital eccentricity if the system satisfies an "apsidal precession resonance".Starting with an eccentric inner BHB, the outer eccentricity can be excited due to the resonance and the enhancement becomes stronger as the inner eccentricity increases.This effect, which can drive the eccentricity of the outer orbits close to unity (even unbound the binary), is overlooked in all previous studies and may have applications for other types of systems, such as the planet around a binary star and star/compact object around a supermassive BH binary with/without a gaseous disc.
Note that the apsidal precession resonance allows momentum exchange to occur efficiently between the inner outer orbits, leading to a transfer of eccentricity.One case that we do not include in our study is the triple systems with e 0 in = 0 and e 0 = 0.9.In this situation, the inner (outer) binary is expected to become (less) eccentric as the resonance occurs.
Figure 6 shows the parameter space where the apsidal precession resonance can play a role, on the basis of ω in (e in ) = ω (e ).Compared to the fiducial cases (left panel), resonance region for the initially stable systems in the right panel is shifted to larger a (or P ).Therefore, the corresponding semimajor axis of the inner BHB (a in ) becomes larger, leading to a much lower GW frequency range (i.e., out of the LISA band).FIG.6: Similar to the top panel of Figure 4, but we include two combinations of initial eccentricities (as labeled).Here, the mass of BHB is set to be m12 = 50 M and the mass ratio is m2/m1 = 0.2.FIG.7: Similar to the middle panel of Figure 1, but we include the change of eccentricities of both inner and outer binaries.
Figure 7 presents the change of eccentricities due to resonance for both inner and outer orbits, taking into account the examples identified in the Figure 6.In the left panel (fiducial case), we see that the inner eccentricity decreases (i.e., ∆e max in 0.1) when the outer eccentricity is excited.In the right panel, we find that the eccentricities evolve in an opposite way.The inner binary can efficiently gain some eccentricity from the outer binary during the resonance (the peak value is about ∆e max in ∼ 0.5), while the the outer eccentricity can decrease by a factor of 0.2.
Since we are interested in the LISA source and the evidence change of e within ∼ 10 years, we do not include this type of system (with e in = 0 and e 0 = 0.9) in our main study.
The formation of compact massive BHB+star systems may be challenging.The progenitor stars of these BHBs usually expand hundreds or thousand of solar radii throughout their evolution, likely dynamically interacting with the stellar tertiary.However, some lowmetallicity massive stellar binaries might remain compact throughout their evolution [12], allowing for dynamicallystable compact triples [101].Alternatively, a BHB can be formed first and eventually capture a long-lived low-mass star.Note that our analysis is not restricted to any specific formation scenarios, and can be adapted to other types of systems by applying scaling relations.
We find that the secular variability of the stellar orbit's apocenter induced by the changing eccentricity is detectable by Gaia for inner BHBs emitting GWs in the LISA frequency range.Assuming that the formation and merger rates of BHBs are in equilibrium, we expect to have hundreds of BHBs in the LISA band in the Galaxy based on the LIGO detection rate [102][103][104][105][106][107].Since our proposed secular variability can only be resolved by Gaia within several kpc, the expected number of sources that LISA and Gaia could see becomes ∼a few. the actual number of LISA sources accompanied by a in our Galaxy is quite uncertain, identifying the secular motion of stellar orbits in current (Gaia Data Release 3; [91]) and future Gaia data is timely.
Our proof-of-concept calculations demonstrate that the long-term evolution of eccentricity of a nearby stellar orbit can serve as a distinctive imprint of such an unseen binary companion.Precise measurements of secular variability are therefore an independent approach to reveal hidden BHBs, in addition to GW detection.We emphasize that the inner binary systems which generate Gaia-detectable variations in the orbits of their stellar tertiary also emit GWs in the LISA band.Therefore, Gaia may provide candidates LISA sources before LISA flies (planned for the 2030's).In this sense, a joint detection with Gaia and LISA [108][109][110] would be a unique multi-messenger tool to understand the evolution, fate and configurations of compact BHBs.

m2 / m1 = 0. 2 ainFIG. 1 :
FIG.1: Apsidal precession resonance in coplanar triple systems where the outer orbital period (P ) is set to be 15 days.All the results are obtained by integrating the SA secular equations including GR effects (but without GW emission).Left panel: Evolution examples of the orbital eccentricities in the inner (black) and outer (red) binaries, with m12 = 50M and initial e 0 in = 0.9 and e 0 = 0. Middle panel: The ratio of apsidal precession rates (Equations 3-4; including the dependence of finite eccentricities, i.e., e 0 in , e > 0) and the maximum change of eccentricity ∆e max of the outer stellar orbit as a function of the semimajor axis ain.The cross-hatched region corresponds to dynamically unstable triple systems[85].The parameters are the same as in the left panel, except ain is relaxed to a range of values.The labeled GW frequency is the peak frequency at pericenter[47].The numerical results are obtained by integrating the SA secular equations over several timescales (as labeled), and the analytical result is given by the energy and angular momentum conservation laws (dashed blue line).Right panel: ∆e max induced by the apsidal precession resonance in the (m2/m1 − ain) plane.We use the same parameters as in the middle panel, taking into account the full range of mass ratios of the BHB.The three black contours (solid, dashed and dot-dashed) specify ∆e max = 0.1, 0.2 and 0.5, respectively.The gray crosses indicate the significant change of e that leads to instability.
FIG.2:The maximum change of eccentricity as a function of the inner binary semimajor axis.We consider the inner BHB with the total mass m12 = 50 M and mass ratio m2/m1 = 0.2 (same as the middle panel of Figure1), and the stellar orbits with different orbital period (as labeled).The results (solid lines) are obtained by the numerical integrations of SA equations for 25 yrs (left panel), 150 yrs (middle panel) and 300 yrs (right panel), and the results (dashed lines) are all from 10 yrs.To compare with the fiducial example (black lines), we fix the initial e 0 (e 0 in ) in the upper (lower) panels (as labeled).

FIG. 3 :
FIG.3: Similar to Figure2, but we set the initial eccentricities as e 0 in = 0.9 and e 0 = 0, and consider different initial inclinations (as labeled; I the inclination angle between Lin and L ).The changes of e and I are shown in the upper and lower panels, respectively.

FIG. 4 :
FIG.4: Top panel: Parameter space in [ain−a or (P )] plane, where the apsidal precession resonance occurs.We consider the BHB with m12 = 50 M , m2/m1 = 0.2 and the initial e 0 in = 0.9, e 0 = 0.The blue region is given by Equations 3-4 including the finite e 0 in and e (as labeled).The dashed lines characterize the frequency of GWs emitted by the inner BHB.Lower panels: Maximum of ∆e max (i.e., ∆e peak ) as a function of a for different m12.The solid and dashed lines are obtained by the numerical integrations for 10, 30 and 100 yrs, respectively.The vertical blue line to a = a Res at e = 0 via fGW = 10 −4 Hz.The cross-hatched region in the bottom panel indicates that the systems eventually become unbound when evolving for a longer timescale.

FIG. 5 :
FIG.5: Detectability of the change of the outer stellar eccentricity in the (P − D) plane with Gaia, which is evaluated by Equations (10)-(11) and ∆e peak is given by the data from Figure4for 10 yrs.The subfigure in the top panel illustrates an example of detectability where the source is at D = 900 pc, and the shaded region corresponds to the systems with θ signal ≥ θGaia.The white dashed lines specify the region with ρ > 2.