A strategy to search for an inner binary black hole from the motion of the tertiary star

TL;DR

This paper proposes a method to detect inner binary black holes in triple star systems through precise radial velocity measurements of the outer star, offering a new way to identify potential gravitational-wave progenitors.

Contribution

It introduces a novel methodology combining approximate and numerical RV analysis to search for inner BBHs in triple systems, including application to real and simulated data.

Findings

Method can constrain parameters of inner BBHs from RV data.

Numerical simulations show feasibility of detecting inner BBHs.

Potential to identify many hidden black hole binaries in the galaxy.

Abstract

There are several on-going projects to search for stars orbiting around an invisible companion. A fraction of such candidates may be a triple, instead of a binary, consisting of an inner binary black hole (BBH) and an outer orbiting star. In this paper, we propose a methodology to search for a signature of such an inner BBH, possibly a progenitor of gravitational-wave sources discovered by {\it LIGO}, from the precise radial velocity (RV) follow-up of the outer star. We first describe a methodology using an existing approximate RV formula for coplanar circular triples. We apply this method and constrain the parameters of a possible inner binary objects in 2M05215658+4359220, which consists of a red giant and an unseen companion. Next we consider co-planar but non-circular triples. We compute numerically the RV variation of a tertiary star orbiting around an inner BBH, generate mock RV curves, and examine the feasibility of the BBH detection for our fiducial models. We conclude that the short-cadence RV monitoring of a star-BH binary provides an interesting and realistic method to constrain and/or search for possible inner BBHs. Indeed a recent discovery of a star--BH binary system LB-1 may imply that there are a large number of such unknown objects in our Galaxy, which are ideal targets for the methodology proposed here.

Tables

Table 1. Table 1: Best-fit parameters for the binary system 2M05215658+4359220 (Thompson et al., 2019 )

parameter	value	meaning
$P_{\mathrm{out}}$	$83.205\pm 0.064\mathrm{days}$	orbital period
$m_{\mathrm{co}}$	${3.3}_{-0.7}^{+2.8}{M}_{\odot }$	mass of an unseen companion
$m_{\mathrm{giant}}$	${3.2}_{-1.0}^{+1.0}{M}_{\odot }$	mass of a red giant
$e_{\mathrm{out}}$	$0.00476\pm 0.00255$	eccentricity
${\omega }_{\mathrm{out}}$	$197.13\pm 32.07\mathrm{deg}$	argument of pericenter
$\sin i$	${0.97}_{-0.12}^{+0.03}$	inclination of the orbital plane
$R_{\mathrm{giant}}$	${30}_{-6}^{+9}{R}_{\odot }$	radius of a red giant

Table 2. Table 2: Initial values of parameters for our fiducial triple.

parameter	initial value
semi-major axis $a_{\mathrm{in}}$	$1.0\mathrm{au}$
semi-major axis $a_{\mathrm{out}}$	$4.0\mathrm{au}$
mass of the primary $m_1$	$10M_{\odot }$
mass of the secondary $m_2$	$10M_{\odot }$
mass of the tertiary $m_{\ast }$	$1M_{\odot }$
inclinations $I_{\mathrm{in}},I_{\mathrm{out}}$	$90\mathrm{deg}$
arguments of pericenter ${\omega }_{\mathrm{in}},{\omega }_{\mathrm{out}}$	$0\mathrm{deg}$
longitudes of ascending node ${\Omega }_{\mathrm{in}},{\Omega }_{\mathrm{out}}$	$0\mathrm{deg}$
true anomaly $f_{\mathrm{in}}$	$30\mathrm{deg}$
true anomaly $f_{\mathrm{out}}$	$120\mathrm{deg}$

Table 3. Table 3: Predicted values for orbital motions calculated from equation ( 1 ) and Table 2 .

parameter	value
orbital period $P_{\mathrm{in}}$	$81.7\mathrm{days}$
orbital period $P_{\mathrm{out}}$	$638\mathrm{days}$
semi-amplitude of Keplerian RV $K_0$	$65.0\mathrm{km}/s$
variation period (2nd term in eq.(1)) $P_{\mathrm{var1}}$	$50.5\mathrm{days}$
variation period (3rd term in eq.(1)) $P_{\mathrm{var2}}$	$43.6\mathrm{days}$
half an inner period $P_{\mathrm{in}}/2$	$40.9\mathrm{days}$
characteristic semi-amplitude of variations $K_{\mathrm{BBH}}$	$130.1m/s$

Table 4. Table 4: Best-fit parameters estimated with RadVel for our fiducial coplanar circular triple. Tildes on each variable indicates that they are fitted from the total RV curve over 800 days.

parameter	value
orbital period ${\tilde{P}}_{\mathrm{out}}$	$619.3041\mathrm{days}$
time of conjunction ${\tilde{T}}_{j,\mathrm{out}}$	$566.12833\mathrm{days}$
eccentricity ${\tilde{e}}_{\mathrm{out}}$	$0.01565794$
pericenter argument ${\tilde{\omega }}_{\mathrm{out}}$	$-1.007588\mathrm{rad}$
the semi-amplitude of RV ${\tilde{K}}_0\sin \tilde{I}$	$65931.314m/s$

Table 5. Table 5: Initial values of parameters for a fiducial binary with an S-type circumbinary planet orbiting around the star.

parameter	initial value
semi-major axis $a_{\mathrm{pl}}$	$0.27\mathrm{au}$
semi-major axis $a_{\mathrm{out}}$	$4.0\mathrm{au}$
mass of the unseen companion $m_{\mathrm{BH}}$	$20M_{\odot }$
mass of the star $m_{\ast }$	$1M_{\odot }$
mass of the planet $m_{\mathrm{pl}}$	$2.7M_J$
inclinations $I_{\mathrm{pl}},I_{\mathrm{out}}$	$90\mathrm{deg}$
arguments of pericenter ${\omega }_{\mathrm{pl}},{\omega }_{\mathrm{out}}$	$0\mathrm{deg}$
longitudes of ascending node ${\Omega }_{\mathrm{pl}},{\Omega }_{\mathrm{out}}$	$0\mathrm{deg}$
true anomaly $f_{\mathrm{pl}}$	$30\mathrm{deg}$
true anomaly $f_{\mathrm{out}}$	$120\mathrm{deg}$

