# Dynamics of Binary System around Supermassive Black Hole

Kei-ichi Maeda  Department of Physics, Waseda University, Shinjuku, Tokyo
169-8555, Japan  Center for Gravitational Physics and Quantum Information,
Yukawa Institute for Theoretical Physics, Kyoto University, 606-8502, Kyoto,
Japan  Priti Gupta  Department of Physics, Indian Institute of Science,
Bangalore 560012, India  Hirotada Okawa  Waseda Institute for Advanced Study
(WIAS), 1-21-1 Nishi Waseda, Shinjuku-ku, Tokyo 169-0051, Japan

###### Abstract

We discuss motion of a binary system around a supermassive black hole. Using
Fermi-Walker transport, we construct a local inertial reference frame and set
up a Newtonian binary system. Assuming a circular geodesic observer around a
Schwarzschild black hole, we write down the equations of motion of a binary.
Introducing a small acceleration of the observer, we remove the interaction
terms between the center of mass (CM) of a binary and its relative
coordinates. The CM follows the observer’s orbit, but its motion deviates from
an exact circular geodesic. We first solve the relative motion of a binary
system, and then find the motion of the CM by the perturbation equations with
the small acceleration.

We show that there appears the Kozai-Lidov (KL) oscillations when a binary is
compact and the initial inclination is larger than a critical angle. In a hard
binary system, KL oscillations are regular, whereas in a soft binary system,
oscillations are irregular both in period and in amplitude, although stable.
We find an orbital flip when the initial inclination is large. As for the
motion of the CM, the radial deviations from a circular orbit become stable
oscillations with very small amplitude.

##  I Introduction

After the remarkable success of the LIGO-Virgo-KAGRA Collaboration Abbott et
al. (2020, 2021), the study of gravitational wave (GW) emission has received a
significant boost. The analysis of data (obtained through the first three
observational runs) produced over a hundred confident detections Collaboration
et al. (2019, 2023) from binary black holes (BH), a binary neutron star (NS)
and BH-NS systems, with more to follow in the next decade. The scientific
insights emerging from the detections have significantly revolutionized our
understanding of the sources. For instance, some notable events revealed
heavier stellar-mass BHs Collaboration and Collaboration (2016), and their
origin is still under discussion. Using the electromagnetic counterpart, we
found that the speed of GWs is very close to the speed of light as predicted
by the general theory of relativity Collaboration and Collaboration (2017).
With the increase in detections, we can probe more fundamental questions like
testing theories of gravity in strong field regimes, finding the redshift
distribution of BHs and their environment, and so on Collaboration and the
Virgo Collaboration (2019, 2020, 2021a, 2021b); Collaboration et al. (2021).

Unlocking the scientific potential of GWs depends on our theoretical
knowledge. In order to filter the GW signal from the detector noise, it is
necessary to model the predicted waveform accurately. Current observations are
from isolated binary systems. It is possible, however, that nature will
provide us with more exotic sources. This paper will examine a three-body
system as one of the likely sources. The environment near supermassive black
holes (SMBHs) in galactic nuclei comprises many stars and compact objects. A
binary system could emerge in these surroundings, which composes a natural
hierarchical triple system Heggie (1975); Hut (1993); Samsing et al. (2014);
Riddle et al. (2015); Antonini et al. (2016a); Stephan et al. (2019). Recent
LIGO events suggest the hierarchical systems as a possible dynamical formation
channel of the heavy merging binary BHs Mapelli et al. (2021); Gayathri et al.
(2020); Gerosa and Fishbach (2021).

In a hierarchical triple, the distance between two bodies (forming an ‘inner’
binary) is much less than the distance to the third body. In 1962, Kozai and
Lidov independently explored the dynamics of restricted hierarchical triples
Kozai (1962); Lidov (1962), revealing a remarkable phenomenon (known as Kozai-
Lidov (KL) resonance) — when the two orbits are inclined relative to each
other, there is a periodic exchange between orbital eccentricity and relative
inclination in secular timescale Shevchenko (Springer, 2017). The orbital
eccentricity can reach extreme values leading to large emission of GWs Randall
and Xianyu (2019a); Hoang and Naoz (2018); Gupta et al. (2020).

There has been extensive work on dynamics of such systems based on Newtonian
or post-Newtonian approximation Naoz et al. (2013b, 2012); Naoz (2016); Naoz
et al. (2019); Teyssandier et al. (2013); Li et al. (2015); Will (2014a, b).
Recently we also find many work focusing on gravitational waves from such
systems Amaro-Seoane et al. (2010); Antonini and Perets (2012); Hoang and Naoz
(2018); Antonini et al. (2016b); Meiron et al. (2017); Robson et al. (2018);
Randall and Xianyu (2019b, a); Hoang et al. (2019); Emami and Loeb (2019);
Gupta et al. (2020); Kuntz and Leyde (2022). It has been shown that KL
resonance leaves an imprint on the waveform and lies in the detectable range
of future space-based detectors like LISA and DECIGO. Indirect observation of
GW from a triple system is also studied by analyzing the cumulative shift of
periastron time of a binary pulsar undergoing KL oscillations Suzuki et al.
(2019, 2020).

In this paper, when we discuss on a binary system near SMBH, we follow another
approach, i.e., a binary system is treated as perturbations of SMBH spacetime.
In the case of a single object in SMBH spacetime, it can be treated a test
particle. But in the case of a binary system, it is not the case because the
self-gravitational mutual interaction is much stronger than the gravitational
tidal force by SMBH. In order to analyze such a hierarchical system, we first
prepare a local inertial frame and set up a binary in this frame. When a
binary is tightly bounded but the mutual gravitational interaction is not so
strong, the binary motion in this frame can be discussed by Newtonian
gravitational dynamics.

Using Fermi normal coordinate system or Fermi-Walker transport, we can
construct a local inertial frame Manasse and Misner (1963); Nesterov (1999);
Delva and Angonin (2012). Using such a technique, there are several
discussions on a tidal force acting on stars near SMBH Banerjee et al. (2019);
Ishii et al. (2005); Cheng and Evans (2013); Kuntz et al. (2021), but only a
few works on a binary system have been discussed Gorbatsievich and Bobrik
(2010); Chen and Zhang (2022); Camilloni et al. (2023). In this paper, we
analyze such a system in detail. Assuming an observer is moving along a
circular geodesic around a Schwarzschild SMBH, we construct a local inertial
frame, and set up a binary system. We then discuss motion of a binary, showing
existence of the KL oscillations when a binary is compact and the initial
inclination angle is larger than a critical value.

The paper is organized as follows: We review how to construct a local inertial
proper reference frame by use of Fermi-Walker transport in Sec. II A. In Sec.
II B, we perform post-Newtonian expansion for a test particle motion in this
frame. In Sec. III, we set up a self-gravitating system in the proper
reference frame and derive the Lagrangian in the Newtonian limit. In Sec. IV,
assuming an observer moving along a circular geodesic in Schwarzschild black
hole, we derive the equations of motion for a binary system. We also discuss
the interaction terms between the center of mass (CM) of a binary and its
relative coordinates. Introducing small acceleration of an observer, we remove
the interaction terms, finding the equations of motion for the CM, which gives
small deviations from a circular geodesic. We then analyze twelve models
numerically and show the properties of binary motions such as the KL
oscillations, chaotic features, and orbital flips in Sec. V. We also discuss
motions of the CM of a binary. Concluding remarks follow in Sec. VI. In the
Appendix A, we provide some numerical and analytic solutions in the coplanar
case. We also present the Lagrange planetary equations of the model and write
down the equations for the orbital parameters of a binary taking averages over
inner and outer binary cycles in Appendix B. We show that this simplified
method recovers numerical results obtained by direct integration of the
equations of motion in the case of a hard binary. It also provides the KL
oscillation time scale and the maximum and minimum values of eccentricity.

Notation used: Greek letters range from 0 to 3, while Roman letters run from 1
to 3; Hatted indices denote tetrad components in a proper reference frame
rotating along an observer; Bar over symbols correspond to quantities in a
static tetrad frame.

##  II Proper reference frame

###  II.1 Proper reference frame of an arbitrary observer in a curved
spacetime

We first discuss how to set up a local inertial frame in a curved spacetime
Misner et al. (1973); Manasse and Misner (1963); Gorbatsievich and Bobrik
(2010). The spacetime metric is given by

$\displaystyle d\bar{s}^{2}=\bar{g}_{\mu\nu}dx^{\mu}dx^{\nu}\,.$ (1)

We then consider an observer, whose orbit is given by a world line $\gamma$
described by

$\displaystyle x^{\mu}=z^{\mu}(\tau)\,,$

where $\tau$ is a proper time of the observer. The 4-velocity is given by

$\displaystyle u^{\mu}(\tau)\equiv{dz^{\mu}\over d\tau}\,.$

We now set up an orthonormal tetrad system $\\{e_{\hat{\alpha}}^{\mu}\\}$
along the observer’s world line $\gamma$, which is defined by the conditions
such that

$\displaystyle
e_{\hat{\alpha}}^{~{}\mu}e_{\hat{\beta}\mu}=\eta_{\hat{\alpha}\hat{\beta}}\,,~{}~{}e_{\hat{0}}^{~{}\mu}=u^{\mu}\,,$

where $\eta_{\hat{\alpha}\hat{\beta}}$ is Minkowski spacetime metric.

For a given 4-velocity $u^{\mu}$, this tetrad system is determined up to
three-dimensional rotations. The tetrad $e_{\hat{\alpha}}^{\mu}$ is
transported along the observer’s world line $\gamma$ as

$\displaystyle{De_{\hat{i}}^{~{}\mu}\over
d\tau}=-\Omega^{\mu\nu}\,e_{\hat{i}\,\nu}\,,$

where

$\displaystyle\Omega^{\mu\nu}\equiv
a^{\mu}u^{\nu}-u^{\mu}a^{\nu}+u_{\alpha}\omega_{\beta}\epsilon^{\alpha\beta\mu\nu}\,.$

Here,

$\displaystyle a^{\mu}$ $\displaystyle\equiv$ $\displaystyle{Du^{\mu}\over
d\tau},$

and

$\displaystyle\omega_{\mu}$ $\displaystyle\equiv$ $\displaystyle{1\over
2}u^{\alpha}\epsilon_{\alpha\mu\rho\sigma}\Omega^{\rho\sigma}={1\over
2}\stackrel{{\scriptstyle(3)}}{{\epsilon}}_{\mu\rho\sigma}\Omega^{\rho\sigma},$

are the acceleration of the observer and the angular velocity of a rotating
spatial basis vector $e_{(a)}^{\mu}$, respectively. A non-rotating tetrad
frame for which $\omega^{\mu}=0$ is called the Fermi-Walker transport. If the
orbit is a geodesic ($a^{\mu}=0$ and $\omega^{\hat{k}}=0$), we find
$\displaystyle{{De_{\hat{a}}^{\mu}\over d\tau}=0}$, which is just a parallel
transport.

Next, we construct a local coordinate system (the observer’s proper reference
system) near the world line $\gamma$, which is described as

$\displaystyle(x^{\hat{\mu}})=(c\tau,x^{\hat{a}})\,,$

where the spatial components $x^{\hat{a}}$ is measured from the point at
$\tau$ on the world line $\gamma$ along the spatial hypersurface
$\Sigma(\tau)$ perpendicular to $\gamma$.

We find that the metric form of this proper reference frame up to the second
order of $x^{\hat{a}}$ is given by

$\displaystyle
g_{\hat{\mu}\hat{\nu}}=\eta_{\hat{\mu}\hat{\nu}}+\varepsilon_{\hat{\mu}\hat{\nu}}+O(|x^{\hat{k}}|^{3}),$
(2)

where

$\displaystyle\varepsilon_{\hat{0}\hat{0}}$ $\displaystyle=$
$\displaystyle-{1\over
c^{2}}\left[2a_{\hat{k}}x^{\hat{k}}+\left(c^{2}\bar{\cal
R}_{\hat{0}\hat{k}\hat{0}\hat{\ell}}-\omega_{\hat{j}\hat{k}}\omega^{\hat{j}}_{~{}\hat{\ell}}\right)x^{\hat{k}}x^{\hat{\ell}}+{\left(a_{\hat{k}}x^{\hat{k}}\right)^{2}\over
c^{2}}\right],$ (3) $\displaystyle\varepsilon_{\hat{0}\hat{j}}$
$\displaystyle=$ $\displaystyle-{1\over
c^{2}}\left[c\,\omega_{\hat{j}\hat{k}}x^{\hat{k}}+{2\over 3}c^{2}\bar{\cal
R}_{\hat{0}\hat{k}\hat{j}\hat{\ell}}x^{\hat{k}}x^{\hat{\ell}}\right],$ (4)
$\displaystyle\varepsilon_{\hat{i}\hat{j}}$ $\displaystyle=$
$\displaystyle-{1\over c^{2}}\left[{1\over 3}c^{2}\bar{\cal
R}_{\hat{i}\hat{k}\hat{j}\hat{\ell}}x^{\hat{k}}x^{\hat{\ell}}\right]\,.$ (5)

where $\bar{\cal R}_{\hat{\mu}\hat{\nu}\hat{\rho}\hat{\sigma}}$ is the tetrad
component of the Riemann curvature of a background spacetime and
$\omega_{\hat{j}\hat{k}}\equiv\epsilon_{\hat{j}\hat{k}\hat{\ell}}\omega^{\hat{\ell}}$
Misner et al. (1973); Manasse and Misner (1963); Gorbatsievich and Bobrik
(2010).

The acceleration and angular frequency in the proper reference frame are
defined by

$\displaystyle a^{\hat{j}}$ $\displaystyle\equiv$ $\displaystyle
e^{\hat{j}}_{~{}\mu}{Du^{\mu}\over d\tau},$ $\displaystyle\omega^{\hat{j}}$
$\displaystyle\equiv$ $\displaystyle{1\over
2}\stackrel{{\scriptstyle(3)}}{{\epsilon}}^{\hat{j}\hat{k}\hat{\ell}}e_{\hat{\ell}\mu}{De_{\hat{k}}^{~{}\mu}\over
d\tau}\,.$

For the Fermi-Walker transport, $\omega^{\hat{k}}=0$. If it is the geodesic
($a^{\hat{j}}=0$ and $\omega^{\hat{k}}=0$), we recover the Fermi normal
coordinates. The tetrad is parallely transformed along the world line.

###  II.2 Test particle motion in a proper reference frame

First, we consider motion of a test particle in the above proper reference
frame. The action for a test particle with mass $m$ in a given spacetime is
given by

$\displaystyle S=-mc\int\sqrt{-ds^{2}}\,.$

For a test particle in the proper reference frame, since the world interval
$ds^{2}$ is given by the metric $g_{\hat{\mu}\hat{\nu}}$, we find the action
for a test particle as

$\displaystyle S=\int d\tau{\cal L}\,,$

where

$\displaystyle{\cal L}$ $\displaystyle\equiv$ $\displaystyle-
mc\sqrt{-\left(\eta_{\hat{\mu}\hat{\nu}}+\varepsilon_{\hat{\mu}\hat{\nu}}\right){dx^{\hat{\mu}}\over
d\tau}{dx^{\hat{\nu}}\over d\tau}},$

is Lagrangian of the test particle and $\varepsilon_{\hat{\mu}\hat{\nu}}$ is a
small deviation from the Minkowski spacetime since
$|\varepsilon_{\hat{\mu}\hat{\nu}}|\ll 1$.

Assuming that the test particle moves slowly in the proper reference frame, we
perform the post-Newtonian expansion in terms of $v^{\hat{j}}/c$ as follows:
First, we expand the square root term in the Lagrangian as

$\displaystyle{\cal L}$ $\displaystyle=$ $\displaystyle-
mc^{2}\sqrt{1-{\\!\\!\mbox{ \boldmath$v$}^{2}\over
c^{2}}-\varepsilon_{\hat{0}\hat{0}}-2\varepsilon_{\hat{0}\hat{j}}{v^{\hat{j}}\over
c}-\varepsilon_{\hat{i}\hat{j}}{v^{\hat{i}}v^{\hat{j}}\over c^{2}}}$
$\displaystyle=$ $\displaystyle-mc^{2}\Big{\\{}1-{\\!\\!\mbox{
\boldmath$v$}^{2}\over 2c^{2}}-{\varepsilon_{\hat{0}\hat{0}}\over
2}-\varepsilon_{\hat{0}\hat{j}}{v^{\hat{j}}\over c}-{1\over
2}\varepsilon_{\hat{i}\hat{j}}{v^{\hat{i}}v^{\hat{j}}\over c^{2}}-{1\over
8}\left[{\\!\\!\mbox{ \boldmath$v$}^{2}\over
c^{2}}+\varepsilon_{\hat{0}\hat{0}}+2\varepsilon_{\hat{0}\hat{j}}{v^{\hat{j}}\over
c}+\varepsilon_{\hat{i}\hat{j}}{v^{\hat{i}}v^{\hat{j}}\over
c^{2}}\right]^{2}+\cdots\Big{\\}}.$

Note that $x^{\hat{0}}=c\tau$. Inserting Eqs. (3), (4), and (5), and expanding
the above Lagrangian in terms of $v^{\hat{j}}/c$, we obtain

$\displaystyle{\cal L}={\cal L}_{0}+{\cal L}_{1/2}+{\cal L}_{1},$

where

$\displaystyle{\cal L}_{0}$ $\displaystyle=$ $\displaystyle{1\over
2}m\\!\\!\mbox{ \boldmath$v$}^{2}-ma_{\hat{k}}x^{\hat{k}}-{1\over
2}m\omega_{\hat{j}\hat{k}}\omega^{\hat{j}}_{~{}\hat{\ell}}x^{\hat{k}}x^{\hat{\ell}}-m\omega_{\hat{j}\hat{k}}v^{\hat{j}}x^{\hat{k}}-{1\over
2}mc^{2}\bar{\cal
R}_{\hat{0}\hat{k}\hat{0}\hat{\ell}}x^{\hat{k}}x^{\hat{\ell}},$
$\displaystyle{\cal L}_{1/2}$ $\displaystyle=$ $\displaystyle-{2\over
3}mc^{2}\bar{\cal
R}_{\hat{0}\hat{k}\hat{j}\hat{\ell}}x^{\hat{k}}x^{\hat{\ell}}\,{v^{\hat{j}}\over
c},$ $\displaystyle{\cal L}_{1}$ $\displaystyle=$ $\displaystyle-{m\over
2}{(a_{\hat{k}}x^{\hat{k}})^{2}\over c^{2}}-{1\over 6}mc^{2}\bar{\cal
R}_{\hat{i}\hat{k}\hat{j}\hat{\ell}}x^{\hat{k}}x^{\hat{\ell}}{v^{\hat{i}}v^{\hat{j}}\over
c^{2}}+{m\over 8c^{2}}\left[\\!\\!\mbox{
\boldmath$v$}^{2}-2a_{\hat{k}}x^{\hat{k}}-\omega_{\hat{j}\hat{k}}\omega^{\hat{j}}_{~{}\hat{l}}x^{\hat{k}}x^{\hat{\ell}}-2\omega_{\hat{j}\hat{k}}x^{\hat{k}}v^{\hat{j}}-c^{2}\bar{\cal
R}_{\hat{0}\hat{k}\hat{0}\hat{\ell}}x^{\hat{k}}x^{\hat{\ell}}\right]^{2}\,,$

are Newtonian, 0.5 PN and 1PN Lagrangian, respectively.

In this expansion, we find 0.5 PN term formally, but it can be an apparent
term which comes from a choice of an observer’s coordinates. In fact, if we
choose an appropriate observer’s acceleration, which appears as 0.5 PN term,
we can remove the above ${\cal L}_{1/2}$. However, if we have multi-particle
system as we will discuss later, this adjustment can be used only for one
particle or the center of mass of the system (See §. III). Hence, we keep the
${\cal L}_{1/2}$ term and discuss “Newtonian” dynamics including such a term.

##  III Self-gravitating Newtonian system in a curved spacetime

###  III.1 Self-gravitating $N$-body system

Now we discuss self-gravitating $N$-body system in a fixed curved background
spacetime, which is given by the metric in Eq. (1) Gorbatsievich and Bobrik
(2010). We are interested in the case where Newtonian dynamics is valid in the
observer’s proper reference frame. The necessary condition is that the typical
scale $\ell_{N\text{-body}}$ of $N$-body system should satisfy

$\displaystyle\ell_{N\text{-body}}\ll{\rm
min}\left[{1\over|a^{\hat{j}}|},{1\over|\omega^{\hat{j}}|},\ell_{\bar{\cal
R}}\right],$

where $\ell_{\bar{\cal R}}$ is the minimum curvature radius defined by

$\displaystyle\ell_{\bar{\cal R}}\equiv{\rm min}\left[|\bar{\cal
R}_{\hat{\mu}\hat{\nu}\hat{\rho}\hat{\sigma}}|^{-{1\over 2}},|\bar{\cal
R}_{\hat{\mu}\hat{\nu}\hat{\rho}\hat{\sigma};\hat{\alpha}}|^{-{1\over
3}},|\bar{\cal
R}_{\hat{\mu}\hat{\nu}\hat{\rho}\hat{\sigma};\hat{\alpha};\hat{\beta}}|^{-{1\over
4}}\right]\,.$

In order to find the metric contributions from $N$-body system, we first focus
on the motion of the $I$th particle, which is gravitating with the other
$(N-1)$ particles. The metric contribution from the those $(N-1)$ particles is
given by

$\displaystyle\varphi^{I}_{\hat{0}\hat{0}}=-{2\Phi_{I}\over c^{2}}\,,$

where $\Phi_{I}$ is the Newtonian potential produced by the $(N-1)$ particles,
which is given by

$\displaystyle\Phi_{I}(x^{\hat{i}})=-G\sum_{J\neq
I}^{N}{m_{J}\over|x^{\hat{i}}-x_{J}^{\hat{i}}|}\,.$

We assume that the other components of $\varphi^{I}_{\hat{\mu}\hat{\nu}}$
vanish because we are interested in Newtonian dynamics in the proper reference
frame.

