Analytic Solutions to Compact Binary Inspirals With Leading Order Spin-Orbit Contribution Using The Dynamical Renormalization Group

TL;DR

This paper derives exact analytic solutions for the motion and spin precession of spinning compact binaries during inspiral, using the dynamical renormalization group to improve accuracy and computational efficiency over traditional methods.

Contribution

It introduces a novel application of the dynamical renormalization group to obtain closed-form solutions for binary inspiral dynamics with spin-orbit effects at leading order.

Findings

Resummed solutions outperform adiabatic approximations in accuracy.

Solutions are computationally faster than numerical integration.

Applicable to arbitrary mass and spin configurations.

Abstract

We calculate the real-space trajectory and spin precession of a generic spinning compact binary inspiral at any time instant using the dynamical renormalization group formalism. This method leads to closed-form analytic solutions to the binary motion through treating radiation reaction as perturbations and resumming the secular growth of perturbative terms. We consider the spin-orbit effects at leading order and the 2.5PN radiation reaction without orbit averaging or precession averaging for arbitrary individual masses and spin magnitudes and orientations. The solutions are written in a moving reference frame, with the orbital angular momentum and binary radial directions aligned along two of the axes. The resummed solutions show improved accuracy compared to adiabatic solutions while also being an order of magnitude faster computationally compared to numerical integration methods.

Equations

$$\mathbf{a}={\mathbf{a}}_{\text{N}}+{\mathbf{a}}_{\text{SO}}+{\mathbf{a}}_{\text{RR}},$$

$${\mathbf{a}}_{\text{N}}$$

$${\mathbf{a}}_{\text{SO}}$$

$${\mathbf{a}}_{\text{RR}}$$

$$\mathbf{S}={\mathbf{S}}_1+{\mathbf{S}}_2,\mathbf{\Sigma }=\frac{{\mathbf{S}}_2}{X_2}-\frac{{\mathbf{S}}_1}{X_1},$$

$${\dot{\mathbf{S}}}_a$$

$$\mathbf{v}=\dot{r}\mathbf{n}+r\omega \mathbf{\lambda }$$

$$\mathbf{a}=(\ddot{r}-r{\omega }^2)\mathbf{n}+(r\dot{\omega }+2\dot{r}\omega )\mathbf{\lambda }+r\varpi \omega \mathbf{l},$$

$$\ddot{r}-r{\omega }^2=$$

$$r\dot{\omega }+2\dot{r}\omega =$$

$$\varpi =$$

$$\frac{\mathrm{d}{S}_n^a}{\mathrm{d}t}$$

$$\frac{\mathrm{d}{S}_{\lambda }^a}{\mathrm{d}t}$$

$$\frac{\mathrm{d}{S}_l^a}{\mathrm{d}t}$$

$$\displaystyle\Omega_{a}\equiv\frac{\nu M\omega}{r}\Bigg{(}2+\frac{3}{2}\frac{m_{b}}{m_{a}}\Bigg{)},$$

$$\frac{\mathrm{d}{S}_n^a}{\mathrm{d}t}$$

$$\frac{\mathrm{d}{S}_{\lambda }^a}{\mathrm{d}t}$$

$$\frac{\mathrm{d}{S}_l^a}{\mathrm{d}t}$$

$$\Omega^{2}_{B}=\frac{M}{R^{3}_{B}}-\frac{\Omega_{B}}{R_{B}^{3}}\big{(}5S_{l}+3\Delta\Sigma_{l}\big{)},$$

$$r(t)=R_B+\delta r(t)+\delta r_S(t),\omega (t)={\Omega }_B+\delta \omega (t)+\delta {\omega }_S(t),$$

$$\delta r\sim v^5{R}_B,\delta \omega \sim v^6/R_B,\delta r_S\sim S{v}^4/R_B,\delta {\omega }_S\sim S{v}^5/R_B^3$$

$$\delta r_S(t)=$$

$$\displaystyle\qquad\qquad\qquad-\sin{(\Omega_{B}(t-t_{0})+\Phi_{B})}\Big{]}+A^{S}_{B}\cos{(\Omega_{B}(t-t_{0})+\Phi_{B})},$$

$$\delta {\omega }_S(t)=$$

$$\displaystyle-\big{(}14S_{l}+6\Delta\Sigma_{l}\big{)}\frac{A_{B}}{R_{B}^{4}}\Omega_{B}(t-t_{0})\cos{(\Omega_{B}(t-t_{0})+\Phi_{B})}-\frac{2A^{S}_{B}\Omega_{B}}{R_{B}}\cos{(\Omega_{B}(t-t_{0})+\Phi_{B})}.$$

$$\delta {\Phi }_S(t)=$$

$$\displaystyle-\Big{(}14S_{l}+6\Delta\Sigma_{l}\Big{)}\frac{A_{B}}{R_{B}^{4}}(t-t_{0})\sin{(\Omega_{B}(t-t_{0})+\Phi_{B})}-\frac{2A^{S}_{B}}{R_{B}}\sin{(\Omega_{B}(t-t_{0})+\Phi_{B})},$$

$$R_B(t_0)$$

$${\Omega }_B(t_0)$$

$${\Phi }_B(t_0)$$

$$A_B^S(t_0)$$

$$r(t)=$$

$$\omega (t)=$$

$$\varphi (t)=$$

$$\displaystyle\frac{64\nu M^{3}}{5}t+\frac{1}{4}R_{R}(t)^{4}+\frac{2\mathscr{S}}{5M^{1/2}}R_{R}(t)^{5/2}+\frac{\mathscr{S}^{2}}{M}R_{R}(t)+\frac{2\mathscr{S}^{8/3}}{\sqrt{3}M^{4/3}}\tan^{-1}\Bigg{(}\frac{1}{\sqrt{3}}+\frac{2M^{1/6}R_{R}(t)^{1/2}}{\sqrt{3}\mathscr{S}^{1/3}}\Bigg{)}$$

$$\displaystyle+\frac{\mathscr{S}^{8/3}}{3M^{4/3}}\ln\Bigg{(}\frac{\Big{(}\mathscr{S}^{1/3}-M^{1/6}R_{R}(t)^{1/2}\Big{)}^{2}}{\mathscr{S}^{2/3}+\mathscr{S}^{1/3}M^{1/6}R_{R}(t)^{1/2}+M^{1/3}R_{R}(t)}\Bigg{)}=\textrm{constant}$$

