Gravitational waves excited by binary stars orbiting around a supermassive black hole

The binary are supposed to revolve around each other along elliptic orbit. We introduce the relative coordinates $\vec{\tilde{r}}=\vec{r}_{1}-\vec{r}_{2}$, with $\tilde{r}=\frac{\tilde{p}}{1+\tilde{e}\cos\tilde{\beta}}$, then the coordinates of the binary with masses $m_{1}$ and $m_{2}$ can be expressed with coordinates of mass center and relative coordinates as $\vec{r}_{1}=\vec{r}_{C}+\vec{\tilde{r}}_{1},\vec{r}_{2}=\vec{r}_{C}+\vec{\tilde{r}}_{2}$, where $\vec{\tilde{r}}_{1}=\frac{m_{2}}{m_{1}+m_{2}}\vec{\tilde{r}},\;\vec{\tilde{r}}_{2}=-\frac{m_{1}}{m_{1}+m_{2}}\vec{\tilde{r}}$. We use the coordinates and parameters with tilde to describe the internal motion of binary stars around each other. The spin axis of the central SBH is along $z_{0}$, as displayed in Fig.1, $\theta_{0}$ is the angle between $z_{0}$ and $z$ which is the normal direction of the surface of binary stars’ internal elliptic orbit. $\phi_{0}$ is the angle between $x_{0}$ and the node line which is the intersection line between the surface of binary stars’ internal orbit and the horizonal surface.

Now let’s discuss the orbit of a mass point in Kerr spacetime, in Boyer-Lindquist coordinates Kerr metric can be expressed as

$$ds^{2}=-\left(1-\frac{2M}{\Sigma}\right)dt^{2}+\frac{\Sigma}{\Delta}dr^{2}-\frac{4aMr\sin^{2}\theta}{\Sigma}dtd\phi+\Sigma d\theta^{2}+\sin^{2}\theta\left(r^{2}+a^{2}+a^{2}\sin^{2}\theta\frac{2Mr}{\Sigma}\right)d\phi^{2},$$

(1)

where $\Delta=r^{2}-2Mr+a^{2}$ and $\Sigma=r^{2}+a^{2}\cos^{2}\theta$. We adopt the Hamilton-Jacobi formulation, the Hamilton-Jacobi equation reads

$$-\frac{\partial S}{\partial\tau}=\frac{1}{2}g^{\alpha\beta}\frac{\partial S}{\partial x^{\alpha}}\frac{\partial S}{\partial x^{\beta}}.$$

(2)

The Kerr BH admits two obvious Killing vectors $\partial_{t}$ and $\partial_{\phi}$, which correspond to two conserved charges, energy $E$ and angular momentum $L_{z}$, so for the mass point with mass $\mu$ ($\mu=m_{1}+m_{2}$ in this paper), the action reads

$$S=\frac{1}{2}\mu^{2}\tau-Et+L_{z}\phi+S_{r}(r)+S_{\theta}(\theta),$$

(3)

which leads to $p_{r}=dS_{r}/dr$ and $p_{\theta}=dS_{\theta}/d\theta$. Inserting (3) into the Hamilton-Jacobi equation (2) and separating variables we have

$$S_{r}(r)=\pm\int dr\frac{1}{\Delta}\sqrt{R(r)},\;S_{\theta}(\theta)=\pm\int d\theta\sqrt{\Theta(\theta)},$$

(4)

with

	$\displaystyle R(r)$	$\displaystyle\equiv[(r^{2}+a^{2})E-aL_{z}]^{2}-\Delta\left[\mu^{2}r^{2}+(L_{z}-aE)^{2}+Q\right],$		(5)
	$\displaystyle\Theta(\theta)$	$\displaystyle\equiv Q+(L_{z}-aE)^{2}-\mu^{2}a^{2}\cos^{2}\theta-\frac{1}{\sin^{2}\theta}\left[L_{z}-aE\sin^{2}\theta\right]^{2},$

where $Q$ is a constant which is related with the separation constant $C$ through $Q=C-(L_{z}-aE)^{2}$. Actually, the variable-separability of Hamilton-Jacobi equation in Kerr spacetime corresponds to a nontrivial symmetry which is described by a Killing tensor

