Eccentric Binary Black Holes with Spin via the Direct Integration of the Post-Newtonian Equations of Motion

TL;DR

This paper develops a method to accurately model gravitational waveforms from generic binary black hole systems with eccentric orbits and spins by directly integrating advanced post-Newtonian equations, enhancing waveform predictions for detectors.

Contribution

It introduces a comprehensive integration of third and a half post-Newtonian equations for spinning, eccentric binary black holes, extending previous models to fully generic configurations.

Findings

Validated waveform models against existing methods for aligned-spin cases.

Analyzed effects of eccentricity and spins on gravitational wave modulations.

Demonstrated relevance for current and future gravitational wave detectors.

Abstract

We integrate the third and a half post-Newtonian equations of motion for a fully generic binary black hole system, allowing both for non-circular orbits, and for one or both of the black holes to spin, in any orientation. Using the second post-Newtonian order expression beyond the leading order quadrupole formula, we study the gravitational waveforms produced from such systems. Our results are validated by comparing to Taylor T4 in the aligned-spin circular cases, and the additional effects and modulations introduced by the eccentricity and the spins are analyzed. We use the framework to evaluate the evolution of eccentricity, and trace its contributions to source terms corresponding to the different definitions. Finally, we discuss how this direct integration equations-of-motion code may be relevant to existing and upcoming gravitational wave detectors, showing fully generic, precessing, eccentric gravitational waveforms from a fiducial binary system with the orbital plane and spin precession, and the eccentricity reduction.

Figures

Tables

Table 1. Table 1: The overlaps between Taylor T4 and EOM methods, documented for different separations, mass ratios ( q 𝑞 q ), and aligned dimensionless spin values ( χ 1 subscript 𝜒 1 \chi_{1} and χ 2 subscript 𝜒 2 \chi_{2} ), at both the consistent T4 flux order and the highest T4 fluxes available. There are several interesting lines to pay attention to: as the spins increase in alignment, there is an inflection point around χ 1 = χ 2 = 0.9 subscript 𝜒 1 subscript 𝜒 2 0.9 \chi_{1}=\chi_{2}=0.9 , where the overlap drops then increases again. In addition, at 50 ​ M 50 𝑀 50M separation, we increase the resolution of the timestep to probe the effects of numerical resolution on the overlap calculations. We find no effect on the overlap of the waveforms due to the timestep resolution.

Separation					Duration			Overlaps $O_{\mathrm{T4}-\mathrm{EOM}}$
$r_{\mathrm{init}}$ $\left[M\right]$	$r_{\mathrm{fin}}$ $\left[M\right]$	$q$	${\chi }_1$	${\chi }_2$	$t$ [M]	orbits	$dt$ [s]	${\mathrm{T4}}_{\mathrm{consistent}}$	${\mathrm{T4}}_{\mathrm{high}}$
100	99	1	0	0	318429.8	$\sim 50$	$1.09061\times {10}^{-3}$	0.99999065	0.999990644
100	99	1	0.3	0.3	318429.8	$\sim 50$	$1.04527\times {10}^{-3}$	0.99993352	0.999933512
100	99	1	0.6	0.6	318429.8	$\sim 50$	$1.04527\times {10}^{-3}$	0.99983625	0.999836244
100	99	1	0.9	0.9	318429.8	$\sim 50$	$1.04527\times {10}^{-3}$	0.99971069	0.999710688
100	99	1	0.3	-0.3	318429.8	$\sim 50$	$1.04527\times {10}^{-3}$	0.99998949	0.999989489
100	99	1	0.6	-0.6	318429.8	$\sim 50$	$1.04527\times {10}^{-3}$	0.9999856287	0.9999856265
100	99	1	0.9	-0.9	318429.8	$\sim 50$	$1.04527\times {10}^{-3}$	0.99997797	0.9999778678
100	99	1	0.9	-0.5	318429.8	$\sim 50$	$1.04527\times {10}^{-3}$	0.9999415165	0.9999415133
100	99.1	2	0	0	318480.2	$\sim 50$	$1.045439\times {10}^{-3}$	0.99999283	0.999992828
100	99.7	10	0	0	318701.62	$\sim 50$	$1.04617\times {10}^{-3}$	0.9999991669	0.9999991668
100	99.95	100	0	0	318833.8	$\sim 50$	$1.0466\times {10}^{-3}$	0.9999999893	0.9999999893
50	47.2	1	0	0	114043.4	$\sim 50$	$3.90593\times {10}^{-4}$	0.994904572	0.994898471
50	47.2	1	0	0	114049.1	$\sim 50$	$2.246257\times {10}^{-4}$	0.994797482	0.994791317
50	47.3	1	0.3	0.3	114043.4	$\sim 50$	$3.74357\times {10}^{-4}$	0.983693663	0.983675238
50	47.3	1	0.6	0.6	114043.4	$\sim 50$	$3.74357\times {10}^{-4}$	0.969583498	0.969553474
50	47.4	1	0.7	0.7	114043.4	$\sim 50$	$3.74357\times {10}^{-4}$	0.981885223	0.981867878
50	47.4	1	0.8	0.8	114043.4	$\sim 50$	$3.74357\times {10}^{-4}$	0.991724021	0.991717488
50	47.4	1	0.85	0.85	114043.4	$\sim 50$	$3.74357\times {10}^{-4}$	0.992601112	0.992597499
50	47.4	1	0.9	0.9	114043.4	$\sim 50$	$3.74357\times {10}^{-4}$	0.990197785	0.990190154
50	47.4	1	0.95	0.95	114043.4	$\sim 50$	$3.74357\times {10}^{-4}$	0.984804119	0.984791162
50	47.2	1	0.3	-0.3	114043.4	$\sim 50$	$3.74357\times {10}^{-4}$	0.995280388	0.995274373
50	47.2	1	0.6	-0.6	114043.4	$\sim 50$	$3.74357\times {10}^{-4}$	0.9971254502	0.997121275
50	47.2	1	0.9	-0.9	114043.4	$\sim 50$	$3.74357\times {10}^{-4}$	0.9988241423	0.998822279
50	47.3	1	0.9	-0.5	114043.4	$\sim 50$	$3.74357\times {10}^{-4}$	0.985182006	0.985164755
50	47.5	2	0	0	114095.5	$\sim 50$	$3.74528\times {10}^{-4}$	0.99467073	0.994667208
50	49.1	10	0	0	114299.9	$\sim 50$	$3.75199\times {10}^{-4}$	0.999489566	0.999489472
50	49.9	100	0	0	114395.1	$\sim 50$	$3.75512\times {10}^{-4}$	0.999994205	0.999994204
20	16	1	0	0	11965	$\sim 20$	$3.92762\times {10}^{-5}$	0.782932526	0.775227367
20	18.2	1	0	0	5982.5	$\sim 10$	$1.963812\times {10}^{-5}$	0.917453743	0.916892613
20	16.5	2	0	0	11979.2	$\sim 20$	$3.93229\times {10}^{-5}$	0.827044247	0.821830277
20	18.8	10	0	0	12050.8	$\sim 20$	$3.95579\times {10}^{-5}$	0.972718171	0.97255266
20	19.9	100	0	0	12087.1	$\sim 20$	$3.96769\times {10}^{-5}$	0.99919055	0.999190336

Equations

$$\vec{S}={\vec{S}}_1+{\vec{S}}_2,$$

$$\frac{\vec{\Sigma }}{M}=\frac{{\vec{S}}_2}{m_2}-\frac{{\vec{S}}_1}{m_1}.$$

$$\frac{\mathrm{d}\vec{v}}{\mathrm{d}t}$$

$$\frac{\mathrm{d}\vec{S}}{\mathrm{d}t}$$

$$\frac{\mathrm{d}\vec{\Sigma }}{\mathrm{d}t}$$

$$h^{ij}$$

$$h_{+,\times }$$

$$\displaystyle\psi=\phi-2M_{\rm ADM}\Omega\ln\,\bigg{(}\frac{\Omega}{\Omega_{0}}\bigg{)}\,.$$

