Analytic Keplerian-type parametrization for general spinning compact binaries with the leading order spin-orbit interaction

TL;DR

This paper presents an analytic Keplerian-type parametrization for the conservative motion of spinning compact binaries, enabling efficient and accurate waveform computation for gravitational wave analysis.

Contribution

It introduces a fully analytic solution for spinning binary dynamics at leading order in spin and post-Newtonian approximation, applicable to arbitrary eccentricities and mass ratios.

Findings

Provides a high-accuracy waveform computation method.

Describes three-dimensional binary motion with arbitrary parameters.

Offers a solution valid up to linear order in spin and leading PN order.

Abstract

We derive a fully analytic Keplerian-type parametrization solution to conservative motion of spinning binary in ADM gauge. This solution is able to describe three dimensional motion of binaries of arbitrary eccentricity, mass ratio and initial configuration of spin angular momentum up to the leading order of post-Newtonian(PN) approximation and a linear order in spin. Based on our results waveforms can be quickly computed with high accuracy.

Equations

$$H=\frac{{\mathbf{p}}^2}{2}-\frac{1}{r}+\frac{1}{c^2{r}^3}\mathbf{L}\cdot {\mathbf{S}}_{\text{eff}},$$

$$\mathbf{L}=\mathbf{r}\times \mathbf{p},$$

$${\mathbf{S}}_{\text{eff}}={\delta }_1{\mathbf{S}}_1+{\delta }_2{\mathbf{S}}_2.$$

$${\delta }_1=2\frac{m_1{m}_2}{M^2}\left(1+\frac{3}{4}\frac{m_2}{m_1}\right),$$

$${\delta }_2=2\frac{m_1{m}_2}{M^2}\left(1+\frac{3}{4}\frac{m_1}{m_2}\right).$$

$$\dot{\mathbf{r}}:=\frac{d\mathbf{r}}{dt}=\{\mathbf{r},H\}.$$

$$\frac{d\mathbf{L}}{dt}$$

$$\frac{d{\mathbf{S}}_1}{dt}$$

$$\frac{d{\mathbf{S}}_2}{dt}$$

$$\{S_i,S_j\}={ϵ}_{ijk}{S}_k.$$

$$\frac{d(\mathbf{L}\cdot \mathbf{L})}{dt}=\frac{2}{c^2{r}^3}\mathbf{L}\cdot ({\mathbf{S}}_{\text{eff}}\times \mathbf{L})=0,$$

$$\frac{d({\mathbf{S}}_1\cdot {\mathbf{S}}_1)}{dt}=\frac{2{\delta }_1}{c^2{r}^3}{\mathbf{S}}_1\cdot (\mathbf{L}\times {\mathbf{S}}_1)=0,$$

$$\frac{d({\mathbf{S}}_2\cdot {\mathbf{S}}_2)}{dt}=\frac{2{\delta }_2}{c^2{r}^3}{\mathbf{S}}_2\cdot (\mathbf{L}\times {\mathbf{S}}_2)=0.$$

$$\displaystyle\left(\begin{array}[]{ccc}\bm{i}\\ \bm{j}\\ \bm{k}\\ \end{array}\right)=\Lambda\,\left(\begin{array}[]{ccc}\bm{e}_{x}\\ \bm{e}_{y}\\ \bm{e}_{z}\\ \end{array}\right)$$

$$\displaystyle\Lambda=\left(\begin{array}[]{ccc}\cos\xi_{1}&\sin\xi_{1}&0\\ -\sin\xi_{1}\,\cos\xi_{2}&\cos\xi_{1}\,\cos\xi_{2}&\sin\xi_{2}\\ \sin\xi_{1}\sin\xi_{2}&-\cos\xi_{1}\sin\xi_{2}&\cos\xi_{2}\\ \end{array}\right)\,.$$

$$\cos \gamma :=\frac{{\mathbf{S}}_1\cdot {\mathbf{S}}_2}{S_1{S}_2},$$

$$\cos {\kappa }_1:=\frac{\mathbf{L}\cdot {\mathbf{S}}_1}{L{S}_1},$$

$$\cos {\kappa }_2:=\frac{\mathbf{L}\cdot {\mathbf{S}}_2}{L{S}_2},$$

$$\frac{d\cos \gamma }{dt}$$

$$\frac{d\cos {\kappa }_1}{dt}$$

$$\frac{d\cos {\kappa }_2}{dt}$$

$${\sigma }_1$$

$${\sigma }_2$$

$$\cos \Delta \psi =\frac{\cos \gamma -\cos {\kappa }_1\cos {\kappa }_2}{\sin {\kappa }_1\sin {\kappa }_2}.$$