Equations

$$V_{\mathrm{RV}}$$

$$K_0\equiv \frac{m_1+m_2}{m_1+m_2+m_{\ast }}{a}_{\mathrm{out}}{\nu }_{\mathrm{out}},$$

$$K_{\mathrm{BBH}}\equiv \frac{m_1{m}_2}{(m_1+m_2{)}^2}\sqrt{\frac{m_1+m_2+m_{\ast }}{m_1+m_2}}{\left(\frac{a_{\mathrm{in}}}{a_{\mathrm{out}}}\right)}^{3.5}{K}_0$$

$$\alpha >{\left(\frac{a_{\mathrm{in}}}{a_{\mathrm{out}}}\right)}_{\mathrm{crit}}\equiv \frac{1-e_{\mathrm{out}}}{2.8}{\left(\frac{(1+m_3/(m_1+m_2))(1+e_{\mathrm{out}})}{\sqrt{1-e_{\mathrm{out}}}}\right)}^{-\frac{2}{5}}\approx \mathrm{0.270.}$$

$$t_{\mathrm{GW}}=\frac{5}{256}\frac{c^5{a}_{\mathrm{in}}^4}{G^3{m}_{12}^2\mu }\approx 1.88\times 10^{11}{\left(\frac{P_{\mathrm{in}}}{\mathrm{day}}\right)}^{8/3}{\left(\frac{m_{12}}{M_{\odot }}\right)}^{-5/3}\mathrm{yrs},$$

$$P_{\mathrm{pl}}\approx 11\mathrm{days}{\left(\frac{a_{\mathrm{pl}}}{30{\mathrm{R}}_{\odot }}\right)}^{3/2}\gtrsim 11\mathrm{days}.$$

$$P_{var1}$$

$$P_{var3}$$

$$\mathcal{H}=\frac{1}{2}\frac{{\mathbf{p}}_{\mathrm{in}}^2}{\mu }+\frac{1}{2}\frac{{\mathbf{p}}_{\mathrm{out}}^2}{{\mu }_{\ast }}-G\frac{m_1{m}_2}{r_{\mathrm{in}}}-G\frac{m_{12}{m}_{\ast }}{r_{\mathrm{out}}}-G\frac{m_1{m}_2{m}_{\ast }}{r_{\mathrm{out}}{m}_{12}}{\rho }^2\frac{1}{2}(3{\cos }^{2}S-1),$$

$$\mu \equiv \frac{m_1{m}_2}{m_{12}}\mathrm{and}{\mu }_{\ast }\equiv \frac{m_{12}{m}_{\ast }}{m_{12}+m_{\ast }}.$$

$${\dot{\mathbf{p}}}_{\mathrm{in}}=-\frac{\partial \mathcal{H}}{\partial {\mathbf{r}}_{\mathrm{in}}}\mathrm{and}{\dot{\mathbf{p}}}_{\mathrm{out}}=-\frac{\partial \mathcal{H}}{\partial {\mathbf{r}}_{\mathrm{out}}}.$$

$${\ddot{\mathbf{r}}}_{\mathrm{in}}=-G\frac{m_{12}}{r_{\mathrm{in}}^3}{\mathbf{r}}_{\mathrm{in}}+G\frac{m_{\ast }}{r_{\mathrm{out}}^3}(3\rho \cos S{\mathbf{r}}_{\mathrm{out}}-{\mathbf{r}}_{\mathrm{in}})$$

$${\ddot{\mathbf{r}}}_{\mathrm{out}}=-G\frac{m_{12}+m_{\ast }}{r_{\mathrm{out}}^3}\left[\left(1+\frac{\mu }{m_{12}}\frac{{\rho }^2}{2}(-3+15{\cos }^{2}S)\right){\mathbf{r}}_{\mathrm{out}}-\frac{\mu }{m_{12}}3\rho \cos S{\mathbf{r}}_{\mathrm{in}}\right],$$

$${\mathbf{r}}_{\mathrm{in}}^{(0)}=a_{\mathrm{in}}\cos S{\hat{\mathbf{r}}}_{\mathrm{out}}+a_{\mathrm{in}}\sin S\hat{\mathbf{u}}\mathrm{and}{\mathbf{r}}_{\mathrm{out}}^{(0)}=a_{\mathrm{out}}{\hat{\mathbf{r}}}_{\mathrm{out}},$$

$${\nu }_{\mathrm{in}}\equiv \sqrt{\frac{G{m}_{12}}{a_{\mathrm{in}}^3}}\mathrm{and}{\nu }_{\mathrm{out}}\equiv \sqrt{\frac{G(m_{12}+m_{\ast })}{a_{\mathrm{out}}^3}}.$$

$${\hat{\mathbf{r}}}_{\mathrm{out}}=\cos \theta \mathbf{i}+\sin \theta \mathbf{j}\mathrm{and}\hat{\mathbf{u}}=-\sin \theta \mathbf{i}+\cos \theta \mathbf{j},$$

$$\ddot{X}=-\frac{G{m}_{12}}{a_{\mathrm{out}}^2}\left[\left(1+\frac{3}{4}{\alpha }^2\frac{\mu }{m_{12}}\right)\cos \theta +\frac{9}{4}{\alpha }^2\frac{\mu }{m_{12}}\cos 2S\cos \theta +\frac{3}{2}{\alpha }^2\frac{\mu }{m_{12}}\sin 2S\sin \theta \right]$$

$$\ddot{Y}=-\frac{G{m}_{12}}{a_{\mathrm{out}}^2}\left[\left(1+\frac{3}{4}{\alpha }^2\frac{\mu }{m_{12}}\right)\sin \theta +\frac{9}{4}{\alpha }^2\frac{\mu }{m_{12}}\cos 2S\sin \theta -\frac{3}{2}{\alpha }^2\frac{\mu }{m_{12}}\sin 2S\cos \theta \right],$$

$$\displaystyle\left(\begin{array}[]{c}X\\ Y\\ \end{array}\right)=\left(\begin{array}[]{c}a_{*}\cos{\theta}\\ a_{*}\sin{\theta}\\ \end{array}\right)+\left(\begin{array}[]{cc}\sin{\theta}&\cos{\theta}\\ -\cos{\theta}&\sin{\theta}\\ \end{array}\right)\left(\begin{array}[]{c}\delta_{X}\sin{2S}\\ \delta_{Y}\cos{2S}\\ \end{array}\right),$$