We then obtain the metric of the observer’s proper reference frame for the
$I$th particle as

$\displaystyle
g^{I}_{\hat{\mu}\hat{\nu}}=\eta_{\hat{\mu}\hat{\nu}}+h^{I}_{\hat{\mu}\hat{\nu}},$
(6)

where

$\displaystyle
h^{I}_{\hat{\mu}\hat{\nu}}=\varepsilon_{\hat{\mu}\hat{\nu}}+\varphi^{I}_{\hat{\mu}\hat{\nu}}\,.$
(7)

The equation of motion for the $I$th particle can be derived by the variation
of the action

$\displaystyle{\cal S}_{I}=\int d\tau{\cal L}_{I},$

where

$\displaystyle{\cal L}_{I}\equiv-
m_{I}c\sqrt{-g^{I}_{\hat{\mu}\hat{\nu}}{dx_{I}^{\hat{\mu}}\over
d\tau}{dx_{I}^{\hat{\nu}}\over d\tau}}\,.$

By use of the metric form (6), we can expand the above Lagrangian up to 0.5 PN
order as

$\displaystyle{\cal L}_{I}$ $\displaystyle=$ $\displaystyle{1\over
2}m_{I}\left({d\\!\\!\mbox{ \boldmath$x$}_{I}\over
d\tau}\right)^{2}-m_{I}\Phi(x_{I})-m_{I}a_{\hat{k}}x_{I}^{\hat{k}}-{1\over
2}m_{I}\omega_{\hat{j}\hat{k}}\omega^{\hat{j}}_{~{}\hat{\ell}}x_{I}^{\hat{k}}x_{I}^{\hat{\ell}}-m_{I}\omega_{\hat{j}\hat{k}}v_{I}^{\hat{j}}x_{I}^{\hat{k}}$
$\displaystyle-{1\over 2}m_{I}c^{2}\bar{\cal
R}_{\hat{0}\hat{k}\hat{0}\hat{\ell}}x_{I}^{\hat{k}}x_{I}^{\hat{\ell}}-{2\over
3}m_{I}c^{2}\bar{\cal
R}_{\hat{0}\hat{k}\hat{j}\hat{\ell}}x_{I}^{\hat{k}}x_{I}^{\hat{\ell}}\,{v_{I}^{\hat{j}}\over
c}.$

The total action of $N$-body system and its Lagrangian are given by summing up
each Lagrangian ${\cal L}_{I}$. We finally obtain

$\displaystyle{\cal S}_{N\mathchar 45\relax{\rm body}}=\int d\tau{\cal
L}_{N\mathchar 45\relax{\rm body}},$

where

$\displaystyle{\cal L}_{N\mathchar 45\relax{\rm body}}={1\over
2}\sum_{I}m_{I}\left({d\\!\\!\mbox{ \boldmath$x$}_{I}\over
d\tau}\right)^{2}+\sum^{N}_{I}\sum^{N}_{J\neq I}{Gm_{I}m_{J}\over
2|x_{I}^{\hat{i}}-x_{J}^{\hat{i}}|}+{\cal L}_{a}+{\cal L}_{\omega}+{\cal
L}_{\bar{\cal R}},$

with

$\displaystyle{\cal L}_{a}$ $\displaystyle=$
$\displaystyle-\sum_{I}^{N}m_{I}a_{\hat{k}}x_{I}^{\hat{k}},$
$\displaystyle{\cal L}_{\omega}$ $\displaystyle=$ $\displaystyle-{1\over
2}\sum_{I}^{N}m_{I}\left[\omega_{\hat{j}\hat{k}}\omega^{\hat{j}}_{~{}\hat{\ell}}x_{I}^{\hat{k}}x_{I}^{\hat{\ell}}+2\omega_{\hat{j}\hat{k}}v_{I}^{\hat{j}}x_{I}^{\hat{k}}\right],$
$\displaystyle{\cal L}_{\bar{\cal R}}$ $\displaystyle=$ $\displaystyle-{1\over
2}\sum_{I}^{N}m_{I}c^{2}\bar{\cal
R}_{\hat{0}\hat{k}\hat{0}\hat{\ell}}x_{I}^{\hat{k}}x_{I}^{\hat{\ell}}-{2\over
3}\sum_{I}m_{I}c^{2}\bar{\cal
R}_{\hat{0}\hat{k}\hat{j}\hat{\ell}}x_{I}^{\hat{k}}x_{I}^{\hat{\ell}}\,{v_{I}^{\hat{j}}\over
c}\,.$

Note that ${\cal L}_{a}$ comes from the inertial force of the accelerated
observer, and ${\cal L}_{\omega}$ originates in the rotation of the observer
(the centrifugal force and the Coriolis force). ${\cal L}_{\bar{\cal R}}$
describes the curvature effect of the third body (the tidal force).

###  III.2 Binary system in a curved spacetime

Next, we discuss a binary system in a fixed curved background. A binary
consists of two point particles with the masses $m_{1}$ and $m_{2}$. The
Lagrangian up to 0.5 PN order is given by

$\displaystyle{\cal L}_{\rm binary}={\cal L}_{\rm N}+{\cal L}_{1/2},$ (8)

where

$\displaystyle{\cal L}_{\rm N}$ $\displaystyle\equiv$ $\displaystyle{1\over
2}\left[m_{1}\left({d\\!\\!\mbox{ \boldmath$x$}_{1}\over
d\tau}\right)^{2}+m_{2}\left({d\\!\\!\mbox{ \boldmath$x$}_{2}\over
d\tau}\right)^{2}\right]+{Gm_{1}m_{2}\over|\\!\\!\mbox{
\boldmath$x$}_{1}-\\!\\!\mbox{ \boldmath$x$}_{2}|}+{\cal L}_{a}+{\cal
L}_{\omega}+{\cal L}_{\bar{\cal R}},$ (9)

with

$\displaystyle{\cal L}_{a}$ $\displaystyle=$
$\displaystyle-\sum_{I=1}^{2}m_{I}a_{\hat{k}}x_{I}^{\hat{k}},$
$\displaystyle{\cal L}_{\omega}$ $\displaystyle=$
$\displaystyle-\sum_{I=1}^{2}m_{I}\left[\epsilon_{\hat{j}\hat{k}\hat{\ell}}\omega^{\hat{\ell}}x_{I}^{\hat{k}}{dx_{I}^{\hat{j}}\over
d\tau}-{1\over 2}\left(\\!\\!\mbox{ \boldmath$\omega$}^{2}\\!\\!\mbox{
\boldmath$x$}_{I}^{2}-(\\!\\!\mbox{ \boldmath$\omega$}\cdot\\!\\!\mbox{
\boldmath$x$}_{I})^{2}\right)\right],$ $\displaystyle{\cal L}_{\bar{\cal R}}$
$\displaystyle=$ $\displaystyle-{1\over 2}\sum_{I=1}^{2}m_{I}\bar{\cal
R}_{\hat{0}\hat{k}\hat{0}\hat{\ell}}x_{I}^{\hat{k}}x_{I}^{\hat{\ell}}\,,$

and

$\displaystyle{\cal L}_{1/2}$ $\displaystyle\equiv$ $\displaystyle-{2\over
3}\sum_{I=1}^{2}m_{I}c^{2}\bar{\cal
R}_{\hat{0}\hat{k}\hat{j}\hat{\ell}}x_{I}^{\hat{k}}x_{I}^{\hat{\ell}}\,{v_{I}^{\hat{j}}\over
c}\,.$ (10)

Introducing the center of mass coordinates and the relative coordinates by

$R$ $\displaystyle=$ $\displaystyle{m_{1}\\!\\!\mbox{
\boldmath$x$}_{1}+m_{2}\\!\\!\mbox{ \boldmath$x$}_{2}\over m_{1}+m_{2}},$ $r$
$\displaystyle=$ $\displaystyle\\!\\!\mbox{ \boldmath$x$}_{2}-\\!\\!\mbox{
\boldmath$x$}_{1}\,,$

we find the Newtonian Lagrangian (Eq. (9)) in terms of  $R$ and  $r$ as

$\displaystyle{\cal L}_{\rm N}={\cal L}_{\rm CM}(\\!\\!\mbox{
\boldmath$R$},\dot{\\!\\!\mbox{ \boldmath$R$}})+{\cal L}_{\rm
rel}(\\!\\!\mbox{ \boldmath$r$},\dot{\\!\\!\mbox{ \boldmath$r$}})\,,$ (11)

where

$\displaystyle{\cal L}_{\rm CM}(\\!\\!\mbox{ \boldmath$R$},\dot{\\!\\!\mbox{
\boldmath$R$}})$ $\displaystyle=$ $\displaystyle{1\over
2}(m_{1}+m_{2})\dot{\\!\\!\mbox{ \boldmath$R$}}^{2}+{\cal L}_{{\rm
CM}\mathchar 45\relax a}(\\!\\!\mbox{ \boldmath$R$},\dot{\\!\\!\mbox{
\boldmath$R$}})+{\cal L}_{{\rm CM}\mathchar 45\relax\omega}(\\!\\!\mbox{
\boldmath$R$},\dot{\\!\\!\mbox{ \boldmath$R$}})+{\cal L}_{{\rm CM}\mathchar
45\relax\bar{\cal R}}(\\!\\!\mbox{ \boldmath$R$},\dot{\\!\\!\mbox{
\boldmath$R$}}),$

with

$\displaystyle{\cal L}_{{\rm CM}\mathchar 45\relax a}$ $\displaystyle=$
$\displaystyle-(m_{1}+m_{2})a_{\hat{k}}x_{I}^{\hat{k}},$ $\displaystyle{\cal
L}_{{\rm CM}\mathchar 45\relax\omega}$ $\displaystyle=$
$\displaystyle-(m_{1}+m_{2})\left[\epsilon_{\hat{j}\hat{k}\hat{\ell}}\omega^{\hat{\ell}}R^{\hat{k}}\dot{R}^{\hat{j}}-{1\over
2}\left(\\!\\!\mbox{ \boldmath$\omega$}^{2}\\!\\!\mbox{
\boldmath$R$}^{2}-\left(\\!\\!\mbox{ \boldmath$\omega$}\cdot\\!\\!\mbox{
\boldmath$R$}\right)^{2}\right)\right],$ $\displaystyle{\cal L}_{{\rm
CM}\mathchar 45\relax\bar{\cal R}}$ $\displaystyle=$ $\displaystyle-{1\over
2}(m_{1}+m_{2})\bar{\cal
R}_{\hat{0}\hat{k}\hat{0}\hat{\ell}}R^{\hat{k}}R^{\hat{\ell}},$

and

$\displaystyle{\cal L}_{\rm rel}(\\!\\!\mbox{ \boldmath$r$},\dot{\\!\\!\mbox{
\boldmath$r$}})$ $\displaystyle=$ $\displaystyle{1\over 2}\mu\dot{\\!\\!\mbox{
\boldmath$r$}}^{2}+{Gm_{1}m_{2}\over r}+{\cal L}_{{\rm rel}\mathchar
45\relax\omega}(\\!\\!\mbox{ \boldmath$r$},\dot{\\!\\!\mbox{
\boldmath$r$}})+{\cal L}_{{\rm rel}\mathchar 45\relax\bar{\cal
R}}(\\!\\!\mbox{ \boldmath$r$},\dot{\\!\\!\mbox{ \boldmath$r$}}),$

with

$\displaystyle{\cal L}_{{\rm rel}\mathchar 45\relax\omega}$ $\displaystyle=$
$\displaystyle-\mu\left[\epsilon_{\hat{j}\hat{k}\hat{\ell}}\omega^{\hat{\ell}}r^{\hat{k}}\dot{r}^{\hat{j}}-{1\over
2}\left(\\!\\!\mbox{ \boldmath$\omega$}^{2}\\!\\!\mbox{
\boldmath$r$}^{2}-\left(\\!\\!\mbox{ \boldmath$\omega$}\cdot\\!\\!\mbox{
\boldmath$r$}\right)^{2}\right)\right]\,,$ $\displaystyle{\cal L}_{{\rm
rel}\mathchar 45\relax\bar{\cal R}}$ $\displaystyle=$ $\displaystyle-{1\over
2}\mu\bar{\cal
R}_{\hat{0}\hat{k}\hat{0}\hat{\ell}}r^{\hat{k}}r^{\hat{\ell}}\,.$

Here, $\mu=m_{1}m_{2}/(m_{1}+m_{2})$ is the reduced mass. When we consider
only ${\cal L}_{\rm N}$, we can separate the variables  $R$ and  $r$. In
particular, when the observer follows the geodesic ($\\!\\!\mbox{
\boldmath$a$}=0$ and $\\!\\!\mbox{ \boldmath$\omega$}=0$), the orbit of
$\\!\\!\mbox{ \boldmath$R$}=0$ is a solution of the equation for  $R$. It
means that the center of mass (CM) follows the observer’s geodesic. We have
only the equation for the relative coordinate  $r$. However, when we include
the 0.5 PN term, it is not the case. The 0.5PN Lagrangian ${\cal L}_{1/2}$ is
rewritten by use of  $R$ and  $r$ as

$\displaystyle{\cal L}_{1/2}={\cal L}_{1/2\mathchar 45\relax{\rm
CM}}(\\!\\!\mbox{ \boldmath$R$},\dot{\\!\\!\mbox{ \boldmath$R$}})+{\cal
L}_{1/2\mathchar 45\relax{\rm rel}}(\\!\\!\mbox{
\boldmath$r$},\dot{\\!\\!\mbox{ \boldmath$r$}})+{\cal L}_{1/2\mathchar
45\relax{\rm int}}(\\!\\!\mbox{ \boldmath$R$},\dot{\\!\\!\mbox{
\boldmath$R$},}\\!\\!\mbox{ \boldmath$r$},\dot{\\!\\!\mbox{ \boldmath$r$}}),$
(12)

where

$\displaystyle{\cal L}_{1/2\mathchar 45\relax{\rm CM}}(\\!\\!\mbox{
\boldmath$R$},\dot{\\!\\!\mbox{ \boldmath$R$}})$ $\displaystyle=$
$\displaystyle-{2\over
3}(m_{1}+m_{2})R_{\hat{0}\hat{k}\hat{j}\hat{\ell}}R^{\hat{k}}R^{\hat{\ell}}\dot{R}^{\hat{j}},$
$\displaystyle{\cal L}_{1/2\mathchar 45\relax{\rm rel}}(\\!\\!\mbox{
\boldmath$r$},\dot{\\!\\!\mbox{ \boldmath$r$}})$ $\displaystyle=$
$\displaystyle-{2\over
3}\mu{(m_{1}-m_{2})\over(m_{1}+m_{2})}R_{\hat{0}\hat{k}\hat{j}\hat{\ell}}r^{\hat{k}}r^{\hat{\ell}}\dot{r}^{\hat{j}},$
$\displaystyle{\cal L}_{1/2\mathchar 45\relax{\rm int}}(\\!\\!\mbox{
\boldmath$R$},\dot{\\!\\!\mbox{ \boldmath$R$},}\\!\\!\mbox{
\boldmath$r$},\dot{\\!\\!\mbox{ \boldmath$r$}})$ $\displaystyle=$
$\displaystyle-{2\over 3}\mu
R_{\hat{0}\hat{k}\hat{j}\hat{\ell}}\left[r^{\hat{k}}r^{\hat{\ell}}\dot{R}^{\hat{j}}+\left(R^{\hat{k}}r^{\hat{\ell}}+r^{\hat{k}}R^{\hat{\ell}}\right)\dot{r}^{\hat{j}}\right].$
(13)

Due of the interaction term ${\cal L}_{1/2\mathchar 45\relax{\rm int}}$, the
orbit of $\\!\\!\mbox{ \boldmath$R$}=0$ is no longer a solution even if the
acceleration vanishes. The motion of the CM ($\\!\\!\mbox{
\boldmath$R$}(\tau$)) couples with the relative motion $(\\!\\!\mbox{
\boldmath$r$}(\tau))$. As a result, not only the orbit of a binary but also
the motion of the CM will become complicated even if the observer’s orbit is a
geodesic.

However, if we introduce an appropriate acceleration  $a$ in 0.5PN order to
cancel the interaction terms, $\\!\\!\mbox{ \boldmath$R$}=0$ will become a
solution, i.e., the CM can follow the observer’s motion as follows:
Integrating by parts the interaction term, we find

$\displaystyle{\cal L}_{1/2\mathchar 45\relax{\rm int}}(\\!\\!\mbox{
\boldmath$R$},\dot{\\!\\!\mbox{ \boldmath$R$},}\\!\\!\mbox{
\boldmath$r$},\dot{\\!\\!\mbox{ \boldmath$r$}})$ $\displaystyle=$
$\displaystyle-{2\over 3}\mu\bar{\cal
R}_{\hat{0}\hat{k}\hat{j}\hat{\ell}}\left[\dot{R}^{\hat{j}}r^{\hat{k}}r^{\hat{\ell}}+\dot{r}^{\hat{j}}\left(R^{\hat{k}}r^{\hat{\ell}}+r^{\hat{k}}R^{\hat{\ell}}\right)\right]$
$\displaystyle\approx$ $\displaystyle 2\mu\left[{1\over 3}{d\bar{\cal
R}_{\hat{0}\hat{k}\hat{j}\hat{\ell}}\over
d\tau}r^{\hat{k}}r^{\hat{\ell}}+\bar{\cal
R}_{\hat{0}\hat{k}\hat{j}\hat{\ell}}r^{\hat{k}}\dot{r}^{\hat{\ell}}\right]R^{\hat{j}}~{}~{}{\rm(integration~{}by~{}part)}\,,$

where the time derivative of the curvature is evaluated along the observer’s
orbit.

If we define the acceleration by

$\displaystyle a_{\hat{j}}={2\mu\over m_{1}+m_{2}}\left[{1\over 3}{d\bar{\cal
R}_{\hat{0}\hat{k}\hat{j}\hat{\ell}}\over
d\tau}r^{\hat{k}}r^{\hat{\ell}}+\bar{\cal
R}_{\hat{0}\hat{k}\hat{j}\hat{\ell}}r^{\hat{k}}\dot{r}^{\hat{\ell}}\right]\,,$

two terms ${\cal L}_{1/2\mathchar 45\relax{\rm int}}$ and ${\cal L}_{{\rm
CM}\mathchar 45\relax a}$ cancel each other. As a result, the Lagrangians for
$R$ and  $r$ are decoupled, and $\\!\\!\mbox{ \boldmath$R$}=0$ becomes an
exact solution of the equation for  $R$, which is derived from the Lagrangian
(${\cal L}_{\rm CM}+{\cal L}_{1/2\mathchar 45\relax{\rm CM}}$). The CM follows
the observer’s orbit and therefore, we obtain the decoupled equation for the
relative coordinate  $r$.

In order to obtain the proper observer’s orbit, which is not a geodesic but
may be close to the geodesic, we have to solve the equation of motion
including small acceleration such that

$\displaystyle{Du_{\rm CM}^{\mu}\over
d\tau}=a^{\mu}=e^{\mu\hat{j}}a_{\hat{j}}={2\mu\over
m_{1}+m_{2}}e^{\mu\hat{j}}\left[{1\over 3}{d\bar{\cal
R}_{\hat{0}\hat{k}\hat{j}\hat{\ell}}\over
d\tau}r^{\hat{k}}r^{\hat{\ell}}+\bar{\cal
R}_{\hat{0}\hat{k}\hat{j}\hat{\ell}}r^{\hat{k}}\dot{r}^{\hat{\ell}}\right]\,.$
(14)

As a result, we first solve the equation for the relative coordinate  $r$,
which is obtained only by the Lagrangian ${\cal L}_{\rm rel}(\\!\\!\mbox{
\boldmath$r$})+{\cal L}_{1/2\mathchar 45\relax{\rm rel}}(\\!\\!\mbox{
\boldmath$r$})$. Note that when $m_{1}=m_{2}$, we have only Newtonian
Lagrangian ${\cal L}_{\rm rel}$ because ${\cal L}_{1/2\mathchar 45\relax{\rm
rel}}=0$ vanishes. After obtaining the solution of $\\!\\!\mbox{
\boldmath$r$}(\tau)$, we find the motion for the CM (or the observer) in the
background spacetime by solving Eq. (14). Using the relative motion
$\\!\\!\mbox{ \boldmath$r$}(\tau)$ with the solution of the CM motion ($x_{\rm
CM}^{\mu}(\tau)$), we will obtain a binary motion in a given curved background
spacetime ($x_{1}^{\mu}(\tau)\,,x_{2}^{\mu}(\tau)$).

##  IV A binary system in a Schwarzschild spacetime

We consider a spherically symmetric supermassive black hole. The background
spacetime is given by Schwarzschild solution as

$\displaystyle d\bar{s}^{2}=-fd\mathfrak{t}^{2}+{1\over
f}d\mathfrak{r}^{2}+\mathfrak{r}^{2}d\Omega^{2}\,,$ (15)

where

$\displaystyle f=1-{\mathfrak{r}_{g}\over\mathfrak{r}}.$ (16)

The gravitational radius $\mathfrak{r}_{g}$ is given by

$\displaystyle\mathfrak{r}_{g}\equiv{2GM\over c^{2}},$

where $M$ is a gravitational mass of the supermassive black hole. In what
follows, we set $G=1$ and $c=1$ for brevity unless specified otherwise.

###  IV.1 Proper reference frame

We first discuss Newtonian dynamics for which Lagrangian ${\cal L}_{\rm N}$ is
given by Eq. (9). In this case, the CM of a binary system follows a geodesic
observer. We consider an observer, which moves along a circular geodesic with
the radius $\mathfrak{r}=\mathfrak{r}_{0}$. The orbit is assumed to be on the
equatorial plane without loss of generality. We introduce a rotating tetrad
system along this geodesic such that