$$\displaystyle\Omega_{R}^{2}(t)R_{R}^{3}(t)+\Omega_{R}(t)\big{(}5S_{l}+3\Delta\Sigma_{l}\big{)}=M$$

$$\displaystyle\Phi_{R}(t)+\frac{1}{32M^{5/2}\nu}R_{R}^{5/2}(t)-\frac{5\Big{(}41S_{l}+15\Delta\Sigma_{l}\Big{)}}{256\nu M^{2}\mathscr{S}^{2}}\Bigg{(}\frac{64\nu M^{3}}{5}t+\frac{1}{4}R_{R}^{4}(t)+\frac{2\mathscr{S}}{5M^{1/2}}R_{R}^{5/2}(t)\Bigg{)}=\textrm{constant}$$

$$A_R(t)=\text{constant}$$

$$\displaystyle A^{S}_{R}(t)-\frac{5A_{R}\Big{(}7S_{l}+3\Delta\Sigma_{l}\Big{)}}{64\nu M^{2}\mathscr{S}^{2}}\Bigg{(}\frac{64\nu M^{3}}{5}t+\frac{1}{4}R_{R}^{4}(t)+\frac{2\mathscr{S}}{5M^{1/2}}R_{R}^{5/2}(t)\Bigg{)}=\textrm{constant}$$