$$K^{\mu\nu}=a^{2}\cos^{2}\theta g^{\mu\nu}+\frac{1}{\sin^{2}\theta}\delta^{\mu}_{\phi}\delta^{\nu}_{\phi}+a^{2}\sin^{2}\theta\delta^{\mu}_{t}\delta^{\nu}_{t}-2a\delta^{\mu}_{t}\delta^{\nu}_{\phi}+\delta^{\mu}_{\theta}\delta^{\nu}_{\theta}.$$

(6)

The Killing tensor satisfies the Killing equation $\nabla_{(\lambda}K_{\mu\nu)}=0$, the bracket denotes symmetrization of the indices. This is a highly nontrivial symmetry which corresponds to the third conserved quantity—Carter constant $Q$—in addition to energy and angular momentum[58]. The explicit integration expression of coordinates can be obtained by setting the partial derivative of Hamilton-Jacobi function with respect to the constants of motion to be zero, that’s to say the motions of particles in Kerr spacetime are integrable.

With the above results, one now is able to give the geodesics of a mass point in Kerr spacetime

	$\displaystyle\Sigma\frac{dr}{d\tau}$	$\displaystyle=\pm\sqrt{V_{r}},$		(7)
	$\displaystyle\Sigma\frac{d\theta}{d\tau}$	$\displaystyle=\pm\sqrt{V_{\theta}},$
	$\displaystyle\Sigma\frac{d\phi}{d\tau}$	$\displaystyle=V_{\phi},$
	$\displaystyle\Sigma\frac{dt}{d\tau}$	$\displaystyle=V_{t},$

with

	$\displaystyle V_{r}$	$\displaystyle=\left[E(r^{2}+a^{2})-L_{z}a\right]^{2}-\Delta\left[\mu^{2}r^{2}+(L_{z}-aE)^{2}+Q\right],$		(8)
	$\displaystyle V_{\theta}$	$\displaystyle=Q-\cos^{2}{\theta}\left[a^{2}(\mu^{2}-E^{2})+\frac{L_{z}^{2}}{\sin^{2}{\theta}}\right],$
	$\displaystyle V_{\phi}$	$\displaystyle=\frac{L_{z}}{\sin^{2}{\theta}}-aE+\frac{a}{\Delta}\left[E(r^{2}+a^{2})-L_{z}a\right],$
	$\displaystyle V_{t}$	$\displaystyle=a\left(L_{z}-aE\sin^{2}{\theta}\right)+\frac{r^{2}+a^{2}}{\Delta}\left[E(r^{2}+a^{2})-L_{z}a\right].$

To solve the equations of geodesic motions (7) numerically, we introduce two angular variables $\psi$ and $\chi$ in place of $r$ and $\theta$ to avoid possible singularities. $\psi$ is defined as

$$r=\frac{p}{1+e\cos\psi}.$$

(9)

It’s evident the periastron and apastron locate at $r_{p}=\frac{p}{1+e}$ and $r_{a}=\frac{p}{1-e}$. $\chi$ is defined as $z=\cos^{2}\theta=z_{-}\cos^{2}\chi$, where $z_{-}$ is given by

$$\beta(z_{+}-z)(z_{-}-z)=\beta z^{2}-z\left[Q+L_{z}^{2}+a^{2}(\mu^{2}-E^{2})\right]+Q,$$

(10)

with $\beta=a^{2}(\mu^{2}-E^{2})$. The parameter $\chi$ ranges from 0 to $2\pi$. As $\chi$ varies from 0 to $2\pi$, $\theta$ goes from one turning point $\theta_{min}$ to the other $\theta_{max}$ and back to $\theta_{min}$. We introduce an “inclination angle” to replace Carter constant $Q=L_{z}^{2}\tan^{2}\iota$. For $\iota=0$ or $\pi$, $Q$ equals to zero, from the geodesic equations (7) and (8) we know mass points remain on the equatorial plane at all times in this case. So $Q=0$ corresponds to the equatorial orbits while $Q\neq 0$ corresponds to non-equatorial orbits. Expanding the radial potential as

$$V_{r}=(\mu^{2}-E^{2})\,(r_{a}-r)\,(r-r_{p})\,(r-r_{3})\,(r-r_{4}),$$