$$\displaystyle\big{|}\bm{L}(t_{0}+dt)-\bm{L}(t_{0})\big{|}^{2}=dt^{\,2}\,\left(\frac{L\,S_{2}\,\delta_{2}\,\sin\kappa_{2}}{c^{2}\,r^{3}}\right)^{2}\,.$$

$$2L^2(1-\cos d{\kappa }_1)=L^{2}d{\kappa }_1^2+O(d{\kappa }_1^3).$$

$$\frac{\mathbf{L}\cdot ({\mathbf{S}}_1\times {\mathbf{S}}_2)}{L{S}_1{S}_2}=\sin {\kappa }_1\sin {\kappa }_2\sin \Delta \psi ,$$

$$=$$

$$\displaystyle+2\,\delta_{2}\,S_{2}\,\big{[}L\,\sigma_{1}\,(\delta_{2}-\delta_{1})+\sigma_{2}\,(\delta_{1}S_{1}+\delta_{2}\sigma_{1}S_{2})\big{]}\,\cos\kappa_{1}-\delta_{2}^{2}\,S_{2}^{2}\,(\sigma_{1}^{2}+\sigma_{2}^{2}-1)\bigg{\}}^{1/2}\,.$$

$$\frac{\mathbf{L}\cdot ({\mathbf{S}}_1\times {\mathbf{S}}_2)}{L{S}_1{S}_2}=\frac{\pm }{{\delta }_2{S}_2}\sqrt{A(x-x_1)(x-x_2)(x-x_3)}.$$

$$\frac{d\mathbf{L}\cdot ({\mathbf{S}}_1\times {\mathbf{S}}_2)}{dt}$$

$$=$$

$$+$$

$$\frac{dx}{dt}$$

$$\frac{d\cos {\kappa }_1}{\sqrt{A(\cos {\kappa }_1-x_1)(\cos {\kappa }_1-x_2)(\cos {\kappa }_1-x_3)}}=\frac{\pm dt}{c^2{r}^3}.$$

$$(r.h.s)=$$

Full text

††preprint: APS/123-QED

I Introduction

Under Newtonian gravity which has only mass as a source, the two body problem can be analytically solved into the well-known Kepler parametrized solution. But under the Einstein gravity which has a stress energy tensor of 10 components as its source, parameters more than mass must be considered. In the classical works Mathisson (2010); Corinaldesi and Papapetrou (1951); Description (1977); Bailey and Israel (1983), multipolar expansion of stress energy tensor on a generic background metric field, has been made under the assumption that matter is highly localized in space. The results, which involve mass as a monopole contribution while spin as a dipole moment, give us understanding on how masses and spin angular momenta of point particles affect the entire dynamics. Unlike Newtonian gravity, spin angular momentum produces extra gravitational effects and thus makes orbit plane of two body and spin axes precess.

Even though dynamics of binary systems in general relativity can be solved by numerical relativity (NR) with high accuracy, it is still worth solving the equation of motion analytically. Actually, the post-Newtonian (PN) approximation which has been well established Blanchet (2014); Bernard et al. (2018), provides spin precession equation Á. Gergely et al. (1998); Gergely (1999); Blanchet et al. (2006); Faye et al. (2006) and the orbital evolution equations up to at least 4PN order Damour et al. (2014); Bernard et al. (2018). Although the PN equations can be integrated numerically with much smaller amount of computer resources than NR but analytical solution would be desirable for quick parameter estimation Veitch et al. (2015); Farr et al. (2016). Until now analytic solutions have been obtained only for some limited cases such as quasi-circular, nearly alignedKlein et al. (2013), slowly spinning Chatziioannou et al. (2013) or almost equal massKönigsdörffer and Gopakumar (2005); Tessmer (2009). General cases with arbitrary spin can be solved up to a certain PN order in terms of elliptic functions Racine (2008) and recent works present a form of the solutions partially in terms of sine-like Jacobi elliptic function with integrating out the oscillating part of orbital motion Marsat et al. (2014); Chatziioannou et al. (2017). In this paper, we present an analytic parametrized solution which involves both spin precession and orbital evolution for the binaries with arbitrary spins, masses and eccentricities of orbit. Namely, our solution parameterizes tilt angle of orbital planes and precession of spin angular momentum with respect to a fixed inertial frame, radial and angular motion of relative displacement of binaries and physical time evolution only via the eccentric anomaly, a Keplerian parameter that serves as an independent variable. Additionally because we express the entire parameterization in closed form witout any approximations, the solution is fully general.

This paper is organized as follows. In §II, dynamics which we are going to solve is presented in terms of Hamiltonian. The Hamiltonian includes Newtonian order of orbital dynamics and the leading order of spin-orbit interaction term (linear in spin). Parameters and reference frames which we choose are also presented. In §III, the ways to get the Keplerian-type parametrizations and their results are given in the following orders : Relative motion of angles between angular momenta (§III.1), the absolute precession of the orbital angular momentum in an inertial frame (§III.2), and finally, the relative position of compact binaries (§III.3). In §IV, we provide an approximant for almost equal mass case and compare it with the previous work Königsdörffer and Gopakumar (2005) as a sanity check. The final secion summarized our results.