$$\displaystyle S^{a}_{+}(t)=\mathcal{S}^{a}_{+R}(t)\exp\Bigg{\{}$$

$$\displaystyle-i\Bigg{[}\frac{2A_{R}(t)}{R_{R}(t)}-\Big{(}19S_{l}+9\Delta\Sigma_{l}\Big{)}\frac{A_{R}(t)}{R_{R}(t)^{4}\Omega_{R}(t)}-\frac{3\nu_{a}}{\nu}\Omega_{R}(t)^{2}A_{R}(t)R_{R}(t)$$

Full text

I Introduction

A worldwide network of gravitational wave (GW) detectors is being developed to monitor the ripples in the fabric of spacetime passing through the earth. This includes the ground-based laser interferometers GEO600, LIGO and VIRGO collaborations currently in operation, and the under-construction spaced-based observatory LISA and cryogenic detector KAGRA in Japan. The successful detection of gravitational waves from inspiraling black holes (BH) and neutron stars (NS) by the LIGO and VIRGO collaborations directly and spectacularly confirmed one of the predictions of Einstein’s Theory of General Relativity. A generic prediction of metric theories of gravity, BH or NS coalescence is a strong source of GWs for interferometric detectors. To successfully identify and analyze the gravitational wave signals, it is necessary to construct a systematic description of the binary black hole dynamics and waveforms during coalescence. A set of expected waveforms portrayed by the intrinsic parameters of the compact binary within the astrophysically interesting region of the parameter space forms a waveform template bank Babak et al. (2006); Ajith et al. (2008); Harry et al. (2016). Using these precise waveform templates, a matched filtering technique is used to try to discover the weak GW signals buried in the detector noise. More accurate templates will help us extract physical information from the observed events in order to gain further knowledge of the black hole or neutron star properties.

The choice of the BH coalescence model is critical for determining the waveform. The last few orbits of the inspiral phase through the merger and ringdown of the BH coalescence have been simulated by Numerical Relativity Baker et al. (2007); Centrella et al. (2010). There have been developments on the analytic understandings for merger and ringdown McWilliams (2019). The slowly-orbiting long inspiral phase can be studied analytically using post-Newtonian (PN) perturbation theory with small velocity and weak field approximations. BH dynamics is described by the Newtonian-like equations of motion in the form of the acceleration of the binary constituents. During the inspiral, the binary slowly loses energy and angular momentum to gravitational radiation starting at 2.5PN Burke and Thorne (1970); Thorne and Hartle (1985). Higher order corrections up to 4PN in the conservative sector have been calculated Damour et al. (2014, 2016); Jaranowski and Schafer (2015); Bernard et al. (2016, 2017); Marchand et al. (2018); Bernard et al. (2018); Foffa and Sturani (2019); Foffa et al. (2019). Solving for the motions is the fundamental step in obtaining the waveforms and deriving the evolution of the theoretical physical measurements in time, such as the GW phase directly measured by the detectors and power loss due to gravitational radiations.

The exact solutions to the motions can be found by numerically integrating these nonlinearly-coupled ordinary differential equations. However, in calculations of template banks, each point in the intrinsic parameter space representing a waveform with different initial conditions requires a new numerical computation. The sample rate of the corresponding waveform directly depends on the precision and step sizes of the solutions of the motions. The discrete nature of the computational solutions also brings the issue of the distance between the templates in the parameter space, which may result in the loss in signal-to-noise ratio due to the mismatch of the template in the match-filter of the signal data. Since the third observing run of LIGO and Virgo are having a weekly rate of observed events, a faster and more accurate way of computation in the signal analysis is critical, with even larger rates expected with future upgrades. A fully analytic waveform solution with continuous parameters would certainly increase calculation efficiency for template-based data analysis.

The adiabatic approximation is often used to find the analytic solutions to the motion, including inspiral radiation reaction effects Drasco et al. (2005); Hughes et al. (2005); Sago et al. (2005). Using the PN expansions of the conserved energy $E$ and flux $\mathcal{F}$, the adiabatic waveforms are obtained by solving the energy balance equation $\mathrm{d}E/\mathrm{d}t=\mathcal{F}$. The balance equation leads to the secular evolutions of the orbital angular frequency $\omega(t)$, from which one can derive the accumulated phase of gravitational waves $\phi(t)=2\int\mathrm{d}\tau\omega(\tau)$. An implicit assumption in the energy balance equation is that $E$ does not change much over an orbital timescale. In other words, the adiabatic solutions are orbit-averaged and thus remove some of the orbital detail. The adiabatic approximation fails to account for secular evolution of some of the orbital elements, which can lead to measurable phasing effects Pound et al. (2005).

When considering spinning black holes, which adds 6 additional degrees of freedom, the binary motions become more complicated. The convention of the PN order counting of the spin here is defined as $|\bm{S}|=\chi m^{2}$, where $m$ is the mass of the object and $\chi$ a dimensionless spin parameter. For a maximally rotating compact BHs, $\chi\sim 1$. The leading contributions from spin-orbit effects enter into the motion at 1.5PN and spin-spin at 2PN, before the leading order radiation reaction force. The major effect of the presence of the spin on the orbital evolution is that a spin component perpendicular to the orbital angular momentum causes the orbital plane to precess. This means the orbital plane will change its orientation when it is not perpendicular to the spin vector. Thus the observed waveform, depending on the orbital orientation with respect to the detector, will modulate due to spin-induced orbital precession. The secular evolutions of the spins themselves are given by the spin precession equations Barker and O’Connell (1979); Thorne and Hartle (1984). With the spin precession equations, it is possible to determine the angular momentum transfer between orbital and spin angular momenta and the total angular momentum loss during the inspiral regime. One of the recent works to construct analytic spin-precessing inspirals is through multiple scale analysis Chatziioannou et al. (2017a, b). This method gives orbit-averaged and precession-averaged closed form solutions by making a clean separation among the orbital time, precession time, and radiation reaction time scales and treating the physical parameters by averaging over the longer time scales to solve for the shorter ones. However, any averaging procedure results in the loss of some of the orbital dynamics.

In order to find analytic solutions to the spinning binary equations of motion and spin precession equations without any averaging procedures, we follow the Dynamical Renormalization Group (DRG) formalism proposed by Galley and Rothstein in Galley and Rothstein (2017). The idea of the DRG method is based on renormalization group theory and the resummation of the singularities for perturbative ordinary differential equation problems Chen et al. (1996). The DRG method applied to binary inspirals starts by treating some of the higher PN order radiation reaction terms as a perturbation to a conservative background orbit. The secular growths of the perturbations are then resummed to preserve the correct power counting of the perturbations. In their work, Galley and Rothstein calculated the resummed solution for a non-spinning binary with leading order radiation up to the second-order corrections and included the PN corrections to the radiation reaction force. In this paper we incorporate spin-orbit effects and the leading-order radiation reaction, using the DRG method to obtain real-time solutions to the generic precessing compact binaries.

The organization of this paper is as follows: In Section II we introduce the PN equations of motion and spin precession equations for compact binary inspirals. We also set up a moving coordinate frame using the radial vector and orbital angular momentum vector in which we will present our solutions. In Section III we summarize the procedures of the DRG method and give the resulting closed-form analytic solutions to the binary motions and spin precession. In Section IV, we compare our DRG resummed solutions to the numerical and adiabatic solutions of the same equations. We also show a rough comparison of the calculation run time between the numerical integration and resummed solution substitution. We conclude in Section V. In Appendix A and B, we present the detailed calculations of the DRG method for orbital motions and spin precession, respectively. In Appendix C we propose a naive transformation of the moving coordinate frame to a fixed observer frame for the purpose of waveform construction.

II Leading Order Spin-Orbit Equations of Motion And Spin Precession Equations

The equations of motion of the compact binaries in the center-of-mass frame, including the Newtonian order, the leading-order spin-orbit contributions at 1.5PN in covariant spin supplementary condition (SSC), and the Burke-Thorne term due to the radiation-reaction force at 2.5PN, are given by Kidder (1995); Burke and Thorne (1970); Blanchet (2014)

[TABLE]

where the terms in the post-Newtonian hierarchy are

[TABLE]

In the expressions above, $\bm{r}$ and $\bm{v}$ are the binary relative center-of-mass separation and velocity, $\hat{\bm{n}}\equiv\bm{r}/r$ and $\dot{r}=\mathrm{d}r/\mathrm{d}t=\hat{\bm{n}}\cdot\bm{v}$. The binary masses are denoted as $m_{1,2}$, the total binary mass $M=m_{1}+m_{2}$, $\nu\equiv m_{1}m_{2}/M^{2}$ and $\Delta\equiv(m_{1}-m_{2})/M$. The combinations of the individual spins are written as

[TABLE]

with $X_{a}=m_{a}/M$. The spin vectors precess due to spin-orbit coupling following the relation of Thorne and Hartle (1985); Kidder (1995)

[TABLE]

where $\{a,b\}$ are the binary labels $\{1,2\}$, and $\bm{L}_{N}=\nu M(\bm{r}\times\bm{v})$ is the Newtonian orbital angular momentum.

In order to obtain the analytic solutions to the inspiral equations of motions and the spin precession equations (1)-(4), we adopt a coordinate frame $\{\bm{\bm{n},\bm{\lambda},\bm{l}}\}$, moving along with the center-of-mass and the orientation of its motion Blanchet (2014); Arun et al. (2009); Marsat et al. (2013), where $\bm{l}=\bm{n}\times\bm{v}/|\bm{n}\times\bm{v}|$ and $\bm{\lambda}=\bm{l}\times\bm{n}$ to complete an orthonormal basis triad. In this moving basis, the relative velocity can be expressed as

[TABLE]

where $\omega$ is the orbital angular frequency of the binary. The relative acceleration $\bm{a}=\mathrm{d}\bm{v}/\mathrm{d}t$ in the moving basis is

[TABLE]

where the orbital plane precession $\varpi$ of the orbit is defined as $\varpi\equiv-\bm{\lambda}\cdot\mathrm{d}\bm{l}/\mathrm{d}t$.

In terms of the moving basis components, the equations of motions (1) are

[TABLE]

where we decompose the spin $\bm{S}=S_{n}\bm{n}+S_{\lambda}\bm{\lambda}+S_{l}\bm{l}$, and similarly for $\bm{\Sigma}$. The spin precession equations (4) become

[TABLE]

where we denote

[TABLE]

which is the norm of the precession vector of the $a$-th spin. The precession frequency $\varpi$, explicitly given by (7c), is of order $\mathcal{O}(S)$. At linear order in spin, the precession equations become

[TABLE]

Thus at order $\mathcal{O}(S)$, the $l$-component of the spin vectors are invariant, which are also the only components that appears in the orbital equations of motion in (7a) and (7b). In solving these two equations by the DRG method, we then are able to treat $S_{l}$ and $\Sigma_{l}$ as time-independent constants. Following Ref Galley and Rothstein (2017), here we ignore the 1PN and 2PN conservative forces, as well as the next-to-leading order spin-orbit effects, which is the same order in the Post-Newtonian expansion as the 2.5PN radiation reaction terms. Instead, we focus on the leading order radiation reaction effects on spinning objects. In order to obtain gravitational wave templates, to be consistent we would need to include at least the 1PN conservative forces.

III DRG Solutions to Dynamics and Spin Precession

The background quasi-circular orbit of a conserved binary with Newtonian and leading spin-orbit effects can be described by

[TABLE]

with constant radius $R_{B}$ and constant angular frequency $\Omega_{B}$. To include the radiation reaction as perturbative effects, we write the orbital solutions as

[TABLE]

where the first time-dependent terms $\delta r(t)$ and $\delta\omega(t)$ are the perturbation due to the 2.5PN radiation reaction force without the spin at a given time $t$. The $\delta r_{S}(t)$ and $\delta\omega_{S}(t)$ represent the perturbations due to the interaction of 1.5PN spin effects and the 2.5PN radiation reaction. The power counting at the initial time $t_{0}$ for each perturbation is given by

[TABLE]

where we keep the spin as a placeholder expansion parameter instead of converting to PN orders for generality. Substituting the perturbed orbital radius and frequency into the equations of motion (7a) and (7b), we find the solutions to the perturbation

[TABLE]

where $\Phi_{B}$, $A_{B}$, and $A_{B}^{S}$ are integration constants. $\{R_{B},\Omega_{B},\Phi_{B},A_{B},A_{B}^{S}\}$ forms a set of bare parameters to be determined by initial conditions. While $e_{B}=A_{B}/R_{B}$ is the small orbital eccentricity of order $\mathcal{O}(v^{5})$ induced by the radiation reaction force, the interaction between the spin and radiation reaction lead to a smaller eccentricity $e^{S}_{B}=A^{S}_{B}/R_{B}\sim\mathcal{O}(Sv^{4})$. The spin-radiation eccentricity deforms the circular orbit out-of-phase compared to the radiation eccentricity, although with a fixed phase difference.

To maintain the power countings of the perturbations, the secularly growing terms in (14) are absorbed into the bare parameters through the relations

[TABLE]

where $\{R_{R},\Omega_{R},\Phi_{R},A_{R}^{S}\}$ are the “renormalized” parameters depending on an arbitrary renormalization scale $\tau$. The quantities $\{\delta_{R}^{v^{5}},\delta_{R}^{S},\delta_{\Omega}^{v^{5}},\delta_{\Omega}^{S},...\}$ are counter-terms, to be determined by renormalizing the perturbation expansions. Introducing the renormalization scale into the perturbation solutions (14) by writing $t-t_{0}=(t-\tau)+(\tau-t_{0})$ and using the counter-terms to cancel all the secular $(\tau-t_{0})$ terms, we find that

[TABLE]

where $r(t)$ and $\omega(t)$ are the orbital radius and frequency defined in the previous section, and $\phi(t)$ is the time integral of $\omega(t)$ representing the orbital phase of the binary inspiral. The renormalized parameters are determined at arbitrary time via the renomalization group equations, determined using the fact that the corresponding bare parameters are independent of the choice of $\tau$. The “beta-functions” of the RG equations are determined by the counter-terms, leading to the first-order equations satisfied by the renormalized parameters. We give the RG solutions in the form of invariance in time as

[TABLE]

where for convenience we have defined $\mathscr{S}\equiv(51+21\Delta\Sigma_{l})/4$. Remember, at this order the $l$-component of the spin vectors are constant. The constants in the equations above can be further determined using the initial conditions by solving (16) at a given time instant. The expressions in (16) combined with (17) give the resummed solution to the 0PN spinning inspiral dynamics valid up to times $(t-t_{0})$ of order $1/(\nu v^{5}(t)\Omega_{R}(t))$. To improve the accuracy, we need to calculate higher order perturbations in the same formalism or include higher PN conservative corrections to the motions.

The background solution to the conserved spin precessions has a constant precession frequency. We renormalize the precession frequency perturbed by the radiation reaction using the same DRG procedure. The resummed solutions to the spin precession equations (4) are

[TABLE]

where $\nu_{a}\equiv(2+\frac{3m_{b}}{2m_{a}})\nu^{2}$ and $S^{a}_{+}\equiv S^{a}_{n}+iS^{a}_{\lambda}$ contains the two precessing components of the spin vector in the moving triad $\{\bm{n},\bm{\lambda},\bm{l}\}$ coordinate system. The exponential preserves the magnitudes of the spin vectors, which is conserved as can easily be seen from Eq. (4). The renormalized parameter $\mathcal{S}^{a}_{+R}(t)$ can be written in terms of invariance over time and other parameters as

[TABLE]

We include the more detailed calculations and renormalization procedures in the appendices for interested readers.

IV Numerical Solution Comparison

To compare our analytic solutions to the orbital equations of motion and the spin precession equations, we solve the sets of equations numerically and compare with the DRG solution solved with the same initial conditions. We choose to compare compact binary systems of total mass $M=1$. The initial conditions for the physical parameters are related to the renormalized parameters through the renormalized solutions (42) and (B.1) at $t_{i}=0$. We choose for our initial conditions

[TABLE]

where the expression for $A_{R}(0)$ comes from

[TABLE]

and $\dot{r}(0)$ is taken to be zero for quasi-circular motion. Meanwhile, we impose a small non-vanishing $\mathcal{O}(v^{5})$ eccentricity $e_{R}=A_{R}/R_{R}(t)$, and a spin-induced eccentricity $e^{S}_{R}=A^{S}_{R}(0)/R_{R}(0)$ at $\mathcal{O}(v^{4}S)$ that runs starting from zero.

For initial spin vectors we consider the compact components maximally rotating, meaning the dimensionless spin parameter $\chi\sim 1$ where for each spin $|\bm{S}_{a}|=\chi_{a}m_{a}^{2}$, with $\chi_{\textrm{max}}=1$ for black holes. In Fig. 1 we compare the resummed solutions to the orbital equations of motion with the numerical and adiabatic solutions Kidder et al. (1993) for two different choices of mass ratio and spins. For the left column, we choose an equal mass binary and anti-aligned spin initial configuration:

[TABLE]

The physical interpretation for the angle of $140^{\circ}$ and $50^{\circ}$ is the angle between the spin vectors and the orbital angular momentum $\bm{L}$. (At linear order in spin, equal mass systems satisfy the spin-orbit resonance orientations Schnittman (2004).) In the right column, we choose a moderate mass ratio ($m_{1}:m_{2}=4$), with a randomly chosen initial spin configuration:

[TABLE]

Specifically, the plots show the orbital radius $r(t)$ and orbital phase $\phi(t)$ for resummed, adiabatic, and numerical solutions to the binary equations of motion. Below each plot of the physical solutions are the fractional errors comparing the numerical results to resummed and adiabatic ones. From these plots, we can see the DRG methods are more accurate compared to the adiabatic solutions, with roughly an order of magnitude improvement in calculating the accumulated orbital phase over most of the inspiral.

We can see the importance of using the DRG method increases as we include higher-order corrections by comparing Fig. 1 to the results in Ref. Galley and Rothstein (2017). In that paper, the authors included the 0PN (i.e. Newtonian) contribution and the leading order radiation reaction term. As can be seen by looking at Fig. 1 of that paper, the DRG and adiabatic results give the same order relative errors.111Note that the authors of Ref. Galley and Rothstein (2017) show how to obtain the result including 1PN contributions, but did not provide any numerical results. They also did the “two-loop” contribution, which includes $\mathcal{O}(v^{10})$ corrections. Including these, the DRG method shows roughly an order of magnitude improvement compared to the adiabatic solution, as can be seen in Fig. 2 of that paper. When including the 1.5PN spin contribution as we did here, there is an order of magnitude improvement, as shown in Fig. 1.

We compare the resummed solutions of the spin precession equations with the numerical solutions to (10) in Fig. 2. The two columns have the same choices for mass ratio and spin configurations as in Fig. 1. In the top two panels, we plot the resummed solutions to the $\bm{n}$- and $\bm{\lambda}$-components, respectively, for the total spin vector (in blue) and the difference between the resummed and numerical solution (in red). We also include an inset plot of the spin precession for the last quarter of the inspiral to illustrate the phase difference. That the error accumulated from the resummed results of the spin precession becomes significant is the consequence of the Post-Newtonian method breaking down for large velocities during the later portion of the inspiral. We expect better accuracy when spin-spin effects and higher PN order terms are incorporated. In the third panel, we plot the angle between the spin vector results from the resummed and numerical solutions.

With the inclusion of radiation reaction, the total angular momentum changes direction and magnitude. In the bottom row of Fig. 2, we show the angular momenta changing throughout the inspiral. The equal mass binary shown in the left panel has a fixed total spin magnitude due to the symmetric form in (4). Both binaries exhibit a rapid loss of orbital and total angular momenta at the end of the inspiral in sync with the drop of the orbital radius in Fig 1.

In Fig. 3 we give a rough comparison of the computational runtime improvement of the DRG methods. The numerical solution for the equations of motion and spin precession was calculated in C++ implementing the ODEINT library Ahnert and Mulansky (2011). We adopt the Dormand-Prince algorithm at fifth order with adaptive step sizes and control the tolerance error to be consistent with the theoretical resummed solution errors. Fig. 3 shows the runtime of the numerical and DRG methods solving the same sets of initial conditions, changing the binary mass ratio Count times in each run. In order to try to have a meaningful comparison, we manipulate the average steps taken per run for the DRG methods to have similar output lengths (i.e., number of time steps for the solution) with the numerical integration. For example, in a total of 50 runs, the numerical integration takes 10 seconds and averagely 11235 steps per run, while the DRG method takes about 1 second and averagely 11438 steps per run. As can be seen, the DRG method is an order of magnitude faster than the numerical solution.

V Conclusion

Using the dynamical renormalization group formalism, we have solved the spinning binary dynamics including the 2.5PN radiation reaction and the leading order spin-orbit effects throughout the inspiral. The solution is obtained by the resummation of the secularly growing perturbations to the physical parameters including orbital radius, angular frequency, orbital phase, and spin precession phases. We solved the resummed solutions to the equations of motion and spin precession equations in a moving triad frame at any time instant. Renormalized parameters defined to describe the resummed solutions are determined using the renormalization group equations and can be written in terms of conserved identities.

The solutions are applicable to arbitrary initial configurations and do not dependent on any specific spin orientations. The comparison of numerical solutions and our analytic solutions shows greater accuracy than the adiabatic solutions and a sizable improvement in computation time compared to the numerical solutions. The use of the DRG method is more important for spinning BHs than for the non-spinning case. However, there is further improvements that can be made. The spin component comparison is not ideal, as shown in Fig. 2 with increasing phase differences. When initial spins are relatively large compare to orbital angular momentum, the discrepancy grows very fast in the early part of the inspiral. This is due to the beginning of the breakdown of the PN expansion. We hope to fix this issue and enhance the accuracy by the inclusion of spin-spin effects and higher-order PN terms into the formulation in the future works Yang .

VI Acknowledgments

We thank Ira Rothstein and Chad Galley for useful discussions. We also thank Shana Li for useful discussions. ZY and AKL were supported in part by the National Science Foundation under Grant No. PHY-1820760.

Appendix A Orbital Equations of Motion Resummation Solutions

We start by investigating the quasi-circular background orbit solution of the conservative spinning binary. In this case, the radius and the orbital angular frequency are constants apart from small non-secular perturbations induced by the presence of spins. The constant radius $R_{B}$ and orbital frequency $\Omega_{B}$ satisfy

[TABLE]

at linear order in spin. Setting the spins to zero, the relation between $R_{B}$ and $\Omega_{B}$ reduces to the usual Newtonian circular motion equation. To solve for $\Omega_{B}$ with a given $R_{B}$, we can either solve the quadratic equation above or substitute for $\Omega_{B}$ iteratively. The analytic solutions to the dynamics of quasi-circular conservative spinning binary systems have been studied Blanchet et al. (2011); Marsat et al. (2013).

A.1 Perturbations of quasi-circular orbits

Next we describe the deviation of the quasi-circular background orbit as a result of the leading order radiation reaction and linear spin-orbit effects by isolating the perturbative corrections $r(t)=R_{B}+\delta r(t)+\delta r_{S}(t)$ and $\omega(t)=\Omega_{B}+\delta\omega(t)+\delta\omega_{S}(t)$. The first time-dependent terms $\delta r(t)$ and $\delta\omega(t)$ are the perturbation that arise due to the 2.5PN radiation reaction force without the spin at a given time $t$, which are given by Galley and Rothstein (2017)

[TABLE]

with bare parameters $\{R_{B},\Omega_{B},A_{B},\Phi_{B}\}$, and $\delta r\sim v^{5}R_{B},\delta\omega\sim v^{6}/R_{B}$ at the initial time $t_{0}$. On the other hand, the terms due to the interaction of 1.5PN spin effects and the 2.5PN radiation reaction start with the power counting of $\delta r_{S}\sim Sv^{4}/R_{B}$ and $\delta\omega_{S}\sim Sv^{5}/R_{B}^{3}$. Expanding the equations of motion (7a) and (7b) to $\mathcal{O}(Sv^{6})$ gives

[TABLE]

with $\delta r(t)$ and $\delta\omega(t)$ the values given in (34a) and (34b). Integrating (35b) with respect to time, solving for $\delta\omega_{S}$, and substituting back into (35a) gives the differential equation for $\delta r_{S}$,

[TABLE]

The differential equation has a solution of the form

[TABLE]

As a result, the angular frequency perturbation $\delta\omega_{S}(t)$ and its time integration $\delta\Phi_{S}(t)$ are given by

[TABLE]

The perturbation $\delta\Phi_{S}(t)$, of order $\mathcal{O}(Sv^{4})$ to the angle $\phi(t)$, is the analog of the orbital phase in planar motion of non-spinning systems. Though it is no longer a physical angle now that the orbital plane precesses due to the spins, it is a combination of the Euler angles, defined in a later section, essential to the time evolution of the moving frame of reference.

We split the perturbation terms into the non-secular terms that remain small permanently, and the secular ones that grows with time. As time progresses, the secular terms gradually become dominant and break down the PN power counting, therefore they need to be resummed.

A.2 Renormalization

The full set of bare solutions to the orbit motion including linear spin-orbit terms and 2.5PN Burke-Thorne terms is given by

[TABLE]

with the corresponding perturbations in (34) and (37). We renormalize these terms by removing the $t_{0}$ dependence with the introduction of counter-terms for the bare parameters. The $\mathcal{O}(v^{5})$ terms were renormalized in Ref. Galley and Rothstein (2017). Thanks to the newly added $\mathcal{O}(Sv^{4})$ perturbations, the bare parameters have to include higher order counter-terms, which means

[TABLE]

In terms of the renormalized initial parameters and the renormalization scale $t-t_{0}=(t-\tau)+(\tau-t_{0})$, the spin-orbit result becomes

[TABLE]

By observation we can write down the counter-terms that cancel the $(\tau-t_{0})$ terms completely as

[TABLE]

Choosing the arbitrary renormalization scale to be $\tau=t_{0}$, the equations of motion is now described by the renormalized quantities $\{R_{R},\Omega_{R},\Phi_{R},A_{R},A_{R}^{S}\}$ as

[TABLE]

The explicit secular terms have been removed thanks to the choice of $\tau$, and the $t_{0}$-dependencies have been absorbed into the counter-terms. The runnings of $\{R_{R},\Omega_{R},\Phi_{R},A_{R},A_{R}^{S}\}$ and their dependence on the initial conditions are then determined by the renormalization group equations.

A.3 Renormalization Group Solutions

Exploiting the fact that the bare quantities $\{R_{B}(t_{0}),\Omega_{B}(t_{0}),\Phi_{B}(t_{0}),A^{S}_{B}(t_{0})\}$ are independent of the arbitrary scale $\tau$, we can write down the renormalization group equations for the renormalized quantities $\{R_{R}(t),\Omega_{R}(t),\Phi_{R}(t),A^{S}_{R}(t)\}$ as

[TABLE]

The right-hand sides of the RG equations, which are called beta functions, includes more iterative time-derivative terms that are of higher orders starting from $\mathcal{O}(v^{11})$ and $\mathcal{O}(S^{2})$. The RG solutions to $\Omega_{R}$, $\Phi_{R}$ and $A^{S}_{R}$ in terms of $R_{R}$ and the initial conditions are

[TABLE]

where $\mathscr{S}\equiv(51S_{l}+21\Delta\Sigma_{l})/4$ is a constant combination of the initial spins, defined for convenience. Substituting in to the radial RG equation, we find

[TABLE]

or, rearranging,

[TABLE]

Integrating both sides gives the exact but implicit relation,

[TABLE]

The parameter $A_{R}$ is unchanged when the spin is added, and from Galley and Rothstein (2017) we learned that $A_{R}$ has a zero $\beta$-function at the order we are working, i.e., $A_{R}$ is a constant, given by initial conditions, proportional to the initial eccentricity $e_{R}(0)=A_{R}(0)/R_{R}(0)\sim\mathcal{O}(v^{5}).$

Using the relation above for $R_{R}(t)$, we can further simply the expressions of $\Phi_{R}(t)$ and $A^{S}_{R}(t)$ in terms of $R_{R}(t)$ and time $t$, eliminating the logarithm and the arctangent terms. Written as an invariant in time, the renormalized quantities with the leading order spin-orbit effect are

[TABLE]

Note that one constraint appears in the RG equations of $R_{R}(t)$, (43a), which indicates the range of effectiveness of the DRG method,

[TABLE]

If $\mathscr{S}=(51S_{l}+21\Delta\Sigma_{l})/4$ is positive, $R_{R}(t)$, which is the dominant part of the binary center-of-mass separation $r(t)$, decreases until $R_{R}(t)=\mathscr{S}^{2/3}M^{-1/3}$. Given a limitation on the smallest value of $R_{R}(t)$ and combining with (A.3), it is possible to determine an approximate end time of the inspiral phase described by the Post-Newtonian equations of motion (1). This could provide useful information to numerical simulations as well.

Appendix B Spin Precession Equations

In this section, we aim to obtain the analytic solutions for the spin precession equations at linear order in spin (10) by applying DRG methods, with the quasi-circular solutions to the equations of motion from the previous section. For a conservative binary system moving in nearly circular motion, solving equations in the form of

[TABLE]

is fairly straightforward for constant radius $R$ and orbital frequency $\Omega$. The solutions are S^{a}_{n}=S^{a}_{\parallel}\sin{\Big{(}(\Omega-\Omega_{a})(t-t_{0})+\Phi\Big{)}} and S^{a}_{\lambda}=S^{a}_{\parallel}\cos{\Big{(}(\Omega-\Omega_{a})(t-t_{0})+\Phi\Big{)}}, where $S^{a}_{\parallel}$ is determined by the initial spin vectors.

With the inclusion of the radiation reaction force and the resulting time-dependence of $r(t)$ and $\omega(t)$, the spin vectors precess in a way entangled with the orbit motion. Defining $S^{a}_{+}\equiv S^{a}_{n}+iS^{a}_{\lambda}$, the precession equations (10) can be combined and written as

[TABLE]

A simple integration with respect to time leads to

[TABLE]

To solve for the integral on the right-hand side, we denote $\nu_{a}\equiv(2+\frac{3m_{b}}{2m_{a}})\nu^{2}$ and recall that $M\sim\Omega_{B}^{2}R_{B}^{3}+\Omega_{B}(5S_{l}+3\Delta\Sigma_{l})$, such that $\Omega_{a}$ in (9) at leading order in spin becomes

[TABLE]

with the 2.5PN radiation perturbation $\{\delta r,\delta\omega\}$ from (34), and the leading order spin-orbit perturbation $\{\delta r_{S},\delta\omega_{S}\}$ from (37).

As a check of self-consistency, notice that we have the choice of substituting $M$ either as a function of the physical values $\{r(t),\omega(t)\}$ using the results from (38), or the bare parameters $\{R_{B},\Omega_{B}\}$, which give the same result after summing up the perturbation expansions.

Substituting the corresponding perturbations back into (53), we obtain the explicit time-dependence of the precession norm $\Omega_{a}$,

[TABLE]

Combined with the expression for $\delta\omega(t)$ in terms of the time-independent bare parameters, we can perform the integration in (52) to write down

[TABLE]

The constant terms $\sin\Phi_{B}$ and $\cos\Phi_{B}$ can be absorbed by redefining the initial condition $i\ln{S}^{a}_{+}(t_{0})$, or via a bare parameter $i\ln\mathcal{S}^{a}_{+B}$,

[TABLE]

The logarithm of the spin components then becomes

[TABLE]

Given the spin vector expansions in terms of the bare parameters $\{R_{B},\Omega_{B},\Phi_{B},A_{B},A^{S}_{B},\mathcal{S}^{a}_{+B}\}$, the next step is to renormalize the spin components by replacing the bare parameters by the renormalized ones plus counter-terms, and splitting $t-t_{0}=(t-\tau)+(\tau-t_{0})$ with a choice of an arbitrary renormalization scale $\tau$.

B.1 Spin Renormalization

To begin, the bare parameter $\ln\mathcal{S}^{a}_{+B}$ is related to the renormalized value $\ln\mathcal{S}^{a}_{+R}$ through

[TABLE]

The renormalization treatment is performed for the natural logarithm of the spin components. As a result, $\mathcal{S}^{a}_{+B}=\mathcal{S}^{a}_{+R}e^{\delta_{\ln S}^{a}}$. The exponential implies that it is the phase of the precession that is renormalized. Dividing the bare parameters into the renormalized parts and the counter-terms and introducing the renormalization scale $\tau$, Eq. (57) then becomes

[TABLE]

The counter-terms are the combined results in Galley and Rothstein (2017) and (41),

[TABLE]

After some algebra, (59) can be simplified to

[TABLE]

Notice that the terms proportional to $(t-\tau)(\tau-t_{0})$ are completely canceled, which was emphasized in Galley and Rothstein (2017) as an important check of self-consistency. Here the cancellation is due to exactly the same set of substitutions we could use to replace $M$ to obtain (53), where the two different choices led to the same expansion result.

To cancel the remaining secular pieces that are proportional to the powers of $(\tau-t_{0})$, the counter-term $\delta_{\ln S}^{a}$ is fixed to be

[TABLE]

Choosing the arbitrary scale $\tau$ to equal $t$, the renormalized solution to $\ln\mathcal{S}^{a}_{+}(t)$ becomes

[TABLE]

or more explicitly in terms of the exponential,

[TABLE]

The renormalized quantities as functions of time have runnings obtained from the RG flow in Section A.3, with only the remaining the spin component bare parameter $\mathcal{S}^{a}_{+R}$ to be done in the next section.

B.2 Spin Component Renormalization Group Solution

The running of the renormalized parameter $\mathcal{S}^{a}_{+R}$ can be determined using (61), which leads to

[TABLE]

It seems to be formally divergent and has the dependence on the cut-off $t_{0}$. However, replacing the derivatives of $R_{R}$ and $\Omega_{R}$ by their RG equations (43a) and (43b), we encounter the non-trivial cancellation and obtain a finite $\beta$-function,

[TABLE]

Notice the similarity in form between the RG equation and (51), the precession equation we start with.

In order to find a solution to the RG equation of the spin component, we can write the relation between the $\tau$-derivative of $i\ln\mathcal{S}^{a}_{+R}$ and the derivative with respect to the renormalized parameter $R_{R}$ as

[TABLE]

Using the RGEs (43a) and (65), we obtain a solution to $\mathcal{S}^{a}_{+R}(\tau)$ in terms of $R_{R}(\tau)$ and initial conditions

[TABLE]

The expressions are not unique in terms of $\Phi_{R}(t)$ and $R_{R}(t)$ due to several RG invariants between them. The invariance over time with spin components can be found from the $\mathcal{S}^{a}_{+R}(\tau)$ solution, which is given by

[TABLE]

Putting all the pieces together, the resummed solution of $S^{a}_{+}(t)$ is given by

[TABLE]

with $\{R_{R}(t),\Omega_{R}(t),\Phi_{R}(t),A_{R}(t),A^{S}_{R}(t),\mathcal{S}^{a}_{+R}(t)\}$ given by (44), (A.3) and (B.2).

The quantity $\mathcal{S}^{a}_{+R}(t_{i})$ depends on the initial conditions of dynamics and spin vectors. For instance, taking the initial input $S^{a}_{n}(t_{i})$ and $S^{a}_{\lambda}(t_{i})$, while getting $A_{R}(t_{i})$, $R_{R}(t_{i})$, $\Omega_{R}(t_{i})$ and $\Phi_{R}(t_{i})$ from numerically solving the initial conditions $r(t_{i})$, $\dot{r}(t_{i})$, $\omega(t_{i})$ and $\phi(t_{i})$ from the dynamics, we can determine the value of $\mathcal{S}^{a}_{+R}(t_{i})$ through

[TABLE]

One immediate validation of the formulation is that the length of the spin vector should be a constant. Thus $\left|S^{a}_{+}(t)\right|=\sqrt{(S^{a}_{n})^{2}+(S^{a}_{\lambda})^{2}}$ should be a constant, since $S^{a}_{l}$ does not change with time. From (69) and (70) we can see that the length is preserved, $\left|S^{a}_{+}(t)\right|=\left|S^{a}_{+R}(t_{i})\right|=\sqrt{(S^{a}_{n}(t_{i}))^{2}+(S^{a}_{\lambda}(t_{i}))^{2}}$ as long as \big{(}M^{1/2}R_{R}^{3/2}(t)-\mathscr{S}\big{)}>0, the same constraint we encounter for the solutions of the orbit equations of motion.

Appendix C The Moving Triad Evolution

In the text, the resummed analytic expressions for the orbital equations of motion and spin precession we obtained are written in terms of the moving triad vectors $\{\bm{n},\bm{\lambda},\bm{l}\}$. To transform the complete results into a fixed frame, we follow the solutions to the evolution equations for the moving triad in Blanchet et al. (2011); Blanchet (2014) for the 1.5PN order conservative dynamics and build the moving triad evolution for the radiative dynamics on quasi-circular orbits.

We start by briefly summarizing the conservative moving triad evolution solution that relies fundamentally on the conservation of the total angular momentum $\bm{J}$. An orthonormal inertial frame $\{\bm{x},\bm{y},\bm{z}\}$ is then introduced with $\bm{J}/J$ as the fixed direction $\bm{z}$. Three Euler angles $\alpha(t),\iota(t),\Phi(t)$ are defined to specify the moving triad within the fixed frame as shown in Figure 4. The azimuth $\alpha$ and the inclination $\iota$ are the standard spherical coordinates of the Newtonian angular momentum direction $\bm{l}$. The angle $\Phi$ is defined to be the angle between $\bm{n}$ and $\bm{x}_{l}$, where

[TABLE]

forming the instantaneous orbital plane and with $\bm{l}$ to complete an auxiliary orthonormal basis {$\bm{x}_{l},\bm{y}_{l},\bm{l}$}.

In terms of the Euler angles, the relation between the moving triad $\{\bm{n}(t),\bm{\lambda}(t),\bm{l}(t)\}$ and the fixed Cartesian frame $\{\bm{x},\bm{y},\bm{z}\}$ can be written as

[TABLE]

The evolution solutions to the Euler angles up to linear order in spin are given by the components of the total angular momentum $\bm{J}=J_{n}(t)\hat{\bm{n}}+J_{\lambda}(t)\hat{\bm{\lambda}}+J_{l}(t)\hat{\bm{l}}$ as

[TABLE]

where $\phi$ is the orbital phase, for which the resummed solution is given by (42c) for the radiative binary orbits.

Finally, expressed in terms of some initial basis $\{\bm{n}_{0},\bm{\lambda}_{0},\bm{l}_{0}\}$ with corresponding Euler angles $\{\alpha_{0},\iota_{0},\Phi_{0}\}$, the moving triad $\{\bm{n}(t),\bm{\lambda}(t),\bm{l}(t)\}$ is given by

[TABLE]

where $\bm{m}\equiv\frac{1}{\sqrt{2}}(\bm{n}+i\bm{\lambda})$ is a complex null vector.

The crucial point of this moving triad solution is the conservation of the total angular momentum and the ability to write out its components in the moving triad for all time, not the physical meaning to $\bm{J}$. In order to apply the triad solutions to a radiative motion where $\bm{J}$ can change, we find such a quantity that satisfies the requirements by observing the calculation of $\mathrm{d}\bm{J}/\mathrm{d}t$ for a conservative quasi-circular orbit. It is conventional to decompose $\bm{J}=\bm{L}+\bm{S}$, where $\bm{S}$ is the total spin specified by the choices of spin variables following Kidder (1995), and $\bm{L}$ is the sum of the non-spinning Newtonian $\bm{L}_{\textrm{N}}$ and the leading order spin-orbit contribution $\bm{L}_{\textrm{SO}}$, given by

[TABLE]

Written in terms of the moving triad components and taking the orbit radius and frequency as constants $R$ and $\Omega$ for the quasi-circular approximation, the spin-orbit momentum becomes

[TABLE]

For a conservative system without radiation, the time derivative to the sum $\bm{J}=\bm{L}_{\textrm{N}}+\bm{L}_{\textrm{SO}}+\bm{S}_{1}+\bm{S}_{2}$ should vanish up to the Newtonian and leading spin order. By carrying out the detail calculation, we find that

[TABLE]

Thus the sum is

[TABLE]

which correspond to 1PN terms to be fixed by including higher-order orbital angular momenta. Notice that the time derivative of the spins in $\dot{\bm{L}}_{\textrm{SO}}$ are completely canceled by $\dot{\bm{S}}$ and $\dot{\bm{L}}_{\textrm{N}}$ at Newtonian order. Therefore we propose that for a radiative quasi-circular binary, the following quantity is conserved:

[TABLE]

Compared with the conservative expressions, we replace the constant orbital radius and frequency by the initial orbital radius and frequency. The conservative spin components are changed into the time-dependent resummed radiative spin component results. The time derivative of this quantity $\mathscr{J}$ is $\sim\mathcal{O}(v^{4}S)$ but we are able to avoid the loss of total angular momentum due to non-spinning radiation at $\mathcal{O}(v^{5})$. Using the substitution with $\mathscr{J}$ instead of $\bm{J}$ into the moving frame solutions (72-75), we can generate 3D-plot of the orbital radius evolution and animations of binary inspiral with spin orientation at every instant.