(11)

we find evolution equations for $\psi$ and $\chi$ to be of the form

	$\displaystyle\frac{d\psi}{dt}$	$\displaystyle=\frac{\sqrt{\mu^{2}-E^{2}}\,\left[(p-r_{3}(1-e))-e(p+r_{3}(1-e)\cos\psi)\right]^{\frac{1}{2}}\,\left[(p-r_{4}(1+e))+e(p-r_{4}(1+e)\cos\psi)\right]^{\frac{1}{2}}}{\left[\gamma+a^{2}\,E\,z(\chi)\right](1-e^{2})},$		(12)
	$\displaystyle\frac{d\chi}{dt}$	$\displaystyle=\frac{\sqrt{\beta\,\left[z_{+}-z(\chi)\right]}}{\gamma+a^{2}\,E\,z(\chi)},$

where $\gamma=E\left[\left(r^{2}+a^{2}\right)^{2}/\Delta-a^{2}\right]-\frac{2MraL_{z}}{\Delta}$. Solving these equations allow us to know the positions of mass center of the binary at any time.

To calculate GWs, we need to calculate the positions of the two compact objects in binary system relative to the mass center. The mass center of the binary system is free-falling around the SBH, thus locally the spacetime around the mass center can be treated as Minkowskian spacetime. Since for binary system the size of their internal orbit is much smaller than that of their external orbit, i.e., $\tilde{p}$ is much smaller than $p$, so the internal motions of binary are dominated by internal Newtonian gravity. For convenience, the two-body problem is reduced to one-body problem, the relative motion is described by the Lagrangian

$$L=\frac{1}{2}m_{r}\dot{\vec{\tilde{r}}}^{2}-V^{(i)}(\tilde{r}),$$

(13)

where $m_{r}=\frac{m_{1}m_{2}}{m_{1}+m_{2}}$ is the reduced mass, and the potential $V^{(i)}(\tilde{r})=-\frac{m_{1}m_{2}}{\tilde{r}}$ describes the internal gravity.

Solving the aforementioned equations that determine the motions of mass center and the internal motions allow us to give the orbit of the binary. In Fig.2 we show the orbits of one of the compact objects in the binary system. From the lower three detail images in Fig.2 one sees that, due to internal motions the orbit have small fluctuations, and thus the orbits fluctuate more frequently as mass increases. The reason roots in that if keeping other parameters fixed, the cycle of orbit becomes shorter as the mass of gravitational source increases. The fluctuations are expected to give rise to higher-frequency GWs.

In the weak-field situation, the metric can be decomposed as $g_{\mu\nu}=\eta_{\mu\nu}+h_{\mu\nu}$, where $\eta_{\mu\nu}$ is the flat metric and $h_{\mu\nu}$ is small perturbation. By introducing the traceless tensor $\bar{h}^{\mu\nu}\equiv h^{\mu\nu}-(1/2)\eta^{\mu\nu}h$ with $h=\eta^{\mu\nu}h_{\mu\nu}$, then the Einstein equations can be recast to the form

$$\Box\bar{h}^{\mu\nu}=-16\pi T^{\mu\nu}.$$

(14)

Here $T^{\mu\nu}$ denotes the energy-momentum tensor of the source. In the slow motion limit, the GW is given by the quadrupole-octupole formula

$\displaystyle\bar{h}^{jk}=\frac{2}{r}\left[\ddot{I}^{jk}-2n_{i}\ddot{S}^{ijk}+n_{i}\dddot{M}^{ijk}\right]_{t^{\prime}=t-r},\$

(15)

with

	$\displaystyle I^{jk}(t^{\prime})$	$\displaystyle=\int x^{\prime j}x^{\prime k}T^{00}(t^{\prime},{\bf x^{\prime}})d^{3}x^{\prime},$		(16)
	$\displaystyle S^{ijk}(t^{\prime})$	$\displaystyle=\int x^{\prime j}x^{\prime k}T^{0i}(t^{\prime},{\bf x^{\prime}})d^{3}x^{\prime},$
	$\displaystyle M^{ijk}(t^{\prime})$	$\displaystyle=\int x^{\prime i}x^{\prime j}x^{\prime k}T^{00}(t^{\prime},{\bf x^{\prime}})d^{3}x^{\prime},$