$$a_{\ast }\equiv \frac{m_{12}}{m_{12}+m_{\ast }}\left(1+\frac{3}{4}{\alpha }^2\frac{\mu }{m_{12}}\right)a_{\mathrm{out}},{\delta }_X\approx \frac{3}{8}\frac{\mu }{m_{12}}{\alpha }^4{a}_{\mathrm{out}}\mathrm{and}{\delta }_Y\approx \frac{9}{16}\frac{\mu }{m_{12}}{\alpha }^4{a}_{\mathrm{out}}.$$

$$V_{\mathrm{RV}}\equiv \dot{z}=\dot{Y}\sin I.$$

$$V_{\mathrm{RV}}$$

$$K_0\equiv \frac{m_{12}}{m_{12}+m_{\ast }}{a}_{\mathrm{out}}{\nu }_{\mathrm{out}}\mathrm{and}{K}_{\mathrm{BBH}}\equiv \frac{m_1{m}_2}{(m_1+m_2{)}^2}\sqrt{\frac{m_{12}+m_{\ast }}{m_{12}}}{\alpha }^{3.5}{K}_0.$$

$$K_{\mathrm{BBH}}=(2\pi G{)}^{1/3}{\left(\frac{m_{123}}{m_{12}}\right)}^{5/3}{P}_{\mathrm{out}}^{-1/3}{m}_{123}^{1/3}{\gamma }^{-2}{\left(\frac{P_{\mathrm{in}}}{P_{\mathrm{out}}}\right)}^{7/3},$$

$$K_{\mathrm{pl}}=(2\pi G{)}^{1/3}{P}_{\mathrm{pl}}^{-1/3}{m}_{\mathrm{pl}}{m}_{\ast }^{-2/3},$$

$$m_{\mathrm{pl}}=2^{7/3}{m}_{\ast }^{2/3}{\left(\frac{m_{12}}{m_{12}+m_{\ast }}\right)}^{5/3}{\gamma }^{-2}(m_{12}+m_{\ast }{)}^{1/3}{\left(\frac{P_{\mathrm{var}}}{P_{\mathrm{out}}}\right)}^{8/3}.$$

$$a_{\mathrm{pl}}>f{\left(\frac{m_{\ast }}{3m_{12}}\right)}^{1/3}{a}_{\mathrm{out}}\left(\mu \equiv \frac{m_{12}}{m_{12}+m_{\ast }}\gtrsim 0.8\right),$$

Full text

techniques: radial velocities - celestial mechanics - (stars:) binaries (including multiple): close

• stars: black holes

††software: Astropy (Astropy Collaboration et al., 2013, 2018), RadVel (Fulton et al., 2018), REBOUND (Rein & Liu, 2012)

1 Introduction

The first direct detection of a gravitational wave (GW) from a binary black-hole (BBH) merger (Abbott et al., 2016) has convincingly established the presence of very compact BBHs in the universe. While the origin, evolution and distribution of such BBHs are not yet understood, several interesting scenarios have been proposed so far. The binary evolution scenario (e.g. Belczynski et al., 2002, 2007, 2012, 2016; Dominik et al., 2012, 2013; Kinugawa et al., 2014, 2016; Spera et al., 2019) assumes that a fraction of massive stars form binary systems, experiences supernovae, and finally evolves into compact binaries including BBHs. The dynamical formation scenario, on the other hand, considers that stars and BHs located in dense clusters experience significant gravitational interactions, and eventually form BBHs ejected from the cluster regions mainly through binary-single, in addition to binary-binary, scattering encounters (e.g. Portegies Zwart & McMillan, 2000; O’Leary et al., 2009; Rodriguez et al., 2016; Tagawa et al., 2016). The primordial BH scenario (e.g. Sasaki et al., 2016, 2018; Bird et al., 2016) proposes that numerous primordial BHs are formed in the very early universe and a tiny fraction of them experiences close encounter and forms tight BBHs through the GW emission.

Regardless of such different formation scenarios, however, there should be a much more abundant population of their progenitor BBHs. Since those BBHs would have wider separations and thus longer orbital periods, they do not generate a detectable GW signal, except for the last few seconds before the final merger. Also they are difficult to be directly imaged unless they are surrounded by appreciable accretion disks. Therefore, the presence of such unseen binaries have to be searched for through their dynamical influence on nearby visible objects, including a radial velocity variation of an outer tertiary star due to the inner BBH. Such a detection is challenging, but if successful, it will not only constrain the formation and evolutionary channel towards the GW emitting BBHs, but also establish a yet unknown class of exciting astrophysical objects.

Indeed there exist a couple of known interesting systems that are relevant for such a strategy. One is a triple system consisting of a white dwarf-pulsar binary and an outer white dwarf orbiting around the inner binary (Ransom et al., 2014). The system was detected with the arrival time analysis of the pulsar. Quite interestingly, the inner and outer orbits of the triple system are near-circular and coplanar; the eccentricities of the inner and outer orbits are $e_{\mathrm{in}}\sim 6.9\times 10^{-4}$ and $e_{\mathrm{out}}\sim 3.5\times 10^{-2}$, and their mutual inclination is $i=(1.20\pm 0.17)\times 10^{-2}~{}\mathrm{deg}$.

The other is a red giant 2M05215658+4359220 with an unseen massive object of $\sim 3M_{\odot}$, possibly a BH or even a neutron star binary (Thompson et al., 2019). The system was discovered from a systematic survey of stars exhibiting anomalous accelerations. The follow-up radial-velocity (RV) observation indicates that the orbit is also near-circular; $e_{\mathrm{out}}=0.00476\pm 0.00255$. These two examples are very encouraging, implying that the dynamical search for unseen companions of visible objects is very rewarding.

A very recent discovery of a star–BH binary system, LB-1, by Liu et al. (2019) is added in the above list. Later, Abdul-Masih et al. (2019) and El-Badry & Quataert (2019) pointed out that the mass of the inner BH would plausibly be $5M_{\odot}\lesssim M\lesssim 20M_{\odot}$, much smaller than the original claim of $\sim 70M_{\odot}$ by Liu et al. (2019). In addition, Shen et al. (2019) put a strong constraint on the presence of an inner BBH for LB-1. Nevertheless, it is very encouraging in a sense that such a new population of star-BH binary systems in a nearly circular orbit really exists, perhaps abundantly in our Galaxy.