$\displaystyle e_{\hat{0}}^{~{}\mu}$ $\displaystyle\equiv$ $\displaystyle
u^{\mu}=\left({\sqrt{1+\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2}}\over\sqrt{f(\mathfrak{r}_{0})}}\,,0\,,0\,,\mathfrak{w}_{0}\right),$
$\displaystyle e_{\hat{1}}^{~{}\mu}$ $\displaystyle=$
$\displaystyle\left(0\,,\sqrt{f(\mathfrak{r}_{0})}\,,0\,,0\right),$
$\displaystyle e_{\hat{2}}^{~{}\mu}$ $\displaystyle=$
$\displaystyle\left({\mathfrak{r}_{0}\mathfrak{w}_{0}\over\sqrt{f(\mathfrak{r}_{0})}}\,,0\,,0\,,{\sqrt{1+\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2}}\over\mathfrak{r}_{0}}\right),$
$\displaystyle e_{\hat{3}}^{~{}\mu}$ $\displaystyle=$
$\displaystyle\left(0\,,0\,,-{1\over\mathfrak{r}_{0}}\,,0\right)\,,$

where the angular frequency $\mathfrak{w}_{0}$ is defined by

$\displaystyle\mathfrak{w}_{0}=u^{3}={d\varphi\over
d\tau}={1\over\mathfrak{r}_{0}(\mathfrak{r}_{0}/M-3)^{1/2}}\,.$ (17)

We then calculate the angular velocity $\omega^{\hat{j}}$ of this rotating
tetrad system as

$\displaystyle\omega^{\hat{j}}={1\over
2}\stackrel{{\scriptstyle(3)}}{{\epsilon}}^{\hat{j}\hat{k}\hat{\ell}}e_{\hat{\ell}\mu}{De_{\hat{k}}^{~{}\mu}\over
d\tau}=\mathfrak{w}_{\rm R}\delta^{\hat{j}}_{~{}\hat{3}}\,,$

where

$\displaystyle\mathfrak{w}_{\rm
R}={1\over\mathfrak{r}_{0}(\mathfrak{r}_{0}/M)^{1/2}},$ (18)

is the angular frequency of the rotating tetrad frame (See Fig. 1). Note that
this angular frequency is different from the observer’s angular frequency
$\mathfrak{w}_{0}$. The difference between two angular frequencies
$\mathfrak{w}_{0}$ and $\mathfrak{w}_{R}$, i.e.,

$\displaystyle\mathfrak{w}_{\rm dS}\equiv\mathfrak{w}_{0}-\mathfrak{w}_{\rm
R}$ (19)

describes the rotation of non-rotating inertial frame moving along a circular
orbit, which may cause de Sitter precession de Sitter (1916).

Figure 1: A local inertial tetrad system $\\{\\!\\!\mbox{
\boldmath$e$}_{\hat{x}}\,,\\!\\!\mbox{ \boldmath$e$}_{\hat{y}}\,,\\!\\!\mbox{
\boldmath$e$}_{\hat{z}}\\}$ rotating with angular frequency $\mathfrak{w}_{\rm
R}$ along a circular orbit.

Next, we calculate the Riemann curvature in the above tetrad system. In the
static tetrad system $\\{e^{\bar{\alpha}}_{~{}\mu}\\}$ with

$\displaystyle
e^{\bar{0}}_{~{}0}=\sqrt{f(\mathfrak{r}_{0})}\,,~{}e^{\bar{r}}_{~{}r}={1\over\sqrt{f(\mathfrak{r}_{0})}}\,,~{}e^{\bar{\theta}}_{~{}\theta}=\mathfrak{r}_{0}\,,~{}e^{\bar{\phi}}_{~{}\phi}=\mathfrak{r}_{0}\,,$

the non-trivial components of the Riemann curvature of the Schwarzschild
spacetime are given by

$\displaystyle\bar{\cal
R}_{\bar{0}\bar{r}\bar{0}\bar{r}}=-{\mathfrak{r}_{g}\over\mathfrak{r}_{0}^{3}}\,,~{}\bar{\cal
R}_{\bar{0}\bar{\theta}\bar{0}\bar{\theta}}=\bar{\cal
R}_{\bar{0}\bar{\phi}\bar{0}\bar{\phi}}={\mathfrak{r}_{g}\over
2\mathfrak{r}_{0}^{3}},$ $\displaystyle\bar{\cal
R}_{\bar{\theta}\bar{\phi}\bar{\theta}\bar{\phi}}={\mathfrak{r}_{g}\over\mathfrak{r}_{0}^{3}}\,,~{}\bar{\cal
R}_{\bar{r}\bar{\theta}\bar{r}\bar{\theta}}=R_{\bar{r}\bar{\phi}\bar{r}\bar{\phi}}=-{\mathfrak{r}_{g}\over
2\mathfrak{r}_{0}^{3}}\,.$

We then introduce the Descartes coordinates $(\bar{x},\bar{y},\bar{z})$ such
that $\bar{x}$-direction is the same as the $r$-direction, but the $\bar{y}$
and $\bar{z}$-directions are $\phi$ and $-\theta$-directions, respectively, as
shown in Fig. 1. Since the tetrad system in this coordinates is given by

$\displaystyle
e^{\bar{0}}_{~{}0}=\sqrt{f(\mathfrak{r}_{0})}\,,~{}e^{\bar{x}}_{~{}r}={1\over\sqrt{f(\mathfrak{r}_{0})}}\,,~{}e^{\bar{y}}_{~{}\phi}=\mathfrak{r}_{0}\,,~{}e^{\bar{z}}_{~{}\theta}=-\mathfrak{r}_{0}\,,$

we find the curvature components in this tetrad system as

$\displaystyle\bar{\cal
R}_{\bar{0}\bar{x}\bar{0}\bar{x}}=-{\mathfrak{r}_{g}\over\mathfrak{r}_{0}^{3}}\,,~{}\bar{\cal
R}_{\bar{0}\bar{y}\bar{0}\bar{y}}=\bar{\cal
R}_{\bar{0}\bar{z}\bar{0}\bar{z}}={\mathfrak{r}_{g}\over
2\mathfrak{r}_{0}^{3}},$ $\displaystyle\bar{\cal
R}_{\bar{y}\bar{z}\bar{y}\bar{z}}={\mathfrak{r}_{g}\over\mathfrak{r}_{0}^{3}}\,,~{}R_{\bar{x}\bar{y}\bar{x}\bar{y}}=\bar{\cal
R}_{\bar{z}\bar{x}\bar{z}\bar{x}}=-{\mathfrak{r}_{g}\over
2\mathfrak{r}_{0}^{3}}\,.$

The transformation matrix between the observer’s tetrad and static tetrad
($\bar{x},\bar{y},\bar{z}$) is given by

$\displaystyle\Lambda_{\hat{0}}^{~{}\bar{\alpha}}$ $\displaystyle=$
$\displaystyle\left(\sqrt{1+\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2}}\,,0\,,\mathfrak{r}_{0}\mathfrak{w}_{0}\,,0\right),$
$\displaystyle\Lambda_{\hat{1}}^{~{}\bar{\alpha}}$ $\displaystyle=$
$\displaystyle\left(0\,,1\,,0\,,0\right),$
$\displaystyle\Lambda_{\hat{2}}^{~{}\bar{\alpha}}$ $\displaystyle=$
$\displaystyle\left(\mathfrak{r}_{0}\mathfrak{w}_{0}\,,0\,,\sqrt{1+\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2}}\,,0\right),$
$\displaystyle\Lambda_{\hat{3}}^{~{}\bar{\alpha}}$ $\displaystyle=$
$\displaystyle\left(0\,,0\,,0\,,1\right).$

It is just the Lorentz boost with velocity

$\displaystyle\mathfrak{v}={\mathfrak{r}_{0}\mathfrak{w}_{0}\over\sqrt{1+\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2}}},$

in $0\,\mathchar 45\relax\,y$ plane. The components in the observer’s proper
reference frame are given by

$\displaystyle\bar{\cal R}_{\hat{0}\hat{x}\hat{0}\hat{x}}=-\bar{\cal
R}_{\hat{y}\hat{z}\hat{y}\hat{z}}=-\left(1+{3\over
2}\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2}\right){\mathfrak{r}_{g}\over\mathfrak{r}_{0}^{3}}\,,~{}$
$\displaystyle\bar{\cal R}_{\hat{0}\hat{y}\hat{0}\hat{y}}=-\bar{\cal
R}_{\hat{z}\hat{x}\hat{z}\hat{x}}={1\over
2}{\mathfrak{r}_{g}\over\mathfrak{r}_{0}^{3}}\,,~{}$ $\displaystyle\bar{\cal
R}_{\hat{0}\hat{z}\hat{0}\hat{z}}=-\bar{\cal
R}_{\hat{x}\hat{y}\hat{x}\hat{y}}=\left({1+3\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2}\over
2}\right){\mathfrak{r}_{g}\over\mathfrak{r}_{0}^{3}},$ $\displaystyle\bar{\cal
R}_{\hat{0}\hat{x}\hat{y}\hat{x}}=-\bar{\cal
R}_{\hat{0}\hat{z}\hat{y}\hat{z}}=-{3\over
2}\epsilon_{1}\mathfrak{r}_{0}\mathfrak{w}_{0}\sqrt{1+\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2}}\,{\mathfrak{r}_{g}\over\mathfrak{r}_{0}^{3}}\,.$

###  IV.2 Equations of motion of a binary

Since the CM of a binary follows the observer’s circular geodesic
($\\!\\!\mbox{ \boldmath$R$}=0$), we have to solve only the equations of
motion for the relative coordinate  $r$. Using
$x=r^{\hat{1}}\,,y=r^{\hat{2}}\,,z=r^{\hat{3}}$, the relative motion of a
binary is given by the Lagrangian

$\displaystyle{\cal L}_{\rm rel}(\\!\\!\mbox{ \boldmath$r$},\dot{\\!\\!\mbox{
\boldmath$r$}})$ $\displaystyle=$ $\displaystyle{1\over 2}\mu\dot{\\!\\!\mbox{
\boldmath$r$}}^{2}+{Gm_{1}m_{2}\over r}+{\cal L}_{{\rm rel}\mathchar
45\relax\mathfrak{w}}(\\!\\!\mbox{ \boldmath$r$},\dot{\\!\\!\mbox{
\boldmath$r$}})+{\cal L}_{{\rm rel}\mathchar 45\relax\bar{\cal
R}}(\\!\\!\mbox{ \boldmath$r$},\dot{\\!\\!\mbox{ \boldmath$r$}}),$ (20)

with

$\displaystyle{\cal L}_{{\rm rel}\mathchar 45\relax\mathfrak{w}}$
$\displaystyle=$
$\displaystyle-\mu\left[\mathfrak{w}_{0}\left(x\dot{y}-y\dot{x}\right)-{\mathfrak{w}_{0}^{2}\over
2}\left(x^{2}+y^{2}\right)\right],$ $\displaystyle{\cal L}_{{\rm rel}\mathchar
45\relax\bar{\cal R}}$ $\displaystyle=$ $\displaystyle-{\mu\over
2}\left(\bar{\cal R}_{\hat{0}\hat{x}\hat{0}\hat{x}}x^{2}+\bar{\cal
R}_{\hat{0}\hat{y}\hat{0}\hat{y}}y^{2}+\bar{\cal
R}_{\hat{0}\hat{z}\hat{0}\hat{z}}z^{2}\right)$ $\displaystyle=$
$\displaystyle-{\mu\mathfrak{r}_{g}\over
4\mathfrak{r}_{0}^{3}}\left[y^{2}-2x^{2}+z^{2}+3\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2}(-x^{2}+z^{2})\right]\,.$

The first and second terms in ${\cal L}_{{\rm rel}\mathchar
45\relax\mathfrak{w}}$ describe the Coriolis force and the centrifugal force,
respectively. The first half terms in ${\cal L}_{{\rm rel}\mathchar
45\relax\bar{\cal R}}$ are the same as those in Newtonian hierarchical triple
system under quadrupole approximation, while the last half terms are
relativistic corrections. Note that in the present approach (approximation up
to the second order of $x^{\hat{a}}$), we cannot go beyond quadrupole
approximation.

In order to analyze the relative motion of a binary, it may be better to go to
a non-rotating initial reference frame. Since the angular frequency of a
rotating tetrad frame is $\mathfrak{w}_{\rm R}$, the position $(x,y,z)$ in the
rotating frame can be replaced by the position
$(\mathsf{x},\mathsf{y},\mathsf{z})$ in a non-rotating Descartes’ coordinate
system by use of the following transformation;

$\displaystyle x$ $\displaystyle=$
$\displaystyle\mathsf{x}\cos\mathfrak{w}_{\rm
R}\tau-\mathsf{y}\sin\mathfrak{w}_{\rm R}\tau,$ $\displaystyle y$
$\displaystyle=$ $\displaystyle\mathsf{x}\sin\mathfrak{w}_{\rm
R}\tau+\mathsf{y}\cos\mathfrak{w}_{\rm R}\tau,$ $\displaystyle z$
$\displaystyle=$ $\displaystyle\mathsf{z}\,.$

The Lagrangian ${\cal L}_{\rm rel}$ in a non-rotating proper reference frame
is given by

$\displaystyle{\cal L}_{\rm rel}={1\over 2}\mu\left({d{\\!\\!\mbox{
\boldmath$\mathsf{r}$}}\over
d\tau}\right)^{2}+{Gm_{1}m_{2}\over\mathsf{r}}+{\cal L}_{\rm rel\mathchar
45\relax dS}(\\!\\!\mbox{ \boldmath$\mathsf{r}$},\dot{\\!\\!\mbox{
\boldmath$\mathsf{r}$}})+{\cal L}_{{\rm rel}\mathchar 45\relax\bar{\cal
R}}(\\!\\!\mbox{ \boldmath$\mathsf{r}$},\tau)\,,$ (21)

where

$\displaystyle{\cal L}_{\rm rel\mathchar 45\relax dS}(\\!\\!\mbox{
\boldmath$\mathsf{r}$},\dot{\\!\\!\mbox{ \boldmath$\mathsf{r}$}})$
$\displaystyle=$ $\displaystyle\mu\mathfrak{w}_{\rm
dS}\left(\dot{\mathsf{x}}\mathsf{y}-\dot{\mathsf{y}}\mathsf{x}\right)+{\mu\over
2}\mathfrak{w}_{\rm dS}^{2}\left(\mathsf{x}^{2}+\mathsf{y}^{2}\right),$
$\displaystyle{\cal L}_{{\rm rel}\mathchar 45\relax\bar{\cal R}}(\\!\\!\mbox{
\boldmath$\mathsf{r}$},\tau)$ $\displaystyle=$
$\displaystyle-{\mu\mathfrak{r}_{g}\over
4\mathfrak{r}_{0}^{3}}\left[\mathsf{x}^{2}+\mathsf{y}^{2}+\mathsf{z}^{2}-3\left(1+\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2}\right)\left(\mathsf{x}\cos\mathfrak{w}_{\rm
R}\tau-\mathsf{y}\sin\mathfrak{w}_{\rm
R}\tau\right)^{2}+3\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2}\mathsf{z}^{2}\right].$

Since the momentum is defined by

$\displaystyle\mathsf{p}_{\mathsf{x}}$ $\displaystyle\equiv$
$\displaystyle{\partial{\cal
L}\over\partial\dot{\mathsf{x}}}=\mu\dot{\mathsf{x}}+\mu\mathfrak{w}_{\rm
dS}\mathsf{y},$ $\displaystyle\mathsf{p}_{\mathsf{y}}$ $\displaystyle\equiv$
$\displaystyle{\partial{\cal
L}\over\partial\dot{\mathsf{y}}}=\mu\dot{\mathsf{y}}-\mu\mathfrak{w}_{\rm
dS}\mathsf{x},$ $\displaystyle\mathsf{p}_{\mathsf{z}}$ $\displaystyle\equiv$
$\displaystyle{\partial{\cal
L}\over\partial\dot{\mathsf{z}}}=\mu\dot{\mathsf{z}}\,,$

we obtain the Hamiltonian as

$\displaystyle{\cal H}_{\rm rel}={\cal H}_{0}+{\cal H}_{1}\,,$ (22)

where

$\displaystyle{\cal H}_{0}$ $\displaystyle=$ $\displaystyle{1\over
2\mu}{\\!\\!\mbox{ \boldmath$\mathsf{p}$}}^{2}-{Gm_{1}m_{2}\over\mathsf{r}},$
$\displaystyle{\cal H}_{1}$ $\displaystyle=$ $\displaystyle{\cal
H}_{1\mathchar 45\relax{\rm dS}}+{\cal H}_{1\mathchar 45\relax\bar{\cal R}},$

with

$\displaystyle{\cal H}_{1\mathchar 45\relax{\rm dS}}$ $\displaystyle=$
$\displaystyle\mathfrak{w}_{\rm
dS}\left(\mathsf{p}_{\mathsf{y}}\mathsf{x}-\mathsf{p}_{\mathsf{x}}\mathsf{y}\right),$
$\displaystyle{\cal H}_{1\mathchar 45\relax\bar{\cal R}}$ $\displaystyle=$
$\displaystyle{\mu\mathfrak{r}_{g}\over
4\mathfrak{r}_{0}^{3}}\left[\mathsf{x}^{2}+\mathsf{y}^{2}+\mathsf{z}^{2}-3\left(1+\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2}\right)\left(\mathsf{x}\cos\mathfrak{w}_{\rm
R}\tau-\mathsf{y}\sin\mathfrak{w}_{\rm
R}\tau\right)^{2}+3\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2}\mathsf{z}^{2}\right].$

The term ${\cal H}_{1\mathchar 45\relax{\rm dS}}$ gives the so-called de
Sitter precession as follows: Let us consider the model with ${\cal
H}^{\prime}={\cal H}_{0}+{\cal H}_{1\mathchar 45\relax{\rm dS}}$. The Hamilton
equations are given as

$\displaystyle\dot{\mathsf{x}}$ $\displaystyle=$ $\displaystyle{\partial{\cal
H}^{\prime}\over\partial\mathsf{p}_{\mathsf{x}}}={\mathsf{p}_{\mathsf{x}}\over\mu}-\mathfrak{w}_{\rm
dS}\mathsf{y},$ $\displaystyle\dot{\mathsf{y}}$ $\displaystyle=$
$\displaystyle{\partial{\cal
H}^{\prime}\over\partial\mathsf{p}_{\mathsf{y}}}={\mathsf{p}_{\mathsf{y}}\over\mu}+\mathfrak{w}_{\rm
dS}\mathsf{x},$ $\displaystyle\dot{\mathsf{z}}$ $\displaystyle=$
$\displaystyle{\partial{\cal
H}^{\prime}\over\partial\mathsf{p}_{\mathsf{z}}}={\mathsf{p}_{z}\over\mu},$

and

$\displaystyle\dot{\mathsf{p}}_{\mathsf{x}}$ $\displaystyle=$
$\displaystyle-{\partial{\cal
H}^{\prime}\over\partial\mathsf{x}}=-{Gm_{1}m_{2}\mathsf{x}\over\mathsf{r}^{3}}-\mathfrak{w}_{\rm
dS}\mathsf{p}_{\mathsf{y}},$ $\displaystyle\dot{\mathsf{p}}_{\mathsf{y}}$
$\displaystyle=$ $\displaystyle-{\partial{\cal
H}^{\prime}\over\partial\mathsf{y}}=-{Gm_{1}m_{2}\mathsf{y}\over\mathsf{r}^{3}}+\mathfrak{w}_{\rm
dS}\mathsf{p}_{\mathsf{x}},$ $\displaystyle\dot{\mathsf{p}}_{\mathsf{z}}$
$\displaystyle=$ $\displaystyle-{\partial{\cal
H}^{\prime}\over\partial\mathsf{z}}=-{Gm_{1}m_{2}\mathsf{z}\over\mathsf{r}^{3}}.$

We then calculate the time evolution of the angular momentum $\\!\\!\mbox{
\boldmath$L$}=(L_{\mathsf{x}},L_{\mathsf{y}},L_{\mathsf{z}})$. Using the
equations of motion, we find

$\displaystyle\dot{L}_{\mathsf{x}}$ $\displaystyle=$ $\displaystyle{d\over
d\tau}(\mathsf{y}\mathsf{p}_{\mathsf{z}}-\mathsf{z}\mathsf{p}_{\mathsf{y}})=\dot{\mathsf{y}}\mathsf{p}_{\mathsf{z}}+\mathsf{y}\dot{\mathsf{p}}_{\mathsf{z}}-\dot{\mathsf{z}}\mathsf{p}_{\mathsf{y}}-\mathsf{z}\dot{\mathsf{p}}_{\mathsf{y}}=-\mathfrak{w}_{\rm
dS}L_{\mathsf{y}},$ $\displaystyle\dot{L}_{\mathsf{y}}$ $\displaystyle=$
$\displaystyle{d\over
d\tau}(\mathsf{z}\mathsf{p}_{\mathsf{x}}-\mathsf{x}\mathsf{p}_{\mathsf{z}})=\dot{\mathsf{z}}\mathsf{p}_{\mathsf{x}}+\mathsf{z}\dot{\mathsf{p}}_{\mathsf{x}}-\dot{\mathsf{x}}\mathsf{p}_{\mathsf{z}}-\mathsf{x}\dot{\mathsf{p}}_{\mathsf{z}}=\mathfrak{w}_{\rm
dS}L_{\mathsf{x}},$ $\displaystyle\dot{L}_{\mathsf{z}}$ $\displaystyle=$
$\displaystyle{d\over
d\tau}(\mathsf{x}\mathsf{p}_{\mathsf{y}}-\mathsf{y}\mathsf{p}_{\mathsf{x}})=\dot{\mathsf{x}}\mathsf{p}_{\mathsf{y}}+\mathsf{x}\dot{\mathsf{p}}_{\mathsf{y}}-\dot{\mathsf{y}}\mathsf{p}_{\mathsf{x}}-\mathsf{y}\dot{\mathsf{p}}_{\mathsf{x}}=0.$

From these equations, we find that the $\mathsf{z}$-component of the angular
momentum is conserved, and $(L_{\mathsf{x}},L_{\mathsf{y}})$ rotates around
the $\mathsf{z}$-axis with the constant angular frequency $\mathfrak{w}_{\rm
dS}$ which is just the de Sitter precession. The full equations are given as