$$\displaystyle M_{\rm ADM}=M\bigg{[}1-\frac{\eta}{2}\gamma+\frac{\eta}{8}(1-\eta)\gamma^{2}+O\bigg{(}\frac{1}{c^{6}}\bigg{)}\bigg{]}\,,$$

$$e=\frac{r_{\max }-r_{\min }}{r_{\max }+r_{\min }},$$

$$\frac{\mathrm{d}\vec{v}}{\mathrm{d}t}=(\ddot{r}-r{\Omega }^2)\hat{n}+(r\dot{\Omega }+2\dot{r}\Omega )\hat{\lambda }+r\varpi \Omega \hat{\ell },$$

$$\frac{\mathrm{d}\vec{v}}{\mathrm{d}t}\simeq -r{\Omega }^2\hat{n}+r\varpi \Omega \hat{\ell }.$$

$$\frac{\mathrm{d}r}{\mathrm{d}t}=v_{\mathrm{insp}}(t)+B\cos (\omega t+\phi ).$$

$${\dot{r}}_{\mathrm{new}}={\dot{r}}_{\mathrm{old}}+\Delta \dot{r},$$

$${\Omega }_{\mathrm{new}}={\Omega }_{\mathrm{old}}+\Delta \Omega ,$$

$$\Delta \dot{r}$$

$$\Delta \Omega$$

$$\mathcal{F}(v)$$

$$\Omega =\frac{|\vec{r}\times \vec{v}|}{r^2},$$

$${\Omega }^2$$

$${\delta }_{aprx1,aprx2}=\left|\frac{{\Omega }_{aprx1}-{\Omega }_{aprx2}}{{\Omega }_{aprx1}}\right|.$$

$${\delta }^{def1,def2}=\left|\frac{{\Omega }^{def1}-{\Omega }^{def2}}{{\Omega }^{def1}}\right|,$$

$$\displaystyle(h_{1}|h_{2})=4\max_{t_{c}}\bigg{|}\int^{f_{\rm high}}_{f_{\rm low}}\frac{h_{1}^{*}(f,t_{c})h_{2}(f)}{S_{n}(f)}\mathrm{d}f\bigg{|}\,,$$

$$O=\frac{(h_1|h_2)}{\mathrm{Norm}(h_1)\mathrm{Norm}(h_2)}.$$

$$\displaystyle\vec{e}=\bigg{(}\frac{v^{2}}{M}-\frac{1}{r}\bigg{)}\vec{r}-\frac{\vec{r}\cdot\vec{v}}{M}\vec{v}\,,$$

$$r(t)$$

$$\dot{r}(t)$$

$$\ddot{r}(t)$$

$$A(t)$$

$$e=\mathrm{Amp}(A)=\mathrm{Amp}\left(\frac{r^2\ddot{r}}{M}\right).$$

$$B(t)$$

$$\sqrt{A^2+B}$$

Full text

I Introduction

The direct detections of gravitational waves (GWs) Abbott et al. (2016a); Abbott et al. (2016b, 2018a) by the LIGO/VIRGO Collaboration Aasi et al. (2015); Acernese et al. (2015); Abbott et al. (2018b) have launched the era of multi-messenger astrophysics, both in providing a new window on events such as binary neutron star (BNS) mergers Abbott et al. (2017a), and by opening for study a completely new field of astrophysical events previously invisible to us, binary black hole (BBH) mergers Abbott et al. (2018a), meriting the Nobel Prize in Physics in 2017.

GWs from eccentric binaries Gopakumar and Iyer (2002); Yunes et al. (2009); Huerta et al. (2014); Hannam et al. (2014); Taracchini et al. (2014); Tanay et al. (2016); Huerta et al. (2017, 2018); Hinder et al. (2018); Lower et al. (2018); Huerta et al. (2019) have become an important topic of study as the next observing runs of the Advanced Laser Interferometer Gravitational-wave Observatory (aLIGO) and VIRGO approach. While it is not expected that LIGO/VIRGO sources in isolated binaries will have significant eccentricity in band, binaries in dense stellar clusters Merritt (2013) may, however, retain significant eccentricity from dynamical interactions prior to entry into the LIGO/VIRGO band. Recent simulations of globular clusters indicate that a distinct population of compact binaries exist (from Refs. Rodriguez et al. (2018, 2016), about 3% of binaries) and enter the LIGO/VIRGO band (canonically taken to be $10$ Hz as the lower bound) with significant eccentricity ($e>0.1$) Rodriguez et al. (2018); Antonini and Rasio (2016). These results add further motivation, for if binaries can form and merge in this way, then there will be a minority (but a distinct minority), that will enter with eccentricity, and which may be missed without templates for the match filtering.

In addition, a large fraction of binaries in globular clusters will have significant eccentricity ($\sim 50\%$ of in-cluster mergers) and will be detectable for the entirety of the LISA band ($\sim 10^{-3}$ – $10^{-1}$ Hz) Breivik et al. (2016); D’Orazio and Samsing (2018). Both 2-body mergers (highly eccentric in cluster mergers in between single binary interactions) and 3-body mergers (when a BBH forms with such high eccentricity that it is essentially a GW capture before it can interact with a third body) will occur in clusters. LISA will be able to measure the 2-body mergers, but not the 3-body mergers Samsing and D’Orazio (2018). Current work is being carried out to see if detection of eccentric sources will be possible with the proposed space-based detector, and current prospects look promising D’Orazio and Samsing (2018); Breivik et al. (2016); Sesana (2016).

In this paper, we develop and calculate the GW waveforms and orbital dynamics from generically spinning, eccentric BBHs. We use the Lagrangian formulation of the post-Newtonian (PN) equations of motion (EOM) in the harmonic gauge for the generation of precessing, eccentric GW signatures Blanchet et al. (1995); Blanchet (1995, 1996); Blanchet et al. (1998); Blanchet (2006); Will (2005); Faye et al. (2006); Blanchet et al. (2006); Will (2011); Marsat et al. (2013); Bohe et al. (2013); Blanchet (2014); Marsat (2015); Bohe et al. (2015); Bernard et al. (2018). Our approach allows us to use any spin values, mass ratios, and eccentricities, without restricting to planar orbits or co-precessing frames, so long as the binary has a large enough separation (around a periapsis passage of $r\sim 10M$), such that the underlying PN theory does not break down Poisson (1995); Mino et al. (1997); Yunes and Berti (2008); Zhang et al. (2011); Sago et al. (2016); Fujita et al. (2018). This approach offers a major step forward as a way to generate eccentric, precessing GW waveforms in a direct and straightforward way, without carrying any of the additional restrictive assumptions of quasi-ellipticity or adiabaticity that many current waveform generators use. Furthermore, we can give more precise waveforms since we do not ignore any timescale effect (see, e.g., Ref. Klein et al. (2018)).

This work is an extension of work first done by Lincoln and Will Lincoln and Will (1990), and of integrations as in Refs. Pati and Will (2000, 2002). In more recent works, Refs. Levin et al. (2011) and Csizmadia et al. (2012) used this framework to 3.5PN to generate eccentric binaries in the PN harmonic gauge, and calculate the relevant waveforms. The research focus of this project is to calculate the orbital quantities and GW waveforms from generic binaries with the requisite accuracy to enable potential future observations from GW detectors. This paper details the methodology, and extends the previous work done by other authors by giving quantitative comparisons to other known PN methods for the first time.

This Lagrangian formulation has larger applications as well. Developing a general EOM that can handle arbitrarily precessing and eccentric BBHs means that we can apply this EOM to the precessing spacetime developed in our previous work Nakano et al. (2016). We must have an accurate EOM for the general precessing BBH that can handle the orbital plane precession of the binary due to spin coupling with the orbit Owen et al. (1998), and also the individual spins precessing in order to evolve the BBH. The Lagrangian formulation excels at all of these, and can directly be used for this evolution.