Indeed, there are several on-going/future projects that search for unseen companions around stars. Gaia was launched at the end of 2013, and performs astrometric survey for a billion of stars in the Galaxy. Since the astrometry of Gaia has great precision especially for bright stars, it can detect a subtle motion of a star around an unseen object. Therefore, in addition to the survey of stellar binaries (e.g. El-Badry & Rix, 2018; Ziegler et al., 2018), there are many proposals to search for star – BH binaries with Gaia (e.g. Breivik et al., 2017; Kawanaka et al., 2017; Mashian & Loeb, 2017; Yamaguchi et al., 2018). Yamaguchi et al. (2018), for instance, estimate that Gaia can detect $200-1000$ binaries in its 5 year operation considering the binary evolution scenario.

Additionally, TESS launched in 2018 is carrying out a photometric survey of nearby stars to search for transit planets. Masuda & Hotokezaka (2019) conclude that TESS will potentially discover $\mathcal{O}(10)$ binaries consisting of a star and an unseen compact object through identifying relativistic effects in their photometric light-curves, for instance.

Thus it is quite likely that Gaia, TESS and other surveys will detect many binary systems with an unseen object. Given the LIGO discovery of very tight BBHs, it is natural to expect that a fraction of those systems are indeed triple systems that host unseen inner BBHs. Therefore it is important to see if one can dynamically distinguish between a single BH and a binary BH in such triple systems. In addition, the stellar evolution of multiple systems is still poorly understood. Future detection of a triple with an inner BBH would shed light on dynamical processes in triple evolution (e.g. Toonen et al., 2016; Hamers et al., 2013; Stephan et al., 2016).

In this paper, we propose a methodology to search for an inner BBH in a triple system from the short-term outer stellar RV variations due to inner BBH perturbation through intensive RV monitoring (see Figure 1). Interestingly, such RV variations were first considered by Schneider & Cabrera (2006) in a context of a possible degeneracy between an inner binary and an S-type circumbinary planet. Later, Morais & Correia (2008) derived an analytic perturbative formula for RV variations in a coplanar circular triple, and pointed out that the inner binary model produces two nearby but distinct frequency modes around twice its orbital frequency. If the two modes can be identified separately, the degeneracy with the planet model, which produces a single frequency mode, can be broken in principle. Their work is further extended for eccentric and inclined triples by (Morais & Correia, 2011).

Such approximate analytic formulae are very useful in understanding the origin and parameter dependence of the RV variations, but they neglect the back-reaction on the inner BBH by the outer orbit. Therefore accurate numerical simulations are needed for broader and more realistic configurations of such triples. This is what we conduct in what follows. We also generate mock RV curves, and discuss the observational feasibility for our fiducial model.

The rest of the paper is organized as follows. In section 2, we first summarize the analytic approximation RV formula for a coplanar circular triple by Morais & Correia (2008), and explain features of the RV variations induced by an inner unseen binary. As a specific application of the formula, we constrain a possible inner binary in a star– unseen-object binary system 2M05215658+4359220. Section 3 introduces a set of fiducial triple configurations. In particular, we focus on a coplanar circular triple, and present a methodology to search for an inner BBH through the RV curves computed numerically. Section 4 applies the methodology for coplanar eccentric triples, and discusses the observational feasibility through mock observations. We discuss how to break the degeneracy of the RV variations due to circumbinary planets and the inner BBHs. Finally, section 6 summarizes and concludes this paper. Appendix presents a detailed derivation of an analytic approximation formula by Morais & Correia (2008) and the parameter degeneracy of the RV variations in the S-type planet and the inner BBH models.

2 Radial-velocity variation for a coplanar circular triple

In this section, we first present an approximation formula for RV variations for a coplanar circular triple with an inner binary, which was derived by Morais & Correia (2008). This analytic formula is useful to show basic feature of RV variations, and also to estimate the expected amplitude of RV variations in designing an RV followup strategy discussed in the next section. As a specific application for that purpose, we put a constraint on a possible inner binary in a binary system 2M05215658+4359220. For definiteness, we summarize a derivation of the approximate formula in appendix A.

2.1 Analytic approximation formula

We consider a coplanar circular hierarchical triple consisting of an inner binary and a star. The RV of the tertiary star orbiting around the inner binary (Morais & Correia, 2008) is given by

[TABLE]

where $I$ is the inclination of the stellar orbit axis with respect to the line of sight, $\nu_{\mathrm{in}}$ and $\nu_{\mathrm{out}}$ are mean motions of the inner and outer orbits, and $S_{0}$ and $\theta_{0}$ are constant phases specified by the initial conditions.

The first term in the right-hand-side of equation (1) corresponds to the Keplerian motion of the star. The semi-amplitude of the unperturbed Keplerian RV for a edge-on view is given by

[TABLE]

where $m_{1}$, $m_{2}$ and $m_{*}$ are the masses of the inner binary and the outer star, and $a_{\mathrm{out}}$ is the semi-major axis of the stellar orbit. For further detail, see appendix A.

The second and third terms express the RV variations induced by an inner binary. The semi-amplitude of the variations is proportional to $K_{\mathrm{BBH}}$ defined as

[TABLE]

with $a_{\mathrm{in}}$ being the semi-major axis of the inner binary.

Equation (1) indicates that the RV variations due to an inner binary consist of two periodic terms with frequencies of $2\nu_{\mathrm{in}}-\nu_{\mathrm{out}}$ and $2\nu_{\mathrm{in}}-3\nu_{\mathrm{out}}$. Since we consider a system with $\nu_{\mathrm{in}}\gg\nu_{\mathrm{out}}$, $K_{\mathrm{BBH}}$ is much smaller than $K_{0}$ by a factor of $(a_{\mathrm{in}}/a_{\mathrm{out}})^{3.5}$, and their frequencies are well approximated by $2\nu_{\mathrm{in}}$. Thus it is not easy to identify the two distinct modes observationally.