II Dynamics with leading order of spin-orbit interaction

II.1 The Hamiltonian

In Arnowitt-Deser-Misner(ADM) type cooridnate Arnowitt et al. (2004) and with the Newton-Wigner-Pryce supplementary spin condition Newton and Wigner (1949); Pryce (1948), the (reduced) Hamiltonian $H$ for compact bianary systems composed of two stars with masses $m_{1}$ and $m_{2}$, is as follows Wex (1995),

[TABLE]

where $c$ is the speed of light, $\bm{r}=\frac{\bm{R}}{G\,(m_{1}+m_{2})}$, $\bm{p}=\frac{\bm{\mathcal{P}}\,(m_{1}+m_{2})^{2}}{m_{1}\,m_{2}}$ with $\bm{R}$ and $\bm{\mathcal{P}}$ being relative position vector and conjugate momentum, respectively. The orbital angular momentum $\bm{L}$ and effective spin $S_{\text{eff}}$ are defined as

[TABLE]

and

[TABLE]

Here $\bm{S}_{1}$ and $\bm{S}_{2}$ are spin vectors of $m_{1}$ and $m_{2}$ and the waiting factors of $\delta_{1}$ and $\delta_{2}$ are,

[TABLE]

For a consistent post-Newtonian approximation, we need to include the first correction of orbital dynamics (1PN) but intentionally exclude it because 1PN correction to orbital motion does not affect the spin precession at the leading order and bringing the 1PN correction into the result is expected to be straightforward Königsdörffer and Gopakumar (2005); Damour and Deruelle (1985).

II.2 The equation of motion

The equation of motion is given by Poisson bracket,

[TABLE]

Likewise the corresponding angular momentum and spin precession equations are given as

[TABLE]

by invoking the fact that angular momentum is a generator of rotations i.e.

[TABLE]

From the evolution of angular momentum and spins, we can find several conserved quantities that can be exploited in deriving the dynamical solutions. First, the magnitudes of the angular momenta are conserved because

[TABLE]

Also in the absence of the angular momentum loss from the binary system, the total angular momentum $\bm{J}=\bm{L}+\bm{S}_{1}+\bm{S}_{2}$ is conserved. Without loss of generality we can define the z-axis of the inertial frame to be aligned with $\bm{J}$ while the the orthogonal axes $x$ and $y$ can be chosen arbitrarily on the plane perpendicular to $z$-axis as shown in Fig. 1.

At the same time, we define an orthonormal non inertial frame $(\bm{i},\bm{j},\bm{k})$,

[TABLE]

where $\Lambda$ is the Euler matrix as shown below

[TABLE]

The geometric meaning of the Euler matrix is presented in Fig. 1.

Next, we define the angles $\gamma$, $\kappa_{1}$ and $\kappa_{2}$ in the range of $[0,\,\pi]$ such that

[TABLE]

as illustrated in Fig. 2.

The rates of changes of these angles are given via Eq.(6) i.e.

[TABLE]

By observing their common factor $\bm{L}\cdot(\bm{S}_{1}\times\bm{S}_{2})$, we can define new conserved quantities $\sigma_{1}$, $\sigma_{2}$ as follows

[TABLE]

Note that if $\delta_{1}-\delta_{2}=0$, or equal mass case, then $\sigma_{1}=\cos\gamma$ implying that the angle $\gamma$ between $\bm{S}_{1}$ and $\bm{S}_{2}$ is conserved.

Through straightforward vector algebra, one can show that $\bm{L}\cdot(\bm{S}_{1}\times\bm{S}_{2})=L\,S_{1}\,S_{2}\,\sin\kappa_{1}\sin\kappa_{2}\sin\Delta\psi$, where $\Delta\psi$ is an angle between the projection of $\bm{S}_{1}$ and $\bm{S}_{2}$ onto x-y plane of the non-inertial frame where $\bm{L}$ is along with $z$-axis (See Fig. 2). It has the following relation with the others,

[TABLE]

When either $\sin\kappa_{1}=0$ or $\sin\kappa_{2}=0$, $\Delta\psi$ is not defined in the geometric sense. But physically we can extend the definition of $\Delta\psi$ properly. For example let us assume that $\kappa_{1}=0$ and $\sin\kappa_{2}\neq 0$ at the initial time $t_{0}$. Then, during an infinitesimal time interval $dt$, $d\bm{S}_{1}=0$ because $|\frac{d\bm{S}_{1}}{dt}|\sim|\bm{L}\times\bm{S}_{1}|\sim\sin\kappa_{1}=0$. On the other hand, $d\bm{L}=\bm{L}(t_{0}+dt)-\bm{L}(t_{0})=-\frac{dt\,\delta_{2}}{c^{2}\,r^{3}}\bm{L}\times\bm{S}_{2}$. Therefore the square of the $d\bm{L}$ is given as