$\displaystyle\dot{\mathsf{x}}$ $\displaystyle=$ $\displaystyle{\partial{\cal
H}\over\partial\mathsf{p}_{\mathsf{x}}}={\mathsf{p}_{\mathsf{x}}\over\mu}-\mathfrak{w}_{\rm
dS}\mathsf{y},$ (23) $\displaystyle\dot{\mathsf{y}}$ $\displaystyle=$
$\displaystyle{\partial{\cal
H}\over\partial\mathsf{p}_{\mathsf{y}}}={\mathsf{p}_{\mathsf{y}}\over\mu}+\mathfrak{w}_{\rm
dS}\mathsf{x},$ (24) $\displaystyle\dot{\mathsf{z}}$ $\displaystyle=$
$\displaystyle{\partial{\cal
H}\over\partial\mathsf{p}_{\mathsf{z}}}={\mathsf{p}_{\mathsf{z}}\over\mu},$
(25)

and

$\displaystyle\dot{\mathsf{p}}_{\mathsf{x}}$ $\displaystyle=$
$\displaystyle-{\partial{\cal
H}\over\partial\mathsf{x}}=-{Gm_{1}m_{2}\over\mathsf{r}^{3}}\mathsf{x}-\mathfrak{w}_{\rm
dS}\mathsf{p}_{\mathsf{y}}-{\mu\mathfrak{r}_{g}\over
2\mathfrak{r}_{0}^{3}}\left\\{\mathsf{x}-3(1+\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2})(\mathsf{x}\cos\mathfrak{w}_{\rm
R}\tau-\mathsf{y}\sin\mathfrak{w}_{\rm R}\tau)\cos\mathfrak{w}_{\rm
R}\tau\right\\},$ (26) $\displaystyle\dot{\mathsf{p}}_{\mathsf{y}}$
$\displaystyle=$ $\displaystyle-{\partial{\cal
H}\over\partial\mathsf{y}}=-{Gm_{1}m_{2}\over\mathsf{r}^{3}}\mathsf{y}+\mathfrak{w}_{\rm
dS}\mathsf{p}_{\mathsf{x}}-{\mu\mathfrak{r}_{g}\over
2\mathfrak{r}_{0}^{3}}\left\\{\mathsf{y}+3(1+\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2})(\mathsf{x}\cos\mathfrak{w}_{\rm
R}\tau-\mathsf{y}\sin\mathfrak{w}_{\rm R}\tau)\sin\mathfrak{w}_{\rm
R}\tau\right\\},$ (27) $\displaystyle\dot{\mathsf{p}}_{\mathsf{z}}$
$\displaystyle=$ $\displaystyle-{\partial{\cal
H}\over\partial\mathsf{z}}=-{Gm_{1}m_{2}\over\mathsf{r}^{3}}\mathsf{z}-{\mu\mathfrak{r}_{g}\over
2\mathfrak{r}_{0}^{3}}(1+3\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2})\mathsf{z}.$
(28)

###  IV.3 Motion with 0.5 PN correction term

Now we consider 0.5 PN terms. As discussed in §. III, we can assume
$\\!\\!\mbox{ \boldmath$R$}=0$ by introduction of the acceleration given by
Eq. (III.2). We first solve the relative coordinates  $r$, and then the motion
of the observer (or the CM).

####  IV.3.1 Equations of motion for relative coordinates

The equation of motion for relative coordinates  $r$ of a binary is now given
by

$\displaystyle\tilde{\cal L}_{\rm rel}(\\!\\!\mbox{
\boldmath$r$},\dot{\\!\\!\mbox{ \boldmath$r$}})$ $\displaystyle=$
$\displaystyle{\cal L}_{\rm rel}(\\!\\!\mbox{ \boldmath$r$},\dot{\\!\\!\mbox{
\boldmath$r$}})+{\cal L}_{1/2\mathchar 45\relax{\rm rel}}(\\!\\!\mbox{
\boldmath$r$},\dot{\\!\\!\mbox{ \boldmath$r$}})\,,$

where ${\cal L}_{\rm rel}$ is given by Eq. (20), while

$\displaystyle{\cal L}_{1/2\mathchar 45\relax{\rm rel}}(\\!\\!\mbox{
\boldmath$r$},\dot{\\!\\!\mbox{ \boldmath$r$}})$ $\displaystyle=$
$\displaystyle-\mu{2(m_{1}-m_{2})\over 3(m_{1}+m_{2})}\left(\bar{\cal
R}_{\hat{0}\hat{x}\hat{y}\hat{x}}x(x\dot{y}-y\dot{x})+\bar{\cal
R}_{\hat{0}\hat{z}\hat{y}\hat{z}}z(z\dot{y}-y\dot{z})\right)$ $\displaystyle=$
$\displaystyle-\mu{\mathfrak{r}_{g}\over\mathfrak{r}_{0}^{3}}{m_{1}-m_{2}\over
m_{1}+m_{2}}\mathfrak{r}_{0}\mathfrak{w}_{0}\sqrt{1+\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2}}\left(-x(x\dot{y}-y\dot{x})+z(z\dot{y}-y\dot{z})\right)\,.$

In non-rotating Fermi-Walker coordinates, we find ${\cal L}_{\rm rel}$ is
given by Eq. (21), while

$\displaystyle{\cal L}_{1/2\mathchar 45\relax{\rm rel}}(\\!\\!\mbox{
\boldmath$\mathsf{r}$},\dot{\\!\\!\mbox{ \boldmath$\mathsf{r}$}})$
$\displaystyle=$
$\displaystyle\mu{\mathfrak{r}_{g}\over\mathfrak{r}_{0}^{3}}{m_{1}-m_{2}\over
m_{1}+m_{2}}\mathfrak{r}_{0}\mathfrak{w}_{0}\sqrt{1+\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2}}\Big{\\{}\cos\mathfrak{w}_{\rm
R}\tau\left[\mathsf{x}\left(\mathsf{x}\dot{\mathsf{y}}-\mathsf{y}\dot{\mathsf{x}}\right)+\mathsf{z}\left(\mathsf{y}\dot{\mathsf{z}}-\mathsf{z}\dot{\mathsf{y}})+\mathfrak{w}_{\rm
R}\mathsf{x}(\mathsf{x}^{2}+\mathsf{y}^{2}-\mathsf{z}^{2}\right)\right]$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}-\sin\mathfrak{w}_{\rm
R}\tau\left[\mathsf{y}\left(\mathsf{x}\dot{\mathsf{y}}-\mathsf{y}\dot{\mathsf{x}}\right)+\mathsf{z}\left(\mathsf{z}\dot{\mathsf{x}}-\mathsf{x}\dot{\mathsf{z}})+\mathfrak{w}_{\rm
R}\mathsf{y}(\mathsf{x}^{2}+\mathsf{y}^{2}-\mathsf{z}^{2}\right)\right]\Big{\\}}.$

The momentum is obtained from the Lagrangian $\tilde{\cal L}_{\rm
rel}(\\!\\!\mbox{ \boldmath$\mathsf{r}$},\dot{\\!\\!\mbox{
\boldmath$\mathsf{r}$}})$ as

$\displaystyle p_{\mathsf{x}}$ $\displaystyle=$
$\displaystyle\mu\dot{\mathsf{x}}+\mu\mathfrak{w}_{\rm
dS}\mathsf{y}+\mu{\mathfrak{r}_{g}\over\mathfrak{r}_{0}^{3}}{m_{1}-m_{2}\over
m_{1}+m_{2}}\mathfrak{r}_{0}\mathfrak{w}_{0}\sqrt{1+\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2}}\left(-\mathsf{x}\mathsf{y}\cos\mathfrak{w}_{\rm
R}\tau+(\mathsf{y}^{2}-\mathsf{z}^{2})\sin\mathfrak{w}_{\rm R}\tau\right),$
$\displaystyle p_{\mathsf{y}}$ $\displaystyle=$
$\displaystyle\mu\dot{\mathsf{y}}-\mu\mathfrak{w}_{\rm
dS}\mathsf{x}+\mu{\mathfrak{r}_{g}\over\mathfrak{r}_{0}^{3}}{m_{1}-m_{2}\over
m_{1}+m_{2}}\mathfrak{r}_{0}\mathfrak{w}_{0}\sqrt{1+\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2}}\left(-\mathsf{x}\mathsf{y}\sin\mathfrak{w}_{\rm
R}\tau+(\mathsf{x}^{2}-\mathsf{z}^{2})\cos\mathfrak{w}_{\rm R}\tau\right),$
$\displaystyle p_{\mathsf{z}}$ $\displaystyle=$
$\displaystyle\mu\dot{\mathsf{z}}+\mu{\mathfrak{r}_{g}\over\mathfrak{r}_{0}^{3}}{m_{1}-m_{2}\over
m_{1}+m_{2}}\mathfrak{r}_{0}\mathfrak{w}_{0}\sqrt{1+\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2}}\mathsf{z}\left(\mathsf{y}\cos\mathfrak{w}_{\rm
R}\tau+\mathsf{x}\sin\mathfrak{w}_{\rm R}\tau\right).$

The Hamiltonian is given by

$\displaystyle\tilde{\cal H}_{\rm rel}(\\!\\!\mbox{
\boldmath$\mathsf{r}$},\\!\\!\mbox{ \boldmath$\mathsf{p}$})$ $\displaystyle=$
$\displaystyle{\cal H}_{\rm rel}(\\!\\!\mbox{
\boldmath$\mathsf{r}$},\\!\\!\mbox{ \boldmath$\mathsf{p}$})+{\cal
H}_{1/2\mathchar 45\relax{\rm rel}}(\\!\\!\mbox{
\boldmath$\mathsf{r}$},\\!\\!\mbox{ \boldmath$\mathsf{p}$})\,,$

where ${\cal H}_{\rm rel}(\\!\\!\mbox{ \boldmath$\mathsf{r}$},\\!\\!\mbox{
\boldmath$\mathsf{p}$})$ is given by Eq. (22), while

$\displaystyle{\cal H}_{1/2\mathchar 45\relax{\rm rel}}(\\!\\!\mbox{
\boldmath$\mathsf{r}$},\\!\\!\mbox{ \boldmath$\mathsf{p}$})$ $\displaystyle=$
$\displaystyle-{1\over
2}\mu\left({\mathfrak{r}_{g}\over\mathfrak{r}_{0}^{3}}\right)^{2}{(m_{1}-m_{2})^{2}\over(m_{1}+m_{2})^{2}}\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2}(1+\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2})\Big{[}\left(-\mathsf{x}\mathsf{y}\cos\mathfrak{w}_{\rm
R}\tau+(\mathsf{y}^{2}-\mathsf{z}^{2})\sin\mathfrak{w}_{\rm R}\tau\right)^{2}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\left(-\mathsf{x}\mathsf{y}\sin\mathfrak{w}_{\rm
R}\tau+(\mathsf{x}^{2}-\mathsf{z}^{2})\cos\mathfrak{w}_{\rm
R}\tau\right)^{2}+\mathsf{z}^{2}\left(\mathsf{y}\cos\mathfrak{w}_{\rm
R}\tau+\mathsf{x}\sin\mathfrak{w}_{\rm R}\tau\right)^{2}\Big{]}\,.$

This Hamiltonian is very complicated, but it should be ignored because it is
beyond quadrupole approximation, although momentum is modified. For an equal
mass binary ($m_{1}=m_{2}$), the 0.5PN correction term vanishes and the
momentum is also the same as the Newtonian one. As a result, the Newtonian
solution is also a solution.

####  IV.3.2 Motion of the CM of a binary and its stability

In order to study stability of the CM of a binary system, we analyze Eq. (14).
Since  $R$ is measured by the circular observer at
$\mathfrak{r}=\mathfrak{r}_{0}$, we can split the 4-velocity $u^{\mu}$ as

$\displaystyle u^{\mu}$ $\displaystyle=$ $\displaystyle
u^{\mu}_{(0)}+u^{\mu}_{(1)},$

where

$\displaystyle u^{\mu}_{(0)}$ $\displaystyle=$
$\displaystyle{d\mathfrak{x}^{\mu}_{(0)}\over
d\tau}=(u^{0}_{(0)}\,,0\,,0\,,u^{3}_{(0)})=\left(\sqrt{1+\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2}\over
f(\mathfrak{r}_{0})}\,,0\,,0\,,\mathfrak{w}_{0}\right),$ $\displaystyle
u^{\mu}_{(1)}$ $\displaystyle=$ $\displaystyle{d\mathfrak{x}^{\mu}_{(1)}\over
d\tau},$

with

$\displaystyle\mathfrak{x}^{\mu}_{(0)}$ $\displaystyle=$
$\displaystyle\left(\sqrt{1+\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2}\over
f(\mathfrak{r}_{0})}\tau\,,\mathfrak{r}_{0}\,,{\pi\over
2}\,,\mathfrak{w}_{0}\tau\right),$ $\displaystyle\mathfrak{x}^{\mu}_{(1)}$
$\displaystyle\equiv$ $\displaystyle e^{\mu}_{~{}\hat{\ell}}R^{\hat{\ell}}.$

The acceleration $a^{\mu}$ is given by the motion of a binary
$x^{\hat{\mu}}(\tau)$ in a rotating frame as

$\displaystyle a^{\mu}=-{3\mu\over
m_{1}+m_{2}}\mathfrak{w}_{0}\sqrt{1+\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2}}{\mathfrak{r}_{g}\over\mathfrak{r}_{0}^{2}}\left[\left(\delta^{\mu}_{0}{\mathfrak{r}_{0}\mathfrak{w}_{0}\over\sqrt{f(\mathfrak{r}_{0})}}+\delta^{\mu}_{3}{\sqrt{1+\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2}}\over\mathfrak{r}_{0}}\right)(\dot{x}x-\dot{z}z)\,-\delta^{\mu}_{1}\sqrt{f(\mathfrak{r}_{0})}\dot{y}x\,-\delta^{\mu}_{2}{1\over\mathfrak{r}_{0}}\dot{y}z\right].$

Here we assume that the deviation from a circular orbit is small, i.e.,
$\mathfrak{x}^{\mu}_{(1)}$ and $u^{\mu}_{(1)}$ are small perturbations.
Ignoring non-linear deviation terms in the equations of motion
$\displaystyle{{Du^{\mu}\over d\tau}=a^{\mu}}$, because the circular orbit
$\mathfrak{x}^{\mu}_{(0)}(\tau)$ is a geodesic, we obtain a linear
differential equation as

$\displaystyle{du^{\mu}_{(1)}\over
d\tau}+2\Gamma^{\mu}_{~{}\rho\sigma}(\mathfrak{r}_{0})u^{\rho}_{(0)}u^{\sigma}_{(1)}+{\partial\Gamma^{\mu}_{~{}\rho\sigma}\over\partial\mathfrak{x}^{\alpha}}(\mathfrak{r}_{0})\,\mathfrak{x}_{(1)}^{\alpha}\,u^{\rho}_{(0)}u^{\sigma}_{(0)}=a^{\mu},$

where $a^{\mu}$ acts as an external force. Describing the deviation as

$\displaystyle\mathfrak{x}^{\mu}_{(1)}=(\mathfrak{t}_{(1)},\mathfrak{r}_{(1)},\mathfrak{\theta}_{(1)},\mathfrak{\varphi}_{(1)})\,,$

we find

$\displaystyle{d^{2}\mathfrak{t}_{(1)}\over
d\tau^{2}}+{\mathfrak{r}_{g}\over\mathfrak{r}_{0}^{2}f(\mathfrak{r}_{0})}\sqrt{1+\mathfrak{r}_{0}^{2}\mathfrak{w}^{2}\over
f(\mathfrak{r}_{0})}{d\mathfrak{r}_{(1)}\over d\tau}=a^{0}=-{3\mu\over
m_{1}+m_{2}}{\mathfrak{w}_{0}^{2}\mathfrak{r}_{g}\over\mathfrak{r}_{0}}{\sqrt{1+\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2}\over
f(\mathfrak{r}_{0})}}(\dot{x}x-\dot{z}z),$ (29)
$\displaystyle{d^{2}\mathfrak{r}_{(1)}\over d\tau^{2}}-{3\mathfrak{r}_{g}\over
2\mathfrak{r}_{0}^{3}}\left(1+\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2}\right)\mathfrak{r}_{(1)}+{\mathfrak{r}_{g}\over\mathfrak{r}_{0}^{2}}f(\mathfrak{r}_{0})\sqrt{1+\mathfrak{r}_{0}^{2}\mathfrak{w}^{2}\over
f(\mathfrak{r}_{0})}{d\mathfrak{t}_{(1)}\over
d\tau}-2\mathfrak{r}_{0}\mathfrak{w}_{0}f(\mathfrak{r}_{0}){d\mathfrak{\varphi}_{(1)}\over
d\tau}$ $\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}=a^{1}={3\mu\over
m_{1}+m_{2}}{\mathfrak{w}_{0}\mathfrak{r}_{g}\over\mathfrak{r}_{0}^{2}}\sqrt{f(\mathfrak{r}_{0})(1+\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2})}\dot{y}x,~{}~{}~{}~{}~{}~{}$
(30) $\displaystyle{d^{2}\mathfrak{\theta}_{(1)}\over
d\tau^{2}}+\mathfrak{w}_{0}^{2}\mathfrak{\theta}_{(1)}=a^{2}={3\mu\over
m_{1}+m_{2}}{\mathfrak{w}_{0}\mathfrak{r}_{g}\over\mathfrak{r}_{0}^{3}}\sqrt{1+\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2}}\,\dot{y}z,$
(31) $\displaystyle{d^{2}\mathfrak{\varphi}_{(1)}\over
d\tau^{2}}+2{\mathfrak{w}_{0}\over\mathfrak{r}_{0}}{d\mathfrak{r}_{(1)}\over
d\tau}=a^{3}=-{3\mu\over
m_{1}+m_{2}}{\mathfrak{w}_{0}\mathfrak{r}_{g}\over\mathfrak{r}_{0}^{3}}\left(1+\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2}\right)(\dot{x}x-\dot{z}z),$
(32)

Integrating Eqs. (29) and (32), we obtain

$\displaystyle{d\mathfrak{t}_{(1)}\over
d\tau}=-{\mathfrak{r}_{g}\over\mathfrak{r}_{0}^{2}f(\mathfrak{r}_{0})}\sqrt{1+\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2}\over
f(\mathfrak{r}_{0})}\,\mathfrak{r}_{(1)}-{3\mu\over
2(m_{1}+m_{2})}{\mathfrak{w}_{0}^{2}\mathfrak{r}_{g}\over\mathfrak{r}_{0}}{\sqrt{1+\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2}\over
f(\mathfrak{r}_{0})}}\left(x^{2}-z^{2}\right),$ (33)
$\displaystyle{d\mathfrak{\varphi}_{(1)}\over
d\tau}=-2{\mathfrak{w}_{0}\over\mathfrak{r}_{0}}\,\mathfrak{r}_{(1)}-{3\mu\over
2(m_{1}+m_{2})}{\mathfrak{w}_{0}\mathfrak{r}_{g}\over\mathfrak{r}_{0}^{3}}\left(1+\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2}\right)\left(x^{2}-z^{2}\right)\,.$
(34)

where we set the integration constants are zero. Plugging Eqs. (33) and (34)
into Eq. (30) with Eq. (17), we obtain the perturbation equation for the
radial coordinates $\mathfrak{r}_{(1)}$ as

$\displaystyle{d^{2}\mathfrak{r}_{(1)}\over
d\tau^{2}}+k^{2}\mathfrak{r}_{(1)}+A\left(x^{2}-z^{2}\right)+B\dot{y}x=0\,,$
(35)

where

$\displaystyle k^{2}$ $\displaystyle=$
$\displaystyle{\mathfrak{r}_{g}(\mathfrak{r}_{0}-3\mathfrak{r}_{g})\over\mathfrak{r}_{0}^{3}(2\mathfrak{r}_{0}-3\mathfrak{r}_{g})},$
(36) $\displaystyle A$ $\displaystyle=$ $\displaystyle{3\mu\over
m_{1}+m_{2}}{\mathfrak{r}_{g}^{2}(\mathfrak{r}_{0}-\mathfrak{r}_{g})\over\mathfrak{r}_{0}^{5}(2\mathfrak{r}_{0}-3\mathfrak{r}_{g})}$
(37) $\displaystyle B$ $\displaystyle=$ $\displaystyle-{3\mu\over
m_{1}+m_{2}}\sqrt{2\mathfrak{r}_{g}^{3}\over\mathfrak{r}_{0}^{7}}{(\mathfrak{r}_{0}-\mathfrak{r}_{g})\over(2\mathfrak{r}_{0}-3\mathfrak{r}_{g})}.$
(38)

We find that $k^{2}>0$ when $\mathfrak{r}_{0}>3\mathfrak{r}_{g}$, while $A>0$
and $B<0$ when $\mathfrak{r}_{0}>3\mathfrak{r}_{g}/2$. The condition for
$k^{2}>0$ is consistent with the fact that the radius of innermost stable
circular orbit (ISCO) is $3\mathfrak{r}_{g}$.

If the binary orbit is bounded ($x(\tau),y(\tau)$ and $z(\tau)$ are finite),
the orbit of the center of mass is also bounded because $k^{2}>0$. We expect
that when $\mathfrak{r}_{0}>3\mathfrak{r}_{g}$ (ISCO radius), a binary system
near SMBH is linearly stable unless a binary is broken.

##  V Numerical Analysis

###  V.1 Validity and Stability

Before showing our numerical results, we discuss validity of the present
approach. The minimum curvature radius at the radius $\mathfrak{r}_{0}$ is
evaluated as