Morais & Correia (2008) pointed out, however, that the separation of the two modes is crucial to distinguish the RV variation of the inner binary against a single periodic modulation induced by a planet orbiting the star. Morais & Correia (2011) extended this consideration for eccentric and inclined cases using a perturbation theory.

The above result shows that a planet around a star in a binary system produces a signal very similar to that expected for an inner binary in a triple system. We discuss the degeneracy between the two models in section 5. In turn, a detection of an inner binary requires a well-organized observational monitor with short-cadence as we discuss below.

2.2 Constraint on a binary system 2M05215658+4359220

Thompson et al. (2019) reported a discovery of a binary system 2M05215658+4359220 that consists of a red giant and an unseen massive object. They first searched for systems exhibiting anomalously large radial accelerations from the Apache Point Observatory Galactic Evolution Experiment (APOGEE) radial velocity data, and selected 200 candidates of such binaries. After checking the photometric variations from the All-Sky Automated Survey for Supernovae (ASAS-SN) data, they identified 2M05215658+4359220 as the most likely binary candidate. Subsequently, they performed the radial velocity follow-up observation with the Tillinghast Reflector Echelle Spectrograph (TRES) on the $1.5~{}\mathrm{m}$ telescope at the Fred Lawrence Whipple Observatory (FLWO). They obtained 11 spectra with the precision of about $0.1~{}\mathrm{kms^{-1}}$ over several months in 2017 and 2018.

Interestingly, the orbital period of the system turned out to be very close to that of the photometric variation of the star, indicating the tidal synchronization. Therefore they assume that the inclination of the rotation axis of the outer red giant $i_{\mathrm{rot}}$ is equal to its orbital inclination $i_{\mathrm{orb}}$ with respect to the line of sight: $i\equiv i_{\mathrm{rot}}=i_{\mathrm{orb}}$. This enabled them to estimate the best-fit parameters of the system (Table 1) from the RV data and the spectroscopic analysis of the red giant.

Thompson et al. (2019) estimated the mass of the unseen companion to be $m_{\mathrm{CO}}=3.3^{+2.8}_{-0.7}~{}\mathrm{M_{\odot}}$. Since the value exceeds a conventionally accepted range of the maximum mass of the neutron star, it could be a single BH, or even a binary neutron star/BH. For simplicity, we assume that the inner binary is near-circular and coplanar with the outer orbit, and set $m_{\mathrm{giant}}=m_{*}$ and $m_{\mathrm{CO}}=m_{1}+m_{2}$; see Table 1. Then we apply the RV approximation formula so as to constrain the viable parameters for the possible inner binary.

Figure 2 plots a color contour map of $K_{\rm\scriptscriptstyle BBH}$ as a function of the mass ratio $m_{2}/m_{1}$ and the orbital period $P_{\mathrm{in}}$ of the possible inner binary for the 2M05215658+4359220 system. We compute the value of $K_{\rm\scriptscriptstyle BBH}$ from equation (3) using the parameters listed in Table 1. The color is coded according to the value of $K_{\rm\scriptscriptstyle BBH}$ that labels the contour curves.

The axis on the right indicates the semi-major axis ratio $(a_{\mathrm{in}}/a_{\mathrm{out}})$ corresponding to the orbital period $P_{\mathrm{in}}$ of the left axis. Note that the approximation formula is degraded towards the upper part of Figure 2. Indeed the three-body system becomes dynamically unstable if it satisfies (Mardling & Aarseth, 2001)

[TABLE]

Thus, the perturbation result completely breaks down there as indicated by the black area in Figure 2. Although the approximation formula is degraded for large $(a_{\mathrm{in}}/a_{\mathrm{out}})$, we can still conservatively estimate the semi-amplitude of RV variations in the allowed region of Figure 2.

Incidentally, the coalescence time of a circular compact binary $t_{\mathrm{GW}}$ via its GW emission is given by Peters (1964):

[TABLE]

where $m_{12}=m_{1}+m_{2}$, $\mu\equiv m_{1}m_{2}/m_{12}$. Equation (5) ensures that we can safely neglect the effect of GW emissions unless $P_{\mathrm{in}}\ll 1~{}\mathrm{days}$.

Considering no detection of anomalous RV variations by TRES beyond $\sim 100~{}\mathrm{m/s}$, Figure 2 shows that an assumed inner coplanar circular binary is constrained to have orbital period less than a couple of weeks if the binary consists of an almost equal mass objects.

Assuming a planet orbiting around the red giant with its stellocentric semi-major axis $a_{\mathrm{pl}}$, the orbital period of planet $P_{\mathrm{pl}}$ is obtained as

[TABLE]

Thus, the period of RV variations by a planet needs to be longer than $11$ days. On the other hand, the period of variations by an inner binary is constrained to be shorter than $P_{\mathrm{in}}/2\sim 13/2$ days from Figure 2. Since the outer orbiting star is a red giant of a radius $\sim 30~{}R_{\odot}$ (see Table 1), we can basically break the degeneracy between an inner binary and an S-type circumbinary planet even if the periodic RV variation is detected.

3 Strategy to search for an inner binary: case for a coplanar circular triple

In this section, we concentrate on coplanar circular triples before considering coplanar eccentric systems in the next section. The main purpose of this section is to describe our methodology to search for a possible inner binary in a star – BH system, and also to check the accuracy of our numerical simulation against the analytic solution. Since there are many parameters characterizing a triple system, we cannot explore an entire parameter space. Instead, we fix a fiducial model, and examine the feasibility of our method for the model.

Table 2 summarizes the initial values of parameters in our fiducial model. As before, we use subscripts “$\mathrm{in}$” and “$\mathrm{out}$” to denote the parameters for inner and outer orbits, respectively. We adopt a $1~{}M_{\odot}$ star as the tertiary for simplicity. Also we assume a relatively large value for $(a_{\mathrm{in}}/a_{\mathrm{out}})$ within the dynamically stable range (Mardling & Aarseth, 2001). In the present section, we consider a circular system, and set inner and outer eccentricities $e_{\mathrm{in}}$, $e_{\mathrm{out}}$ to $10^{-5}$. We vary this value when we consider eccentric systems later in section 4. As a reference, Table 3 lists the predicted values for orbital motions calculated from equation (1).