being the mass quadrupole moment, current quadrupole moment and mass octupole moment respectively, and $r^{2}=\bf{x}\cdot\bf{x}$, ${\bf n}=\bf{x}/r$, where $\bf{x}$ is the location of the observer.

In the NK prescription, a flat-space trajectory which is “equivalent” to a geodesic in Kerr spacetime is reconstructed by projecting the Boyer-Lindquist coordinates onto a fictitious spherical polar grid. A particle moving along geodesic path in Kerr spacetime can be viewed as the particle moving along curved path in flat-space, so it’s convenient for us to calculate GWs with Cartesian coordinates. In Cartesian coordinates the positions of one of the compact objects is given by

	$\displaystyle x_{1}$	$\displaystyle=r\sin\theta\cos\phi-\frac{m_{2}}{m_{1}+m_{2}}\frac{\tilde{p}}{1+\tilde{e}\cos\tilde{\beta}}\left(\sin(\tilde{\beta}+\gamma)\cos\phi_{0}+\cos(\tilde{\beta}+\gamma)\cos\theta_{0}\sin\phi_{0}\right),$		(17)
	$\displaystyle y_{1}$	$\displaystyle=r\sin\theta\sin\phi+\frac{m_{2}}{m_{1}+m_{2}}\frac{\tilde{p}}{1+\tilde{e}\cos\tilde{\beta}}\left(-\sin(\tilde{\beta}+\gamma)\sin\phi_{0}+\cos(\tilde{\beta}+\gamma)\cos\theta_{0}\cos\phi_{0}\right),$
	$\displaystyle z_{1}$	$\displaystyle=r\cos\theta+\frac{m_{2}}{m_{1}+m_{2}}\frac{\tilde{p}}{1+\tilde{e}\cos\tilde{\beta}}\cos(\tilde{\beta}+\gamma)\sin\theta_{0},$

where $\gamma$ is the angle between the line connecting mass center and one of the compact objects and the line that lies in the plane of binary internal orbit and is perpendicular to the node line. For the position of the other compact object one only has to replace $\frac{m_{2}}{m_{1}+m_{2}}$ with $-\frac{m_{1}}{m_{1}+m_{2}}$ in the above expression.

The expression (15) are valid for a general extended source in flat-space. For a single mass point with mass $m_{1}$ the energy-momentum tensor in flat spacetime is given by

$$T^{\mu\nu}(t^{\prime},{\bf x^{\prime}})=\mu\int_{-\infty}^{\infty}\frac{\rm{d}x^{\prime\mu}_{p}}{\rm{d}\tau}\frac{\rm{d}x^{\prime\nu}_{p}}{\rm{d}\tau}\delta^{4}\left(x^{\prime}-x^{\prime}_{p}(\tau)\right)\;\;\rm{d}\tau=\mu\left(\frac{\rm{d}\tau}{\rm{d}t^{\prime}_{p}}\right)^{2}\frac{\rm{d}x^{\prime\mu}_{p}}{\rm{d}\tau}\frac{\rm{d}x^{\prime\nu}_{p}}{\rm{d}\tau}\delta^{3}\left({\bf x}^{\prime}-{\bf x}_{p}^{\prime}(t^{\prime})\right).$$

(18)

For binary stars we know the energy-momentum tensor is the sum of two such terms with $\mu=m_{1}$ and $m_{2}$ respectively. The momenta (16) now can be simplified to

	$\displaystyle I^{jk}$	$\displaystyle=\mu x^{\prime j}_{p}x^{\prime k}_{p}\;,$		(19)
	$\displaystyle S^{ijk}$	$\displaystyle=v^{i}\ I^{jk}\;,$
	$\displaystyle M^{ijk}$	$\displaystyle=x^{\prime i}_{p}\ I^{jk}\;.$

Substituting (19) into (15) one then obtains the GW generated by a mass point. For binary system the GW is the superposition of that generated by two compact objects.

In the standard transverse-traceless (TT) gauge, the waveform is given by the TT projection of (15). The orthonormal basis of spherical coordinates is defined as