The regime of final inspiral to plunge and merger can only be modeled by numerical relativity Pretorius (2005); Campanelli et al. (2006a); Baker et al. (2006); Campanelli et al. (2006b); Campanelli et al. (2006); Campanelli et al. (2007a, b); Lousto and Healy (2015); Lousto et al. (2016). This does, however, provide us a unique opportunity to use the numerical relativity regime; we can stitch our PN evolution onto the beginning of NR simulations, and thus create a full waveform model for the binary using hybridization of waveforms Campanelli et al. (2009); Ajith et al. (2012).

This paper is organized as follows. Sec. II gives an overview of the PN EOMs for generic spins and eccentricities. Sec. III shows the power and flexibility of the EOM code, and demonstrates our main results with orbital and waveform quantities for some fiducial systems with eccentricity and spins. Sec. IV quantifies the comparisons of the EOM code with the Taylor T4 method, the initial conditions for quasi-circular orbits, including orbital frequency comparison, and waveform overlaps. In Sec. V, we discuss a simple test of orbital eccentricity by using the EOM. Sec. VI contains discussion, conclusions, and future work.

Throughout this paper, we use the geometric unit system, where $c=G=1$, with the useful conversion factor $1M_{\odot}=1.477\;{\rm{km}}=4.926\times 10^{-6}\;{\rm{s}}$, although we will keep some $1/c$ factors to count PN orders.

II post-Newtonian Formulation

II.1 Orbital Position and Spin Equations of Motion

The general case of unaligned spins (i.e., spins that are neither aligned nor anti-aligned with the orbital angular momentum) leads to precession of the individual spins as well as precession of the orbital angular momentum vector, and therefore precession of the orbital plane. To describe this complicated motion requires two sets of evolution equations: one set describing the positions (trajectories) of the point masses, and one set describing the evolution of the spin of each point mass. We will review the position equations first, then the precession equations. There are several PN papers that go through this in detail. See Ref. Blanchet (2014) for the non-spinning terms, and Refs. Marsat et al. (2013); Bohe et al. (2013) for the spin dependent terms.

The system comprises of two bodies ($a=1,2$), each described by a mass $m_{a}$, a spin $\vec{S}_{a}$, a spin precession vector $\vec{\Omega}_{a}$, and position and velocity vectors $\vec{r}_{a}$, $\vec{v}_{a}$. We also define the total combinations $M=m_{1}+m_{2}$ for the total mass, $\eta=m_{1}m_{2}/M^{2}$ for the symmetric mass ratio, $\mu=M\eta$ for the reduced mass, and

[TABLE]

for the total spin, as well as the relative combinations $\vec{r}=\vec{r}_{1}-\vec{r}_{2}=r\hat{n}$ (where $r=|\vec{r}|$) denoting the relative separation vector, $\vec{v}=\vec{v}_{1}-\vec{v}_{2}$ for the relative velocity, and the relative spin difference $\vec{\Sigma}$, given by

[TABLE]

The schematic EOM for the orbital positions, in terms of relative variables and in the center-of-mass frame, then takes the form:

[TABLE]

with the non-spinning components of these $\mathcal{A}$ and $\mathcal{B}$ terms defined in Ref. Blanchet (2014), and the spin terms and expressions for ${\mathcal{C}}_{1}$, ${\mathcal{C}}_{2}$, ${\mathcal{C}}_{3}$, and ${\mathcal{C}}_{4}$ defined in Refs. Faye et al. (2006); Blanchet et al. (2006); Marsat et al. (2013); Bohe et al. (2013).

The corresponding precession equations from Ref. Bohe et al. (2013) take the form:

[TABLE]

Following Refs. Will (2005); Blanchet et al. (2006); Blanchet (2014), we use the Tulczyjew spin supplementary condition (SSC) to define a spin vector with conserved Euclidean norm. For this SSC, the higher order spin-orbit terms have been derived in Refs. Faye et al. (2006); Blanchet et al. (2006); Marsat et al. (2013); Bohe et al. (2013); Blanchet (2014). The next-to-leading order spin-spin terms were derived in Ref. Bohe et al. (2015), and the leading order cubic-in-spin terms were derived in Ref. Marsat (2015). To date, leading order quartic- and quintic-in-spin contributions to the EOM have not been derived in PN harmonic coordinates with this SSC. See Appendix B for some details.

Up to 3.5PN 111A PN order $N$ is said to be a term of order $(v/c)^{2N}$. order in this formalism, and for maximal spins ($|\vec{\chi}_{\mathrm{a}}|\simeq 1$ where $\vec{\chi}_{\mathrm{a}}=\vec{S}_{\mathrm{a}}/m_{\mathrm{a}}^{2}$), spin-orbit effects contribute to the EOM at 1.5PN, 2.5PN, and 3.5PN. Spin-spin effects contribute at 2PN and 3PN. Cubic-in-spin effects start from 3.5PN. Quartic- and higher order in-spin effects are beyond 3.5PN.

Equations (5a) and (5b), along with the acceleration equation (4) comprise the EOMs of the two-body system. For a given set of initial conditions, this system of ODE’s is easily solved numerically 222We have used Mathematica’s built-in RK4 solver to solve the EOMs..

We close this section with some brief comments about radiation reaction. Through 3.5PN order (inclusive), radiation reaction effects arise exclusively in the non-spin part of $\mathcal{B}$. Radiation reaction effects first appear in the non-spin part of $\mathcal{A}$ at 6PN. Spin effects enter radiation reaction at 4PN, and radiation reaction effects enter spin contributions at 4PN. Furthermore, for quasi-circular orbits, the non-spin part of $\mathcal{B}$ contains only radiation reaction effects. In that case, radiation reaction is completely controlled to 3.5PN by turning on/off the non-spin part of $\mathcal{B}$.

II.2 Generation of Gravitational Radiation

To calculate the polarization of GWs, $h_{+}$ and $h_{\times}$, we need to define the principle axes for the GWs. In particular, we define $\hat{N}$ which is the radial direction to the observer, $\hat{p}$, which lies along the line of nodes (for our purposes we can set to be along the $y$ axis), and $\hat{q}$ which is orthogonal to $\hat{N}$ and $\hat{p}$ Will and Wiseman (1996). Note that $\hat{N}\neq\hat{p}\times\hat{q}$.

This allows us to define the inclination $i$, and the phase $\phi$, with respect to the Newtonian angular momentum $\vec{L}=\mu\vec{r}\times\vec{v}$ (see Fig. 1).

For the GW calculation, we use the 2PN order formula beyond the leading order quadrupole approximation by following Ref. Will and Wiseman (1996):

[TABLE]

where $R$ denotes the distance between the observer and the binary, and the individual PN terms are broken up into non-spinning, spin-orbit (SO), and spin-spin (SS) contributions. For example, $P^{3/2}Q_{\rm SO}^{ij}$ is the 1.5PN SO contribution, and TT means the transverse-traceless part. We have implemented all of the non-spinning contributions to the quadrupole moment up to 2PN, in order to be in agreement with the GW calculations for Taylor T4.

For the Taylor T4 GW waveforms, we follow Ref. Blanchet (2014), which expands the waveforms $h_{+}$ and $h_{\times}$ into PN orders as powers of the frequency variable $x=(M\Omega)^{2/3}$,

[TABLE]

to the desired PN order. The phase variable $\psi$ is related to the binary orbital phase $\phi$ by the auxiliary phase variable Blanchet (2014):

[TABLE]

The constant frequency $\Omega_{0}$ can be chosen at will, and for our analysis is set to the one related to $10$ Hz (chosen naively as the entry frequency into the LIGO/VIRGO band). The mass variable in the GW calculation is the ADM mass of the binary:

[TABLE]

where $\gamma=M/r$ is an expression of the PN parameter of the system. The $H_{(p/2)+,\times}(\psi,\cos i,\sin i;\ln x)$ terms are the specific expansion coefficients, and are dependent on the auxiliary phase variable $\psi$, and the inclination of the binary $i$. The log terms of the frequency variable first occur at 3PN in the GW waveform. Explicit expressions for $H_{(p/2)+,\times}$ appear in Eqs. (322) and (323) of Ref. Blanchet (2014), which can be easily computed from the output of our EOM integrator.