Then, we perform N-body simulations for the coplanar circular triple with initial conditions in Table 2, using a public N-body package REBOUND (Rein & Liu, 2012). We use WHFast integrator (Rein & Tamayo, 2015), which is one of the fast and accurate symplectic integrators, and set the calculation time interval to be $10^{-6}~{}\mathrm{yr}/2\pi$. We output the RV data in $0.1$-day cadence over 800 days in total.

Since we are looking for a tiny RV variation relative to a much larger Keplerian component, we need to determine the base-line Keplerian motion accurately and subtract it from the total RV. Although it may seem a trivial and straightforward procedure, it is not the case in reality.

Figure 3 illustrates an example of RV residuals. First we compute the total RV signal of the tertiary star numerically. Since we know the initial values of the input parameters (Table 2), we subtract the first-term in the right-hand-side of equation (1) assuming $P_{\mathrm{out}}=638$ days and $e_{\mathrm{out}}=0$ (Table 3) from the total RV signal. The result is plotted as a black dashed line, which exhibits a large-amplitude modulation of a period close to $P_{\mathrm{out}}$. This discrepancy comes from the fact that the orbital period of the outer star is indeed affected by the motion of the inner binary. In other words, $P_{\mathrm{out}}$ is not a constant unlike as assumed in the approximate formula, equation (1).

Therefore, we next fit the total RV signal using a RV fitting code RadVel (Fulton et al., 2018), and determine the best-fit Keplerian motion locally over the over 800 days. In doing so, we fix a constant motion ($\gamma$ velocity) and neglect the RV jitters. The resulting best-fit values are summarized in Table 4. The blue line in Figure 3 is the RV residual after removing the locally fitted Keplerian component with $\tilde{P}_{\mathrm{out}}\approx 619.3$ days and $\tilde{e}_{\mathrm{out}}\approx 0.016$. It is in good agreement with the RV variation induced by an inner binary, the second and third terms in equation (1), which is plotted as a red line in Figure 3.

The above result implies that the removal of the accurate Keplerian component from the total RV signal is crucial to extract the RV variations due to an inner binary. The good agreement between our numerical result and the analytic formula for the RV variation simultaneously justifies the robustness of our methodology proposed here. There is a tiny discrepancy between the blue and red lines, but it is likely due to the limitation of the perturbative approximation, because our fiducial set of parameters are not fully in the perturbative regime.

Therefore, in what follows, we use fitting with RadVel to determine the base-line Keplerian component, and define the residuals as RV variations.

4 Mock observation of a coplanar eccentric triple

This section considers effects of eccentricities for coplanar triples. Although Morais & Correia (2011) provide analytic perturbative RV formulae even for non-coplanar eccentric cases, the applicability of the result is limited by the validity of the perturbative approach. Also the formulas are fairly complicated and it is not easy to understand the overall behavior of RV variations. Therefore, we first provide examples for the predicted RV variations from numerical simulations, and then perform mock observations including noise and finite sampling. Finally we discuss the feasibility of the detection of an inner binary based on our methodology.

4.1 RV variations for coplanar eccentric triples

For definiteness, we consider four cases of different eccentricities: $(e_{\mathrm{in}},e_{\mathrm{out}})=(10^{-5},10^{-5})$, $(10^{-5},0.2)$, $(0.1,10^{-5})$, and $(0.1,0.2)$. The other parameters are fixed as those listed in Table 2. As in section 3, we generate the RV data in 0.1-day cadence over 800 days using REBOUND, and extract the base-line Keplerian motion fitted with RadVel

Figure 4 shows the RV variations (left) and the corresponding Lomb-Scargle (LS) periodogram (right). The LS periodograms are calculated using a community-developed core Python package for Astronomy, Astropy (Astropy Collaboration et al., 2013, 2018).

The RV variation for a coplanar circular triple exhibits a clear periodic signal as shown in the top-left panel of Figure 4. According to equation (1), the variation should have two distinct frequencies $\nu_{\mathrm{var1}}\equiv 2\nu_{\mathrm{in}}-\nu_{\mathrm{out}}$, and $\nu_{\mathrm{var3}}\equiv 2\nu_{\mathrm{in}}-3\nu_{\mathrm{out}}$, corresponding to

[TABLE]

and

[TABLE]

respectively. The LS periodograms in the top-right panel of Figure 4 indicates that the RV variation that we compute numerically is dominated by the mode with $\sim 1/P_{\mathrm{var3}}\approx 0.02~{}{\rm day}^{-1}$, and the mode with $\sim 1/P_{\mathrm{var1}}\approx 0.023~{}{\rm day}^{-1}$ is barely visible. For comparison, we also computed the LS periodogram directly for the second and third terms of equation (1) with the same parameters and 0.1-day cadence over 800 days. We found that the secondary peak at $\nu_{\mathrm{var1}}$ is not visible either, but correctly reproduced only when the total base-line duration is significantly larger than the 800 days we adopted here. Therefore, the identification of the two distinct modes due to an inner binary is very challenging in practice.

The LS periodograms also imply that non-zero eccentricities generate a variety of additional frequency modes in the LS periodograms, and that the value of the frequency $\nu_{\mathrm{var3}}$ is fairly insensitive to $e_{\mathrm{in}}$, but decreases as $e_{\mathrm{out}}$ increases. The sensitivity of the RV variation curve in time domain, instead of the LS periodogram, on $e_{\mathrm{in}}$ and $e_{\mathrm{out}}$ may be useful to estimate those eccentricities of the system. Thus we conclude that the presence of the inner binary itself can be inferred robustly. Also we find that the eccentricity is useful in distinguishing between the inner binary and the planetary signals as we discuss in section 5.

4.2 Mock observations

In order to examine the feasibility of our methodology to detect an inner binary, we perform mock observations of the RV signal including observational noise and finite sampling effect. If the tertiary star is a solar-type star located at $\sim 1~{}\mathrm{kpc}$, the apparent magnitude is around $15$. According to Plavchan et al. (2015), the statistical RV error is at the level of a few tens of m/s with an ideal 10m telescope with $\sim 1000$ second exposure. While the noise is crucially dependent on the nature of the star, we neglect any systematic/non-Gaussian noise, and add $20$ m/s Gaussian noise into the RV signal.