[TABLE]

Because $\bm{L}(t_{0})$ and $\bm{S}_{1}(t_{0})$ are aligned and $\bm{S}_{1}(t_{0}+dt)=\bm{S}_{1}(t_{0})$, the infinitesimal change of $\kappa_{1}$ is the angle between $\bm{L}(t_{0}+dt)$ and $\bm{L}(t_{0})$ . Thus the right hand side of Eq.(24) is

[TABLE]

Since $\kappa_{1}$ can only increase during $dt$, we get $\frac{d\kappa_{1}}{dt}=\frac{\delta_{2}\,S_{2}\,\sin\kappa_{2}}{c^{2}\,r^{3}}$. Thus by setting $\sin\Delta\psi=-1$, replacing the vectorial expression $\bm{L}\cdot(\bm{S}_{1}\times\bm{S}_{2})$ in Eqs.(21) to the algebraic expression $L\,S_{1}\,S_{2}\,\sin\kappa_{1}\sin\kappa_{2}\sin\Delta\psi$ becomes valid even in the case of $\kappa_{1}=0$. Likewise we arrive at same extension $\sin\Delta\psi=-1$ when $\kappa_{2}=\pi$ (with $\sin\kappa_{1}\neq 0$), while $\sin\Delta\psi=1$ in the case of $\kappa_{1}=\pi$ ($\sin\kappa_{2}\neq 0$) and $\kappa_{2}=0$ ($\sin\kappa_{1}\neq 0$). Mathematically, this kind of extension is nothing but taking the limit of $\sin\kappa_{1,2}\rightarrow 0$ along the physically preferred paths. If $\sin\kappa_{1}=0$ and $\sin\kappa_{2}=0$ are both satisfied at some point, since every angle becomes constant, we do not need to define $\Delta\psi$.

We have three time dependent variables $\kappa_{1}(t)$, $\kappa_{2}(t)$ and $\gamma(t)$ to solve and two conserved quantities, $\sigma_{1}$ and $\sigma_{2}$. So, if we determine a numerical value of any one angle (say $\kappa_{1}$), the others are determined automatically by the initially fixed values of $\sigma_{1}$ and $\sigma_{2}$. This fact allows the relations of Eqs.(21), where $\kappa_{1}$, $\kappa_{2}$ and $\gamma$ are non-linearly entangled, to be expressed in a single variable equation, for example $\kappa_{1}$. This specific choice of $\kappa_{1}$ is followed from the fact that $\gamma$ is not a good parameter. As mentioned earlier, $\gamma$ fails to parametrize the dynamics when mass ratio becomes close to 1. Choosing $\gamma$ as a paramter was used in Kesden et al. (2015), and we can see that there are divergences as $q\rightarrow 1$ in the equations (7a) and (7b) in Kesden et al. (2015).

III The Keplerian-type parametrization

III.1 The Keplerian-type parametrization of the Angles

In this section, we solve the Eq.(21)s. We reduce the common factor $\bm{L}\cdot(\bm{S}_{1}\times\bm{S}_{2})$, using $\sigma_{1}$ and $\sigma_{2}$, into single variable ($\cos\kappa_{1}$) dependent expression,

[TABLE]

The sign of $\bm{L}\cdot(\bm{S}_{1}\times\bm{S}_{2})$ is determined by the sign of $\sin\Delta\psi$. For simplicity, let us denote $\cos\kappa_{1}=x$ and factorize the expression inside the square root in terms of its 3rd polynomial roots of $x=(x_{1},x_{2},x_{3})$ and introduce a coefficient $A:=2\,L\,S_{1}\,\delta_{1}\,(\delta_{2}-\delta_{1})$ to simplify the above equation as

[TABLE]

