Post-Newtonian Quasicircular Initial Orbits for Numerical Relativity

TL;DR

This paper develops a method using post-Newtonian approximations to set initial conditions for black hole binary simulations that result in low-eccentricity inspirals, improving the accuracy of numerical relativity models.

Contribution

It introduces a new approach combining 3PN Hamiltonian and high-order radiation reaction expressions to generate initial data with minimal eccentricity for black hole binaries.

Findings

Using 3PN Hamiltonian and high-order radiation reaction yields the lowest eccentricity.

The method surpasses previous inspiral data accuracy based on 3PN and 4PN evolutions.

Applicable to spinning, nonprecessing, unequal mass binaries.

Abstract

We use post-Newtonian (PN) approximations to determine the initial orbital and spin parameters of black hole binaries that lead to low-eccentricity inspirals when evolved with numerical relativity techniques. In particular, we seek initial configurations that lead to very small eccentricities at small separations, as is expected for astrophysical systems. We consider three cases: (i) quasicircular orbits with no radial velocity, (ii) quasicircular orbits with an initial radial velocity determined by radiation reaction, and (iii) parameters obtained form evolution of the PN equations of motion from much larger separations. We study seven cases of spinning, nonprecessing, unequal mass binaries. We then use several definitions of the eccentricity, based on orbital separations and waveform phase and amplitude, and find that using the complete 3PN Hamiltonian for quasicircular orbits to obtain the tangential orbital momentum, and using the highest-known-order radiation reaction expressions to obtain the radial momentum, leads to the lowest eccentricity. The accuracy of this method even exceeds that of inspiral data based on 3PN and 4PN evolutions.

Tables