$\displaystyle\ell_{\bar{\cal R}}$ $\displaystyle=$ $\displaystyle{\rm
min}\left[|\bar{\cal R}_{\hat{\mu}\hat{\nu}\hat{\rho}\hat{\sigma}}|^{-{1\over
2}},|\bar{\cal
R}_{\hat{\mu}\hat{\nu}\hat{\rho}\hat{\sigma};\hat{\alpha}}|^{-{1\over
3}},|\bar{\cal
R}_{\hat{\mu}\hat{\nu}\hat{\rho}\hat{\sigma};\hat{\alpha};\hat{\beta}}|^{-{1\over
4}}\right]$ $\displaystyle\sim$ $\displaystyle{\rm
min}\left[\left({\mathfrak{r}_{g}\over\mathfrak{r}_{0}^{3}}\right)^{-{1\over
2}}\,,\left({\mathfrak{r}_{g}\over\mathfrak{r}_{0}^{4}}\right)^{-{1\over
3}}\,,\left({\mathfrak{r}_{g}\over\mathfrak{r}_{0}^{5}}\right)^{-{1\over
4}}\right]$ $\displaystyle\sim$
$\displaystyle\mathfrak{r}_{0}\left(\mathfrak{r}_{0}\over{\mathfrak{r}_{g}}\right)^{1/4}\geq
3\sqrt[4]{3}\mathfrak{r}_{g}\sim 8{\rm AU}\left({M\over
10^{8}M_{\odot}}\right)\,.$

The equality is held at the ISCO radius $\mathfrak{r}_{0}=3\mathfrak{r}_{g}$.

When we put a binary at $\mathfrak{r}=\mathfrak{r}_{0}$, the binary size
$\ell_{\rm binary}$ should satisfy

$\displaystyle\ell_{\rm binary}\ll\ell_{\bar{\cal R}}$

As for stability of a binary, the mutual gravitational interaction between a
binary should be much larger than the tidal force by a third body. The
condition is given by

$\displaystyle{Gm_{1}m_{2}\over
r^{2}}\gg{\mu\mathfrak{r}_{g}\over\mathfrak{r}_{0}^{3}}\,r\,,$

which gives the constraint on a binary size $\ell_{\rm binary}$ as

$\displaystyle\ell_{\rm binary}$ $\displaystyle\ll$
$\displaystyle\left({m_{1}+m_{2}\over 2M}\right)^{1\over 3}\,\mathfrak{r}_{0}$
$\displaystyle\hskip-42.67912pt\approx 4.64\times
10^{-3}\left({(m_{1}+m_{2})\over 20M_{\odot}}\right)^{1\over 3}\left({M\over
10^{8}M_{\odot}}\right)^{-{1\over 3}}\,\mathfrak{r}_{0}.$

When we are interested in the orbit near the ISCO radius
$\mathfrak{r}_{0}=3\mathfrak{r}_{g}$, we find

$\displaystyle\ell_{\rm binary}\ll 3\times 10^{-2}{\rm
AU}\left({(m_{1}+m_{2})\over 20M_{\odot}}\right)^{1\over 3}\left({M\over
10^{8}M_{\odot}}\right)^{{2\over 3}}.$

As for the relativistic effect in a binary becomes important when

$\displaystyle\ell_{\rm binary}$ $\displaystyle\leq$
$\displaystyle{G(m_{1}+m_{2})\over c^{2}}$ $\displaystyle\approx$
$\displaystyle 2\times 10^{-7}{\rm AU}\left({(m_{1}+m_{2})\over
20M_{\odot}}\right)\,.$

Hence, for a binary with the size of

$\displaystyle 2\times 10^{-7}{\rm AU}\left({(m_{1}+m_{2})\over
20M_{\odot}}\right)\ll\ell_{\rm binary}$ $\displaystyle~{}~{}~{}~{}\ll 3\times
10^{-2}{\rm AU}\left({(m_{1}+m_{2})\over 20M_{\odot}}\right)^{1\over
3}\left({M\over 10^{8}M_{\odot}}\right)^{{2\over 3}}\,,$

we may apply the present Newtonian approach to the ISCO radius.

###  V.2 Normalization and Initial Data

We show some numerical examples for an equal-mass binary ($m_{1}=m_{2}$).
Hence, we have to solve Eqs. (23) -(28). In order to solve these basic
equations, we shall introduce dimensionless variables as follows: The length
scale of a binary is normalized by an initial semi-major axis $a_{0}$, while
time scale is normalized by an initial binary mean motion $n_{0}$, which is
defined by

$\displaystyle n_{0}\equiv\bigg{(}{G(m_{1}+m_{2})\over
a_{0}^{3}}\bigg{)}^{1/2}\,.$

Note that the initial binary period is given by $P_{\rm in}=2\pi/n_{0}$.
Introducing

$\displaystyle\tilde{\tau}=n_{0}\tau,$
$\displaystyle\tilde{\mathsf{x}}={\mathsf{x}\over
a_{0}}\,,~{}\tilde{\mathsf{y}}={\mathsf{y}\over
a_{0}}\,,~{}\tilde{\mathsf{z}}={\mathsf{z}\over
a_{0}}\,,~{}\tilde{\mathsf{r}}={\mathsf{r}\over a_{0}},$
$\displaystyle\tilde{p}_{\mathsf{x}}={p_{\mathsf{x}}\over\mu
a_{0}n_{0}}\,,~{}\tilde{p}_{\mathsf{y}}={p_{\mathsf{y}}\over\mu
a_{0}n_{0}}\,,~{}\tilde{p}_{\mathsf{z}}={p_{\mathsf{z}}\over\mu
a_{0}n_{0}}\,,$

we find the dimensionless equations of motion as

$\displaystyle{d\tilde{\mathsf{x}}\over d\tilde{\tau}}$ $\displaystyle=$
$\displaystyle\tilde{p}_{\mathsf{x}}-\tilde{\mathfrak{w}}_{\rm
dS}\tilde{\mathsf{y}},$ (40) $\displaystyle{d\tilde{\mathsf{y}}\over
d\tilde{\tau}}$ $\displaystyle=$
$\displaystyle\tilde{p}_{\mathsf{y}}+\tilde{\mathfrak{w}}_{\rm
dS}\tilde{\mathsf{x}},$ (41) $\displaystyle{d\tilde{\mathsf{z}}\over
d\tilde{\tau}}$ $\displaystyle=$ $\displaystyle\tilde{p}_{\mathsf{z}},$ (42)

and

$\displaystyle{d\tilde{\mathsf{p}}_{\mathsf{x}}\over d\tilde{\tau}}$
$\displaystyle=$
$\displaystyle-{\tilde{\mathsf{x}}\over\tilde{\mathsf{r}}^{3}}-\tilde{\mathfrak{w}}_{\rm
dS}\tilde{\mathsf{p}}_{\mathsf{y}}-{\mathfrak{r}_{g}\over
2\mathfrak{r}_{0}}\epsilon^{2}\left\\{\tilde{\mathsf{x}}-3\left(1+\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2}\right)(\tilde{\mathsf{x}}\cos\tilde{\mathfrak{w}}_{\rm
R}\tilde{\tau}-\tilde{\mathsf{y}}\sin\tilde{\mathfrak{w}}_{\rm
R}\tilde{\tau})\cos\tilde{\mathfrak{w}}_{\rm R}\tilde{\tau}\right\\},$ (43)
$\displaystyle{d\tilde{\mathsf{p}}_{\mathsf{y}}\over d\tilde{\tau}}$
$\displaystyle=$
$\displaystyle-{\tilde{\mathsf{y}}\over\tilde{\mathsf{r}}^{3}}+\tilde{\mathfrak{w}}_{\rm
dS}\tilde{\mathsf{p}}_{\mathsf{x}}-{\mathfrak{r}_{g}\over
2\mathfrak{r}_{0}}\epsilon^{2}\left\\{\tilde{\mathsf{y}}+3\left(1+\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2}\right)(\tilde{\mathsf{x}}\cos\tilde{\mathfrak{w}}_{\rm
R}\tilde{\tau}-\tilde{\mathsf{y}}\sin\tilde{\mathfrak{w}}_{\rm
R}\tilde{\tau})\sin\tilde{\mathfrak{w}}_{\rm R}\tilde{\tau}\right\\},$ (44)
$\displaystyle{d\tilde{\mathsf{p}}_{\mathsf{z}}\over d\tilde{\tau}}$
$\displaystyle=$
$\displaystyle-{\tilde{\mathsf{z}}\over\tilde{\mathsf{r}}^{3}}-{\mathfrak{r}_{g}\over
2\mathfrak{r}_{0}}\epsilon^{2}\left(1+3\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2}\right)\tilde{\mathsf{z}}\,,$
(45)

where

$\displaystyle\tilde{\mathfrak{w}}_{\rm dS}$ $\displaystyle\equiv$
$\displaystyle{\mathfrak{w}_{\rm dS}\over
n_{0}}\,,~{}~{}\tilde{\mathfrak{w}}_{\rm R}\,\equiv\,{\mathfrak{w}_{\rm
R}\over n_{0}},$ $\displaystyle\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2}$
$\displaystyle=$
$\displaystyle\left({2\mathfrak{r}_{0}\over\mathfrak{r}_{g}}-3\right)^{-1},$

and

$\displaystyle\epsilon^{2}\equiv{1\over\mathfrak{r}_{0}^{2}n_{0}^{2}}=\left({a_{0}\over\mathfrak{r}_{0}}\right)^{3}\left({\mathfrak{r}_{0}\over\mathfrak{r}_{g}}\right)\left({2M\over
m_{1}+m_{2}}\right).$ (46)

This $\epsilon$ corresponds to the initial semi-major axis $a_{0}$ as

$\displaystyle a_{0}=\epsilon^{2\over
3}\,\mathfrak{r}_{0}\,\left({m_{1}+m_{2}\over
2M}\right)^{1/3}\left({\mathfrak{r}_{0}\over\mathfrak{r}_{g}}\right)^{-1/3}\,.$

Using $\epsilon$, we find

$\displaystyle\tilde{\mathfrak{w}}_{\rm dS}$ $\displaystyle=$
$\displaystyle\epsilon\,\left[\left({\mathfrak{r}_{g}\over
2\mathfrak{r}_{0}-3\mathfrak{r}_{g}}\right)^{1\over
2}-\left({\mathfrak{r}_{g}\over 2\mathfrak{r}_{0}}\right)^{1\over 2}\right],$
$\displaystyle\tilde{\mathfrak{w}}_{\rm R}$ $\displaystyle=$
$\displaystyle\epsilon\,\left({\mathfrak{r}_{g}\over
2\mathfrak{r}_{0}}\right)^{1\over 2}.$

As we discuss in §. V.1, the necessary conditions for a stable “Newtonian”
binary is given by the length scale of a binary. If we set $\ell_{\rm
binary}\sim a_{0}$, we find the condition for $\epsilon$ as

$\displaystyle{m_{1}+m_{2}\over
2M}\left({\mathfrak{r}_{0}\over\mathfrak{r}_{g}}\right)^{-1}\ll\epsilon\ll\sqrt{\mathfrak{r}_{0}\over\mathfrak{r}_{g}}.$

When a binary locates near the ISCO radius, this condition is

$\displaystyle O(10^{-7})\ll\epsilon\ll O(1),$

for $m_{1}=m_{2}=10M_{\odot},M=10^{8}M_{\odot}$.

In order to solve Eqs. (23) -(28), we first give masses $m_{1},m_{2},M$, the
radius of the circular orbit $\mathfrak{r}_{0}$, and $\epsilon$, which
corresponds to the initial semi-major axis of a binary $a_{0}$. We then
provide the initial data of a binary, i.e.,
$\tilde{\mathsf{x}}(0)\,,\tilde{\mathsf{y}}(0)\,,\tilde{\mathsf{z}}(0)$ and
$\tilde{\mathsf{p}}_{\mathsf{x}}(0)\,,\tilde{\mathsf{p}}_{\mathsf{y}}(0)\,,\tilde{\mathsf{p}}_{\mathsf{z}}(0)$.

Since the motion of a binary can be approximated by an elliptic orbit, we
shall fix the initial values by assuming an elliptic orbit given by

$\displaystyle\mathsf{r}={a(1-e^{2})\over 1+e\cos f}\,,$

where $a$ is a semi-major axis, $e$ is the eccentricity, and $f$ is true
anomaly. Since the orbital plane is not, in general, $\mathsf{z}=0$, we have
to introduce three angular variables; the argument of periapsis $\omega$, the
ascending node $\Omega$ and the inclination angle $I$. We have the relations
between the position $\bf{\mathsf{r}}=(\mathsf{x},\mathsf{y},\mathsf{z})$ of
the component of a binary and the orbital parameters
$(\omega\,,\Omega\,,a\,,e\,,I\,,f)$ as

$\displaystyle\begin{pmatrix}\mathsf{x}\\\ \mathsf{y}\\\ \mathsf{z}\\\
\end{pmatrix}$ $\displaystyle=$
$\displaystyle\mathsf{r}\begin{pmatrix}\cos\Omega\cos(\omega+f)-\sin\Omega\sin(\omega+f)\cos
I\\\ \sin\Omega\cos(\omega+f)+\cos\Omega\sin(\omega+f)\cos I\\\
\sin(\omega+f)\sin I\\\ \end{pmatrix}$ (47)

the initial position of a binary can be fixed by the orbital parameters. As
for the initial velocity, we have the relation between the mean anomaly
$\mathfrak{l}$ and the true anomaly $f$ as

$\displaystyle d\mathfrak{l}=(1-e^{2})^{{3\over 2}}\left(1+e\cos
f\right)^{-2}df.$ (48)

In Newtonian dynamics, $\mathfrak{l}=n(\tau-\tau_{0})$, where $n$ is the mean
motion. Hence, in this approximation, the time derivative ($d/d\tau$) is given
by the derivative with respect to the true anomaly ($d/df$).

Assuming $f=0$ at $\tau=0$, we find

$\displaystyle\tilde{\mathsf{x}}(0)$ $\displaystyle=$
$\displaystyle(1-e_{0})\cos\Omega_{0}\cos\omega_{0}-\sin\Omega_{0}\sin\omega_{0}\cos
I_{0},$ $\displaystyle\tilde{\mathsf{y}}(0)$ $\displaystyle=$
$\displaystyle(1-e_{0})\sin\Omega_{0}\cos\omega_{0}+\cos\Omega_{0}\sin\omega_{0}\cos
I_{0},$ $\displaystyle\tilde{\mathsf{z}}(0)$ $\displaystyle=$
$\displaystyle(1-e_{0})\sin\omega_{0}\sin I_{0},$

and

$\displaystyle{d\tilde{\mathsf{x}}\over d\tilde{\tau}}(0)$ $\displaystyle=$
$\displaystyle-\sqrt{1+e_{0}\over
1-e_{0}}\,\left[\cos\Omega_{0}\sin\omega_{0}+\sin\Omega_{0}\cos\omega_{0}\cos
I_{0}\right],$ $\displaystyle{d\tilde{\mathsf{y}}\over d\tilde{\tau}}(0)$
$\displaystyle=$ $\displaystyle\sqrt{1+e_{0}\over
1-e_{0}}\,\left[-\sin\Omega_{0}\sin\omega_{0}+\cos\Omega_{0}\cos\omega_{0}\cos
I_{0}\right],$ $\displaystyle{d\tilde{\mathsf{z}}\over d\tilde{\tau}}(0)$
$\displaystyle=$ $\displaystyle\sqrt{1+e_{0}\over 1-e_{0}}\,\cos\omega_{0}\sin
I_{0}.$

Hence, when we prepare the initial orbital parameters
($e_{0}\,,I_{0}\,,\omega_{0}\,,\Omega_{0}$), we can provide the initial data
for the normalized evolution equations (23) -(28).

###  V.3 Binary motion near ISCO radius

In a hierarchical triple system, there are several important features. One is
the so-called Kozai-Lidov (KL) oscillations. If the system is inclined more
than some critical angle, there appears an oscillation between the
eccentricity and inclination angle. The second interesting feature is an
orbital flip, which may appear when the inclination angle evolves into near
$90^{\circ}$. The last one which we show is a chaotic feature in the long-time
evolution.

Here, we show our numerical results. In order to discuss the properties of a
binary orbit, it is more convenient to use the orbital parameters assuming
that the binary motion is close to an elliptic orbit.

In order to extract the orbital parameters from the orbit given by the
Cartesian coordinates, one can use the osculating orbit when the orbit is
close to an ellipse. However, one must be careful with the definitions of
orbital elements when using the osculating method. For instance, the magnitude
of the normalized Laplace-Runge-Lenz vector, which is defined by

$\displaystyle\\!\\!\mbox{ \boldmath$e$}\equiv\tilde{\\!\\!\mbox{
\boldmath$\mathsf{p}$}}\times(\tilde{\\!\\!\mbox{
\boldmath$\mathsf{r}$}}\times\tilde{\\!\\!\mbox{
\boldmath$\mathsf{p}$}})-{\tilde{\\!\\!\mbox{
\boldmath$\mathsf{r}$}}\over\tilde{\mathsf{r}}}\,,$ (49)

is commonly used for a measure of orbital eccentricity, but it is not always
appropriate. It may show an “apparent” rise in eccentricity or unphysical
rapid oscillations especially when the eccentricity is very small Will (2019).
We take caution and extract the elements’ information from physical orbit. In
that case, it may be better to define the eccentricity by the averaged one
over one cycle as

$\displaystyle\langle e\rangle\equiv{\mathsf{r}_{\rm max}-\mathsf{r}_{\rm
min}\over\mathsf{r}_{\rm max}+\mathsf{r}_{\rm min}}.$

where $\mathsf{r}_{\rm max}$ and $\mathsf{r}_{\rm min}$ correspond to orbital
separation at adjacent turning points of an eccentric orbit.

The inclination angle $I$ is defined as mutual inclination between angular
momenta of the inner and outer binary. In the present case, since the outer
binary is just a circular motion on the equatorial plane, the inclination is
given by

$\displaystyle I=\cos^{-1}\left({L_{\mathsf{z}}\over|\\!\\!\mbox{
\boldmath$L$}|}\right)\,,$

where $\\!\\!\mbox{ \boldmath$L$}\equiv\\!\\!\mbox{
\boldmath$\mathsf{r}$}\times\\!\\!\mbox{ \boldmath$\mathsf{p}$}$ is the
angular momentum of a binary.

The other two essential angles $\Omega$ and $\omega$ governing the orientation
of the orbital plane and the orbit are also computed in the post-process. The
line that marks the intersection of the orbital plane with the reference plane
(the equatorial plane in the present case) is called the node line, and the
point on the node line where the orbit passes above the reference plane from
below is called the ascending node. The angle between the reference axis (say
$\mathsf{x}$-axis) and node line vector  $N$ is the longitude of ascending
node $\Omega$. First, node line is defined as,

$\displaystyle\\!\\!\mbox{ \boldmath$N$}=\\!\\!\mbox{
\boldmath$e$}_{\mathsf{z}}\times\\!\\!\mbox{ \boldmath$h$}\,.$

where $\\!\\!\mbox{ \boldmath$e$}_{\mathsf{z}}$ is normal to the reference
plane (the unit vector in the $\mathsf{z}$ direction) and $\\!\\!\mbox{
\boldmath$h$}=\\!\\!\mbox{ \boldmath$L$}/\mu$ is specific angular momentum
vector of a binary. Thus, $\Omega$ is computed as,

$\displaystyle\Omega=\cos^{-1}(N_{\mathsf{x}}/N)\,.$

The argument of periapsis $\omega$ is the angle between node line and
periapsis measured in the direction of motion. Therefore,

$\displaystyle\omega=\cos^{-1}\bigg{(}\frac{\\!\\!\mbox{
\boldmath$N$}\cdot\\!\\!\mbox{ \boldmath$e$}}{N\,e}\bigg{)}\,.$

When the orbit can be approximated well by the osculating one,  $e$ is given
by the normalized Laplace-Runge-Lenz vector (49). Otherwise, we define the
averaged eccentricity vector by

$\displaystyle\langle\\!\\!\mbox{
\boldmath$e$}\rangle\equiv-{\left(\\!\\!\mbox{ \boldmath$\mathsf{r}$}_{\rm
min}+\\!\\!\mbox{ \boldmath$\mathsf{r}$}_{\rm
max}\right)\over\left(\mathsf{r}_{\rm min}+\mathsf{r}_{\rm max}\right)}$

pointing towards the periapsis, where $\\!\\!\mbox{
\boldmath$\mathsf{r}$}_{\rm max}$ and $\\!\\!\mbox{
\boldmath$\mathsf{r}$}_{\rm min}$ are numerical data of position vector. We
have used both definitions and found that most results agree well.

####  V.3.1 Firmness of a binary and stability

As we discuss in §. V.1, a binary near SMBH may be broken when the tidal force
by SMBH is stronger than the mutual gravitational attractive force of a
binary. We introduce a firmness parameter of a binary $\mathfrak{f}$ defined
by

$\displaystyle\mathfrak{f}$ $\displaystyle\equiv$ $\displaystyle{{\rm
gravitational~{}force}\over{\rm
tidal~{}force~{}by~{}SMBH}}={{Gm_{1}m_{2}/\ell_{\rm
binary}^{2}}\over{\mu\mathfrak{r}_{g}\,\ell_{\rm
binary}/\mathfrak{r}_{0}^{3}}}$ $\displaystyle=$
$\displaystyle\left({m_{1}+m_{2}\over
2M}\right)\left({\mathfrak{r}_{0}\over\ell_{\rm
binary}}\right)^{3}\approx{1\over\epsilon^{2}}\left({\mathfrak{r}_{0}\over\mathfrak{r}_{g}}\right)\,.$

If the firmness parameter $\mathfrak{f}$ is smaller than $O(1)$, we expect
that the tidal force will break a binary. Hence, $\mathfrak{f}>1$ is a
necessary condition for stability of a binary. When we are interested in a
binary motion near SMBH, this stability condition gives $\epsilon\leq O(1)$.
We can confirm numerically that our present system is really unstable when
$\epsilon\geq 0.5$. In what follows, we numerically analyze a binary system
under the conditions of $\epsilon\leq 0.4$.

####  V.3.2 Kozai-Lidov oscillation

Here, we show some numerical examples, which shows the KL oscillations in the
long-time evolution.