Before proceeding integrations, one needs to clarify the existence of the real roots of $x$. Here is a brief proof: Let us assume that there is no instance such that $\kappa_{1}(t)$, $\kappa_{2}(t)$ and $\gamma(t)$ satisfy $\bm{L}\cdot(\bm{S}_{1}\times\bm{S}_{2})=0$ during evolutions. This assumption says that if $\bm{L}\cdot(\bm{S}_{1}\times\bm{S}_{2})>0$ initially this inequalitiy always holds by virtue of continuity. Also the followings hold forever, $\frac{d\cos\kappa_{1}}{dt}>0$, $\frac{d\cos\gamma}{dt}>0$ (if $\delta_{1}>\delta_{2}$) and $\frac{d\cos\kappa_{2}}{dt}<0$. Then we are able to find a real value $T$ such that $\cos\kappa_{1}>0$, $\cos\gamma>0$ and $\cos\kappa_{2}<0$ at times $t>T$ unless $\frac{d\cos\kappa_{1}}{dt}\rightarrow 0+$, $\frac{d\cos\gamma}{dt}\rightarrow 0+$ and $\frac{d\cos\kappa_{2}}{dt}\rightarrow 0-$ fast enough as $t\rightarrow\infty$. But we exclude this dissipative possibility because governing dynamics is symmetric under the reflection of $t\rightarrow-t$.

Meanwhile the time derivative of $\bm{L}\cdot(\bm{S}_{1}\times\bm{S}_{2})$ is

[TABLE]

It is obvious that $\frac{d\bm{L}\cdot(\bm{S}_{1}\times\bm{S}_{2})}{dt}<0$ when $t>T$. Similarly we can estimate its second time derivative and we can decompose it into two components i.e. $\frac{d^{2}\bm{L}\cdot(\bm{S}_{1}\times\bm{S}_{2})}{dt^{2}}=\hat{A}(t)+\hat{B}(t)\,\dot{r}$ where $\hat{A}(t)<0$, $\hat{B}(t)>0$ but $\dot{r}$ oscillates sinusoidally between positive and negative values with zero averaged value i.e. $\frac{1}{2\pi}\int_{0}^{2\pi}du\,\dot{r}=0$. Thus in much longer time scale than one orbital period, the absolute value of $\frac{d\bm{L}\cdot(\bm{S}_{1}\times\bm{S}_{2})}{dt}$ increases, i.e. decreasing of $\bm{L}\cdot(\bm{S}_{1}\times\bm{S}_{2})$ accelerates. This suggests that at some point $\bm{L}\cdot(\bm{S}_{1}\times\bm{S}_{2})$ will pass zero, or $\bm{L}\cdot(\bm{S}_{1}\times\bm{S}_{2})>0$ cannot hold forever. This is contradictory to the assumption. Likewise another assumption that $\bm{L}\cdot(\bm{S}_{1}\times\bm{S}_{2})<0$ also encounters the same contradiction. Finally we come to the conclusion that whatever the initial condition is, binary systems eventually meet the configuration which corresponds to $\bm{L}\cdot(\bm{S}_{1}\times\bm{S}_{2})=0$. Actually this happens twice, once when $\bm{L}\cdot(\bm{S}_{1}\times\bm{S}_{2})$ is increasing or once when it is decreasing. Let us match cosine value of these two $\kappa_{1}$ values to $x_{2}$ and $x_{3}$. Then it is obvious that they are real and the absolute values are less than or equal to $1$. Additionally, from the fact that cubic polynomicals of which all coefficients are real, cannot have two real roots and a single complex root, we can conclude that $x_{1}$ is also real. Note that our assumption does not include the equal mass case, i.e., $\delta_{1}=\delta_{2}$ in the above proof. Since there is no reason for any discontinuity toward the equal mass case, we can extend the validity of our formalism to $\delta_{1}=\delta_{2}$ although we need a special caution in applying the formalism presented below when two masses become very close. The case of nearly equal mass is treated separately in §IV.

Without loss of generality, we assume that $A>0$ and $x_{2}\leq\cos\kappa_{1}\leq x_{3}$, because it must be non negative within the square root as illustrated in Fig. 3.

Eq.(21b) becomes,

[TABLE]

or,

[TABLE]

It is enough that $r$ is given in the Newtonian order, because we restrict the spin orbit interaction only up to the leading order. In Newtonian order we can use the well known expressions of $r=a\,(1-e\,\cos u)$ and $n\,(t-t_{0})=u-e\,\sin u$, where $u$ is the eccentric anomaly, $n$ is the mean angular speed and $t$ is physical time with an initial value $t_{0}$ when $u=0$, $a$ is the semi-major axis of the ellipse, and $e$ is the eccentricity. Thus the right hand side of Eq.(30) can be integrated to

[TABLE]

where we use the Newtonian expression of the magnitude of the orbital angular momentum $L=n^{\frac{1}{3}}\,a\,(1-e^{2})^{\frac{1}{2}}+O(\frac{1}{c})$ with $\nu=2\,\arctan(\sqrt{\frac{1+e}{1-e}}\,\tan\frac{u}{2})$ and $\alpha$ refers to proper initial value which will be defined later. Note that the given expression of Eq.(III.1) is somewhat deceiving. Because the $\pm$ sign in front of $\nu$ dependent term changes as $\cos\kappa_{1}$ passes its maximum $x_{3}$ or minimum $x_{2}$. For example, if $\cos\kappa_{1}$ arrives at $x_{3}$ when $\nu=\nu_{3}$, then when $\nu>\nu_{3}$ it would be better if right hand side of Eq.(III.1) is written as