$${\bf e}_{r}=\frac{\partial}{\partial r},\quad{\bf e}_{\Theta}=\frac{1}{r}\frac{\partial}{\partial\Theta},\quad{\bf e}_{\Phi}=\frac{1}{r\sin\Theta}\frac{\partial}{\partial\Phi},$$

(20)

where the angles $\{\Theta,\Phi\}$ denote the observation point’s latitude and azimuth respectively. After TT projection the waveform is given by

$$h^{jk}_{TT}=\frac{1}{2}\left(\begin{array}[]{ccc}0&0&0\\ 0&h^{\Theta\Theta}-h^{\Phi\Phi}&2h^{\Theta\Phi}\\ 0&2h^{\Theta\Phi}&h^{\Phi\Phi}-h^{\Theta\Theta}\end{array}\right),$$

(21)

with

	$\displaystyle h^{\Theta\Theta}$	$\displaystyle=\cos^{2}{\Theta}\left[h^{xx}\cos^{2}{\Phi}+h^{xy}\sin{2\Phi}+h^{yy}\sin^{2}{\Phi}\right]+h^{zz}\sin^{2}{\Theta}-\sin{2\Theta}\left[h^{xz}\cos{\Phi}+h^{yz}\sin{\Phi}\right],$		(22)
	$\displaystyle h^{\Phi\Theta}$	$\displaystyle=\cos{\Theta}\left[-\frac{1}{2}h^{xx}\sin{2\Phi}+h^{xy}\cos{2\Phi}+\frac{1}{2}h^{yy}\sin{2\Phi}\right]+\sin{\Theta}\left[h^{xz}\sin{\Phi}-h^{yz}\cos{\Phi}\right],$
	$\displaystyle h^{\Phi\Phi}$	$\displaystyle=h^{xx}\sin^{2}{\Phi}-h^{xy}\sin{2\Phi}+h^{yy}\cos^{2}{\Phi}.$

The usual “plus” and “cross” waveform polarizations are given by $h^{\Theta\Theta}-h^{\Phi\Phi}$ and $2h^{\Theta\Phi}$ respectively.

Due to the GW emissions, the compact objects will lose energy and angular momentum, although the amount of energy and angular momentum carried away by GWs is very small, the trajectory will alter accordingly, for a more precise study we consider backreaction of energy loss to the orbit. We study generic inclined-eccentric orbits of the center of mass of binary stars rather than just orbits on the equatorial plane. As we will see, since the amplitude of GW generated by internal motions of binary stars is small compared to that generated by external motions, we only consider the radiation reaction caused by the external motion of binary stars around the SBH while that caused by the internal motion is neglected. We consider the flux of generic orbits to 2nd order post-Newtonian approximation[57]

	$\displaystyle(\dot{E})_{\rm 2PN}$	$\displaystyle=-\frac{32}{5}\frac{\mu^{2}}{M^{2}}\left(\frac{M}{p}\right)^{5}(1-e^{2})^{3/2}\left[g_{1}(e)-q\left(\frac{M}{p}\right)^{3/2}g_{2}(e)\cos\iota-\left(\frac{M}{p}\right)g_{3}(e)+\pi\left(\frac{M}{p}\right)^{3/2}g_{4}(e)\right.$		(23)
		$\displaystyle\left.-\left(\frac{M}{p}\right)^{2}g_{5}(e)+q^{2}\left(\frac{M}{p}\right)^{2}g_{6}(e)-\frac{527}{96}q^{2}\left(\frac{M}{p}\right)^{2}\sin^{2}\iota\right],\$
	$\displaystyle(\dot{L}_{z})_{\rm 2PN}$	$\displaystyle=-\frac{32}{5}\frac{\mu^{2}}{M}\left(\frac{M}{p}\right)^{7/2}(1-e^{2})^{3/2}\left[g_{9}(e)\cos\iota+q\left(\frac{M}{p}\right)^{3/2}\{g_{10}^{a}(e)-\cos^{2}\iota g_{10}^{b}(e)\}-\left(\frac{M}{p}\right)g_{11}(e)\cos\iota\right.$
		$\displaystyle\left.+\pi\left(\frac{M}{p}\right)^{3/2}g_{12}(e)\cos\iota-\left(\frac{M}{p}\right)^{2}g_{13}(e)\cos\iota+q^{2}\left(\frac{M}{p}\right)^{2}\cos\iota\left(g_{14}(e)-\frac{45}{8}\,\sin^{2}\iota\right)\right],$
	$\displaystyle(\dot{Q})_{\rm 2PN}$	$\displaystyle=-\frac{64}{5}\frac{\mu^{2}}{M}\left(\frac{M}{p}\right)^{7/2}\,\sqrt{Q}\,\sin\iota\,(1-e^{2})^{3/2}\,\left[g_{9}(e)-q\left(\frac{M}{p}\right)^{3/2}\cos\iota g_{10}^{b}(e)-\left(\frac{M}{p}\right)g_{11}(e)\right.$
		$\displaystyle\left.+\pi\left(\frac{M}{p}\right)^{3/2}g_{12}(e)-\left(\frac{M}{p}\right)^{2}g_{13}(e)+q^{2}\left(\frac{M}{p}\right)^{2}\,\left(g_{14}(e)-\frac{45}{8}\,\sin^{2}\iota\right)\right],$