We expect the KL oscillations to occur when the inclination angle is larger
than the critical value. Under the quadrupole approximation in a Newtonian
hierarchical system, the critical inclination angle is given by $I_{\rm
crit}^{\rm(N)}=\sin^{-1}\sqrt{2/5}\approx 39.23^{\circ}$. In the present
model, it can be obtained by the double-averaging analysis of the Lagrange
planetary equations, which shows that the critical value is slightly larger as
the radius $\mathfrak{r}_{0}$ gets smaller, and it increases to $41.6^{\circ}$
near the ISCO radius (See Appendix B.2). It seems hard to obtain the exact
critical value by numerical simulations, although we have found consistent
results.

The typical time scale of the KL oscillations is given by Suzuki et al.
(2019); Antognini (2015); Shevchenko (Springer, 2017)

$\displaystyle T_{\rm KL}\sim{P_{\rm out}^{2}\over P_{\rm in}}\,,$ (50)

where $P_{\rm in}$ and $P_{\rm out}$ are the periods of an inner binary and of
an outer binary, respectively. For a hierarchical triple system, since $P_{\rm
in}\ll P_{\rm out}$, $T_{\rm KL}\gg P_{\rm out}$, which means that the KL
oscillation is a secular effect.

In the first model (Model Ic), we choose
$\mathfrak{r}_{0}=3.5\mathfrak{r}_{g}$ and $\epsilon=0.1$. It corresponds to
$a_{0}=0.0023\mathfrak{r}_{g}\approx 0.0045{\rm AU}$ for
$m_{1}=m_{2}=10M_{\odot}$ and $M=10^{8}M_{\odot}$. Since
$G(m_{1}+m_{2})/c^{2}\approx 30{\rm km}\ll a_{0}$, the binary motion can be
described by a Newtonian orbit in a local inertial frame. We then adopt the
initial conditions as the eccentricity $e_{0}=0.01$ and inclination
$I_{0}=60^{\circ}$. The numerical results are given in Figs. 2 and 3.

Figure 2: This figure shows Kozai-Lidov oscillation between orbital
eccentricity $e$ (the red curve) and relative inclination $I$ (the blue
curve). For this particular model (Model Ic), we choose
$\mathfrak{r}_{0}=3.5\mathfrak{r}_{g}$ and $\epsilon=0.1$. The initial data is
$e_{0}=0.01$, $I_{0}=60^{\circ}$, $\omega_{0}=60^{\circ}$ and
$\Omega_{0}=30^{\circ}$.  Figure 3: Typical orbits in the same model as Fig.
2. The red, blue and magenta curves denote a few cycles near $n_{0}\tau=2120$,
$4150$, and $6385$ respectively. The inclination angles of red, blue and
magenta curves are $I\approx 38.5^{\circ}$, $61.1^{\circ}$, and
$38.3^{\circ}$, while the eccentricities of those curves are $e\approx 0.80$,
$0.0017$ and $0.80$, respectively.

The KL oscillation period is obtained numerically as $n_{0}T_{\rm KL}\sim
4000$. Since $n_{0}P_{\rm in}\sim 2\pi$ and $\omega_{0}P_{\rm out}\sim 2\pi$,
we can evaluate it by Eq. (50) as $n_{0}{P_{\rm out}^{2}/P_{\rm in}}\sim
2\pi(n_{0}/\omega_{0})^{2}\sim 2500$, which is consistent with the above
numerical value.

####  V.3.3 Chaotic feature

One of the well-known features of three body system is a chaotic behaviour in
a binary motion. The model given in the previous subsection shows a stable KL
oscillation. The KL oscillation period is regular. The chaotic feature is not
seen. This is just because a binary is extremely compact and very hard. The
firmness parameter is $\mathfrak{f}\sim 350$.

Figure 4: Chaotic Kozai-Lidov oscillation. The red and blue curves denote the
eccentricity and the inclination. For this model (Model IIc), we choose
$\mathfrak{r}_{0}=3.5\mathfrak{r}_{g}$ and $\epsilon=0.4$. The initial data is
$e_{0}=0.01$, $I_{0}=60^{\circ}$, $\omega_{0}=60^{\circ}$ and
$\Omega_{0}=30^{\circ}$.

However, if a binary is not so highly compact, we find some chaotic features.
In Fig. 4, we show one example (Model IIc). The model parameters are given by
$\epsilon=0.4$ and $\mathfrak{r}_{0}=3.5\mathfrak{r}_{g}$, and the initial
parameters are chosen as $e_{0}=0.01$, $I_{0}=60^{\circ}$,
$\omega_{0}=60^{\circ}$ and $\Omega_{0}=30^{\circ}$. The larger value of
$\epsilon$ corresponds to a larger-scale binary, i.e., the initial semi-major
axis is $a_{0}=0.0058\mathfrak{r}_{g}\approx 0.0115{\rm AU}$. We can see
clearly the KL oscillation, but the period is not strictly regular. Since the
firmness parameter is $\mathfrak{f}\sim 22>O(1)$, the system is still stable,
but shows some irregular behaviours in the KL oscillations. The maximum values
of the eccentricity is also random as shown in Fig. 4.

####  V.3.4 Orbital flip

Another interesting feature is an orbital flip, i.e., an inclination angle
goes beyond $90^{\circ}$. It may occur when the initial inclination angle is
near $90^{\circ}$.

Figure 5: Orbital flip in KL oscillation (Model IIa). We choose
$\mathfrak{r}_{0}=3.5\mathfrak{r}_{g}$ and $\epsilon=0.4$. The initial data is
$e_{0}=0.01$, $I_{0}=85^{\circ}$, $\omega_{0}=60^{\circ}$ and
$\Omega_{0}=30^{\circ}$. The red and blue curves denote the eccentricity and
the inclination. The inclination angle evolves beyond $90^{\circ}$ several
times.

In Fig. 5, we find that the orbital flip accompanying KL oscillations occurs
periodically. The model parameters are given by $\epsilon=0.4$ and
$\mathfrak{r}_{0}=3.5\mathfrak{r}_{g}$, and the initial parameters are chosen
as $e_{0}=0.01$, $I_{0}=85^{\circ}$, $\omega_{0}=60^{\circ}$ and
$\Omega_{0}=30^{\circ}$ (Model IIa). The inclination blue curve in Fig. 5
evolves beyond $90^{\circ}$ several times, but the time period is irregular.

One interesting observation is there exists a strong correlation between an
orbital flip and large eccentricity. When an orbital flip occurs, the
eccentricity becomes very close to unity.

In Fig. 6, we show the time evolution of the angular momentum vector, which is
defined by $\\!\\!\mbox{ \boldmath$L$}=\\!\\!\mbox{
\boldmath$\mathsf{r}$}\times\\!\\!\mbox{ \boldmath$\mathsf{p}$}$. The
$z$-component of the angular momentum becomes negative near $n_{0}\tau\approx
60-62\,,73-80\,,$ and $92-97$. The corresponding vectors are shown by the
colored arrows $\\!\\!\mbox{ \boldmath$L$}_{1}\,,\\!\\!\mbox{
\boldmath$L$}_{2}\,,$ and $\\!\\!\mbox{ \boldmath$L$}_{3}$, respectively.

Figure 6: Time evolution of angular momentum. The parameters and initial
conditions are the same as those in Fig. 5. The red, light-red, orange,
yellow, green, magenta, and blue curves denote the angular momentum vectors in
the periods $n_{0}\tau=40-50,50-60,60-70,70-80,80-90,90-100$ and $100-110$,
respectively. The $z$-components of $\\!\\!\mbox{
\boldmath$L$}_{1}\,,\\!\\!\mbox{ \boldmath$L$}_{2}$ , and $\\!\\!\mbox{
\boldmath$L$}_{3}$ become negative.

####  V.3.5 Summary of various models

We summarize our numerical results in Table 1. We have simulated three types
of models (I, II and III). The parameters of Model I, Model II and Model III
are ($\mathfrak{r}_{0}=3.5\mathfrak{r}_{g}$, $\epsilon=0.1$),
($\mathfrak{r}_{0}=3.5\mathfrak{r}_{g}$, $\epsilon=0.4$), and
($\mathfrak{r}_{0}=7\mathfrak{r}_{g}$, $\epsilon=0.1$), respectively. For each
model, we choose the initial data as (a) $e_{0}=0.01\,,I_{0}=85^{\circ}$, (b)
$e_{0}=0.9\,,I_{0}=85^{\circ}$, (c) $e_{0}=0.01\,,I_{0}=60^{\circ}$, and (d)
$e_{0}=0.9\,,I_{0}=60^{\circ}$. The initial argument of periapsis and
ascending node are chosen as $\omega_{0}=60^{\circ}$ and
$\Omega_{0}=30^{\circ}$ for all models. We also performed the simulation with
different values of those two parameters, the results do not change so much.

In Model I, since a binary is very compact (the firmness parameter
$\mathfrak{f}\sim 350$), it is very stable. As shown in Fig. 2, we find the
regular KL oscillation, although the oscillation period $T_{\rm KL}$ is not
strictly constant but slightly disperse within $1-9\%$. When the initial
eccentricity is large (Model Ib and Id), the oscillation period gets smaller
and minimum eccentricity becomes larger. Except for Model Ic, the KL
oscillation type is the so-called libration, which shows the oscillation of
the argument of periapsis around $90^{\circ}$. In Model Ic, the KL oscillation
seems to be the rotation type, which denotes the argument of periapsis
increases monotonically. However in the present model, it does increase on
average but not monotonically (sometimes going back and forth).

Model | $\mathfrak{r}_{0}$ | $\epsilon$ | $a_{0}$ | $e_{0}$ | $I_{0}$ | $\omega_{0}$ | $\Omega_{0}$ | KL oscillation type | $n_{0}T_{\rm KL}$ | $e_{\rm min}/e_{\rm max}$ | $I_{\rm min}/I_{\rm max}$ | orbital flip
---|---|---|---|---|---|---|---|---|---|---|---|---
Ia |  |  |  | $0.01$ | $85^{\circ}$ | $60^{\circ}$ | $30^{\circ}$ |  $\mathbb{L}$ [$36^{\circ}\leq\omega\leq 144^{\circ}$] | $2396-2430$ | $0.01/0.995$ | $39.1^{\circ}/85.5^{\circ}$ | No
Ib |  |  |  | $0.9$ | $85^{\circ}$ | $60^{\circ}$ | $30^{\circ}$ |  $\mathbb{L}$ [$55^{\circ}\leq\omega\leq 125^{\circ}$] | $580-600$ | $0.672/0.998$ | $54.8^{\circ}/87.5^{\circ}$ | No
Ic | $3.5\mathfrak{r}_{g}$ | $0.1$ | $0.0046\mathfrak{r}_{g}$ | $0.01$ | $60^{\circ}$ | $60^{\circ}$ | $30^{\circ}$ | $\mathbb{R}^{*}$ | $3600-4300$ | $0.003/0.8$ | $38.5^{\circ}/60.5^{\circ}$ | No
Id |  |  |  | $0.9$ | $60^{\circ}$ | $60^{\circ}$ | $30^{\circ}$ |  $\mathbb{L}$ [$56^{\circ}\leq\omega\leq 124^{\circ}$] | $753-770$ | $0.49/0.95$ | $47.3^{\circ}/74.3^{\circ}$ | No
IIa |  |  |  | $0.01$ | $85^{\circ}$ | $60^{\circ}$ | $30^{\circ}$ | $\mathbb{R}^{*}$ | $50-78$ | ($0.15-0.34$)/ | ($33.6^{\circ}-42.0^{\circ}$)/ | Yes
|  |  |  |  |  |  |  |  |  | $(0.97-1.00)$ | $(137.9^{\circ}-143.9^{\circ})$ | (irregular)
IIb |  |  |  | $0.9$ | $85^{\circ}$ | $60^{\circ}$ | $30^{\circ}$ |  $\mathbb{L}$ [$(50^{\circ}-53^{\circ})\leq$ | $39-44$ | $(0.533-0.589)/$ | $(50.2^{\circ}-53.7^{\circ})/$ | Yes
|  |  |  |  |  |  |  |  $\omega\leq(126^{\circ}-131^{\circ})$] |  | $(0.978-0.999)$ | $(126.4^{\circ}-128.6^{\circ})$ | (irregular)
IIc | $3.5\mathfrak{r}_{g}$ | $0.4$ | $0.0116\mathfrak{r}_{g}$ | $0.01$ | $60^{\circ}$ | $60^{\circ}$ | $30^{\circ}$ | $\mathbb{R}^{*}$ | $116-250$ | $(0.04-0.08)/$ | $(36.5^{\circ}-37.7^{\circ})/$ | No
|  |  |  |  |  |  |  |  |  | $(0.82-0.90)$ | $(59^{\circ}-65^{\circ})$ |
IId |  |  |  | $0.9$ | $60^{\circ}$ | $60^{\circ}$ | $30^{\circ}$ |  $\mathbb{L}$ [$(50^{\circ}-56^{\circ})\leq$ | $54-60$ | $(0.41-0.49)/$ | $(43^{\circ}-46^{\circ})/$ | No
|  |  |  |  |  |  |  |  $\omega\leq(124^{\circ}-130^{\circ})$] |  | $(0.92-0.97)$ | $(69^{\circ}-75^{\circ})$ |
IIIa |  |  |  | $0.01$ | $85^{\circ}$ | $60^{\circ}$ | $30^{\circ}$ |  $\mathbb{L}$ [$38^{\circ}\leq\omega\leq 142^{\circ}$] | $6776-6928$ | $0.009/0.997$ | $39.1^{\circ}/85.2^{\circ}$ | No
IIIb |  |  |  | $0.9$ | $85^{\circ}$ | $60^{\circ}$ | $30^{\circ}$ |  $\mathbb{L}$ [$56^{\circ}\leq\omega\leq 124^{\circ}$] | $1585-1628$ | $0.675/0.998$ | $55.1^{\circ}/87.2^{\circ}$ | No
IIIc | $7\mathfrak{r}_{g}$ | $0.1$ | $0.0073\mathfrak{r}_{g}$ | $0.01$ | $60^{\circ}$ | $60^{\circ}$ | $30^{\circ}$ |  $\mathbb{L}$ [$44^{\circ}\leq\omega\leq 136^{\circ}$] | $9450-10600$ | $0.004/0.792$ | $38.3^{\circ}/60.3^{\circ}$ | No
IIId |  |  |  | $0.9$ | $60^{\circ}$ | $60^{\circ}$ | $30^{\circ}$ |  $\mathbb{L}$ [$56^{\circ}\leq\omega\leq 124^{\circ}$] | $2045-2074$ | $0.493/0.949$ | $47.5/74.5$ | No

Table 1: The properties of the KL oscillations. The oscillations in Models I
and III are regular, although there are small amount of dispersion in the KL
oscillation time scale. For Model II, the oscillations are irregular and there
appears dispersion in the maximum and minimum values of the eccentricity and
inclination, whose values are shown in the bracket (). $\mathbb{L}$ and
$\mathbb{R}$ denote libration and rotation type, respectively.
$\mathbb{R}^{*}$ means that the argument of periapsis does increase on average
but not monotonically, and sometimes going back and forth.

For Model II, a binary is slightly less compact (the firmness parameter
$\mathfrak{f}\sim 22$). It is still stable but becomes irregular both in the
oscillation period and in the amplitude as shown in Fig. 4. There appear
dispersion in the maximum and minimum values of the eccentricity and
inclination, whose values are shown in the bracket (). When the initial
inclination angle is large (Models IIa and IIb), we find irregular orbital
flips. The KL oscillation is either libration type (Models IIb and IId), in
which the maximum and minimum angles of the argument of the periapsis disperse
as shown by the brackets (), or irregular rotation type (Models IIa and IIc),
in which the inclination angle increases on average but not monotonically
(sometimes going back and forth).

Model III is the case with a hard binary as Model I, but its location is a
little far from SMBH ($\mathfrak{r}_{0}=7\mathfrak{r}_{g}$). Hence the
relativistic effect as well as the tidal force are smaller than those in Model
I. As a result, the KL oscillation becomes more regular and stable. The KL
oscillation timescale is larger than that in Model I because the outer orbital
period $P_{\rm out}$ is larger.

###  V.4 Motion of the CM of a binary

Next we show the motion of the CM of a binary. In order to know how much the
motion deviates from a circular motion with the radius $\mathfrak{r}_{0}$, we
solve the radial perturbation equation (35). Using the previous numerical
solutions, we integrate Eq. (35) with the initial conditions such that
$\mathfrak{r}_{(1)}(0)=0$ and $\dot{\mathfrak{r}}_{(1)}(0)=0$. Here, we
discuss some typical cases.

####  V.4.1 Model Ia

In this case, we find stable and regular oscillations with the period $T_{\rm
osc}\equiv{2\pi/k}\approx 2.65P_{\rm out}$, where $k$ is defined by Eq. (36)
and $P_{\rm out}={2\pi/\mathfrak{w}_{0}}$ is the period of the circular
motion. The oscillation center is given by
$\langle\mathfrak{r}_{1}\rangle\approx 1.60\times 10^{-4}\mathfrak{r}_{0}$
with the amplitude $\Delta\mathfrak{r}_{1}\approx 1.71\times
10^{-4}\mathfrak{r}_{0}$, but the center increases to $2.18\times
10^{-4}\mathfrak{r}_{0}$ when the eccentricity becomes close unity, keeping
the oscillation amplitude $\Delta\mathfrak{r}_{1}$ to be constant (See Fig. 7
(left)).

There exists a good correlation between the shift of the oscillation center
and the eccentricity. Since the oscillation amplitude is very small compared
with the circular radius $\mathfrak{r}_{0}$, the deviations from the circular
motion can be treated as perturbations, which confirms our approach.

Figure 7: The radial deviation $\mathfrak{r}_{(1)}$ normalized by the radius
of the circular orbit $\mathfrak{r}_{0}$ are shown by blue curves for Model Ia
(left) and Model IIa (right). The deviation from a circular orbit in Model Ia
is regular oscillations around $\langle\mathfrak{r}_{1}\rangle$, which
increases when the eccentricity increases (the red curve), while that in Model
IIa becomes irregular oscillations without correlation with the eccentricity
(the red curve).

####  V.4.2 Model IIa

In this case, the KL oscillation is stable but irregular. As a result, we find
irregular oscillations of the CM as shown by a blue curve in Fig 7 (right).
There is no correlation with the eccentricity $e$. Although the oscillations
are irregular, the amplitude is very small and then the deviation from the
circular motion can be treated as stable perturbations.

####  V.4.3 Other Models

We summarize the results in Table 2 for the Models given in Table 1. We find
that Models (I) and (III) give stable and regular oscillations of the CM. It
is just because those binary motions also show stable and regular KL
oscillations. On the other hand, Model (II) shows irregular oscillations
because of the irregular KL oscillations in those binary motions.

Model | $\mathfrak{r}_{0}({\rm AU})$ | $a_{0}({\rm AU})$ | $\epsilon$ | $\mathfrak{f}$ | $e_{0}$ | $I_{0}$ | $[\langle\mathfrak{r}_{1}\rangle\pm\Delta\mathfrak{r}_{1}]/\mathfrak{r}_{0}$ | $T_{\rm ocs}/P_{\rm out}$ | oscillation property
---|---|---|---|---|---|---|---|---|---
Ia |  |  |  |  | $0.01$ | $85^{\circ}$ | $[(1.60-2.18)\pm 1.71]\times 10^{-4}$ |  | regular (good correlation with $e$)
Ib |  |  |  |  | $0.9$ | $85^{\circ}$ | $[0.18\pm 3.51]\times 10^{-3}$ |  | regular
Ic | $6.91$ | $0.0091$ | $0.1$ | $350$ | $0.01$ | $60^{\circ}$ | $[(7.01-7.37)\pm 7.04]\times 10^{-4}$ | $2.65$ | regular (good correlation with $e$)
Id |  |  |  |  | $0.9$ | $60^{\circ}$ | $[(4.09-4.40)\pm 4.51]\times 10^{-4}$ |  | regular (good correlation with $e$)
IIa |  |  |  |  | $0.01$ | $85^{\circ}$ | $[(1.60-2.91)\pm 1.88]\times 10^{-4}$ |  | irregular (no correlation with $e$)
IIb |  |  |  |  | $0.9$ | $85^{\circ}$ | $[(2.89-3.18)\pm 2.85]\times 10^{-4}$ |  | regular (no correlation with $e$)
IIc | $6.91$ | $0.0229$ | $0.4$ | $21.9$ | $0.01$ | $60^{\circ}$ | $[(4.30-5.77)\pm 6.09]\times 10^{-4}$ | $2.65$ | irregular (no correlation with $e$)
IId |  |  |  |  | $0.9$ | $60^{\circ}$ | $[(4.05-4.31)\pm 4.41]\times 10^{-4}$ |  | irregular (no correlation with $e$)
IIIa |  |  |  |  | $0.01$ | $85^{\circ}$ | $[(1.23-1.56)\pm 1.28]\times 10^{-5}$ |  | regular (good correlation with $e$)
IIIb |  |  |  |  | $0.9$ | $85^{\circ}$ | $[(1.22-1.43)\pm 1.35]\times 10^{-5}$ |  | regular (good correlation with $e$)
IIIc | $13.8$ | $0.0144$ | $0.1$ | $700$ | $0.01$ | $60^{\circ}$ | $[(5.78-5.98)\pm 5.83]\times 10^{-5}$ | $1.32$ | regular (good correlation with $e$)
IIId |  |  |  |  | $0.9$ | $60^{\circ}$ | $[(3.14-3.33)\pm 3.38]\times 10^{-5}$ |  | regular (good correlation with $e$)

Table 2: The oscillations of the CM, which are the radial deviations from a
circular geodesic motion. $\langle\mathfrak{r}_{1}\rangle$ denotes the center
of the oscillations, while $\Delta\mathfrak{r}_{1}$ gives the amplitude of the
oscillations (or dispersion for the irregular oscillations (Model II)). The
typical oscillation period is given by $T_{\rm osc}\equiv 2\pi/k$, where $k$
is defined by Eq. (36).

###  V.5 Comparison with double-averaging approach