III Eccentric, Precessing Binaries

We now demonstrate the power of our framework for the most generic systems of BBHs, ones with non-unity mass ratios, spin orientations, and eccentricities. Since this parameter space is very large, we use a fiducial example for a demonstration of the orbital and waveform quantities, and know that the formalism is generic for any set of mass, spin, and eccentricity parameters.

For the following example, we follow Ref. Rodriguez et al. (2018) to draw a fiducial binary from the histograms that globular cluster simulations produce. Other than picking physically relevant parameters, we do not restrict ourselves to any particular parameters.

The system that we chose to evaluate is a binary with an initial periapsis passage of $35M$, an initial value eccentricity picked to be $0.4$, a mass ratio $q=3/2$, and dimensionless spin parameter values of $\chi_{1}=(-0.3,0,0)$ and $\chi_{2}=(0,0,0.3)$. The initial periapsis passage and initial eccentricity were selected such that a binary with masses comparable to detected LIGO sources (we use a total mass of 25 solar masses in the analysis of the GW waveform later to dimensionalize the units) would fall just into the LIGO band at periapsis. The chosen spin parameters have low spin to lie in the physically relevant globular cluster results, and are initially strongly misaligned to give lots of spin-spin and spin-orbit precession.

The extrinsic parameters for waveform production that we chose for this fiducial system are optimized for ease. The distance from the source we set to be $500$ Mpc (a redshift of $\sim 0.1$, small enough not to need to take cosmological effects into account), with the orientation set to $(i,\varpi,\Omega)=(0,0,0)$ where $i$, $\varpi$ and $\Omega$ are the extrinsic orbital parameters, and denote the inclination angle, argument of periapsis and longitude of the ascending node, respectively. This is referred to as the optimal orientation, where the binary is face on (inclination equal to zero), the argument of periapsis is set to zero, and the longitude of the ascending node is equal to zero (the binary is not tilted with respect to the observer). These 13 parameters are enough to completely specify the binary that we are describing (though it will not give sky localization, as this set of parameters is only for one detector, and a second detector would be needed for localization purposes). For a real source, it is a straightforward generalization to put in the detector antenna patterns.

The final thing that we will need for the analysis is the eccentricity definition which we use for generic eccentricities. This is given as:

[TABLE]

which is a workable definition of eccentricity for the orbits Levin et al. (2011); Csizmadia et al. (2012) in the situation of non-negligible eccentricity. It is noted that this definition relies on the position of the binary at different points in the orbit; therefore, it is not an instantaneous definition at a point, but is instead averaged over the orbit.

We start by evolving this binary and plotting the orbital trajectories and spin vectors, shown in Figs. 2, 3, and 4. The eccentricity of the orbit is measured following Eq. (13), from which we calculate the eccentricity in the beginning of the evolution as $e_{\rm init}=0.449$ and the final one as $e_{\rm fin}=0.074$. The discrepancy between the initial eccentricity in the code and the eccentricity that we calculate is due to the setup of the initial conditions, as the initial inputs to the code are calculated as Newtonian initial parameters given a PN $\dot{r}$ and $\omega$, which still does not fully incorporate the PN effects.

The eccentricity reduction can be observed fully in the orbital separation as a function of time (see Fig. 2), and efficiently radiates the eccentricity away over the evolution.

With the orbital quantities, we can calculate the waveforms as well, using the GW waveform prescription in Ref. Will and Wiseman (1996), using the quadrupole formula up to 2PN corrections in Eq. (8). The plus and cross polarizations are tabulated in Figs. 5 and 6.

We calculate the initial and final periapsis and apoapsis GW frequencies directly from the waveform as opposed to an approximate GW frequency from the orbital frequency in order to be as accurate as possible in the calculation of the GW frequency for the entry into a relevant detector (such as LIGO). For this binary, the measured initial frequencies are $3.5$ Hz at apoapsis and $10.5$ Hz at periapsis, and the end frequencies are $55.5$ Hz at apoapsis and $66.2$ Hz at periapsis.

The spin precession frequency is calculated in the same way as the periapsis and apoapsis frequencies by the peak to peak calculation on the waveform. This is not an exact calculation for this frequency, as the spin precession has some ambiguity in the waveform. Since one needs to be able to define the peak in the precessional modulation, this may off by an orbit or two. For the purposes of these calculations, we use the peak amplitude in the waveform for each larger spin precession modulation, getting a value of $0.1$ Hz.

An interesting note with these calculations is that the peak periapsis frequency actually drops with time. This is an artifact of averaging over half an orbit, and when calculating the frequency through $\mathrm{d}\phi/\mathrm{d}t$, this is not an issue. The period in the later stages of the orbit are cleanly broken up into periapsis and apoapsis frequencies. We observe that the sharp periapsis smooths out as the binary circularizes, as expected.

IV Comparisons and Validation

We now turn to comparing the methodology of the direct integration to our implementation of the Taylor T4 method (see Appendix C based on, e.g., Ref. Boyle et al. (2007)), restricted to the case when both should be valid.

Taylor T4 is a PN approximant that assumes quasi-circular, non-precessing orbits. So to compare to Taylor T4, we restrict ourselves to black holes (BHs) that are sufficiently separated, have vanishing orbital eccentricity, and with both spins only along the axis of the orbital angular momentum $\vec{L}$ (each could be either aligned, anti-aligned, or zero).

To begin, we will define a consistent comparison for both Taylor T4 and the direct integration EOMs, and outline the steps needed to reduce to this comparison. We will then compare the orbital frequency and other orbital quantities, and the GW waveforms generated by these systems.

IV.1 Quasi-Circular Orbits in the Direct Integration

Until now we have made no assumptions about eccentricity; all of the expressions to this point are valid for general orbits. Here, we reduce the EOMs to the case of quasi-circular orbits, which is important for our comparisons to other quasi-circular approximants.

To do this, following Ref. Faye et al. (2006), we introduce a moving orthonormal triad $\{\hat{n},\hat{\lambda},\hat{\ell}\}$, where $\hat{n}$ is the same as above, $\hat{\ell}=\hat{n}\times\vec{v}/|\hat{n}\times\vec{v}|$, and $\hat{\lambda}=\hat{\ell}\times\hat{n}$ (see Appendix A). Notice that $\hat{n}$ and $\hat{\lambda}$ span the orbital plane, while $\hat{\ell}$ is normal to it. In this basis, general kinematic considerations lead to the acceleration equation

[TABLE]

where the dot denotes $\mathrm{d}/\mathrm{d}t$, $\Omega$ is the orbital frequency, and $\varpi$ is the orbital plane precession frequency defined by $\varpi=-\hat{\lambda}\cdot\mathrm{d}\hat{\ell}/\mathrm{d}t$. Quasi-circular orbits $\ddot{r}\simeq\dot{r}\simeq\dot{\Omega}\simeq 0$, the acceleration equation reduces to

[TABLE]

By identifying Eqs. (15) and (38), one obtains PN expansions for $\Omega^{2}$ and $\varpi$.

Introducing the frequency-related parameter $x\equiv(M\Omega)^{2/3}$ gives a relation $x(\gamma)$, which may then be inverted order-by-order to obtain a relation $\gamma(x)$. With $\gamma(x)$ in hand, one may re-express any function of the coordinate-related parameter $\gamma$ in terms of the frequency-related parameter $x$. This is often preferable, since expressions in terms of the frequency-related parameter are gauge-invariant, whereas expressions in terms of the coordinate-related parameter are not. Examples of quantities that are useful to express as functions of $x$ are the energy, $E(x)$, the flux, $\mathcal{F}(x)$, and the orbital phase, $\phi(x)$. Expressions for these are given in Ref. Ajith et al. (2007) in terms of phase variable $v$ (see Appendix C), and is the jumping off point for the Taylor T4 approximant (see, e.g., Ref. Boyle et al. (2007)).