The observational cadence also affects the detectability of the inner binary. We consider two different cases; 10 and 30 data points randomly selected from one-day cadence data over 100 days. In both cases, we assume that the same cadence data are available for 600 days. Then we fit the data over 600 days with RadVel to extract the base-line Keplerian component, and obtain the RV residuals.

The results are shown in Figure 6 for the 30 percent sampling rate, and in Figure 6 for the 10 percent sampling rate. Blue lines in the left panels of Figures 6 and 6 correspond to the original RV variation curves of Figure 4. We note that we add the noise on the total RV curve, instead of the RV variation curve (blue lines). Thus the extracted Keplerian component is not exactly the same depending on the Gaussian noise. This is why the data points in the left panels of the two figures do not necessarily distribute around the blue lines. This clearly implies that an accurate determination and extraction of the base-line Keplerian component is crucial in our methodology.

Except for the limitation, Figures 6 and 6 indicate that the RV variations can clearly reveal the presence of the inner binary if the relevant short-cadence RV follow-up is performed. The LS periodograms in our mock data are very similar to the idealized cases shown in Figure 4; the peaks around the frequency $\nu_{\mathrm{var3}}$ are fair robust even in the (modest) eccentricities we assumed here.

While the inner binary orbital period can be inferred relatively easily, it would be difficult to determine the mass ratio $m_{2}/m_{1}$ of the binary. As equation (3) indicates, the RV variations are more insensitive to the mass ratio than to the inner binary separation. Even a small uncertainty of the amplitude due to noise would affect the estimate of the mass ratio.

5 Degeneracy with an S-type circumbinary planet

As mentioned in sections 1 and 2, an S-type circumbinary planet (see Figure 7) produces a RV variation similar to that due to an inner binary. We examine the degree of a possible degeneracy by performing mock observations presented in section 4.

We consider the configuration of a star– BH binary system with a planet around a star. Following Appendix B, we choose a set of parameters for a planet (Table 5) so as to mimic the RV variations produced by our fiducial triple model (Table 2). We also consider non-zero eccentricities of the planetary and stellar orbits: $(e_{\mathrm{pl}},e_{\mathrm{out}})=(10^{-5},10^{-5})$, $(10^{-5},0.2)$, $(0.1,10^{-5})$, and $(0.1,0.2)$. While we choose the same set of values for $(e_{\mathrm{in}},e_{\mathrm{out}})$ in subsection 4.2, $e_{\mathrm{in}}$ and $e_{\mathrm{pl}}$ do not have to be the same value in the two different pictures. We then perform the mock observations similar in section 4. We here consider the 30 percent sampling rate and examine if we can distinguish between the two pictures observationally.

Figure 8 shows the resulting RV variations (left) and the LS periodograms (right) for star–BH binaries with a planet orbiting around the star. As expected, a circular inner binary (top panels in Figure 6) is difficult to be distinguished from a star with a planet. The RV variation curves due to a planet are not so much affected by $e_{\mathrm{pl}}$ nor $e_{\mathrm{out}}$, in contrast to their sensitivity due to an inner binary on $e_{\mathrm{in}}$ and $e_{\mathrm{out}}$ (Figure 6). Thus eccentricities of the inner binary is indeed helpful to break the degeneracy between the two models. In appendix B, we provide the parameter correspondence between the two possibilities for coplanar and circular systems.

6 Summary

The LIGO-VIRGO detection of gravitational waves implies that the universe harbors a number of massive BBHs that have been dismissed in traditional astronomy. Since the GW emission from such BBHs is detectable only for the last few seconds before the final merger, their progenitors should be even more numerous and remain undetected as a population of unseen binaries with longer orbital periods.

We propose a methodology to detect such wider-separation BBHs from the motion of a star orbiting around them. The radial velocity of the tertiary star in a star – BBH triple system exhibits a periodic variation due to the orbital motion of the inner BBH. The variation signal is small compared with the Keplerian component of the stellar motion. Nevertheless it is comparable to, or even larger than, a typical amplitude of the radial velocities for observed exoplanetary systems. Thus such a velocity variation signal is indeed detectable if such star-BBH triples are identified.

There are several on-going projects that search for star – compact-object binaries based on astrometry, transit and radial-velocity surveys. They will soon discover a number of star-BH binary candidates, and a fraction of such binaries may be star-BBH triples, instead of mere binaries. Indeed a couple of binary systems, 2M05215658+4359220 (Thompson et al., 2019) and LB-1 (Liu et al., 2019), provide particularly interesting examples.

We have generated a series of mock radial velocity curves for coplanar systems, and showed that a short-period variation signal due to the inner BBH is detectable relatively easily as long as the star is sufficiently bright to allow for high-resolution spectroscopy, but that the mass ratio of the inner BBHs is difficult to estimate. The eccentricity of the tertiary star is useful in breaking the degeneracy between a planet orbiting the star and an inner BBH.

We find that a short-period variation in the 2M05215658+4359220 binary, if detected at all, should be due to the inner unseen binaries (neutron stars and/or BHs), since the red giant in the system cannot accommodate a planet with short orbital periods.

The discovery of LB-1 strongly indicated a presence of abundant star-BH binary systems in our Galaxy. LB-1 was originally claimed to contain a $68^{+11}_{-13}M_{\odot}$ unseen companion, which is intriguingly close to the total mass of several BBHs discovered by LIGO-VIRGO. Although subsequent studies (e.g. Abdul-Masih et al., 2019; El-Badry & Quataert, 2019) suggested that the mass of the unseen companion is much smaller, the system still remains a promising target of the future RV follow-up. The current analysis focused on coplanar triples, but the inner and outer orbits are more likely to be mutually inclined. The extension of our methodology to LB-1 in a more general configuration will be examined and reported elsewhere.

Finally, we would like to emphasize that the methodology presented in this paper is no more a mere theoretical idea. It will find a realistic application in the near future that might unveil a hitherto unknown population of astronomical objects in a complementary manner to gravitational-wave observations.

We would like to thank Alexandre C. M. Correia for letting us know of their important work on the degeneracy between the inner binary and the star with a planet, and Zheng Zheng for useful correspondences on LB-1. Simulations and analyses in this paper made use of REBOUND, RadVel, and Astropy. We gratefully acknowledge the support from Grants-in-Aid for Scientific Research by Japan Society for Promotion of Science (JSPS) No.18H01247 and No.19H01947, and from JSPS Core-to-core Program “International Network of Planetary Sciences”.