In Appendix B, we present one of standard approaches on hierarchical triple
system, which is the Lagrange planetary equations for the orbital parameters.
Since we are interested in the long-time behaviour such as the Kozai-Lidov
mechanism, taking averages of the Hamiltonian over two periods of the inner
and outer binaries, we can analyze the simplified doubly-averaged planetary
equations.

Here, we compare our numerical results with those obtained in the double-
averaging (DA) approach. We show the evolution of the eccentricity for Models
Ia, IIa and IIIa in the left, center and right panels of Fig. 8, respectively.

For Model Ia, two results agree very well although the oscillation period of
the DA approach is slightly longer than that of the direct integration (DI)
method. For Model IIa, the DA approach gives a regular periodic oscillations,
but the DI method shows irregular chaotic oscillations. Two results do not
agree well, although the maximum values of the eccentricity are almost the
same. For Model IIIa, two results agree almost completely. For other models,
we also find the similar results.

We conclude that the DA approach for Models I and III may give good results
although the period of oscillations deviates slightly. On the other hand, for
Model II, which shows chaotic feature in the KL oscillations, the DA approach
does not give correct results.

Figure 8: The evolution of the eccentricity for Models Ia (left), IIa (center)
and IIIa (right). The red curves denote the results by the double-averaging
(DA) approach, while the blue ones by the direct integration (DI) method. For
Models Ia and IIIa, two results agree well, but for Model IIa, they do not
agree well.

##  VI Concluding Remarks

In this paper, we discuss motion of a binary system near SMBH. Using Fermi-
normal coordinates, we set up Newtonian self-gravitating system in the local
proper reference frame. Assuming a circular geodesic observer around a
Schwarzschild SMBH, we write down the equations of motion of a binary. To
remove the interaction terms between the CM of a binary and its relative
coordinates, we introduce a small acceleration of the observer. As a result,
the CM follows the observer’s orbit, but its motion deviates from an exact
geodesic. Since the relative motion is decoupled from the system, we first
solve it, and then find the motion of the CM by the perturbation equations
with the small acceleration, which is given by the relative motion.

We show that there appears the KL oscillations when a binary is compact enough
and the inclination angle is larger than the critical value. If the firmness
parameter $\mathfrak{f}$ is larger than a few hundreds, the oscillations are
regular. However when $\mathfrak{f}$ is around a few tens, although we find
the stable KL oscillations, but the oscillations become irregular both in the
oscillation period and in the amplitude. Especially, if the initial
inclination is large, we find an orbital flip.

One of the most interesting and important subjects of a binary system near
SMBH is the gravitational waves. When the KL oscillations is found in a binary
motion, we expect a large amount of the GW emission because the eccentricity
becomes large. The large eccentricity also provides much higher frequencies
than that from a circular binary Randall and Xianyu (2019a); Hoang et al.
(2019); Gupta et al. (2020). Another interesting point on the GWs from the KL
oscillations is that the large amount of the GW emission repeats periodically
with the KL oscillation time scale. It is a good advantage in the observations
because we have a certain preparation time for next observations.

There are two GW sources in a hierarchical triple system: One is the GWs from
an inner binary and the other is those from the outer binary. The time scale
of the emission of gravitational waves from a circular binary is evaluated as
Maggiore (2007)

$\displaystyle\tau_{\rm GW}={5\over 256}{c^{5}R_{0}^{4}\over
G^{3}m^{2}\mu}\,,$

where $m=m_{1}+m_{2}$ and $\mu=m_{1}m_{2}/(m_{1}+m_{2})$ and $R_{0}$ is the
initial distance. Hence, the ratio of the time scale of an outer binary to
that of an inner binary is

$\displaystyle{\tau_{\rm outer}\over\tau_{\rm inner}}={m_{1}m_{2}\over
M^{2}}\left({\mathfrak{r}_{0}\over\ell_{\rm binary}}\right)^{4}.$

In the present model, as we discuss in §. V.1, we have some constraints. From
stability condition of a binary (V.1), we have

$\displaystyle{\ell_{\rm
binary}\over\mathfrak{r}_{0}}\,\mbox{\raisebox{-4.30554pt}{$\stackrel{{\scriptstyle\textstyle<}}{{\textstyle\sim}}$}}\,\left({m_{1}+m_{2}\over
2M}\right)^{1/3}\,,$

while from the validity of Newtonian dynamics, we find

$\displaystyle{\ell_{\rm
binary}\over\mathfrak{r}_{0}}\,\mbox{\raisebox{-4.30554pt}{$\stackrel{{\scriptstyle\textstyle>}}{{\textstyle\sim}}$}}\,\left({m_{1}+m_{2}\over
2M}\right){\mathfrak{r}_{g}\over\mathfrak{r}_{0}}\,.$

If a binary exists near the ISCO radius ($\mathfrak{r}_{0}\sim
3\mathfrak{r}_{g}$), we find

$\displaystyle 2\times
10^{-5}\,\mbox{\raisebox{-4.30554pt}{$\stackrel{{\scriptstyle\textstyle<}}{{\textstyle\sim}}$}}\,{\tau_{\rm
outer}\over\tau_{\rm
inner}}\,\mbox{\raisebox{-4.30554pt}{$\stackrel{{\scriptstyle\textstyle<}}{{\textstyle\sim}}$}}\,8\times
10^{16}\,.$

Hence, in most cases, gravitational waves from the outer orbit are less
effective compared with those from the inner binary. However when the binary
is close to instability range, it is not the case. In fact, in our examples
discussed in this paper, if we assume a circular binary, we find

$\displaystyle{\tau_{\rm outer}\over\tau_{\rm
inner}}\sim\left\\{\begin{array}[]{cl}3\times 10^{-3}&{\rm(Model~{}I)}\\\
8\times 10^{-5}&{\rm(Model~{}II)}\\\ 8\times 10^{-3}&{\rm(Model~{}III)}\\\
\end{array}\right.$

The GWs from the outer binary become much larger than those from the inner
binary. However if there exists the Kozai-Lidov oscillation, the emission time
scale is reduced by the factor Peters and Mathews (1963)

$\displaystyle F(e_{\rm in})\approx{768\over 429}\left(1-e_{\rm
in}^{2}\right)^{7/2},$

when $e_{\rm in}\approx 1$. As a result, the GWs from the inner binary may
become larger that those from the outer binary.

In recent years, three-body systems and the emission of GWs from them have
received significant attention Amaro-Seoane et al. (2010); Antonini and Perets
(2012); Hoang and Naoz (2018); Antonini et al. (2016b); Meiron et al. (2017);
Robson et al. (2018); Randall and Xianyu (2019b, a); Hoang et al. (2019);
Emami and Loeb (2019); Gupta et al. (2020). Our future work will involve
evaluating the GWs from the present hierarchical triple setting using the
black hole perturbation approach, since near the ISCO radius the quadrupole
formula may not be valid Cardoso et al. (2021).

In this paper, we assume that the CM of a binary moves along a circular orbit,
but an eccentric orbit is interesting to be studied since the Kozai-Lidov
oscillation may be modulated on a longer timescaleLithwick and Naoz (2011);
Katz et al. (2011); Naoz et al. (2013); Li et al. (2014); Liu et al. (2015).
However, for such a highly eccentric orbit, the present proper reference frame
expanded up to the second order of the spatial coordinates $x^{\hat{a}}$ may
not be sufficient. We may need higher-order terms in the metric, where the
derivatives of the Riemann curvature appearNesterov (1999); Delva and Angonin
(2012). Although the basic equations are very complicated, such an extension
is straightforward.

Another natural direction would be an extension to a rotating SMBH that may
allow us to study the precession of the binary orbit around the Kerr black
hole. Such systems can reveal the impact of spin on GWs emitted from a nearby
binary. Recent research Camilloni et al. (2023) considers the secular dynamics
of the binary system distorted by a much larger Kerr black hole’s tidal
forces. This is done by deriving the magnetic and electric tidal moment at
quadrupole orders.

###### Acknowledgements.

We thank Luc Blanchet, Vitor Cardoso, Eric Gourgoulhon and Haruka Suzuki for
useful discussions. This work was supported in part by JSPS KAKENHI Grant
Numbers JP17H06359, JP19K03857 (K.M.), and by JP20K03953 (H.O.). P.G.
acknowledges support from C.V. Raman fellowship grant. K.M. would like to
acknowledges the Yukawa Institute for Theoretical Physics at Kyoto University,
where the present work was begun during the Visitors Program of FY2022. He
would also thank Niels Bohr Institute/Niels Bohr International Academy, The
Max Planck Institute for Gravitational Physics (Albert Einstein Institute),
The laboratory AstroParticle and Cosmology (APC), and CENTRA/Instituto
Superior Técnico, where most of this work was performed.

##  Appendix A Coplanar Binary

In this Appendix, we analyze motion of a coplanar binary, that is
$\mathsf{z}=0\,,{\mathsf{p}}_{z}=0$. It is an exact solution for Eqs. (25) and
(28).

###  A.1 Coplanar motion: Numerical results

In the case of coplanar motion of a binary, the relative inclination angle $I$
is always zero. We then have the coupled equations (23), (24), (26), and (27)
for $\mathsf{x}$ and $\mathsf{y}$. We first show numerical results for Model I
($\mathfrak{r}_{0}=3.5\mathfrak{r}_{g}\,,\epsilon=0.1$) and II
($\mathfrak{r}_{0}=3.5\mathfrak{r}_{g}\,,\epsilon=0.4$) in Fig. 9. We choose
the initial conditions as $e_{0}=0.9$ and $\omega_{0}=60^{\circ}$. For Model
I, as shown in Fig. 9 (left), the orbit is approximately elliptic, but the
periapsis is rotating because of relativistic effect and the shift is quite
regular. On the other hand, for Model II, we find some irregular behaviours in
the orbit as shown in Fig. 9(right). As discussed in the text, this model is
not tightly bounded (the firmness parameter $\mathfrak{f}\sim 22$) and the
effects of the tidal force by SMBH is not so small. As a result, the system
shows some chaotic feature. In the present coplanar case, the orbital shape is
deformed from an ellipse.

Figure 9: The coplanar orbits for Model I (left) and for Model II (right). We
choose the initial eccentricity $e_{0}=0.9$. The orbits start from the red
point, and end at the blue point [($n_{0}\tau=200$ (Model I) and
$n_{0}\tau=50$ (Model II)].

It is further confirmed by the time evolution of the eccentricity, which is
given in Fig. 10. For Model I, the eccentricity is almost constant ($e\sim
0.9$). On the other hand, for Model II, it oscillates irregularly between
$e\sim 0.7$ and $0.9$. Although the orbital shape is not well-approximated by
an ellipse, we evaluate the eccentricity by the osculating orbit. If we use
the averages eccentricity over one cycle, the eccentricity oscillates between
$\langle e\rangle\sim 0.75$ and $0.87$ , which is slightly different from the
values in Fig. 10.

Figure 10: Time evolution of the eccentricity for Model I (the blue curve) and
for Model II (the red curve). The eccentricity for Model I is almost constant
($e\approx 0.9$), while that for Model II oscillates irregularly between
$e\sim 0.7$ and $0.9$.

###  A.2 Circular motion

There exists an exact circular motion of a binary as follows: Assuming a
circular solution as

$\displaystyle\zeta\equiv\mathsf{x}+i\mathsf{y}=\rho_{0}\exp\left[i\theta(\tau)\right],~{}~{}~{}~{}{\rm
with}~{}~{}\theta(\tau)=\bar{\theta}(\tau)+\mathfrak{w}_{\rm dS}\tau,$

where the radius $\rho_{0}$ is constant, we find two equations for
$\bar{\theta}$ as

$\displaystyle-\dot{\bar{\theta}}^{2}+{G(m_{1}+m_{2})\over\rho_{0}^{3}}-{\mathfrak{r}_{g}\over
4\mathfrak{r}_{0}^{3}}(1+3\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2})-{3\mathfrak{r}_{g}\over
4\mathfrak{r}_{0}^{3}}(1+\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2})\cos(2(\bar{\theta}+2\mathfrak{w}_{\rm
R}\tau))=0,$ (52) $\displaystyle\ddot{\bar{\theta}}+{3\mathfrak{r}_{g}\over
4\mathfrak{r}_{0}^{3}}(1+\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2})\sin(2(\bar{\theta}+2\mathfrak{w}_{\rm
R}\tau))=0.$ (53)

Since the derivative of Eq. (52) with respect to $\tau$ gives Eq. (53) unless
$\bar{\theta}=0$, we first solve Eq. (53). Setting $\eta\equiv
2(\bar{\theta}+2\mathfrak{w}_{\rm R}\tau)$, we find

$\displaystyle\ddot{\eta}+{3\mathfrak{r}_{g}\over
2\mathfrak{r}_{0}^{3}}(1+\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2})\sin\eta=0\,,$

which can be integrated as

$\displaystyle\dot{\eta}^{2}-{3\mathfrak{r}_{g}\over\mathfrak{r}_{0}^{3}}(1+\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2})\cos\eta={\rm
constant}\,.$

By use of this equation, we can eliminate the term with $\cos\eta$ in Eq.(52),
resulting in $\dot{\eta}$ and $\dot{\bar{\theta}}$ to be constant. It follows
that $\cos\eta=\cos(2(\bar{\theta}+2\mathfrak{w}_{R}\tau))$ is constant. We
then obtain a solution such that

$\displaystyle\bar{\theta}+2\mathfrak{w}_{\rm R}\tau={n\over 2}\pi\,,$

where $n$ is an integer. Eq. (52) becomes

$\displaystyle-4\mathfrak{w}_{\rm
R}^{2}=-{G(m_{1}+m_{2})\over\rho_{0}^{3}}+{\mathfrak{r}_{g}\over
4\mathfrak{r}_{0}^{3}}(1+3\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2})+{3\mathfrak{r}_{g}\over
4\mathfrak{r}_{0}^{3}}(1+\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2})(-1)^{n}\,.$

Since $\mathfrak{w}_{\rm R}=\sqrt{\mathfrak{r}_{g}/2\mathfrak{r}_{0}^{3}}$, we
obtain two analytic solutions as follows:
When $n$ is even,

$\displaystyle\rho_{0}$ $\displaystyle=$
$\displaystyle\rho_{0}^{(+)}\equiv\sqrt[3]{(m_{1}+m_{2})\over
2M(1+\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2})}~{}\mathfrak{r}_{0},$ (54)
$\displaystyle\theta$ $\displaystyle=$
$\displaystyle\theta^{(+)}\equiv\bar{\theta}^{(+)}+\mathfrak{w}_{\rm
dS}\tau=(\mathfrak{w}_{\rm dS}-2\mathfrak{w}_{\rm
R})\tau+m\pi~{}~{}(m\in\mathbb{Z})\,,$ (55)

while when $n$ is odd,

$\displaystyle\rho_{0}$ $\displaystyle=$
$\displaystyle\rho_{0}^{(-)}\equiv\sqrt[3]{(m_{1}+m_{2})\over
M(1+\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2})}~{}\mathfrak{r}_{0},$ (56)
$\displaystyle\theta$ $\displaystyle=$
$\displaystyle\theta^{(-)}\equiv\bar{\theta}^{(-)}+\mathfrak{w}_{\rm
dS}\tau=(\mathfrak{w}_{\rm dS}-2\mathfrak{w}_{\rm R})\tau+\left(m+{1\over
2}\right)\pi~{}~{}(m\in\mathbb{Z}).$ (57)

Using these solutions, we can find the analytic solution for the motion of the
CM of a binary. The relative coordinate $(x,y,z)$ of a binary in the rotating
proper frame are given by

$\displaystyle x$ $\displaystyle=$
$\displaystyle\mathsf{x}\cos(\mathfrak{w}_{\rm
R}\tau)-\mathsf{y}\sin(\mathfrak{w}_{\rm
R}\tau)=\rho_{0}^{(\pm)}\cos[(\mathfrak{w}_{\rm dS}-\mathfrak{w}_{\rm
R})\tau+\phi_{0}^{(\pm)}],$ $\displaystyle y$ $\displaystyle=$
$\displaystyle\mathsf{x}\sin(\mathfrak{w}_{\rm
R}\tau)+\mathsf{y}\cos(\mathfrak{w}_{\rm
R}\tau)=\rho_{0}^{(\pm)}\sin[(\mathfrak{w}_{\rm dS}-\mathfrak{w}_{\rm
R})\tau+\phi_{0}^{(\pm)}],$ $\displaystyle z$ $\displaystyle=$
$\displaystyle\mathsf{z}=0,$

where $\rho_{0}^{(\pm)}$ are given by the previous solutions, and
$\phi_{0}^{(+)}=m\pi$ and $\phi_{0}^{(-)}=\left(m+{1\over 2}\right)\pi$. We
then find the perturbation equation for $\mathfrak{r}_{(1)}$ as

$\displaystyle{d^{2}\mathfrak{r}_{(1)}\over
d\tau^{2}}+k\mathfrak{r}_{(1)}+C^{(\pm)}\cos^{2}[(\mathfrak{w}_{\rm
dS}-\mathfrak{w}_{\rm R})\tau+\phi_{0}^{(+)}]=0\,,$ (58)

where

$\displaystyle C^{(\pm)}\equiv\left[A+(\mathfrak{w}_{\rm dS}-\mathfrak{w}_{\rm
R})B\right]\left(\rho_{0}^{(+)}\right)^{2}\,.$

Since two independent solutions of the homogeneous equation are $\sin k\tau$
and $\cos k\tau$, introducing two unknown functions $u(\tau)$ and $v(\tau)$,
we set

$\displaystyle\mathfrak{r}_{(1)}=u(\tau)\sin k\tau+v(\tau)\cos k\tau\,.$

Inserting this into Eq. (58), we find

$\displaystyle\ddot{u}\sin k\tau+\ddot{v}\cos k\tau+2k\left(\dot{u}\cos
k\tau-\dot{v}\sin k\tau\right)+C^{(\pm)}\cos^{2}[(\mathfrak{w}_{\rm
dS}-\mathfrak{w}_{\rm R})\tau+\phi_{0}^{(\pm)}]=0\,.$

We assume one constraint equation such that

$\displaystyle\dot{u}\sin k\tau+\dot{v}\cos k\tau=0\,,$ (59)

which yields

$\displaystyle\ddot{u}\sin k\tau+\ddot{v}\cos k\tau+k\left(\dot{u}\cos
k\tau-\dot{v}\sin k\tau\right)=0.$

We then find

$\displaystyle k\left(\dot{u}\cos k\tau-\dot{v}\sin
k\tau\right)+C^{(\pm)}\cos^{2}[(\mathfrak{w}_{\rm dS}-\mathfrak{w}_{\rm
R})\tau+\phi_{0}^{(\pm)}]=0.$ (60)

From Eqs. (59) and (60), we obtain

$\displaystyle\dot{u}$ $\displaystyle=$ $\displaystyle-{C^{(\pm)}\over k}\cos
k\tau\cos^{2}[(\mathfrak{w}_{\rm dS}-\mathfrak{w}_{\rm
R})\tau+\phi_{0}^{(\pm)}],$ $\displaystyle\dot{v}$ $\displaystyle=$
$\displaystyle{C^{(\pm)}\over k}\sin k\tau\cos^{2}[(\mathfrak{w}_{\rm
dS}-\mathfrak{w}_{\rm R})\tau+\phi_{0}^{(\pm)}].$

We can integrate these equations as

$\displaystyle u$ $\displaystyle=$ $\displaystyle u_{0}-{C^{(\pm)}\over
2k}\left[{1\over k}\sin k\tau+{1\over 2}\left({\sin[(k+2(\mathfrak{w}_{\rm
dS}-\mathfrak{w}_{\rm R}))\tau+2\phi_{0}^{(\pm)}]\over k+2(\mathfrak{w}_{\rm
dS}-\mathfrak{w}_{\rm R})}+{\sin[(k-2(\mathfrak{w}_{\rm dS}-\mathfrak{w}_{\rm
R}))\tau-2\phi_{0}^{(\pm)}]\over k-2(\mathfrak{w}_{\rm dS}-\mathfrak{w}_{\rm
R})}\right)\right],$ $\displaystyle v$ $\displaystyle=$ $\displaystyle
v_{0}-{C^{(\pm)}\over 2k}\left[{1\over k}\cos k\tau+{1\over
2}\left({\cos[(k+2(\mathfrak{w}_{\rm dS}-\mathfrak{w}_{\rm
R}))\tau+2\phi_{0}^{(\pm)}]\over k+2(\mathfrak{w}_{\rm dS}-\mathfrak{w}_{\rm
R})}+{\cos[(k-2(\mathfrak{w}_{\rm dS}-\mathfrak{w}_{\rm
R}))\tau-2\phi_{0}^{(\pm)}]\over k-2(\mathfrak{w}_{\rm dS}-\mathfrak{w}_{\rm
R})}\right)\right],$

where $u_{0}$ and $v_{0}$ are integration constants. As a result, we obtain
general solution as

$\displaystyle\mathfrak{r}_{(1)}=-{C^{(\pm)}\over 2k^{2}}+u_{0}\sin
k\tau+v_{0}\cos k\tau-{C^{(\pm)}\over 2[k^{2}-4(\mathfrak{w}_{\rm
dS}-\mathfrak{w}_{\rm R})^{2}]}\cos\left[2(\mathfrak{w}_{\rm
dS}-\mathfrak{w}_{\rm R})\tau+2\phi_{0}^{(\pm)}\right].$

The initial conditions determines the integration constants $u_{0}$ and
$v_{0}$. For example, if we assume $\mathfrak{r}_{(1)}(0)=0$ and
$\dot{\mathfrak{r}}_{(1)}(0)=0$, we find

$\displaystyle\mathfrak{r}_{(1)}=-{C^{(\pm)}\over k^{2}}\sin^{2}{k\tau\over
2}+{(-1)^{n}C^{(\pm)}\over k^{2}-4(\mathfrak{w}_{\rm dS}-\mathfrak{w}_{\rm
R})^{2}}\sin\left[\left(\mathfrak{w}_{\rm dS}-\mathfrak{w}_{\rm R}+{k\over
2}\right)\tau\right]\sin\left[\left(\mathfrak{w}_{\rm dS}-\mathfrak{w}_{\rm
R}-{k\over 2}\right)\tau\right],\,$ (61)