where $q=a/M$. The angle $\iota$ shows the inclination of internal orbit with respect to the external orbit through the Carter constant,

$$Q=L_{z}^{2}\tan^{2}\iota.$$

(24)

We defer the expressions of the coefficients $g_{1}(e)\sim g_{14}(e)$ to the Appendix. This expression of fluxes are the results of combining Tagoshi’s[46], Shibata’s [53] and Ryan’s fluxes [55, 56]. Numerical analysis indicates that these results ameliorate the unphysical inspiral properties for nearly circular and polar orbits, such as the inspiral moves rapidly away from circularity for the trajectories near zero eccentricity $e=0$, and there is discontinuity between trajectories with $\iota$ slightly less than $90^{\circ}$ and those slightly above $90^{\circ}$, etc., thus the fluxes (23) produce more physical reasonable results than previous ones. From the fluxes (23), we obtain the evolution of eccentricity and semi-latus rectum

	$\displaystyle\dot{e}$	$\displaystyle=\frac{\partial e}{\partial E}\dot{E}+\frac{\partial e}{\partial Lz}\dot{Lz}+\frac{\partial e}{\partial Q}\dot{Q},$		(25)
	$\displaystyle\dot{p}$	$\displaystyle=\frac{\partial p}{\partial E}\dot{E}+\frac{\partial p}{\partial Lz}\dot{Lz}+\frac{\partial p}{\partial Q}\dot{Q},$

we obtain

	$\displaystyle\dot{e}=$	$\displaystyle\frac{1-e^{2}}{2p}\left((1+e)\frac{N_{1}(r_{p})}{N(r_{p})}-(1-e)\frac{N_{1}(r_{a})}{N(r_{a})}\right)\dot{E}+\frac{1-e^{2}}{2p}\left((1+e)\frac{N_{2}(r_{p})}{N(r_{p})}-(1-e)\frac{N_{2}(r_{a})}{N(r_{a})}\right)\dot{L}$		(26)
		$\displaystyle+\frac{1-e^{2}}{2p}\left((1+e)\frac{N_{3}(r_{p})}{N(r_{p})}-(1-e)\frac{N_{3}(r_{a})}{N(r_{a})}\right)\dot{\iota},$
	$\displaystyle\dot{p}=$	$\displaystyle\left(-\frac{(1+e)^{2}}{2}\frac{N_{1}(r_{p})}{N(r_{p})}-\frac{(1-e)^{2}}{2}\frac{N_{1}(r_{a})}{N(r_{a})}\right)\dot{E}+\left(-\frac{(1+e)^{2}}{2}\frac{N_{2}(r_{p})}{N(r_{p})}-\frac{(1-e)^{2}}{2}\frac{N_{2}(r_{a})}{N(r_{a})}\right)\dot{L}$
		$\displaystyle+\left(-\frac{(1+e)^{2}}{2}\frac{N_{3}(r_{p})}{N(r_{p})}-\frac{(1-e)^{2}}{2}\frac{N_{3}(r_{a})}{N(r_{a})}\right)\dot{\iota},$