Appendix A The derivation of the RV variations induced by an inner binary

We here derive the equation (1) following Morais & Correia (2008). We consider a coplanar circular triple consisting of an inner binary ($\rho\equiv|\bm{r}_{\mathrm{in}}|/|\bm{r}_{\mathrm{out}}|\ll 1$) and a star orbiting around the binary (see the right panel of Fig.9). The Hamiltonian of the system to the second order of $\rho$ is written as

[TABLE]

where $S\equiv\angle(\bm{r}_{\mathrm{in}}\bm{r}_{\mathrm{out}})$, $m_{12}\equiv m_{1}+m_{2}$, $\bm{p}_{\mathrm{in}}$ and $\bm{p}_{\mathrm{out}}$ are the momenta corresponding to $\bm{r}_{\mathrm{in}}$ and $\bm{r}_{\mathrm{out}}$ in terms of Jacobi coordinates, respectively. The reduced masses $\mu$ and $\mu_{*}$ are defined as

[TABLE]

Using the Hamiltonian $\mathcal{H}$, the equations of motion are obtained from the canonical equations:

[TABLE]

Therefore, the explicit forms of equations of motion for the inner binary and outer star are

[TABLE]

and

[TABLE]

respectively.

The 0-th order solutions of equations (A4) and (A5) for a circular system are as follows:

[TABLE]

where $S=S_{0}+(\nu_{\mathrm{in}}-\nu_{\mathrm{out}})t$, and $a_{\mathrm{in}}$ and $a_{\mathrm{out}}$ are the semi-major axes of the inner and outer orbits, respectively, with $S_{0}$ being a constant phase determined by the initial positions. The mean motions $\nu_{\mathrm{in}}$ and $\nu_{\mathrm{out}}$ are defined as

[TABLE]

In terms of the Cartesian coordinate on the orbital plane, unit vectors $\hat{\bm{r}}_{\mathrm{out}}$ and $\hat{\bm{u}}$ are

[TABLE]

where $\theta\equiv\theta_{0}+\nu_{\mathrm{out}}t$, and $\bm{i}$ and $\bm{j}$ are constant base vectors (see Figure 9), with $\theta_{0}$ being a constant phase determined by the initial positions. Substituting equations (A6) - (A8) into the equation (A5), we obtain

[TABLE]

and

[TABLE]

where $X$ and $Y$ are defined as $\bm{r}_{\mathrm{out}}=X\bm{i}+Y\bm{j}$.

The solutions of these equations are written as

[TABLE]

where

[TABLE]

The above expressions for $\delta_{X}$ and $\delta_{Y}$ are derived on the assumption of $\nu_{\mathrm{out}}/\nu_{\mathrm{in}}\ll 1$.

The radial velocity $V_{\mathrm{RV}}$ is defined as

[TABLE]

Therefore, we can write down the following RV approximation formula:

[TABLE]

where

[TABLE]

Appendix B Correspondence of parameters between a triple and

S-type circumbinary planet

In section 5, we find that for a coplanar circular triple the RV variations are almost degenerate with those produced by an S-type circumbinary planet. If a star is a giant, we can sometimes rule out the degeneracy using the same way as for the system 2M05215658+4359220. However, this is not always applicable. Thus, we here consider a parameter correspondence for a coplanar circular case when we cannot distinguish these two possibilities.

The RV variations are basically characterized by its semi-amplitude $K_{\rm var}$ and period $P_{\rm var}$. The latter is equal to $P_{\rm pl}$ and $P_{\mathrm{in}}/2$ in the planet-star and inner BBH interpretations, respectively.

The semi-amplitude of RV variations induced by an inner binary is estimated from equation (3):

[TABLE]

where $P_{\mathrm{out}}$ is the orbital period of the outer star, $m_{12}=m_{1}+m_{2}$, $m_{123}=m_{12}+m_{*}$, and $\gamma\equiv\left(\sqrt{m_{2}/m_{1}}+\sqrt{m_{1}/m_{2}}\right)$.

On the other hand, the semi-amplitude of variations induced by a coplanar circular planet around the outer star is

[TABLE]

where $m_{\mathrm{pl}}$ is the mass of planet, and we assume that $m_{\mathrm{pl}}\ll m_{*}$.

We cannot distinguish these two interpretations from observation if both of them induce the RV variations with similar amplitude and period. For this case, using equations (B1) and (B2), we can obtain the parameter correspondence considering $K_{\mathrm{var}}=K_{\rm\scriptscriptstyle BBH}=K_{\mathrm{pl}}$, and $P_{\mathrm{var}}=P_{\mathrm{pl}}=P_{\mathrm{in}}/2$:

[TABLE]

Equation (B3) relates the parameters in the two different interpretations.

Figures 10 shows the parameter correspondence between $m_{\mathrm{pl}}$ in the planet-star interpretation and $m_{12}$ in the inner BBH interpretation for $P_{\mathrm{out}}=100~{}\mathrm{days}$. We assume that $m_{2}/m_{1}=1$ and $m_{*}=1~{}M_{\odot}$. Given the values of the RV variation semi-amplitude and period, $K_{\rm var}$ (color-coded) and $P_{\rm var}$ (black contour), the corresponding values of $m_{\mathrm{pl}}$ and $m_{\mathrm{12}}$ can be read off from the figures.

The dotted area above the blue curve is excluded in the inner BBH picture from the instability of the triple. We use the dynamical stability criterion (Mardling & Aarseth, 2001) (see equation (4)) to compute the region. The tiny black region is excluded by the instability of the planet orbiting the star. To compute the region, we adopt the following Hill instability criterion (e.g. Holman & Wiegert, 1999; Barnes & O’Brien, 2002):

[TABLE]

where $f=0.36$ is derived from numerical simulations, and $a_{\mathrm{pl}}$ is the stellocentric semi-major axis of planet. Note that $m_{12}$ is interpreted as the mass of a single BH in a planet-star picture. Figure 10 shows that most of regions accept both interpretations, indicating that a Hot Jupiter around a star in a star-BH binary system may mimic the RV variation induced by an inner binary in the system.