[TABLE]

hence it has a periodic structure behind and more complicated feature than presented in Eq.(III.1). But as we can see later, the presented expression of Eq.(III.1) will be justified.

Now let us solve the left hand side of Eq.(30) keeping in mind that $x_{2}\leq x\leq x_{3}<x_{1}$,

[TABLE]

Since $0\leq\frac{x-x_{2}}{x_{3}-x_{2}}\leq 1$, we can introduce a new variable $y$ defined through

[TABLE]

which is restricted to $0\leq y\leq\frac{\pi}{2}$. Then the last line of the Eq.(33) becomes

[TABLE]

Identifying this integral as the Elliptic Integral of the first kind $F$,

[TABLE]

the Eq.(35) can be written in the following compact form

[TABLE]

There could be other possible forms instead of Eq.(37) but, under the assumption of $x_{2}\leq\cos\kappa_{1}\leq x_{3}$, the above form is well defined.

Finally by integrating both sides of Eq.(30), Keplerian parametrization of $\cos\kappa_{1}$ can be obtained as,

[TABLE]

where $\Upsilon=\frac{\sqrt{A(x_{1}-x_{2})}}{2}(\alpha+\frac{\nu+e\sin\nu}{c^{2}\,L^{3}})$, $\beta=\sqrt{\frac{x_{3}-x_{2}}{x_{1}-x_{2}}}$ and $\text{sn},\text{cn}$ and dn are the sine, cosine and tangent like Jacobi elliptic functions Whittaker and Watson (2009).

It is worth pointing out that $\pm$ sign which was supposed to be inside the definition of $\Upsilon$ according to Eq.(III.1), is removed. Because of symmetric property, $\text{sn}^{2}(X,b)=\text{sn}^{2}(-X,b)$, we have a freedom to choose definition of $\Upsilon$ among four cases of $\pm\frac{\sqrt{A(x_{1}-x_{2})}}{2}(\alpha\pm\frac{\nu+e\sin\nu}{c^{2}\,L^{3}})$. Firstly, we have chosen $\frac{\sqrt{A(x_{1}-x_{2})}}{2}(\pm\alpha+\frac{\nu+e\sin\nu}{c^{2}\,L^{3}})$ so that the sign of $\nu$ dependent term is positive. And we make $\pm$ sign in front of $\alpha$ be absorbed into the definition of $\alpha$. That is,

[TABLE]

where $\cos\kappa_{1(0)}$ is an initial value of $\cos\kappa_{1}$ when $t=t_{0}$, or $u=0$. The $\pm$ sign is determined as : $+$ for the case that $\frac{d\cos\kappa_{1}}{dt}>0$ initially while $-$ for $\frac{d\cos\kappa_{1}}{dt}<0$ initially. Lastly, as we have seen, the expression of Eq.(III.1) could not capture its oscillating feature. But we do not need to be concerned of it since the final form of our solution Eq.(38) trivially captures the feature. Taking Eq.(32) as an example again, the reflection symmetry of sn, $\text{sn}(X+b,d)=\text{sn}(X-b,d)$ when sn has its maximum or minimum at $X$, allows (here we intentionally omit $\frac{\sqrt{A(x_{1}-x_{2})}}{2}$ and $\beta$)

[TABLE]

so that our solution is justified to be valid independent of duration.

In the case of $A<0$ i.e. $m_{1}<m_{2}$, the form of our solution is still valid with rearranging the order of the roots as $x_{1}<x_{3}\leq x_{2}$. Also in the case of no dynamics i.e. $\frac{d\cos\kappa_{1}}{dt}=0$ all the time, our solution is valid with $x_{2}=x_{3}$. In the case of $A=0$ i.e. $m_{1}=m_{2}$, our solution is still valid with $x_{1}\rightarrow\infty$ and thus $\beta=0$.

Additionally, because the sine-like Jacobian elliptic function $\text{sn}(\Upsilon,\beta)$ is periodic, we can get the period of the precessional motion. Let $N$ be the number of orbital cycles. During $N$ cycles, $\Upsilon$ increases as

[TABLE]

while $\text{sn}^{2}(\Upsilon,\beta)$ has a following period Whittaker and Watson (2009),

[TABLE]

So $\Delta\Upsilon_{P}=\Delta\Upsilon_{N}$ yields the number of cycles $N=N_{\text{prec}}$ of the binary system during one period of precession,

[TABLE]

where $\mathcal{F}_{2,1}$ is a hypergeometric function.

Furthermore, invoking the definitions of $\sigma_{1}$ and $\sigma_{2}$, we can get the relationship between $\gamma$, $\kappa_{1}$ and $\kappa_{2}$