IV.2 Eccentricity Reduction

For comparing to adiabatic methods such as Taylor T4, we need to be able to accurately and reliably give quasi-circular initial conditions to the direct EOMs numerically. This is more challenging than at first glance, because if we were to simply give Newtonian (or even PN Healy et al. (2017); Ramos-Buades et al. (2019)) initial conditions (detailed above), there would be an error on the order of the neglected PN terms in the initial trajectories. This would manifest as a spurious eccentricity and add undesirable dynamics into the simulation.

We follow a simple procedure to remove this unwanted eccentricity. This procedure has been developed to set up low eccentricity numerical relativity initial data Pfeiffer et al. (2007); Buonanno et al. (2011). We begin by modeling the inspiral as a superposition of two effects, the (real) inspiral, which is a smooth decrease in the orbital separation as a function of time, and the (unphysical) oscillation due to the spurious eccentricity.

We start with a simple assumption for the inspiral part and oscillatory part, namely:

[TABLE]

We take the inspiral model $v_{\mathrm{insp}}(t)$ to be a simple polynomial, which we can fit for $v_{\mathrm{insp}}(t)=v_{0}+v_{1}t+v_{2}t^{2}$ with coefficients $v_{0}$, $v_{1}$ and $v_{2}$, and the oscillatory piece $B\cos(\omega\,t+\varphi)$ with the amplitude $B$, frequency $\omega$, and initial phase $\varphi$. With this model, we run our EOM code with the quasi-circular initial conditions detailed above, and fit the data with this model. We can then subtract out the oscillatory piece, and iterate on this model as many times as we need to attain the quasi-circular initial conditions. This is illustrated in Fig. 7.

Mathematically, this is taking the initial conditions $\dot{r}$ and $\Omega$, and after fitting the inspiral, we update

[TABLE]

and

[TABLE]

with $\Delta\dot{r}$ and $\Delta\Omega$ given by:

[TABLE]

where $r_{0}$ is the initial separation of the binary, and $\Omega_{0}$ is the initial orbital frequency.

IV.3 Consistent PN Order in Taylor T4

The flux formula in Taylor T4 is higher order than what we can calculate in the EOM formalism. To combat this, we need to tailor the T4 fluxes to a consistent PN order. From, e.g., Eq. (A.13) of Ref. Ajith et al. (2007), we can see that the leading order of $v$ is a 2.5PN term. We can also see the series expansion in the flux is out to $v^{7}$, which is 3.5PN beyond leading order. The absolute PN order for the fluxes is then 6PN, which is far beyond the highest order terms in the EOM formalism. To be consistent with our non-adiabatic EOM formalism above, we must truncate the flux terms above next to leading order so that the total re-expanded rational fraction is consistently 3.5PN, i.e., for non-spinning binaries:

[TABLE]

This is what we mean by a “consistent” PN order when we compare to the EOMs in the following sections. When quantifying comparisons, we will use the consistent PN order and the high PN order (keeping the flux terms to 6PN) to track the effects of PN orders.

IV.4 Orbital Frequency Comparison: No Spin

With the orbital quantities from solving Eq. (51), and the direct EOM orbital quantities from solving Eq. (4), we can explore how to compare these to one another.

First, we need to give both methodologies the same initial conditions, to ensure an apples to apples comparison. This is achieved by using the orbital frequency calculated for the EOM orbit initially (see Eq. (IV.4) below), and setting the initial frequency of the T4 code to the same initial value. This is simple using the relation $\Omega=(v^{3}/M)$, and solving for $\Omega$. We also define the binary phase to start in the same location on the orbit, i.e., the binary starts on the $x$-axis and the angular momentum is defined in the usual right handed sense.

The first and most direct comparisons that we can make are with the orbital frequencies themselves. The orbital frequency via the Taylor T4 method is given by the second integration equation, $\mathrm{d}\phi/\mathrm{d}t=v^{3}/M$, as $\mathrm{d}\phi/\mathrm{d}t=\Omega$.

When we talk about the direct EOM method, however, we need to be careful in how we define the orbital phase. Since EOM method directly outputs the trajectories and velocities, we can calculate the orbital frequency in the Newtonian sense:

[TABLE]

which specify the orbital frequency given the orbital trajectories and velocities, $\vec{r}$ and $\vec{v}$. To contrast this, we can use an alternate definition from PN given in Ref. Blanchet (2014) for no spins as

[TABLE]

where $r_{0}^{\prime}$ is a gauge constant that we set to $10M$, which is also used in the direct EOM method, and $\gamma$ is the PN parameter $M/r$. This allows us to quantify a PN order by dropping higher order terms and plotting the differences.

For the following analysis, we are restricting to a non-spinning binary in a quasi-circular orbit, that starts in the same initial position for the evolutions of the binary with both Taylor T4 and our direct EOM integration.

We give the results of these different PN orders in Fig. 8 by calculating

[TABLE]

We can see immediately that they follow a strict hierarchy of decreasing relative difference as the PN orders are increased, i.e., the differences $\delta_{\rm 3PN,iPN}$ decrease with increasing PN order $i$, which is a good initial sanity check.

We similarly define the difference stemming from using different definitions of the frequency itself,

[TABLE]

and note that the difference $\delta^{\rm 3PN,n\times v/r}$ accounts for roughly the same order error as the 3PN term (which is reassuring, since the orbital frequency is given to 3PN order in the alternate PN definition stated above). There is a curious thing happening at the initial separation $r=50M$, where this difference $\delta^{\rm 3PN,n\times v/r}$ is larger than the 3PN error. This is explained in Fig. 9, where it becomes apparent that the 3PN term has a zero crossing, and thus dips lower than the difference from the frequency definition. This zero crossing is directly dependent on the gauge term in Eq. (IV.4), and can be shifted at will for different choices of $r_{0}^{\prime}$.

The relative difference of separate PN orders gives a measure of the error in terms of the PN order. We can use this to quantify any comparison in terms of the orbital frequency PN order, and to check the PN scaling of the discrepancies of Taylor T4 and the EOM methods, $\delta_{\rm 3PN,T4}$. What we find is a strong scaling as a function of separation. When the binary is closely separated ($r\leq 12M$), the discrepancy is greater than $1\%$ after an initial discrepancy below the 3PN line. As the binary inspirals, the discrepancy grows, crossing the 2PN line after a short amount of evolution ($t\sim 3000M$), and crosses the 1PN line after ($t\sim 10800M$). When we look at larger separations of $50M$, we see again this same trend of a sharply scaling function of separation, but in this case, it stays below the 2PN line. When we look at the $100M$ case, we find that the discrepancy between the Taylor T4 and EOM codes again stays below the 2PN line, but with a lower relative difference.

Of course, to demonstrate that this is indeed a PN scaling and not another effect, we need to show that the relative differences between the PN orders scale as the proper powers of $r^{s}$ when we look over a large range of separations. We do this in Fig. 9. This figure shows the relative differences scaling as a function of separation. In each case, we fit a power law to the line, and record the slope. The slopes that we find scale with orbital separation $r^{s}$, where the slope $s$ is measured for the individual relative errors. We find that the 1PN error scales with separation as $s\approx-0.976$, the 2PN error scales with separation as $s\approx-1.9995$, the 3PN error scales with separation as $s\approx-4.002$ at far separations ($r\approx 1000M$), and $s\approx-2.7577$ at small separations ($r\sim 12M$). The 3PN error does not scale as $s\approx-3.0$ because of the gauge dependent logarithm term. When looking at separations that are even farther out than $r\approx 100M$, like $r\sim 1000000M$, we find that the slope approaches the correct PN scaling of $-3$. The scaling of the different PN definitions is also tracked here, and while there is not a clear power law trend on this plot, we do note that the definitional discrepancy falls at or below the 3PN line (except for the special case around $50M$), so we can say from this that either definition is acceptable in the EOM code. For our purposes, we will use the geometric definition of $\Omega=|\hat{n}\times\vec{v}|/r$. When we measure the scaling of the Taylor T4–EOM relative error in Fig. 9, we find it is $s\approx-4.545$. We do not attribute this discrepancy to the PN order since this does not obey any obvious PN scaling.