where we have

	$\displaystyle N_{1}(r)$	$\displaystyle=Er^{4}+a^{2}Er^{2}-2aM(L_{z}-aE)r,$		(27)
	$\displaystyle N_{2}(r)$	$\displaystyle=-L_{z}\sec^{2}\iota\;r^{2}+2M(L_{z}\sec^{2}\iota-aE)r-a^{2}L_{z}\tan^{2}\iota,$
	$\displaystyle N_{3}(r)$	$\displaystyle=L_{z}^{2}\tan\iota\sec^{2}\iota(r(2M-r)-a^{2}),$
	$\displaystyle N(r)$	$\displaystyle=-2(\mu^{2}-E^{2})r^{3}+3\mu^{2}Mr^{2}-(a^{2}(\mu^{2}-E^{2})+L_{z}^{2}\sec^{2}\iota)r+L_{z}^{2}M\tan^{2}\iota+(L_{z}-aE)^{2}M.$

Using the relation between inclination angle and Carter constant $Q=L_{z}^{2}\tan^{2}\iota$, one has

$$\dot{\iota}=\frac{1}{2L_{z}^{2}\tan(\iota)\sec^{2}(\iota)}\dot{Q}-\frac{\sin(\iota)\cos(\iota)}{L_{z}}\dot{L_{z}}.$$

(28)

With all the results above, we are able to calculate the GW numerically with the NK method. In Figs.3-6, we exhibit the $h_{+}$ waveforms of EMRIs and B-EMRIs for different parameter selections. The orientation of the internal elliptic orbit plane is determined by the angles $\theta_{0}$ and $\phi_{0}$, $\phi_{0}$ is fixed to be $\frac{\pi}{6}$. The observer is supposed to be located at $r=D$, $\Theta=\frac{\pi}{3}$, $\Phi=0$. From the figure we see that, the waveforms of B-EMRIs have the identical profile as that of EMRIs, but the waveforms of B-EMRIs have oscillations which are caused by the internal motion of the binary. We will analyze the waveforms in detail in the next section.

The internal orbit is calculated by Newtonian dynamics in Eq. 13. It is interesting and important to consider the higher order corrections e.g. the Gravito-electromagnetic force effect [60]. Unlike single free-falling compact object which always locates at the origin of the free-fall frame (FFF), every compact objects of the binary system will not locate at the origin of FFF due to the internal interaction of binary. According to equivalence principle only at the origin of FFF the spacetime is flat, if there is departure from the origin gravity will play a role. So for the binary we take into account the gravitational effects caused by the departure from the origin of FFF. For simplicity, we consider a circular orbit on the equatorial plane in the following. The natural choice of FFF is Riemann normal coordinates

$$ds^{2}=-\left(1+R_{0i0j}x^{i}x^{j}\right)c^{2}d\tau^{2}-\frac{4}{3}R_{0jik}x^{j}x^{k}d\tau dx^{i}+\left(\delta_{ij}-\frac{1}{3}R_{ikjl}x^{k}x^{l}\right)dx^{i}dx^{j}.$$

(29)

If $R_{ikjl}$, which is small compared to $R_{0i0j}$ and $R_{0jik}$, is neglected, the metric above can be rewritten as

$$ds^{2}=-\left(1-2\frac{\Phi}{c^{2}}\right)c^{2}d\tau^{2}-\frac{4}{c}(\vec{A}\cdot d\vec{x})cd\tau+\delta_{ij}dx^{i}dx^{j},$$

(30)

where

	$\displaystyle\Phi$	$\displaystyle=-\frac{1}{2}R_{0i0j}x^{i}x^{j},$		(31)
	$\displaystyle A_{i}$	$\displaystyle=\frac{1}{3}R_{0jik}x^{j}x^{k},$

are the scalar and vector potentials of GEM fields. Analogous to the electromagnetic fields, the GEM force can be introduced

$$\vec{F}=-m\vec{E}-2m\frac{\vec{v}}{c}\times\vec{B},$$

(32)

where

	$\displaystyle E_{i}$	$\displaystyle=R_{0i0j}x^{j},$		(33)
	$\displaystyle B_{i}$	$\displaystyle=-\frac{1}{2}\epsilon_{ijk}R_{0l}^{\;\;\;jk}x^{l}$