Table 1. Table 1: Initial data parameters for the quasi-circular configurations with a smaller mass black hole (labeled 1), and a larger mass spinning black hole (labeled 2). The 40 configurations are split into 8 families (labeled #1 – #8) of fixed mass ratio, spins, and initial separations. Within each family, 5 different choices for the QC momentum parameters are used (id0 – id4). For a given family, the parameters that remain fixed are given for id0 only. The mass ratio q = m H ​ 1 / m H ​ 2 𝑞 subscript 𝑚 𝐻 1 subscript 𝑚 𝐻 2 q=m_{H1}/m_{H2} is given in terms of the Christodoulou masses of each BH. The parameters m 1 , 2 p subscript superscript 𝑚 𝑝 1 2 m^{p}_{1,2} are chosen such that m H ​ 1 subscript 𝑚 𝐻 1 m_{H1} and m H ​ 2 subscript 𝑚 𝐻 2 m_{H2} agree with the PN parameters m 1 , 2 subscript 𝑚 1 2 m_{1,2} to at least 5 significant digits (this also means that m H ​ 1 / M + m H ​ 2 / M = 1 + 𝒪 ​ ( 10 − 5 ) subscript 𝑚 𝐻 1 𝑀 subscript 𝑚 𝐻 2 𝑀 1 𝒪 superscript 10 5 m_{H1}/M+m_{H2}/M=1+{\cal O}(10^{-5}) ). The dimensionless spins χ 1 subscript 𝜒 1 \chi_{1} and χ 2 subscript 𝜒 2 \chi_{2} were obtained using the IH formalism. The orbital frequency Ω ​ orb . Ω orb \Omega{\rm orb.} was measured directly from the numerical orbital trajectory at t = 200 ​ M 𝑡 200 𝑀 t=200M . The initial puncture locations are r → 1 = ( x 1 , 0 , 0 ) subscript → 𝑟 1 subscript 𝑥 1 0 0 \vec{r}_{1}=(x_{1},0,0) and r → 2 = ( x 2 , 0 , 0 ) subscript → 𝑟 2 subscript 𝑥 2 0 0 \vec{r}_{2}=(x_{2},0,0) with mass parameters m p / M superscript 𝑚 𝑝 𝑀 m^{p}/M , momentum P → / M = ∓ ( P r , − P t , 0 ) → 𝑃 𝑀 minus-or-plus subscript 𝑃 𝑟 subscript 𝑃 𝑡 0 \vec{P}/M=\mp(P_{r},-P_{t},0) , spins S → i = ( 0 , 0 , S i ) subscript → 𝑆 𝑖 0 0 subscript 𝑆 𝑖 \vec{S}_{i}=(0,0,S_{i}) . The last column gives the ADM mass of each configuration. Note that in the table P r subscript 𝑃 𝑟 P_{r} is multiplied by 10 3 superscript 10 3 10^{3} and P t subscript 𝑃 𝑡 P_{t} by 10.

Run ID	$q$	${\chi }_1$	${\chi }_2$	$M{\Omega }_{\mathrm{orb}.}$	$x_1/M$	$x_2/M$	$S_1/M^2$	$S_2/M^2$	$m_1^p/M$	$m_2^p/M$	${10}^3\cdot P_r/M$	$10\cdot P_t/M$	$M_{\mathrm{ADM}}$
#1id0	0.333	0.8000	-0.5000	0.0275	-8.2500	2.7500	0.0500	-0.2812	0.1491	0.6580	0.0000	0.6841	0.993076
#1id1	$\cdots$	$\cdots$	$\cdots$	0.0258	$\cdots$	$\cdots$	$\cdots$	$\cdots$	0.1491	0.6580	0.4448	0.6880	0.993182
#1id2	$\cdots$	$\cdots$	$\cdots$	0.0259	$\cdots$	$\cdots$	$\cdots$	$\cdots$	0.1491	0.6580	0.0000	0.6880	0.993182
#1id3	$\cdots$	$\cdots$	$\cdots$	0.0259	$\cdots$	$\cdots$	$\cdots$	$\cdots$	0.1491	0.6580	0.4490	0.6879	0.993179
#1id4	$\cdots$	$\cdots$	$\cdots$	0.0257	$\cdots$	$\cdots$	$\cdots$	$\cdots$	0.1491	0.6580	0.4551	0.6882	0.993188
#2id0	0.333	-0.8000	0.5000	0.0277	-7.8750	2.6250	-0.0500	0.2812	0.1489	0.6576	0.0000	0.6869	0.992637
#2id1	$\cdots$	$\cdots$	$\cdots$	0.0270	$\cdots$	$\cdots$	$\cdots$	$\cdots$	0.1489	0.6576	0.4467	0.6888	0.992690
#2id2	$\cdots$	$\cdots$	$\cdots$	0.0272	$\cdots$	$\cdots$	$\cdots$	$\cdots$	0.1489	0.6576	0.0000	0.6888	0.992689
#2id3	$\cdots$	$\cdots$	$\cdots$	0.0268	$\cdots$	$\cdots$	$\cdots$	$\cdots$	0.1489	0.6576	0.4521	0.6893	0.992705
#2id4	$\cdots$	$\cdots$	$\cdots$	0.0254	$\cdots$	$\cdots$	$\cdots$	$\cdots$	0.1489	0.6576	0.4831	0.6933	0.992818
#3id0	0.250	0.8000	-0.8000	0.0316	-8.4880	2.1220	0.0320	-0.5120	0.1186	0.4928	0.0000	0.6114	0.994076
#3id1	$\cdots$	$\cdots$	$\cdots$	0.0277	$\cdots$	$\cdots$	$\cdots$	$\cdots$	0.1186	0.4927	0.4196	0.6176	0.994249
#3id2	$\cdots$	$\cdots$	$\cdots$	0.0278	$\cdots$	$\cdots$	$\cdots$	$\cdots$	0.1186	0.4927	0.0000	0.6176	0.994249
#3id3	$\cdots$	$\cdots$	$\cdots$	0.0283	$\cdots$	$\cdots$	$\cdots$	$\cdots$	0.1186	0.4927	0.4302	0.6165	0.994220
#3id4	$\cdots$	$\cdots$	$\cdots$	0.0281	$\cdots$	$\cdots$	$\cdots$	$\cdots$	0.1186	0.4927	0.4365	0.6169	0.994230
#4id0	1.000	-0.8000	-0.8000	0.0314	-5.5000	5.5000	-0.2000	-0.2000	0.3029	0.3029	0.0000	0.9411	0.991047
#4id1	$\cdots$	$\cdots$	$\cdots$	0.0269	$\cdots$	$\cdots$	$\cdots$	$\cdots$	0.3029	0.3029	0.9437	0.9519	0.991351
#4id2	$\cdots$	$\cdots$	$\cdots$	0.0268	$\cdots$	$\cdots$	$\cdots$	$\cdots$	0.3029	0.3029	0.0000	0.9519	0.991349
#4id3	$\cdots$	$\cdots$	$\cdots$	0.0299	$\cdots$	$\cdots$	$\cdots$	$\cdots$	0.3029	0.3029	0.8941	0.9444	0.991140
#4id4	$\cdots$	$\cdots$	$\cdots$	0.0256	$\cdots$	$\cdots$	$\cdots$	$\cdots$	0.3028	0.3028	1.0522	0.9554	0.991449
#5id0	0.750	-0.8500	0.6375	0.0270	-6.2857	4.7143	-0.1561	0.2082	0.2192	0.4479	0.0000	0.8779	0.990803
#5id1	$\cdots$	$\cdots$	$\cdots$	0.0255	$\cdots$	$\cdots$	$\cdots$	$\cdots$	0.2192	0.4479	0.6879	0.8824	0.990926
#5id2	$\cdots$	$\cdots$	$\cdots$	0.0257	$\cdots$	$\cdots$	$\cdots$	$\cdots$	0.2192	0.4479	0.0000	0.8824	0.990926
#5id3	$\cdots$	$\cdots$	$\cdots$	0.0258	$\cdots$	$\cdots$	$\cdots$	$\cdots$	0.2192	0.4479	0.6788	0.8815	0.990902
#5id4	$\cdots$	$\cdots$	$\cdots$	0.0244	$\cdots$	$\cdots$	$\cdots$	$\cdots$	0.2192	0.4479	0.7284	0.8863	0.991034
#6id0	0.700	0.0000	0.0000	0.0234	-7.0588	4.9412	0.0000	0.0000	0.4002	0.5771	0.0000	0.8219	0.991345
#6id1	$\cdots$	$\cdots$	$\cdots$	0.0226	$\cdots$	$\cdots$	$\cdots$	$\cdots$	0.4002	0.5771	0.5066	0.8246	0.991418
#6id2	$\cdots$	$\cdots$	$\cdots$	0.0226	$\cdots$	$\cdots$	$\cdots$	$\cdots$	0.4002	0.5771	0.0000	0.8246	0.991417
#6id3	$\cdots$	$\cdots$	$\cdots$	0.0228	$\cdots$	$\cdots$	$\cdots$	$\cdots$	0.4002	0.5771	0.5031	0.8241	0.991405
#6id4	$\cdots$	$\cdots$	$\cdots$	0.0221	$\cdots$	$\cdots$	$\cdots$	$\cdots$	0.4002	0.5771	0.5227	0.8266	0.991471
#7id0	0.250	0.0000	0.0000	0.0277	-8.4000	2.1000	0.0000	0.0000	0.1909	0.7920	0.0000	0.5959	0.993687
#7id1	$\cdots$	$\cdots$	$\cdots$	0.0268	$\cdots$	$\cdots$	$\cdots$	$\cdots$	0.1909	0.7919	0.3541	0.5977	0.993737
#7id2	$\cdots$	$\cdots$	$\cdots$	0.0269	$\cdots$	$\cdots$	$\cdots$	$\cdots$	0.1909	0.7919	0.0000	0.5977	0.993736
#7id3	$\cdots$	$\cdots$	$\cdots$	0.0265	$\cdots$	$\cdots$	$\cdots$	$\cdots$	0.1909	0.7919	0.3626	0.5983	0.993754
#7id4	$\cdots$	$\cdots$	$\cdots$	0.0261	$\cdots$	$\cdots$	$\cdots$	$\cdots$	0.1909	0.7919	0.3705	0.5994	0.993783
#8id0	0.333	0.8000	0.8000	0.0306	-7.1250	2.3750	0.0500	0.4500	0.1483	0.4601	0.0000	0.7049	0.991784
#8id1	$\cdots$	$\cdots$	$\cdots$	0.0302	$\cdots$	$\cdots$	$\cdots$	$\cdots$	0.1483	0.4601	0.5315	0.7047	0.991779
#8id2	$\cdots$	$\cdots$	$\cdots$	0.0306	$\cdots$	$\cdots$	$\cdots$	$\cdots$	0.1483	0.4601	0.0000	0.7047	0.991778
#8id3	$\cdots$	$\cdots$	$\cdots$	0.0293	$\cdots$	$\cdots$	$\cdots$	$\cdots$	0.1483	0.4601	0.5616	0.7082	0.991880
#8id4	$\cdots$	$\cdots$	$\cdots$	0.0297	$\cdots$	$\cdots$	$\cdots$	$\cdots$	0.1483	0.4601	0.5333	0.7068	0.991841

Table 2. Table 2: The measured eccentricities using e D ​ 2 subscript 𝑒 𝐷 2 e_{D2} and e A subscript 𝑒 𝐴 e_{A} for all 40 runs, as well as the relative difference between these two measures. Note that the relative differences are largest for small eccentricities due the waveform and trajectory noise dominating the small eccentricity effects. The eccentricities are given in the form ( e D ​ 2 , e A ) subscript 𝑒 𝐷 2 subscript 𝑒 𝐴 (e_{D2},e_{A}) [rel. difference] in the table below.

run#	id0	id1	id2	id3	id4
1	(0.0120, 0.0103) [7%]	(0.0007, 0.0007) [5%]	(0.0051, 0.0048) [3%]	(0.0007, 0.0000) [100%]	(0.0008, 0.0008) [1%]
2	(0.0086, 0.0082) [2%]	(0.0017, 0.0015) [7%]	(0.0072, 0.0063) [7%]	(0.0028, 0.0026) [4%]	(0.0200, 0.0147) [15%]
3	(0.0227, 0.0262) [7%]	(0.0032, 0.0051) [23%]	(0.0060, 0.0052) [7%]	(0.0062, 0.0096) [22%]	(0.0051, 0.0080) [22%]
4	(0.0250, 0.0200) [11%]	(0.0004, 0.0000) [100%]	(0.0081, 0.0090) [5%]	(0.0170, 0.0151) [6%]	(0.0090, 0.0062) [19%]
5	(0.0133, 0.0116) [7%]	(0.0006, 0.0009) [19%]	(0.0069, 0.0053) [13%]	(0.0026, 0.0027) [1%]	(0.0093, 0.0076) [10%]
6	(0.0100, 0.0091) [5%]	(0.0035, 0.0020) [28%]	(0.0086, 0.0084) [1%]	(0.0036, 0.0035) [2%]	(0.0062, 0.0058) [3%]
7	(0.0072, 0.0081) [6%]	(0.0006, 0.0007) [9%]	(0.0042, 0.0043) [1%]	(0.0021, 0.0022) [3%]	(0.0053, 0.0059) [6%]
8	(0.0031, 0.0030) [2%]	(0.0033, 0.0034) [2%]	(0.0029, 0.0027) [4%]	(0.0100, 0.0107) [4%]	(0.0069, 0.0074) [4%]

Equations

$$H$$

$$P_r=0,\frac{\partial H}{\partial r}=0.$$

$$\Omega =\left(\frac{\partial H}{\partial P_{\varphi }}\right),$$

$$P_t=\frac{P_{\varphi }}{r},$$

$$M_{\mathrm{ADM}}=M+H,$$

$$M\Omega$$

$$\frac{P_t}{M}$$

$$\frac{M_{\mathrm{ADM}}}{M}$$

$$\frac{x_1}{M}$$

$$\frac{y_1}{M}=\frac{y_2}{M}$$

$$\frac{z_1}{M}=\frac{z_2}{M}$$

$$M_{\mathrm{ADM}}=M+E_{\mathrm{Orb}},$$

$$\frac{d{M}_{\mathrm{ADM}}}{dt}=\frac{dM}{dt}+\frac{d{E}_{\mathrm{Orb}}}{dt}.$$

$$-\frac{d{E}_{\mathrm{GW}}}{dt}-\frac{dM}{dt}=\frac{d{E}_{\mathrm{Orb}}}{dt},$$

$$\frac{dr}{dt}=-\left(\frac{d{E}_{\mathrm{GW}}}{dt}+\frac{dM}{dt}\right){\left(\frac{d{E}_{\mathrm{Orb}}}{dr}\right)}^{-1}.$$

$$\frac{dr}{dt}=\frac{\partial H}{\partial P_r},$$

$$\frac{dr}{dt}$$

$$H$$

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

$$r(t)$$

$$\Phi (t)$$

$$e_{D1}=\mathrm{Amp}\left(\frac{D(t)-D_{\sec }(t)}{D(t)}\right),$$

$$e_{D2}=\mathrm{Amp}\left(\frac{D(t{)}^2\ddot{D}(t)}{M_{\mathrm{ADM}}}\right),$$

$$h^{\mathrm{TT}}=\frac{2}{R}{\frac{d^2{J}_{ij}^{\mathrm{TT}}}{d{t}^2}|}_{\mathrm{ret}}{\bar{m}}^i{\bar{m}}^j,$$

$$m^j=\frac{\partial x^j}{\partial \theta }+\frac{i}{\sin \theta }\frac{\partial x^j}{\partial \varphi },$$

$${\psi }_{4(22)}=A_{22}(t)\exp [i{\phi }_{22}(t)],$$

$$A_{22}(t)$$

$${\phi }_{22}(t)$$

$${\omega }_{22}$$

$$K=\frac{64\mu m^2}{R{A}^4}\sqrt{\frac{\pi }{5}}$$

$$e_A=\left(\frac{8}{39}\right)\mathrm{Amp}\left(\frac{A_{22}(t)-A_{22\sec }(t)}{A_{22}(t)}\right),$$

$$e_{\omega }=\left(\frac{8}{21}\right)\mathrm{Amp}\left(\frac{{\omega }_{22}(t)-{\omega }_{22\sec }(t)}{{\omega }_{22}(t)}\right),$$

$$e_{\varphi }=\left(\frac{4}{21}\right)\mathrm{Amp}\left[{\phi }_{22}(t)-{\phi }_{22\sec }(t)\right].$$

$$h_{(22)}=B_{22}(t)e^{i{\vartheta }_{22}(t)},$$

Full text

I Introduction

With the 2005 breakthroughs Pretorius (2005); Campanelli et al. (2006a); Baker et al. (2006), numerical relativity techniques are now routinely used to simulate the late inspiral and merger of black-hole binary systems (BHB), and its detailed predictions of the gravitational waves (GWs) produced by those systems have been recently observed by LIGO Abbott et al. (2016a, b, c). Direct numerical solutions of Einstein’s equations were used in the validation and interpretation of these discoveries Abbott et al. (2016a, d, e, f, c, b, g); Lovelace et al. (2016).

Because LIGO is only sensitive to the last few orbits of stellar-mass BHB mergers, and because of the large computational resources needed to produce these simulations, BHB simulations start at relatively small initial separations in astrophysical terms. On the other hand, by the time a BHB enters the LIGO sensitivity band we expect that the eccentricity will be very small because any eccentricity of astrophysical origin would have been most likely radiated away well before the late-inspiral at a rate proportional to $d^{19/12}$ Peters (1964), with $d$, the separation of the binary (see, for instance, Fig. 6 of Ref. Mroue et al. (2010) or Fig. 9 in Ref. Lousto et al. (2016)).

Given the current limitations of numerical relativity to carry out BHB simulations down to merger (see though Ref. Lousto and Zlochower (2013)) from initial separations much larger than $\sim 25M$ Lousto and Healy (2015); Szilagyi et al. (2015), it is crucial to give initial orbital parameters corresponding to an eccentricity sufficiently low that its effect on the waveforms is below the relevant accuracy requirements.

The problem of finding low-eccentricity parameters was recognized early on. One powerful technique to generate these parameters uses post-Newtonian (PN) quasicircular approximations to generate orbital parameters. The first systematic study (up to third order, i.e., 3PN) for subsequent use in numerical evolutions was performed in Ref. Baker et al. (2002). The technique that was used there to provide initial data is for the first time described in detail in this current paper, as well as its completion to 3PN and extension to incorporate 4PN corrections, in Sec. II.1.

Husa et al. Husa et al. (2008) pioneered an extension to this method. Rather than using quasicircular parameters, they used the PN equations of motion to evolve non-spinning binaries from large separations ($\sim 100M$) down to the separation of the start of the subsequent numerical evolution ($\sim 10M$). They then used the evolved PN orbital parameters to construct initial data. With this technique, they were able to generate reasonably low eccentricity values.

Our group used a similar technique to provide low eccentricity initial data for spinning and precessing binaries in Ref. Campanelli et al. (2009) and in subsequent papers.

The most accurate way to produce low-eccentricity data is to use an iterative procedure to correct the initial orbital parameters until sufficiently low eccentricities are obtained Pfeiffer et al. (2007); Buonanno et al. (2011); Purrer et al. (2012); Buchman et al. (2012). Each step in this iterative procedure requires the evolution of a binary for a few orbits. Thus this requires several full numerical evolutions before starting the definitive one.

In this paper, we revisit the scenario of quasicircular orbits as defined by the PN approximation. We update our previous approach from incomplete 3PN to complete 3PN, and include the possibility of an initial radial velocity based on instantaneous radiation reaction terms. Finally, we also consider the direct evolution of the spinning 3PN and 4PN equations of motion from large separations down to the desired starting separation to provide initial parameters for the numerical evolutions.

The paper is organized as follows. In Sec. II, we describe how to compute quasicircular orbits up to 3PN and how to calculate the radial inspiral velocity based in the radiation reaction terms. We also describe the 4PN equations of motion and how a PN inspiral simulation can be used to provide the alternative initial parameters. In Sec. III, we explore eight different spinning (but non-precessing) BHB configurations with initial orbital parameters given by incomplete 3PN and complete 3PN quasicircular parameters (with and without an initial radial velocity), and 3PN and 4PN inspiral parameters (a total of 40 individual simulations). We conclude with a discussion in Sec. IV of the benefits of each method to obtain a first and computationally simple set of initial configurations leading to small eccentricity for applications not requiring initial eccentricities below $\sim 10^{-3}$.

Throughout this paper we use geometric units where $G=c=1$. The vacuum general relativity field equations are scale invariant. In the case of a black hole binary, rescaling the total mass while keeping the mass ratio fixed, keeping the dimensionless spins fixed, and rescaling the momenta by the same factor as the masses leads to an equivalent solution. When reporting quantities with dimension, we rescale each by an appropriate power of an arbitrary positive constant $M$ (which has dimensions of mass).

II Post-Newtonian orbits

II.1 Quasicircular Initial Parameters

Our construction of quasicircular orbital parameters is based on PN dynamics in the Arnowitt-Deser-Misner transverse traceless (ADM-TT) gauge. This gauge is closely related (up to 2PN) to isotropic gauge used in puncture initial data Kelly et al. (2007); Tichy et al. (2003). From the Hamiltonian, i.e., the conservative part of the PN dynamics, we get the ADM mass $M_{\rm ADM}$, orbital separation $r$, and tangential linear momentum $P_{t}$ all in terms of the orbital frequency $\Omega$, or, alternatively, $M_{\rm ADM}$, $P_{t}$, and $\Omega$ can be obtained as functions of $r$. The mass, separation, and linear momenta are all calculated in the center of mass. To get the locations of the two BHs, we also need expressions for the position of the center of mass. Finally, the dissipative part, i.e., the radiation reaction gives the radial momentum $P_{r}$.

In an actual numerical simulation, we use the ADM-TT positions of the two BHs, their tangential and radial momenta, and their dimensionless spins to construct the initial data. That is, we take these ADM-TT position and momentum parameters, as well as the parameters $\vec{S}_{i}=m_{i}^{2}\vec{\chi}_{i}$, where $m_{i}$ is the PN mass parameter and $\chi_{i}$ is the PN dimensionless spin of particle $i$, and plug them directly into the Brandt-Brügmann initial puncture formalism Brandt and Brügmann (1997). The puncture masses of the two BHs is set by demanding that the numerical horizon (Christodoulou) masses match the ADM-TT mass parameters.

II.1.1 From the Hamiltonian

We start from the 3PN ADM-TT Hamiltonian,

[TABLE]

The terms in Eq. (3) were taken from Ref. Buonanno et al. (2006) augmented with the next-to-leading order (NLO) spin-orbit coupling, $H_{\rm SO,2.5PN}$, derived by Ref. Damour et al. (2008) and the NLO $S_{1}$-$S_{2}$ coupling, $H_{\rm S_{1}S_{2},3PN}$, derived by Refs. Steinhoff et al. (2008a, b) that was used in Ref. Campanelli et al. (2009). The above Hamiltonian also has the NLO spin-squared ($S_{1}^{2}$ and $S_{2}^{2}$) terms, $H_{\rm S^{2},3PN}$, presented in Ref. Steinhoff et al. (2008c).

In standard spherical coordinates, $\{r,\,\theta,\,\phi\}$, the quasicircular conditions in absence of radiation reaction are

[TABLE]

With this information, we can solve for $P_{\phi}$ as a function of $r$ and can then obtain the initial orbital frequency

[TABLE]

as a function of $r$, as well. Note that since it is not necessary to solve for $\dot{P}_{\phi}$, we do not need to evaluate $\partial H/\partial\phi$. The tangential linear momentum is given by

[TABLE]

where $P_{t}=P_{1y}=-P_{2y}$ for the individual BHs (the BHs are assumed to lie on the $x$-axis initially). The ADM mass is given by

[TABLE]

where $M=m_{1}+m_{2}$ is the total mass, and the Hamiltonian $H$ is calculated with $P_{r}=0$ and the solution of $P_{\phi}$. Finally, the total ADM angular momentum of this system is $\vec{J}=\{S_{1x}+S_{2x},\,S_{1y}+S_{2y},\,L_{z}+S_{1z}+S_{2z}\}$ with $L_{z}=P_{\phi}$, and the BH’s spins, $\vec{S}_{1}$ and $\vec{S}_{2}$.

In the following, we summarize some quantities derived from the above analysis. Here, we use the mass ratio, $q=m_{1}/m_{2}$, the nondimensional spin, $\chi_{1x}=S_{1x}/m_{1}^{2}$, $\chi_{1y}=S_{1y}/m_{1}^{2}$, $\chi_{1z}=S_{1z}/m_{1}^{2}$, $\chi_{2x}=S_{2x}/m_{2}^{2}$, $\chi_{2y}=S_{2y}/m_{2}^{2}$, and $\chi_{2z}=S_{2z}/m_{2}^{2}$. The symmetric mass ratio is given by $\eta=q/(1+q)^{2}$. First, all quantities are written in terms of the orbital separation $r$. The orbital frequency, tangential linear momentum, and the ADM mass are given by

[TABLE]

where $P_{1y}=-P_{2y}=P_{t}$.

In Appendix A, we present the ADM mass, orbital separation, and the tangential linear momentum in terms of the orbital frequency $\Omega$.

II.1.2 From the center-of-mass vector

In the above analysis, we used the center-of-mass Hamiltonian, i.e., the center of mass was located at the origin and the vector $\vec{r}$ is the displacement $\vec{r}=\vec{x}_{1}-\vec{x}_{2}$. To obtain the positions of each BH, we need to know the position of the BHs relative to the center of mass. To do this, we use expressions for the center-of-mass (given BH positions $\vec{x}_{1}$ and $\vec{x}_{2}$) based on the nonspinning terms in Ref. Damour et al. (2000), the spin-orbit terms in Ref. Damour et al. (2008), the $S_{1}$-$S_{2}$ terms in Ref. Steinhoff et al. (2008b), and the spin-squared ($S^{2}$) terms in Ref. Hergt and Schaefer (2008).

Then, given that we want $\vec{x}_{1}-\vec{x}_{2}$ to lie along the $x$-axis, and the center of mass to lie on the origin, we solve for $\vec{x}_{1}$ and $\vec{x}_{2}$. The resulting expressions for the BH positions are given by

[TABLE]

where $x_{2}=r-x_{1}$. Note that for the precessing case, the two BHs are displaced from the $x$ axis. However, as required $y_{1}-y_{2}=z_{1}-z_{2}=0$.

II.1.3 From radiation reaction

The radial momentum (velocity), which is driven by the emission of gravitational waves and tidal heating drives the inspiral. To get low-eccentricity orbital parameters, the corresponding radial momentum generally needs to be included.

To obtain the radial momentum, we begin with the time derivative of the ADM mass

[TABLE]

where $E_{\rm Orb}$ is the orbital energy for quasicircular orbits. Both the orbital energy and the BH masses change during the inspiral (the change in the latter is due to tidal heating effects). Thus, we have

[TABLE]

The ADM mass loss produces the flux of gravitational wave energy leaving the binary. Hence

[TABLE]

and since $dE_{\rm Orb}/dt=(dr/dt)(dE_{\rm Orb}/dr)$, we have

[TABLE]

For nonspinning and nonprecessing cases, we use the orbital energy and GW energy flux summarized in Ref. Ajith et al. (2007), and for the precessing case, we can use the formulas summarized in Appendix A of Ref. Ossokine et al. (2015). Ref. Ajith et al. (2007) uses the work of Alvi Alvi (2001) to calculate $dM/dt$. While higher-order correction to $dM/dt$ in the nonprecessing case are known (see, e.g., Ref. Chatziioannou et al. (2013)), additional correction for the precessing case are not known. Thus we use the formula given in Ref. Ajith et al. (2007) for both the non-precessing and precessing cases.

Given $E_{\rm Orb}$, $dE_{\rm GW}/dt$, and $dM/dt$, it is straightforward to calculate $dr/dt$. However, what we need is a formula for $P_{r}$, rather than $dr/dt$ itself. To obtain this, we note that

[TABLE]

In the right hand side of the above equation, we pick up only $O(P_{r}^{0})$ and $O(P_{r}^{1})$ terms because $O(P_{r}^{2})$ and higher order terms can be ignored in the PN approximation. In practice, we have

[TABLE]

and we may solve this equation with respect to $P_{r}$. This treatment is basically the same as the one presented in Ref. Buonanno et al. (2006) if we ignore the $O(P_{r}^{0})$ term. Note that we have ignored effects due to the change in mass of the two BHs, i.e., $dM/dt$. This is because the effects of $dM/dt$ are comparable in magnitude to the 4PN corrections to $dE_{\rm GW}/dt$. A complete self-consistent formulation of these 4PN terms has not not been completely worked out yet, and is thus also ignored here.

II.2 Orbital Evolution

Given the extent of some of the PN expressions used in the integration of the equations of motion for the orbital evolution, we do not provide them explicitly here, but give a detailed set of references where these expressions can be found.

Our original implementation of the 3PN orbital equations of motion in the ADM-TT gauge (that for the sake of simplicity we will label this as “3PNevolution”) was first described in Ref. Campanelli et al. (2009). That code was based on the formulation developed in Ref. Buonanno et al. (2006) in which the Hamiltonian has the form

[TABLE]

The equations of motion were obtained from Eqs. (2.23)-(2.25) and (3.1) in Ref. Buonanno et al. (2006), and the radiation reaction force from Eq. (3.27) there, as well. In addition, we added higher-order PN terms derived in Refs. Damour et al. (2008); Steinhoff et al. (2008a, c) to the above Hamiltonian, and higher-order corrections to the radiation reaction derived in Refs. Arun et al. (2009); Blanchet et al. (2011). The Hamiltonian is then given by Eq. (3).

The new code that, for the sake of simplicity, we will call the “4PNevolution” code includes, in addition to the above expressions, higher-order terms we list below.

1. (i)

New higher-order radiation reaction force terms from Ref. Ossokine et al. (2015). We removed spin terms in Eq. (3.27) of Ref. Buonanno et al. (2006) because the orbital averaging of this force becomes zero. In the EOB approach, the same treatment has been done (see, e.g., Ref. Pan et al. (2010)). 2. (ii)

Spin $S^{3}$ and $S^{4}$ terms in the Hamiltonian Levi and Steinhoff (2015). 3. (iii)

4PN nonspinning local term from Eq. (5.13) in Ref. Damour et al. (2014), as well as nonlocal terms from Eq. (7.9) (or (7.12a)) of Ref. Damour et al. (2015). Note that we assume “quasicircular” orbits for this nonlocal term. The consistent result for the 4PN nonspinning term was recently derived in Ref. Bernard et al. (2016). 4. (iv)

3.5PN spin-orbit coupling (NNLO SO) from Eq. (140a) of Ref. Hartung et al. (2013) which is in in the center-of-mass frame. This has been confirmed in Ref. Levi and Steinhoff (2016a). 5. (v)

4PN S1S2 coupling terms (NNLO S1S2) from Eq. (140b) of Ref. Hartung et al. (2013). This has been confirmed by Ref. Levi and Steinhoff (2014) (note that there was a typo in Eq. (140b)).

We note that we do not include any corrections for:

1. (i)

The change in mass of the BHs during the evolution. 2. (ii)

4PN spin-squared (NNLO $S^{2}$) in the EFT gauge. It was derived in Ref. Levi and Steinhoff (2016b), but it has not yet been confirmed in the ADM-TT gauge.

III Full Numerical Evolutions

III.1 Methods

Due to large computational expense of numerical relativity simulations of merging BHBs, in order to make systematic studies and build a data bank of full numerical simulations, it is crucial to develop efficient numerical algorithms. To this end, we evolve the following BHB data sets using the LazEv Zlochower et al. (2005) implementation of the moving puncture approach Campanelli et al. (2006a); Baker et al. (2006) with the conformal function $W=\sqrt{\chi}=\exp(-2\phi)$ suggested by Ref. Marronetti et al. (2008). For the runs presented here, we use centered, sixth-order finite differencing in space Lousto and Zlochower (2008) and a fourth-order Runge Kutta time integrator. This sixth-order spatial finite differencing speeds up the code by a factor of $4/3$ compared to an eighth-order implementation (mostly, this is due to the reduction in the number of ghostzones). We also used a Courant factor (CFL) of $1/3$ instead of the previous CFL of $1/4$ Zlochower et al. (2012), gaining another speedup factor of 4/3. We verified that, for the runs presented here, increasing the CFL and reducing the finite difference order still lead to acceptable conservation of the BH masses and spins during the evolution, as well as an acceptable gravitational wave phase error during the entire simulation (below $10^{-5}$).

This plus the use of the new Xsede supercomputer Comet at SDSC 111https://portal.xsede.org/sdsc-comet led to typical evolution speeds of $250M/{\rm day}$ on 16 nodes. Note that our previous Lousto and Zlochower (2013); Lousto et al. (2016) comparable simulations averaged $\sim 100M/{\rm day}$.

Our code uses the EinsteinToolkit Löffler et al. (2012); ein / Cactus cac / Carpet Schnetter et al. (2004) infrastructure. The Carpet mesh refinement driver provides a “moving boxes” style of mesh refinement. In this approach, refined grids of fixed size are arranged about the coordinate centers of both holes. The Carpet code then moves these fine grids about the computational domain by following the trajectories of the two BHs.

We use AHFinderDirect Thornburg (2004) to locate apparent horizons. We measure the magnitude of the horizon spin using the isolated horizon (IH) algorithm detailed in Ref. Dreyer et al. (2003) and as implemented in Ref. Campanelli et al. (2007). Note that once we have the horizon spin, we can calculate the horizon mass via the Christodoulou formula ${m_{H}}=\sqrt{m_{\rm irr}^{2}+S_{H}^{2}/(4m_{\rm irr}^{2})}\,,$ where $m_{\rm irr}=\sqrt{A/(16\pi)}$, $A$ is the surface area of the horizon, and $S_{H}$ is the spin angular momentum of the BH (in units of $M^{2}$).

We use the Antenna code Campanelli et al. (2006b) to calculate the gravitational waveform via the Weyl scalar $\psi_{4}$. Here, we decompose $\psi_{4}$ into $(\ell,m)$ modes. For the present work, we analyze the $(\ell=2,m=\pm 2)$ modes, in particular, because of the relatively simple way the eccentricity of (non-precessing) binaries can be extracted from them.

III.2 Initial Orbital Parameters

To compute the numerical initial data, we use the puncture approach Brandt and Brügmann (1997) along with the TwoPunctures Ansorg et al. (2004) thorn. In this approach the 3-metric on the initial slice has the form $\gamma_{ab}=(\psi_{BL}+u)^{4}\delta_{ab}$, where $\psi_{BL}$ is the Brill-Lindquist conformal factor, $\delta_{ab}$ is the Euclidean metric, and $u$ is (at least) $C^{2}$ on the punctures. The Brill-Lindquist conformal factor is given by $\psi_{BL}=1+\sum_{i=1}^{n}m_{i}^{p}/(2|\vec{r}-\vec{r}_{i}|),$ where $n$ is the total number of ‘punctures’, $m_{i}^{p}$ is the mass parameter of puncture $i$ ($m_{i}^{p}$ is not the horizon mass associated with puncture $i$), and $\vec{r}_{i}$ is the coordinate location of puncture $i$. For the initial (conformal) extrinsic curvature we take the analytic form $\hat{K}_{ij}^{BY}$ given by Bowen and York Bowen and York (1980). In the puncture formalism, there are 15 non-trivial free parameters. These are the initial coordinate separation of the two BHs, the three components of the linear momentum and spin of each BH, and finally, the mass parameter of each BH. The momentum, spin and separation parameters are obtained directly from the various PN approximations described above. The mass parameters, however, have to be set by demanding that the total ADM mass matches the PN prediction and the mass ratio matches the desired value.

In this work, we evolve eight sets of simulations of spinning, nonprecessing BHBs spanning a range of mass ratios $1/3\leq q\leq 1$, including nonspinning cases. For each configuration, we consider five different approximations to generate low-eccentricity data. These are:

• id0

Quasicircular data using an incomplete 3PN Hamiltonian.

• id1

Quasicircular data with radiation reaction-driven radial momentum using all 3PN terms.

• id2

Quasicircular data with no radial momentum using the full 3PN Hamiltonian.

• id3

Inspiral parameters from a 3PN evolution from large separations.

• id4

Inspiral parameters from a 4PN evolution from large separations.

The initial data parameters for the runs presented here are given in Table 1 (note that we consider 5 momentum variations of the same 8 basic configurations). In this paper, we chose to provide initial data by specifying an initial orbital separation rather than an initial orbital frequency (this choice is convenient when obtaining the parameters based from PN inspiral evolutions), although the formalism allows for specifying the initial orbital frequency instead.

It is interesting to note how the various approximations change the momentum parameters. In Fig. 1, we show the tangential momentum versus radial momentum for all 40 configurations.

Finally, it is important to note that we are testing these variations of choices of the initial data parameters in the context of Brandt-Brügmann puncture data. While Brandt-Brügmann data are quite popular, the effects of the different parameters choices on the eccentricities produced by other methods for solving binary data, such as the extended conformal thin sandwich methods used by the Spectral Einstein Code (SpEC) spec_web for both conformally flat Cook and Pfeiffer (2004); Caudill et al. (2006) and conformally Kerr Lovelace et al. (2008) backgrounds, as well as on the newly-developed puncture-based conformally Kerr data Ruchlin et al. (2016), are not determined here.

III.3 Results

We measure eccentricity, $e$, using several orbital and waveform-based methods. For our orbital based methods, we use the time variation of the proper distance of the coordinate line joining the two BHs, which we refer to as the simple proper distance, or SPD. More specifically, the SPD is the proper length of the part of the coordinate ray joining the two centers that is outside both apparent horizons. 222 We found in Ref. Lousto and Zlochower (2013), Figure 9, that the coordinate separation shows substantial gauge effects during the first few orbits. While the SPD varies consistently on the orbital frequency, the coordinate separation can vary at twice this frequency. We observe the same effect in the runs in this paper, hence we choose to work with the SPD for eccentricity computations. In the Newtonian limit the separation vector ($\vec{r}$) between the two BHs (i.e., Newtonian point particles) is given by

[TABLE]

We therefore define two distance eccentricities as

[TABLE]

and

[TABLE]

where $D(t)$ is the numerical SPD, a dot indicates a coordinate time derivative, $D_{\rm sec}(t)$ is the non-oscillatory part of $D(t)$, and ${\rm Amp}$ indicates the amplitude of a sinusoidal function. In practice, $D_{\rm sec}$ is constructed by fitting $D(t)$ to a low-order polynomial in $\sqrt{t}$ plus a sinusoidal function. The polynomial part of the fit is $D_{\rm sec}$. The difference $D(t)-D_{\rm sec}(t)$ is then dominated by the sinusoidal part of $D(t)$. Both $e_{D1}$ and $e_{D2}$ are only accurate to linear order in the eccentricity for Newtonian orbits. Note that $e_{D1}$ was introduced in Ref. Pfeiffer et al. (2007) and $e_{D2}$ was introduced in Ref. Campanelli et al. (2009).

Following Ref. Mroue et al. (2010), we similarly define three eccentricities based on the amplitude, frequency, and phase of the $(\ell=2,\,m=2)$ mode of the Weyl scalar, $\psi_{4(22)}$. To find how the eccentricity affects the waveform, we start with the quadrapole formula (see, e.g., Ref. Carroll (2004))

[TABLE]

where $J_{ij}^{\rm TT}$ is transverse and traceless part of the quadrapole moment (evaluated at retarded time), the complex null-vector $m^{j}$ is given by

[TABLE]

$\theta$ and $\phi$ are the standard spherical-polar coordinates, and $R$ is the distance from the binary to the observer. To obtain $\psi_{4}$, we take the second time derivative of the $h^{\rm TT}$. For two point particles with separation vector $\vec{D}(t)$ given by the low-eccentricity Newtonian trajectory for two point particles, i.e., Eq. (58), $\psi_{4(22)}$ is given by

[TABLE]

where

[TABLE]

$K$ is an overall normalization factor (given by

[TABLE]

in the Newtonian limit, where $m$ is the total mass of the binary, $\mu$ is its reduced mass, $A$ is the average separation of the binary, and $R$ is the distance from the binary to the observer) and $\omega_{22}=d\varphi_{22}/dt$. Hence we define three waveform-based eccentricity measures

[TABLE]

[TABLE]

and

[TABLE]

We choose to define eccentricities using $\psi_{4}$ rather than $h^{\rm TT}$ because $\psi_{4}$ is directly calculated from the simulation results whereas $h^{\rm TT}$ requires a double time integral. We compare the performance of each of these measures of the eccentricity in Fig. 2. Note how all are roughly equivalent for run #1. Note that $e_{\phi}$ and $e_{\omega}$ were introduced in Ref. Mroue et al. (2010), but with slightly different coefficients than presented here. Note that due to the need for fitting the secular part, all eccentricity measures, except $e_{D2}$, can suffer from a significant biases due to the form of the fitting function and chosen interval for the fit. However, $e_{D2}$ can be calculated directly from simulation data and, while with some more residual gauge dependence than $e_{A}$, $e_{\phi}$, and $e_{\omega}$, $e_{D2}$ provides a critical sanity check for those other measures.

For reference, we also include eccentricity measures based on $h^{\rm TT}$. The $(\ell=2,m=2)$ mode of $h^{TT}$ [see Eq. (61)] is given by

[TABLE]

where

[TABLE]

and $\varpi_{22}=d\vartheta_{22}/dt$. The corresponding eccentricity measures are

[TABLE]

[TABLE]

and

[TABLE]

The main result of this work is the comparison of the performance of the 5 choices for low-eccentricity data. In most cases, id0 produced the most eccentric data, and id1 produced the best. Similarly, in most cases, id4 (4PN inspiral parameters) were higher eccentricity than id3 (3PN inspiral parameters). We show the sinusoidal dependence of $e_{A}$ for all five initial data variations for run #7 in Fig. 3. The eccentricities are given by the amplitudes of these oscillations.

In Table 2, we provide the eccentricities $e_{D2}$ and $e_{A}$ (representative of using trajectories and waveforms respectively) for all configurations. As both of the these are approximations, the two are not supposed to agree exactly. We see about a $\lesssim 20\%$ difference between the two measures for most configurations. The largest relative differences occur at low-eccentricities. This is due to waveform and trajectory noise interfering with our ability to determine the secular and oscillating parts of the amplitude of $\psi_{4}$ and the second-time-derivative of $D$.

In Fig. 4, we show the average of $e_{D2}$ and $e_{A}$ for all configurations. From the plot, it is clear that id1 is the best overall choice for initial parameters. On average, the next best choice is id3, which was proposed in Ref. Husa et al. (2008). Perhaps surprisingly, id4 is, on average, worse than id3. This is remarkable because id4 is generated with higher-order PN terms than id3 (both are based on PN inspirals). This suggest that completing all 3PN terms (id2) and including the radial momentum (id1) consistently leads to the best results regarding reduction of eccentricities. Notably, the evolution from large separations may lead to some initial data with higher eccentricities (run #4 for 3PN and run #2 for 4PN).

The observation that our 4PN evolution does not give the lowest initial eccentricity may be a consequence of the poor convergence of the PN series. Although one cannot make a strong statement concerning this for the comparable mass case, we note that the 3PN calculation has a larger region of validity than 4PN for extreme mass ratio inspirals (see, e.g., Refs. Yunes and Berti (2008); Zhang et al. (2011); Sago et al. (2016)).

IV Conclusions and Discussion

In this paper, we measured the eccentricity resulting from fully non-linear numerical evolutions of BHBs using orbital parameters obtained from various post-Newtonian approximations. We find that using the full 3PN Hamiltonian to generate quasicircular orbital parameters and the addition of a radiation reaction-driven inspiral momentum leads to significantly lower eccentricities than the other methods we tried. This includes inspiral parameters based on 3PN and 4PN evolutions from large separations. Interestingly, we find that the (incomplete) 4PN inspiral parameters are sometimes more eccentric than the 3PN parameters.

These new quasicircular (plus radial) parameters can also serve as the initial seed for iterative methods for further reducing the eccentricity (such as in Refs. Pfeiffer et al. (2007); Buonanno et al. (2011); Purrer et al. (2012)). Choosing good starting parameters will reduce the number of iterations required and thus reduce the overall computational cost considerably.

An interesting application of these measurements of eccentricity from the waveforms, or more precisely, the analogs which directly use the strain, i.e., $e_{\varpi}$, $e_{\vartheta}$, and $e_{B}$, is that they can be used to give upper bounds to the binary’s eccentricity based solely on the observed waveform. Our measures here do not include precession effects, which can mimic true eccentricity, so these measurements would provide an upper bound.

An open question for follow-up work will be to find similarly performing initial parameters for highly-precessing binaries, as well Buonanno et al. (2011). The quasicircular data given in Sec. II.1, and 3PN and 4PN evolutions described in Sec. II.2 are applicable to precessing cases.

Acknowledgements.

The authors gratefully acknowledge the NSF for financial support from Grants No. PHY-1607520, No. PHY-1707946 , No. ACI-1550436 , No. AST-1516150, No. ACI-1516125. Computational resources were provided by XSEDE allocation TG-PHY060027N, and by NewHorizons and BlueSky Clusters at Rochester Institute of Technology, which were supported by NSF grant No. PHY-0722703, No. DMS-0820923, No. AST-1028087, and No. PHY-1229173. HN is also supported by MEXT Grant-in-Aid for Scientific Research on Innovative Areas, “New Developments in Astrophysics Through Multi-Messenger Observations of Gravitational Wave Sources”, No. 24103006, and JSPS Grant-in-Aid for Scientific Research (C), No. 16K05347.

Appendix A Initial parameters in terms of the orbital frequency

The relationship $\Omega(r)$ in Eq. (15) can be inverted to solve for $r(\Omega)$, and hence all quasicircular orbital parameters can be expressed as a function of the orbital frequency. In terms of the orbital frequency $\Omega$, the ADM mass, orbital separation, and the tangential linear momentum are given by

[TABLE]