This begs the question of what is causing this rapid drop in relative difference as the separation is increased.

To investigate this, we need to consider what kinds of effects could be at play: this could be a numeric effect (e.g., one of the codes is not calculating to the requisite precision, which is causing a discrepancy at close separations where time steps are smaller); an eccentricity that is deviating the orbit from quasi-circularity in the EOM code; or an effect of the adiabatic approximation breaking down.

To test whether or not this is a numeric effect, we doubled the precision of both codes and re-ran the test at $20M$ again, with the same results. Therefore, we conclude that this is not a numeric issue. To test whether this is an eccentric effect, we ran the EOM code with different levels of the eccentricity remover (which should, in principle, solve the quasi-circularity problem if indeed it is one), and checked the relative differences. We find preliminarily that the residual eccentricity does scale the relative difference of the orbital frequencies, not quite one-to-one, and is removed by the iterative eccentricity removal procedure.

The final plausible possibility is the adiabatic approximation breaking down at close separations, which leads to a high discrepancy between the two methods. This is the most challenging possibility to eliminate. The best method for tracing this is to orbit average all of the EOM code, essentially making it adiabatic. We have done this for the lowest order Peters-Mathews test (0PN (Newtonian) conservative terms, with a 2.5PN radiation reaction term added), and find that this does not seem to affect the evolution at lowest order. Of course, to show that this is indeed unaffected in the general sense, we need to go beyond the leading order evolutions and demonstrate this for 3.5PN.

IV.5 Waveform Comparisons for Aligned Spins and different Mass Ratios

With the orbital frequency analysis concluded, we now turn to calculating the GW waveform overlaps between Taylor T4 and the EOM methods. We pick several fiducial separations, mass ratios, and spins (both aligned and anti-aligned, denoted by $\chi_{1}$ and $\chi_{2}$). The initial separations for which we choose to calculate the overlaps are $100M$, $50M$, and $20M$, with mass ratios $q=m_{2}/m_{1}$ ($m_{2}\geq m_{1}$) of $1$, $2$, $10$, and $100$. The spin parameters that we pick are $\chi_{1}=\{0,0.3,0.6,0.9\}$ and $\chi_{2}=\{0,0.3,0.6,0.9,-0.3,-0.5,-0.6,-0.9\}$, at separations of $50M$ and $100M$.

We calculate the overlap and maximize over time and phase (e.g., Ref. Buonanno et al. (2009)), by calculating the inner product in the frequency domain

[TABLE]

where $S_{n}$ is a noise power spectrum density of a detector 333Here, we used LIGO’s target sensitivity curve, Zero-Detuned High-Power (v2) Shoemaker (2009); it has since been superseded by v5 Barsotti et al. (2018).. The maximization over time is handled by maximizing over $t_{c}$, and the phase maximization is handled by the shifting of $h_{1}$ in the frequency domain. The low frequency cutoff on the integration is set to a reasonable frequency for a detector (for this analysis we set the low frequency cutoff to $10$ Hz, which is a reasonable if a bit ambitious lower frequency bound for LIGO). The high frequency cutoff is not set, to capture the maximum overlap of the waveform if the endpoint is not exactly set to the same frequency.

We then normalize (using the euclidean norm) over $h_{1}$ and $h_{2}$ to obtain the overlap:

[TABLE]

The results are tabulated in Table 1 (see also Figs. 10 and 11). We keep all of the parameters that we used to calculate the overlaps in the table: the initial and final separations, the mass ratio, the aligned dimensionless spin values $\chi_{1}$ and $\chi_{2}$, the total simulation time in units of $M$, the number of orbits the waveform spanned, the time step of the overlap calculation, and finally the maximized overlap for both the Taylor T4 to EOM comparison at a consistent PN order (3.5PN), and also at the highest T4 order (6PN).

We see that the overlap is a strong function of the orbital separation: as the separation increases, the overlap increases from a bad overlap at $20M$ of only $O\sim 0.8$ to an overlap of $O\sim 0.99999$ at $100M$. In addition, as the mass ratio increases, the overlap also increases. For example, at a separation of $50M$, holding the spins of the individual BHs to zero, we increase in overlap from $O\sim 0.98$ at $q=1$ to $O\sim 0.99999$ at $q=100$.

When we move to explore the spin parameter space, we hold the mass ratio fixed and compute the overlaps for separations of $100M$ and $50M$. As the dimensionless spin parameter increases in value to higher positive $\chi$ effective ($\chi$ effective is the spin values projected along the orbital angular momentum), the overlap goes down, at $50M$, from $O\sim 0.995$ at spins of zero, to $O\sim 0.99$ at a high positive $\chi$ effective.

When the spins are anti aligned with each other, keeping the effective spin zero, the overlap stays fairly constant, but increases slightly, with the $\chi_{1}=0.9,\,\chi_{2}=-0.9$ having an overlap of $O\sim 0.999$.

A final discussion point is to highlight the effects of PN at higher order with Taylor T4, specifically when we add higher order radiation reaction flux terms. The high order T4 generically differs in overlap from the consistent order T4 by $10^{-8}$ at $100M$, $10^{-5}$ at $50M$, and $10^{-3}$ at $20M$. Though the overlap difference increases as the separation decreases, the PN effects at 4PN do not account for the discrepancy of the overlaps between Taylor T4 and the direct integration EOM.

In addition to the results we tabulate in Table 1, we perform stability tests on the overlap between Taylor T4 and the direct integration EOM by inputting a small $x$-component perturbation to the two spin vectors, which will cause a small amount of spin precession, hence marked ${\rm EOM_{P}}$. Specifically, we give the dimensionless spin vectors $\chi_{1}=(10^{-4},0,0.3)$, and $\chi_{2}=(-10^{-4},0,0.3)$, and we run the same overlaps with the consistent T4 method, and also run the overlap with the EOM code with no $x$-component perturbation to the spins. The overlaps that we obtain are $O_{{\rm T4}-{\rm EOM}}=0.983693663044$, which is exactly the overlap that we obtain when running without the $x$ perturbation. This is easily verifiable by redoing the overlap analysis with the EOM code with and without the perturbation. We obtain an overlap of $O_{{\rm EOM}-{\rm EOM_{P}}}=1.0$, which clearly shows that a perturbation to the spin directions do not affect the overlaps.

V Generic Testbed for Other Approximants: Eccentricity Definitions

This direct EOM that we present is generic and can be applied broadly to many different physical systems, and can be used as a benchmark test for other formalisms. To illustrate this, we apply to a simple test of orbital eccentricity put forward recently Loutrel et al. (2019).

Eccentricity in general relativity is difficult and unintuitive to define Loutrel et al. (2019); Chandrasekhar (1992); Yunes et al. (2009); Memmesheimer et al. (2004). As such, most Newtonian definitions are not sufficient, and can give wildly different results. To illustrate this, we are going to take two different parameterizations of the Newtonian eccentricity, the Runge-Lenz vector and a low eccentricity definition of $e$ used in NR.

The Newtonian Runge-Lenz vector is defined as Poisson and Will (2014):

[TABLE]

which, for a Newtonian orbit, is a constant of the motion, and implies that the eccentricity itself is another constant of the motion. To obtain the scalar eccentricity, we simply take the magnitude of this vector. This measure is shown in Fig. 12, in a non-spinning binary system at a low PN order at an initial separation of $20M$, as it evolves down to $10M$.

To contrast this, we will use another definition of eccentricity taken from numerical relativity initial data, using the radial acceleration $\ddot{r}$, following Refs. Campanelli et al. (2009) and Husa et al. (2008). Under the double assumptions of low eccentricity initial data and adiabatic frequency $\Omega$, we can estimate $r(t)$ and its derivatives (to leading order in $e$) by