where $n=2m$ or $2m+1$. As a result, $\mathfrak{r}_{(1)}$ oscillates around
zero. As for the other variables
$\mathfrak{t}_{(1)}\,,\mathfrak{\theta}_{(1)}\,,\mathfrak{\varphi}_{(1)}$,
although some of them may diverge as $\tau\rightarrow\infty$, no singularity
appear in the evolution equations. Hence, we conclude that the coplanar
circular orbit is linearly stable.

##  Appendix B Planetary equations for a binary system in Newtonian limit

In order to understand our numerical results, it may be better to introduce
the Lagrange planetary equations, which give time evolution of the orbital
parameters such as the semi-major axis, eccentricity and inclination. To
derive the planetary equations, we treat the proper Hamiltonian with unit mass
$\mu=1$, which is given by

$\displaystyle\bar{{\cal H}}=\bar{{\cal H}}_{0}+\bar{{\cal H}}_{1},$

where

$\displaystyle\bar{{\cal H}}_{0}$ $\displaystyle=$ $\displaystyle{1\over
2}\bar{\\!\\!\mbox{
\boldmath$\mathsf{p}$}}^{2}-{G(m_{1}+m_{2})\over\mathsf{r}},$
$\displaystyle\bar{{\cal H}}_{1}$ $\displaystyle=$ $\displaystyle\bar{{\cal
H}}_{1\mathchar 45\relax{\rm dS}}+\bar{{\cal H}}_{1\mathchar 45\relax\bar{\cal
R}},$

with

$\displaystyle\bar{{\cal H}}_{1\mathchar 45\relax{\rm dS}}$ $\displaystyle=$
$\displaystyle\mathfrak{w}_{\rm
dS}\left(\bar{\mathsf{p}}_{\mathsf{y}}\mathsf{x}-\bar{\mathsf{p}}_{\mathsf{x}}\mathsf{y}\right),$
$\displaystyle\bar{{\cal H}}_{1\mathchar 45\relax\bar{\cal R}}$
$\displaystyle=$ $\displaystyle{\mathfrak{r}_{g}\over
4\mathfrak{r}_{0}^{3}}\left[\mathsf{x}^{2}+\mathsf{y}^{2}+\mathsf{z}^{2}-3\left(1+\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2}\right)\left(\mathsf{x}\cos\mathfrak{w}_{\rm
R}\tau-\mathsf{y}\sin\mathfrak{w}_{\rm
R}\tau\right)^{2}+3\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2}\mathsf{z}^{2}\right].$

The position $\bf{\mathsf{r}}=(\mathsf{x},\mathsf{y},\mathsf{z})$ of a binary
should be described in the non-rotating proper reference frame.

For the unperturbed Hamiltonian $\bar{{\cal H}}_{0}$, it is just the same as
that of a binary in Newtonian dynamics. We find an elliptic orbit, which is
described by

$\displaystyle\mathsf{r}={a(1-e^{2})\over 1+e\cos f}\,,$ (62)

where $\mathsf{r},a,e$ and $f$ are the radial distance from the center of
mass, the semi-major axis, the eccentricity, and the true anomaly. This
orbital plane is inclined with the inclination angle $I$ from the equatorial
plane in the proper reference frame.

Hence, the relative position
$\bf{\mathsf{r}}=(\mathsf{x},\mathsf{y},\mathsf{z})$ of a binary is given by
the orbital parameters $(\omega\,,\Omega\,,a\,,e\,,I\,,f)$ as Eq. (47) with
Eq. (62). Introducing the Delaunay variables as

$\displaystyle\left\\{\begin{array}[]{l}\mathfrak{l}=n(t-t_{0})\\\
\mathfrak{g}=\omega\\\ \mathfrak{h}=\Omega\\\ \end{array}\right.~{}~{}{\rm
and}\hskip
14.22636pt\left\\{\begin{array}[]{l}\mathfrak{L}=\sqrt{G(m_{1}+m_{2})a}\\\
\mathfrak{G}=\sqrt{G(m_{1}+m_{2})a(1-e^{2})}\\\
\mathfrak{H}=\sqrt{G(m_{1}+m_{2})a(1-e^{2})}\cos I\\\ \end{array}\right.\,,$

where

$\displaystyle n\equiv{2\pi\over P}=\sqrt{{G(m_{1}+m_{2})\over a^{3}}},$

is the mean motion, we find new unperturbed Hamiltonian as

$\displaystyle\widetilde{\bar{\cal H}}_{0}=-{G^{2}(m_{1}+m_{2})^{2}\over
2\mathfrak{L}^{2}}.$

Including the perturbations $\bar{\cal H}_{1}$, we obtain the Hamiltonian the
Delaunay variables as

$\displaystyle\widetilde{\bar{\cal H}}=\widetilde{\bar{\cal H}}_{0}+\bar{\cal
H}_{1}\,.$

After some calculations, the Hamilton equations are reduced to

$\displaystyle\dot{a}$ $\displaystyle=$ $\displaystyle-{2\over
na}{\partial\bar{\cal H}_{1}\over\partial\mathfrak{l}},$ (64)
$\displaystyle\dot{e}$ $\displaystyle=$ $\displaystyle{\sqrt{1-e^{2}}\over
na^{2}e}{\partial\bar{\cal H}_{1}\over\partial\omega}-{1-e^{2}\over
na^{2}e}{\partial\bar{\cal H}_{1}\over\partial\mathfrak{l}},$ (65)
$\displaystyle\dot{I}$ $\displaystyle=$ $\displaystyle{1\over na^{2}\sin
I\sqrt{1-e^{2}}}{\partial\bar{\cal H}_{1}\over\partial\Omega}-{\cos I\over
na^{2}\sin I\sqrt{1-e^{2}}}{\partial\bar{\cal H}_{1}\over\partial\omega},$
(66) $\displaystyle\dot{\mathfrak{l}}$ $\displaystyle=$ $\displaystyle
n+{2\over na}{\partial\bar{\cal H}_{1}\over\partial a}+{1-e^{2}\over
na^{2}e}{\partial\bar{\cal H}_{1}\over\partial e},$ (67)
$\displaystyle\dot{\omega}$ $\displaystyle=$
$\displaystyle-{\sqrt{1-e^{2}}\over na^{2}e}{\partial\bar{\cal
H}_{1}\over\partial e}+{\cos I\over na^{2}\sin
I\sqrt{1-e^{2}}}{\partial\bar{\cal H}_{1}\over\partial I},$ (68)
$\displaystyle\dot{\Omega}$ $\displaystyle=$ $\displaystyle-{1\over na^{2}\sin
I\sqrt{1-e^{2}}}{\partial\bar{\cal H}_{1}\over\partial I}.$ (69)

The partial derivative $\partial/\partial\mathfrak{l}$ can be replaced by that
of the true anomaly as

$\displaystyle{\partial\over\partial\mathfrak{l}}=(1-e^{2})^{-{3\over
2}}\left(1+e\cos f\right)^{2}{\partial\over\partial f}.$

Hence, once we find the perturbation Hamiltoninan $\bar{\cal H}_{1}$ in terms
of the orbital parameters, we obtain the planetary equations.

The proper Hamiltonian is described by the orbital parameters by inserting the
relation given in Eq.(47) with Eq. (62). We then find the perturbed
Hamiltonian as

$\displaystyle\bar{\cal H}_{1}=\bar{{\cal H}}_{1\mathchar 45\relax{\rm
dS}}+\bar{{\cal H}}_{1\mathchar 45\relax\bar{\cal R}}\,,$

where

$\displaystyle\bar{{\cal H}}_{1\mathchar 45\relax{\rm dS}}$ $\displaystyle=$
$\displaystyle\mathfrak{w}_{\rm dS}\mathsf{r}^{2}(a,e,f)\left\\{n\cos
I(1-e^{2})^{-3/2}(1+e\cos f)^{2}-\mathfrak{w}_{\rm
dS}\left(\cos^{2}(\omega+f)+\sin^{2}(\omega+f)\cos^{2}I\right)\right\\}~{}~{}~{}~{}~{}~{}~{}$
(70) $\displaystyle\bar{{\cal H}}_{1\mathchar 45\relax\bar{\cal R}}$
$\displaystyle=$ $\displaystyle{\mathfrak{r}_{g}\over
4\mathfrak{r}_{0}^{3}}\mathsf{r}^{2}(a,e,f)\Big{\\{}1-3(1+\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2})\left[\cos(\Omega+\mathfrak{w}_{\rm
R}\tau)\cos(\omega+f)-\sin(\Omega+\mathfrak{w}_{\rm R}\tau)\sin(\omega+f)\cos
I\right]^{2}$ (71)
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+3\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2}\sin^{2}(\omega+f)\sin^{2}I\Big{\\}}\,.$

We then obtain the planetary equations for the present hierarchical triple
system from Eqs. (64) - (69).

###  B.1 Double-averaging approach

Here, instead of solving the Lagrange planetary equations themselves, which is
equivalent to our numerical methods in the text, we take average of the
perturbed Hamiltonian over two periods, the inner and outer orbital periods,
and then analyze the simplified equations, because we are interested in the
long-time behaviour of the present system such as the Kozai-Lidov mechanism.

The double-averaged Hamiltonian is defined by

$\displaystyle\langle\langle\bar{\cal H}_{1}\rangle\rangle\equiv{1\over
2\pi}\int_{0}^{2\pi}d\mathfrak{l}_{\rm out}\left({1\over
2\pi}\int_{0}^{2\pi}d\mathfrak{l}\,\bar{\cal H}_{1}\right)$

Since the outer orbit is circular, we find that $\mathfrak{l}_{\rm out}=f_{\rm
out}=\mathfrak{w}_{0}\tau$. We also have

$\displaystyle d\mathfrak{l}={1\over\sqrt{1-e^{2}}}\left({\mathsf{r}\over
a}\right)^{2}df\,.$

Inserting Eqs. (70) and (71) into the above integrals, we find the doubly-
averaged Hamiltonian as

$\displaystyle\langle\langle\bar{\cal H}_{1}\rangle\rangle$ $\displaystyle=$
$\displaystyle\mathfrak{w}_{\rm dS}na^{2}\sqrt{1-e^{2}}\cos I-{a^{2}\over
8}\Big{\\{}(2+3e^{2})\Big{[}\mathfrak{w}_{\rm dS}^{2}(3+\cos
2I)+{\mathfrak{r}_{g}\over
8\mathfrak{r}_{0}^{3}}(1+3\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2})(1+3\cos
2I)\Big{]}$ (72)
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+10e^{2}\sin^{2}I\cos
2\omega\left[\mathfrak{w}_{\rm dS}^{2}+{3\mathfrak{r}_{g}\over
8\mathfrak{r}_{0}^{3}}(1+3\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2})\right]\Big{\\}}.$

Using the doubly-averaged Hamiltonian Eq.(72), we obtain the Lagrange
planetary equations as

$\displaystyle\dot{e}$ $\displaystyle=$ $\displaystyle{5\over
4}\left(\mathfrak{w}_{\rm dS}^{2}+{3\mathfrak{r}_{g}\over
8\mathfrak{r}_{0}^{3}}(1+3\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2})\right){e\sqrt{1-e^{2}}\over
n}\left(1-\cos 2I\right)\sin(2\omega),$ (73) $\displaystyle\dot{I}$
$\displaystyle=$ $\displaystyle-{5\over 4}\left(\mathfrak{w}_{\rm
dS}^{2}+{3\mathfrak{r}_{g}\over
8\mathfrak{r}_{0}^{3}}(1+3\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2})\right){e^{2}\over
n\sqrt{1-e^{2}}}\sin 2I\sin(2\omega),$ (74) $\displaystyle\dot{\omega}$
$\displaystyle=$ $\displaystyle{1\over 4n}\left(\mathfrak{w}_{\rm
dS}^{2}+{3\mathfrak{r}_{g}\over
8\mathfrak{r}_{0}^{3}}(1+3\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2})\right)\Big{[}\sqrt{1-e^{2}}[3+5\cos
2I+5\left(1-\cos 2I\right)\cos 2\omega]$ (75)
$\displaystyle+{5e^{2}\over\sqrt{1-e^{2}}}(1+\cos 2I)(1-\cos
2\omega)\Big{]}+{2\sqrt{1-e^{2}}\over n}\mathfrak{w}_{\rm dS}^{2},$
$\displaystyle\dot{\Omega}$ $\displaystyle=$ $\displaystyle\mathfrak{w}_{\rm
dS}+{\cos I\over 2n\sqrt{1-e^{2}}}\left(\mathfrak{w}_{\rm dS}^{2}+{3\over
8}{\mathfrak{r}_{g}\over\mathfrak{r}_{0}^{3}}(1+3\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2})\right)\left[-(2+3e^{2})+5e^{2}\cos(2\omega)\right].$
(76)

The semi-major axis $a$ is constant in the present approximation. Also, we can
easily check from Eqs. (73) and (74) that

$\displaystyle{d\over d\tau}\left(\sqrt{1-e^{2}}\cos I\right)=0\,,$

which corresponds to conservation of the $z$-component of the angular
momentum.

###  B.2 KL oscillations

Introducing a “potential” by $V_{S}\equiv-\langle\langle\bar{\cal
H}_{1}\rangle\rangle$, we rewrite the above planetary equations as

$\displaystyle\dot{e}$ $\displaystyle=$ $\displaystyle-{\sqrt{1-e^{2}}\over
na^{2}e}{\partial{V_{S}}\over\partial\omega},$ (77) $\displaystyle\dot{I}$
$\displaystyle=$ $\displaystyle{\cos I\over na^{2}\sin
I\sqrt{1-e^{2}}}{\partial{V_{S}}\over\partial\omega},$ (78)
$\displaystyle\dot{\omega}$ $\displaystyle=$
$\displaystyle{\sqrt{1-e^{2}}\over na^{2}e}{\partial{V_{S}}\over\partial
e}-{\cos I\over na^{2}\sin I\sqrt{1-e^{2}}}{\partial{V_{S}}\over\partial I},$
(79) $\displaystyle\dot{\Omega}$ $\displaystyle=$ $\displaystyle{1\over
na^{2}\sin I\sqrt{1-e^{2}}}{\partial{V_{S}}\over\partial I}\,.$ (80)

We obtain the closed form of a set of the differential equations for
$e\,,I\,,$ and $\omega$ by Eqs. (77), (78) and (79). It gives several
properties of KL oscillations such as the oscillation amplitude of the
eccentricity and the oscillation time scale as analyzed in Newtonian and 1PN
hierarchical triple system Suzuki et al. (2020).

The potential is written by use of $\eta\equiv\sqrt{1-e^{2}}$ and
$\mu_{I}\equiv\cos I$ as

$\displaystyle V_{S}$ $\displaystyle\equiv$
$\displaystyle-\langle\langle\bar{\cal
H}_{1}\rangle\rangle={a^{2}\mathfrak{r}_{g}(1+3\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2})\over
32\mathfrak{r}_{0}^{3}}v_{S}(\eta,\mu_{I}),$

where

$\displaystyle v_{S}(\eta,\mu_{I})$ $\displaystyle\equiv$ $\displaystyle
2(-1+3\mu_{I}^{2}\eta^{2})(1+\alpha_{\rm dS})+12C_{KL}+4\alpha_{\rm
dS}\left(2-{3n\over\mathfrak{w}_{\rm dS}}\mu_{I}\eta\right),$

with

$\displaystyle\alpha_{\rm dS}$ $\displaystyle\equiv$
$\displaystyle{8\mathfrak{w}_{\rm dS}^{2}\mathfrak{r}_{0}^{3}\over
3\mathfrak{r}_{g}(1+3\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2})},$
$\displaystyle C_{\rm KL}$ $\displaystyle\equiv$
$\displaystyle(1-\eta^{2})\left[(1+2\alpha_{\rm dS})-{5\over 2}(1+\alpha_{\rm
dS})(1-\mu_{I}^{2})\sin^{2}\omega\right].$

Note that when $\alpha_{\rm dS}=0$, we find the same equations for Newtoninan
hierarchical triple system with quadrupole approximation. The terms with
$\alpha_{\rm dS}$ give relativistic corrections.

Introducing the normalized time $\tilde{\tau}$, which is defined by

$\displaystyle\tilde{\tau}\equiv{\tau\over\tau_{\rm KL}},$

with the typical KL time scale

$\displaystyle\tau_{\rm
KL}\equiv{32n\mathfrak{r}_{0}^{3}\over\mathfrak{r}_{g}(1+3\mathfrak{r}_{0}^{2}\mathfrak{w}_{0}^{2})}\left(\sim{P_{\rm
out}^{2}\over P_{\rm in}}\right)\,,$

the above planetary equation is rewritten as

$\displaystyle{d\eta\over d\tilde{\tau}}$ $\displaystyle=$
$\displaystyle{\partial{v_{S}}\over\partial\omega},$
$\displaystyle{1\over\mu_{I}}{d\mu_{I}\over d\tilde{\tau}}$ $\displaystyle=$
$\displaystyle-{1\over\eta}{\partial{v_{S}}\over\partial\omega},$
$\displaystyle{d\omega\over d\tilde{\tau}}$ $\displaystyle=$
$\displaystyle-{\partial{v_{S}}\over\partial\eta}+{\mu_{I}\over\eta}{\partial{v_{S}}\over\partial\mu_{I}}.$

From these equations, we can easily show that

$\displaystyle{d(\mu_{I}\eta)\over d\tilde{\tau}}=0~{}\,,~{}~{}~{}{dv_{S}\over
d\tilde{\tau}}=0\,,$

which means there exist two conserved quantities $\vartheta\equiv\mu_{I}\eta$
and $C_{\rm KL}$ just as the Newtoninan and 1PN hierarchical triple system
under dipole approximation. Using these two conserved quantities, we obtain a
single equation for $\eta$ as

$\displaystyle{d\eta^{2}\over
d\tilde{\tau}}=-24\sqrt{2}\sqrt{f(\eta^{2})g(\eta^{2})},$

with

$\displaystyle f(\eta^{2})$ $\displaystyle\equiv$
$\displaystyle(1+2\alpha_{\rm dS})(1-\eta^{2})-C_{\rm KL},$ $\displaystyle
g(\eta^{2})$ $\displaystyle\equiv$ $\displaystyle-5(1+\alpha_{\rm
dS})\vartheta^{2}+\left[5(1+\alpha_{\rm dS})\vartheta^{2}+3+\alpha_{\rm
dS}+2C_{\rm KL}\right]\eta^{2}-(3+\alpha_{\rm dS})\eta^{4}.$

Setting $\xi=\eta^{2}$, we find

$\displaystyle{d\xi\over d\tilde{\tau}}=-24\sqrt{2(1+2\alpha_{\rm
dS})(3+\alpha_{\rm dS})}\sqrt{(\xi-\xi_{0})(\xi-\xi_{+})(\xi-\xi_{-})},$

where

$\displaystyle\xi_{0}$ $\displaystyle=$ $\displaystyle 1-{C_{\rm KL}\over
1+2\alpha_{\rm dS}},$ $\displaystyle\xi_{\pm}$ $\displaystyle=$
$\displaystyle{1\over 2}\Big{[}\left(1+{5(1+\alpha_{\rm dS})\over
3+\alpha_{\rm dS}}\vartheta^{2}+{2\over 3+\alpha_{\rm dS}}C_{\rm
KL}\right)\pm\sqrt{\left(1+{5(1+\alpha_{\rm dS})\over 3+\alpha_{\rm
dS}}\vartheta^{2}+{2\over 3+\alpha_{\rm dS}}C_{\rm
KL}\right)^{2}-{20(1+\alpha_{\rm dS})\over 3+\alpha_{\rm
dS}}\vartheta^{2}}\Big{]},$

are the solutions of $f(\xi)=0$ and $g(\xi)=0$, respectively.

We can find the relativistic corrections with $\alpha_{\rm dS}$, which is
evaluated as

$\displaystyle\alpha_{\rm dS}$ $\displaystyle=$
$\displaystyle{3\mathfrak{r}_{g}^{2}\over\mathfrak{r}_{0}^{2}\left(1+\sqrt{1-{3\mathfrak{r}_{g}\over
2\mathfrak{r}_{0}}}\right)^{2}}\leq{1\over
3\left(1+{1\over\sqrt{2}}\right)^{2}}\approx 0.114382.$

The equality is found in the case of the ISCO radius
($\mathfrak{r}_{0}=3\mathfrak{r}_{g})$.

Analyzing the above equation, we find that there exists KL oscillation in this
system just the same as in Newtonian hierarchical triple system, and we can
classify the KL solutions by the sign of $C_{\rm KL}$ into two cases: (1)
$C_{\rm KL}>0$ (rotation) and (2) $C_{\rm KL}<0$ (libration).

####  B.2.1 $C_{\rm KL}>0$ (rotation)

In this case, $0<\xi_{-}<1<\xi_{+}$ and $0<\xi_{0}<1$. This is possible if

$\displaystyle 0<C_{\rm KL}<1+2\alpha_{\rm dS}\,.$

Hence we find the maximum and minimum values of the eccentricity as

$\displaystyle e_{\rm max}=\sqrt{1-\xi_{-}}~{}~{},~{}~{}e_{\rm
min}=\sqrt{1-\xi_{0}}.$

The KL time scale is given by

$\displaystyle T_{\rm KL}=\tau_{\rm KL}\,\mathfrak{T}_{\rm KL}^{\rm(rot)},$
(81)

where

$\displaystyle\mathfrak{T}_{\rm KL}^{\rm(rot)}\equiv{1\over
12\sqrt{2(1+2\alpha_{\rm dS})(3+\alpha_{\rm
dS})}}\times\int_{\xi_{-}}^{\xi_{0}}d\xi{1\over\sqrt{(\xi-\xi_{0})(\xi-\xi_{+})(\xi-\xi_{-})}}\,.~{}~{}~{}~{}~{}~{}~{}$
(82)

####  B.2.2 $C_{\rm KL}<0$ (libration)