[TABLE]

Now, in the non-inertial frame where the orbital angular momentum is fixed as $\bm{L}=(0,0,L)$, we have a parametrization on how the spin angular momentums evolve around $\bm{L}$.

In the non inertial frame, each angular momentum is expressed as,

[TABLE]

where $\Lambda^{-1,T}$ is the inverse matrix of the transpose of $\Lambda$ and the subscripts $n$ and $i$ refer to non-inertial and inertial frame respectively, where components of vector are described. Since $\bm{J}=\bm{L}+\bm{S}_{1}+\bm{S}_{2}$ , we get the following three relations

[TABLE]

Fom Eq.(46c), we get $\cos\xi_{2}$ in terms of $\cos\kappa_{1}$,

[TABLE]

On the other hands, from Eq.(46a) and Eq.(46b), we get the expression for $\xi_{3}$

[TABLE]

with Eq.(23),

[TABLE]

This completes parametrization of all relative angles in Kepelerian type.

III.2 The Keplerian-type Parametrization of the Orbital Angular Momentum

Now, we are in the position to solve a Keplerian parametrization of $\xi_{1}$. The easiest way is to get components of $\frac{d\bm{L}}{dt}$ in the non inertial frame, and solve Eq.(6a), also in the non inertial frame.

Starting from the inertial-expression of $\frac{d\bm{L}}{dt}$, we can get components of $\frac{d\bm{L}}{dt}$ in the non-inertial frame ,

[TABLE]

From the componential expression of Eq.(6a) in the non-inertial frame, we get

[TABLE]

The first row in the equation can be written in terms of $\cos\kappa_{1}$ using the relations Eqs.(44) and (46),

[TABLE]

where $\alpha_{1}$, $\alpha_{2}$, $\beta_{1}$ and $\beta_{2}$ are defined as follows

[TABLE]

The integration can be written as, again with $x=\cos\kappa_{1}$,

[TABLE]

At this point, we encounter similar issue discussed around Eq.(III.1). The sign of $\frac{dx}{dt}$ oscillates during evolution and this makes the integration messy. To resolve this, we use $y$ introduced in Eq.(34) again but imposing that $y$ increases as $t$ increases instead of restricting the range of $y$. (In fact, by comparing with Eq.(38), $y=\text{am}(\Upsilon,\beta)$.) Then, the integration can be written without the $\pm$ sign as

[TABLE]

The above integration is identified as the Elliptic integral of the third kind Whittaker and Watson (2009), denoted as $\Pi$,

[TABLE]

Therefore

[TABLE]

where $\xi_{1(0)}$ is an initially determined parameter and am denotes the Jacobi amplitude.

Up to now Keplerian parametrization for evolution of angular momentum has been derived in the inertial frame. Especially, the orbital angular momentum $\bm{L}$ dynamics is totally determined by Eqs.(75) and (47) up to leading order precession. The last step is a parametrization of relative motion of binaries up to the leading order of spin-orbit interaction.

III.3 Keplerian-type Parametrization of the relative motion of binaries

From the Hamiltonian Eq.(1), the equation governing $\bm{r}$ is

[TABLE]

or,

[TABLE]

As can be seen from Fig. 1, $\bm{r}$ can be written,

[TABLE]

To get an explicit expression for $\dot{\bm{r}}$, we must take time derivative of $\bm{r}$ in the inertial frame, and bring it back to the non-inertial frame using the Euler matrix,

[TABLE]

Now we have all ingredients (Eqs.(45) and (85)) to express the vector identity $\bm{L}=\bm{r}\times\bm{p}$ in non-inertial frame components,

[TABLE]

where

[TABLE]

Note that $D$ trivially vanishes as Eq.(92) requires as Eq.(69) is satisfied. The third component of $\bm{L}$ reads to the following equation

[TABLE]

For perturbatively consistent integration of $\frac{d\varphi}{dt}$, we need 1.5PN accurate $r$ expression. It is because $\frac{d\varphi}{dt}$ depends on $r$ from Newtonian order while we want to get $\varphi$ with an accuracy of 1.5PN. For a moment, we need to find the Keplerian parametrization of $r$. Using $\bm{p}^{2}=\dot{r}^{2}+\frac{L^{2}}{r^{2}}$ and $H=E$,

[TABLE]

Note that $\bm{L}\cdot\bm{S}_{\text{eff}}$ is a constant with the following simplification

[TABLE]

Eq.(95) is easily integrated Damour and Deruelle (1985) and the radial motion is described, up to leading order of spin-orbit coupling, as

[TABLE]

[TABLE]

Now let us go back to Eq.(94). We’ve already known that the integration of the second term of right hand side of Eq.(94) is the Elliptic integral of the thrid kind and the integration of the first term is easily done via Eq.(97),