[TABLE]

and then define

[TABLE]

so that the eccentricity is given by

[TABLE]

$A(t)$ is also shown in Fig. 12. We note that both measures, as shown in Fig. 12, oscillate – and for $A(t)$, we expect the amplitude of these oscillations to be the actual eccentricity. We thus try to complement $A(t)\sim\sin(\Omega t)$ with the corresponding $\cos(\Omega t)$, using

[TABLE]

and then should have

[TABLE]

However, in Fig. 12 we see that

The curve for $A(t)$ drifts downwards, and at late times oscillates entirely below 0. 2. 2.

The curve for $\sqrt{A^{2}+B}$ coincides (to $\sim 1\%$) with the Runge-Lenz eccentricity measure. 3. 3.

The curve for $\sqrt{A^{2}+B}$, which is supposed to serve as the non-oscillatory envelope of $A(t)$, does oscillate just like $A(t)$; furthermore, it both starts off above and drifts upwards faster than $A(t)$. 4. 4.

The amplitudes of oscillations themselves shrink for both curves (which is what we would have expected the full eccentricity to do).

As all “physical” eccentricities remain very low ($e\sim 10^{-3}$), this suggests that the adiabatic approximation is breaking down; we shall analyze just how:

We no longer take $\Omega$ as a constant, but rather as $\Omega(t)$; we model its leading order behavior as $\Omega(t)\propto\left(t_{c}\!-\!t\right)^{-3/8}$, with $t_{c}$ the projected leading order time of coalescence Abbott et al. (2017b). Then $\dot{\Omega}\sim\Omega/\left(t_{c}\!-\!t\right)$ and $\ddot{\Omega}\sim\Omega/\left(t_{c}\!-\!t\right)^{2}$. We note that for the evolution shown in Fig. 12, $\left(t_{c}\!-\!t\right)$ decreases from $\sim 10^{4}$ to $\sim 10^{3}$, while $\Omega$ changes over a few $\sim 10^{-2}$.

We then see that the derivatives of

[TABLE]

more accurately behave like

[TABLE]

with numerical factors of order unity, $\alpha$ and $\beta$.

For $\ddot{r}$, when $\left(t_{c}\!-\!t\right)\sim 10^{4}$ the drift term is $\sim\left(t_{c}\!-\!t\right)^{-2}\Omega^{-2}\sim(10^{4})^{-2}(10^{-2})^{-2}=10^{-4}$, and so is smaller than the oscillating eccentricity term $e$; however when $\left(t_{c}\!-\!t\right)\sim 10^{3}$, the drift term grows to $\sim(10^{3})^{-2}(10^{-2})^{-2}=10^{-2}$, and so at late times dominates over the oscillations. This explains why the term $A(t)$ oscillates with the amplitude of the eccentricity early on, but is later dominated by the drift.

For $\dot{r}$ the situation is different: even at early times, $\left(t_{c}\!-\!t\right)\sim 10^{4}$, the drift term, $\sim\left(t_{c}\!-\!t\right)^{-1}\Omega^{-1}\sim(10^{4})^{-1}(10^{-2})^{-1}=10^{-2}$, is already larger than the eccentricity $e$, and so the behavior of $B(t)$ is governed by the drift, rather than the eccentricity oscillations.

As the Runge-Lenz definition Eq. (27) is essentially given by the velocity components, it adopts these features of $\dot{r}$, and measures the drift associated to non-adiabaticity, rather than the eccentricity itself; the eccentricity is given only by the oscillations.

VI Discussion

We have developed and constructed a direct integration of the PN EOM in the harmonic gauge for the construction of eccentric BBHs with arbitrary spins. This work represents a step forward for the modeling of these systems by extending known methods (such as Taylor T4, and the current LIGO/VIRGO collaboration methods) to binaries that evolve the spin procession equations and orbital motion equations for a fully generic GW waveform. This formalism is not limited by eccentricity as the post-Keplerian eccentricity expansion models (e.g., Refs. Gopakumar and Iyer (2002); Yunes et al. (2009); Huerta et al. (2014); Tanay et al. (2016); Hinder et al. (2018)), and is capable of handling eccentricity along with spin precession.

We test the validity of this method by comparing to known results from the Taylor T4 method. In particular, we look at the orbital frequencies produced by systems with the same initial conditions, and compare. The results we find is that the relative difference in the orbital frequencies is a strong function of separation, that does not scale with a high order PN effect, and scales as a function of separation as $r^{s}$ with $s\sim-4.5$.

This effect is not overly worrying, even though we cannot ascribe exactly the cause of this scaling. This is most likely the fundamental difference between Taylor T4 and our method. The various PN dynamics methods do not produce the same results even among the various Taylor T approximants Buonanno et al. (2009). We rule out a numeric issue, and an eccentric effect as the sole cause. We have confirmed that the adiabatic approximation does not affect the evolution of the binary at low order, and are planning on verifying this for the full 3.5PN dynamics in future work.

This begs the question of what is causing this rapid drop in relative difference as the separation is increased.

To investigate this, we need to consider what kinds of effects could be at play: this could be a numeric effect (e.g., one of the codes is not calculating to the requisite precision, which is causing a discrepancy at close separations where time steps are smaller); an eccentricity that is deviating the orbit from quasi-circularity in the EOM code; or an effect of the adiabatic approximation breaking down.

To test whether or not this is a numeric effect, we doubled the precision of both codes and re-ran the test at $20M$ again, with the same results. Therefore, we conclude that this is not a numeric issue. This is an eccentric effect, we ran the EOM code with different levels of the eccentricity remover (which should, in principle, solve the quasi-circularity problem if indeed it is one), and checked the relative differences. We find preliminarily that the residual eccentricity does scale the relative difference of the orbital frequencies, not quite one-to-one, and is removed by the iterative eccentricity removal procedure.

The final plausible possibility is the adiabatic approximation breaking down at close separations, which leads to a high discrepancy between the two methods. This is the most challenging possibility to eliminate. The best method for tracing this is to orbit average all of the EOM code, essentially making it adiabatic. We have done this for the lowest order Peters-Mathews test (0PN conservative terms, with a 2.5PN radiation reaction term added), and find that this does not seem to affect the evolution at lowest order. Of course, to show that this is indeed unaffected in the general sense, we need to go beyond the leading order evolutions and demonstrate this for 3.5PN.

We also compare the GW waveforms of Taylor T4 to the EOM methods, maximizing over time and phase, to give a quantification of the overlap. The results that we find are consistent with the findings of the orbital frequency analysis: the overlap is a strong function of the orbital separation from $20M$ to $100M$. The overlaps increase strongly as a function of mass ratio, with the best overlaps of the waveforms we ran being at a mass ratio $q=100$. In addition, we performed spinning waveform overlaps to test the validity of the spins in the EOM code. We find that the spins do not modify the overlap of the waveform when the effective spin of the binary remains zero, and drops the overlap by about a percent when the spins are aligned and low spins ($\chi_{1}=\chi_{2}=0.3$), to about two percent when the spins are moderately spinning ($\chi_{1}=\chi_{2}=0.6$). The overlap then gets better as the spin gets stronger. We also tested the stability of the code to small perturbations in spin by giving the EOM code a small non-zero spin unaligned with the orbital angular momentum, and found no effect on the overlap.

Finally, we elucidated a fiducial binary system, with realistic parameters drawn from galactic binary simulations, to demonstrate the flexibility and power of this direct integration EOM code. We used initial conditions such that the binary would be in the very low frequency end of the LIGO band at periapsis, and let the binary evolve for 200 orbits. We then output the binary orbital trajectories, spins, and velocities. We recover the orbital plane precession of the binary due to spin-orbit coupling, the spin-spin precession of the individual spins and spin totals, the eccentricity reduction in the binary, and calculate the initial and final eccentricity of the binary using our geometric definition. We take the orbital quantities and calculate the waveform, recovering the eccentric signal imprinted on the outgoing gravitational radiation, estimate the periapsis and apoapsis frequencies of this radiation, and show the spin precession modulation that is also imparted on the binary. We do this for an optimally oriented binary, but leave the code generic so that any binary orientation can be used, providing us with a fully generic, precessing, eccentric binary GW waveform.