In FFF the equations of motion of the compact object in binary system can be wriiten as

$$m_{i}\frac{d^{2}\vec{x}_{i}}{d\tau^{2}}=-m_{i}m_{j}\frac{\vec{x}_{i}-\vec{x}_{j}}{|\vec{x}_{i}-\vec{x}_{j}|^{3}}+\vec{F}_{i}$$

(34)

where $i=1,2$ label the compact objects which constitute the binary.

Although the equations of motion can be written in a simple form in FFF, it is not straightforward to calculate GEM force in FFF since Riemann tensor is conventionally derived in local non-rotating frame (LNRF). The transformations between the coordinates of LNRF ($X^{t},X^{r},X^{\theta},X^{\phi}$) and BL coordinates are given by

	$\displaystyle dX^{t}$	$\displaystyle=(\Sigma\Delta/A)^{1/2}dt,$		(35)
	$\displaystyle dX^{r}$	$\displaystyle=(\Sigma/\Delta)^{1/2}dr,$
	$\displaystyle dX^{\theta}$	$\displaystyle=\Sigma^{1/2}d\theta,$
	$\displaystyle dX^{\phi}$	$\displaystyle=-\frac{2Mar\sin\theta}{(\Sigma A)^{1/2}}+(A/\Sigma)^{1/2}\sin\theta d\phi,$

where $A=(r^{2}+a^{2})^{2}-a^{2}\Delta\sin^{2}\theta$. It’s shown Riemann tensor in LNRF takes very simple form[61, 62]. For circular orbit on the equatorial plane, the free-fall observer is moving in the $X^{\phi}$ direction with the speed[61]

$$u=\frac{\pm M^{1/2}(r\mp 2aM^{1/2}r^{1/2}+a^{2})}{\Delta^{1/2}(r^{3/2}\pm aM^{1/2})}$$

(36)

where the upper signs refer to prograde orbits and the lower ones refer to retrograde orbits. The LNRF and the local inertial frame (LIF) are connected by a Lorentz transformation

$$\Lambda^{\mu}_{\mu^{\prime}}=\left(\begin{array}[]{cccc}\gamma&0&0&-\gamma\beta\\ 0&1&0&0\\ 0&0&1&0\\ -\gamma\beta&0&0&\gamma\\ \end{array}\right),$$

(37)

where $\beta=u/c,\gamma=\sqrt{1-\beta^{2}}$. Through the Lorentz transformation (37) Riemann tensor in LNRF can be transformed to the ones in LIF $R_{\mu^{\prime}\nu^{\prime}\rho^{\prime}\sigma^{\prime}}=\Lambda^{\mu}_{\mu^{\prime}}\Lambda^{\nu}_{\nu^{\prime}}\Lambda^{\rho}_{\rho^{\prime}}\Lambda^{\sigma}_{\sigma^{\prime}}R_{\mu\nu\rho\sigma}$. If a free gyro placed at the origin of FFF, its spin axis is fixed in FFF, but in general the gyro will precess relative to the LNRF and LIF, so there exist rotation between LIF and FFF, the angular velocity is $\omega=-\sqrt{M/r^{3}}$[63]. The tensors in LIF and the that in FFF are connected through a rotation transformation. The components of GEM force in FFF can be given in terms of that in LIF as

	$\displaystyle F_{x}$	$\displaystyle=F_{r^{\prime}}\cos(\omega\tau)+F_{\phi^{\prime}}\sin(\omega\tau),$		(38)
	$\displaystyle F_{y}$	$\displaystyle=F_{\theta^{\prime}},$
	$\displaystyle F_{z}$	$\displaystyle=-F_{r^{\prime}}\sin(\omega\tau)+F_{\phi^{\prime}}\cos(\omega\tau).$

Substituting the GEM force (38) into the equations of motion (34) and solving the equations allows to obtain the positions of compact objects in FFF. Considering the orbit of origin of FFF is geodesic in external SBH spacetime, it’s straightforward to obtain the positions of compact objects of the binary in SBH spacetime, based on this one is able to calculate the GWs with the NK method. In Fig.7 we present the waveform of GW with and without taking into account GEM force.