[TABLE]

and $\varphi_{0}$ is an initially given value.

In principle, we have all parametrizations for describing the Newtonian orbital motion with the leading order of spin-orbit coupling. In the inertial frame, it is written as

[TABLE]

IV Almost equal mass Approximation

In this section, we restrict our parametrizations in the case of $\delta_{2}-\delta_{1}\sim 0$, and provide the simpler version trucated at the order of $O((\delta_{2}-\delta_{1})^{2})$. Such a limiting case was considered by Tessmer (2009) but the result was not provided in a closed form. Here we provide it only by elementary functions. In the exactly equal mass case, the solution was derived by Königsdörffer and Gopakumar (2005) analytically so that we have a chance to compare directly.

For the purpose of expansion for nearly same mass binaries, we choose a dimensionless parameter $\delta$,

[TABLE]

which is assumed to be non-negative (i.e. $m_{1}\geq m_{2}$). The targeted approximant is achieved straightforwardly by changing of variable $\delta_{2}\rightarrow\delta+\delta_{1}$ and expanding it in terms of $\delta$ up to the first order. We should be careful about the fact that, since $A\,(x-x_{1})(x-x_{2})(x-x_{3})$ becomes the 2nd order polynomial in the case of equal mass, $x_{1}$ should blow up as $\delta\rightarrow 0$ and indeed by concrete calculation, it turns out that $x_{1}=\frac{1}{\delta}\,\frac{\delta_{1}\,\left(S_{1}^{2}+2\sigma_{1}\,S_{1}\,S_{2}+S_{2}^{2}\right)}{2\,L\,S_{1}}+O(\delta^{0})$. And we further find that

[TABLE]

In order to get the Keplerian parametrizations for $\xi_{1}$ and $\varphi$ expanded up to the first order in $\delta$, it is essential to get the first order approximant of the type $\Pi(\mathcal{A},\,\text{am}(\Upsilon,\beta),\,\beta)$. With $\mathcal{A}\sim O(\delta^{1})$, $\beta\sim O(\delta^{1})$ and $\Upsilon\sim O(\delta^{0})$, we have

[TABLE]

Then it staightforwardly follows that

[TABLE]

where

[TABLE]

and $\hat{\Upsilon}$ is the truncated one,

[TABLE]

Similarly, Eq.(98a) is expanded to

[TABLE]

where

[TABLE]

with $\hat{\nu}=2\,\arctan\left(\sqrt{\frac{1+\hat{e}}{1-\hat{e}}}\,\tan\frac{u}{2}\right)$. Note that $\hat{k}$, $\hat{e}$ are just approximants of $k$ and $e_{\varphi}$ in Eqs.(98) truncated at the order of $O(\delta^{2})$.

IV.1 Equal mass binaries limit

From this we can easily compare the leading terms with the results of Königsdörffer and Gopakumar (2005), which deals with equal mass case. By setting $\delta=0$, the expression for $\xi_{1}$ becomes,

[TABLE]

Similarly,

[TABLE]

Since the followings hold perturbatively,

[TABLE]

and

[TABLE]

we can make some part of $\varphi$ more compact by introducing

[TABLE]

Now we arrive the compact form in the accuracy of $O(\frac{1}{c^{4}})$ such that

[TABLE]

where

[TABLE]

These results are coincident with the results of Königsdörffer and Gopakumar (2005).

V Conclusion

We have solved the spin precession equation and the Newtonian orbital motion in ADM gauge for binary system with arbitrary spins, mass ratio and eccentricity up to leading order of spin orbit coupling. We arrived at quasi-Keplerian parametric solutions in a simple closed form. The elliptic functions are essential in our parametrizations and they have been thoroughly speculated by mathematicians. So our solutions are expected to give us systematic and mathematically deep understandings on how rotating bodies move in gravitational field. On the other hand, numerical computations for the elliptic functions are rather expensive compared to the elementary functions. However, our rough numercial estimations suggest that a couple of the constants defined in this paper such as $\beta$, $\frac{x_{2}-x_{3}}{\alpha_{1}+x_{2}}$ tend to be very small in almost all initial configurations, implying that we would be able to re-express the elliptic functions in terms of the elementary functions without significant numercial errors. This could be our future work.

We also expect that the solutions presented here can be directly used to get efficient ready-to-use time domain gravitational waveform templates modulated by spin precession and swings of orbital plane, which would be useful in analyzing the gravitational wave data. Furthermore, the result of this paper will provide useful knowledges to complete phenomenological models such as effective-one-body which requires both numerical calibrations and analytic solutions.

Acknowledgement

Gihyuk Cho acknowledges financial support through BK21 plus program of the Ministry of Education, Government of korea.