We hope to apply this EOM code to parameter estimation for future gravitational wave detections. We currently have the code implemented in Mathematica, and the analysis takes approximately a minute to complete. This is far too long for parameter estimation purposes currently, and we are planning on focusing on this in the future.

There are several relevant detectors for this fiducial source. The periapsis frequency passage is at the threshold of detectability for LIGO/VIRGO at design sensitivity Moore et al. (2015); Martynov et al. (2016); Barsotti et al. (2018); it will fall into the band for future LIGO upgrades such as A+; and third generation ground-based detectors will have both the periapsis and apoapsis frequencies in band.

The precessional frequency will be detectable by planned space based detectors such as the LISA mission Amaro-Seoane et al. (2017), though the source outlined above will be too weak for detection, these frequencies are in band if a nearby binary happens to have these parameters.

We also note that there is a planned Chinese space-borne GW detector in the millihertz frequencies, TianQin Luo et al. (2016), and a planned Japanese space-born detector, DECihertz laser Interferometer Gravitational wave Observatory (DECIGO/B-DECIGO) Seto et al. (2001); Sato et al. (2017) in between the frequency ranges of ground based and LISA detectors. These will make it possible to detect a binary’s dynamics by multiband GW astronomy.

Acknowledgments

B. I. and M. C. received support from National Science Foundation (NSF) grants AST-1028087, AST-1516150, PHY-1305730, and PHY-1707946. O. B. is supported through the Frontiers in Gravitational Wave Astrophysics Initiative at the RIT’s Center for Computational Relativity and Gravitation. He also acknowledges support from NSF grant PHY-1607520. H.N. acknowledges support from JSPS KAKENHI Grant No. JP16K05347 and No. JP17H06358. Computational resources were provided by the BlueSky Cluster at Rochester Institute of Technology. The BlueSky cluster was supported by NSF grants AST-1028087, PHY-0722703 and PHY-1229173.

Appendix A Another expression for PN EOM

Equation (4) can be reformulated into the alternate form:

[TABLE]

We note that while the unit vector $\hat{\ell}$ is orthogonal to $\hat{n}$ and $\vec{v}$, the triad $\{\hat{n},\vec{v},\hat{\ell}\}$ does not form an orthogonal basis, because $\hat{n}$ and $\vec{v}$ will, in general, not be orthogonal to one another (although they will be approximately orthogonal in the special case of quasi-circular orbits). This equation resembles Eq. (129) in Ref. Blanchet (2014), except for the additional $\tilde{\mathcal{C}}$ coefficient. It is noted that the term with coefficient $\tilde{\mathcal{C}}$ represents a component of the acceleration directed out of the instantaneous orbital plane spanned by $\hat{n}$ and $\vec{v}$. It therefore gives rise to precession of the orbital plane and correspondingly ought to vanish identically when spins are aligned or anti-aligned with the orbital angular momentum.

To obtain Eq. (38), we may expand each cross-product in Eq. (4) by using

[TABLE]

where $\Omega$ is the instantaneous orbital frequency, $\hat{\lambda}\equiv\hat{\ell}\times\hat{n}$, and we have introduced the notation

[TABLE]

and

[TABLE]

The PN expansion for $\Omega$ is only calculated in the case of quasi-circular orbits. To avoid assuming quasi-circular orbits, we may replace $r\Omega$ in the above identities with the (exact) identity

[TABLE]

After using Eqs. (39)–(42) to eliminate the cross-products in Eq. (4) and recollecting terms, one obtains an acceleration equation in the form of Eq. (38), with

[TABLE]

It is noted that this last equation is used in the case of quasi-circular orbits to compute the orbital plane precession frequency $\varpi$. Also we note that in the case when spins are aligned or anti-aligned, $\mathcal{C}$ vanishes identically (since in that case $S_{\mathrm{n}}=\tilde{S}_{\lambda}=0$ and $\Sigma_{\mathrm{n}}=\tilde{\Sigma}_{\lambda}=0$).

Appendix B Spin Precession Equations

We overview the spin precession equations for two bodies with spins that are unaligned with the orbital angular momentum. As mentioned previously, these equations are needed to complete the EOM.

There is no unique definition of center of mass in a relativistic theory. The notion of spin inherits this ambiguity, which gives rise to non-physical degrees of freedom that need to be controlled. This is done by imposing a “spin supplementary condition” (or SSC) that eliminates the non-physical degrees of freedom. The SSC is imposed on the spin tensor, out of which are constructed a spin 4-vector and eventually a spin 3-vector, with just three physical degrees of freedom. Following Refs. Faye et al. (2006); Blanchet et al. (2006); Marsat et al. (2013); Bohe et al. (2013), we use the Tulczyjew SSC Kidder (1995).

[TABLE]

We define spin 3-vectors with conserved norm by $\vec{S}_{\mathrm{a}}$, where $\mathrm{a}=1,\,2$ is the label of the particles. In terms of these conserved norm spin vectors, the precession equations are written as

[TABLE]

where $\vec{\Omega}_{\mathrm{a}}$ with $\mathrm{a}=1,\,2$ are the precession vectors for mass 1 and 2, respectively. Each precession vector can be decomposed as

[TABLE]

The precession vectors are only expanded to quadratic-in-spin order, $O{(SS)}$, because they get multiplied by a spin vector in the precession equation. Terms in the precession vectors at $O{(SSS)}$ would contribute at $O{(SSSS)}$ in the precession equations, which is beyond the scope of this work. Each contribution to the precession vector is then decomposed into a PN expansion of the form,

[TABLE]

Explicit expressions for $\vec{\Omega}_{1}$ and $\vec{\Omega}_{2}$ are given in Ref. Bohe et al. (2013).

Appendix C The Taylor T4 Method

Taylor T4 is an adiabatic approximation 444An adiabatic approximation means that we do not consider the change in any quantity that is smaller than one orbit, i.e., the inspiral of the BHs will not affect the orbit. As such, since we consider spin precession on a different timescale than the orbital timescale, these approximations do not have spin precession built into them, and are only for aligned spins. that evolves the phase and frequency of the BBH. There are two parts to this method. The first is the phase evolution of the binary, which we will detail here, and the second is the calculation of the gravitational radiation from the orbital phasing of the binary.

All of the formulas of energy and flux can be written in terms of the frequency variable $v\equiv(M\mathrm{d}\phi/\mathrm{d}t)^{1/3}$, and we start with the basic definition of energy conservation, namely:

[TABLE]

where the time rate of change of the energy is equal to the flux leaving the system plus the mass rate of change due to GW absorption. In practice, since this BH absorption is a relatively tiny contribution Isoyama and Nakano (2018) (although 2.5PN order relative to the leading flux), we neglect it for the rest of this analysis, bringing our conservation statement to

[TABLE]

We start with the energy balance equation (50), integrate it to find $v(t)$, and therefore $\Omega(t)$, since $\mathrm{d}\phi/\mathrm{d}t=\Omega$. Taylor T4 explicitly takes the rational fraction,

[TABLE]

along with

[TABLE]

re-expands the fraction as a consistent PN series, then numerically integrates to obtain $v(t)$ and $\phi(t)$.

Full expressions of $E(v)$, ${\cal F}(v)$, and $\mathrm{d}E/\mathrm{d}v$ that we use are given in Ref. Ajith et al. (2007). These equations assume a quasi-circular orbit for the trajectory. Taylor T4 approximants integrate the energy balance equations and have been thoroughly explored and developed. Further approximations can be constructed in the frequency domain (the Taylor F-approximants). These methods have been compared, see Ref. Buonanno et al. (2009).

Appendix D Tabulated results on overlap calculations