{'fraction_non_alphanumeric': 0.0324455078484474, 'fraction_numerical': 0.012857255541240502, 'mean_word_length': 4.770596267637687, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Measuring orbital eccentricity and periastron advance in quasi-circular black hole simulations

We compare different methods of computing the orbital eccentricity of quasi-circular binary black hole systems using the orbital variables and gravitational wave phase and frequency. For eccentricities of about a per cent, most methods work satisfactorily. For small eccentricity, however, the gravitational wave phase allows a particularly clean and reliable measurement of the eccentricity. Furthermore, we measure the decay of the orbital eccentricity during the inspiral and find reasonable agreement with post-Newtonian results. Finally, we measure the periastron advance of non-spinning binary black holes, and we compare them to post-Newtonian approximations. With the low uncertainty in the measurement of the periastron advance, we positively detect deviations between fully numerical simulations and post-Newtonian calculations.


I. INTRODUCTION
The inspiral and merger of binary black holes or neutron stars is one of the most promising sources for current and future generations of gravitational wave detectors such as LIGO and VIRGO. The late stage of the inspiral, corresponding to the final few orbits and merger of the binary, is highly dynamical and involves strong gravitational fields, and it must be handled by numerical relativity. Breakthroughs in numerical relativity have allowed a system of two inspiraling black holes to be evolved through merger and the ringdown of the remnant black hole [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15].
During the inspiral of an isolated binary, the orbit circularizes via the emission of gravitational waves [16,17]. As a result, even binaries starting with some eccentricity at the beginning of their stellar evolution are expected to have negligible eccentricity by the time the frequency of the emitted gravitational radiation enters the frequency band of ground based detectors.
However, different physical scenarios [18][19][20][21][22][23][24][25][26] suggest that binaries could approach merger with a significant eccentricity without being circularized by radiation reaction. This implies that eccentric binaries are a potential gravitational wave source for ground based interferometers. For example, in globular clusters, the Kozai mechanism [18] could increase the eccentricity of an inner binary's orbit through a secular resonance caused by a third perturbing black hole on an outer orbit [20]. Many-body encounters of black holes in globular clusters could also result in the merger of highly eccentric binaries [19]. Ref. [24] predicted that 30% of the hierarchical triple black hole systems formed in a globular cluster will possess eccentricities greater than 0.1 when their emitted gravitational waves pass through a frequency of 10Hz.
For these reasons considerable attention has been paid to eccentric binaries. Analytical waveform templates have been constructed for the gravitational wave signal emitted by compact binaries moving in inspiraling eccentric orbits [27][28][29]. In this case, orbits involve three different time scales: orbital period, periastron advance and radiation reaction time scales. By combining these three time scales, one computes "postadiabatic" short-period contributions to the orbital phasing and gravitational wave polarizations. These gravitational wave polarizations are needed for astrophysical measurements with gravitational wave interferometers. Refs. [30][31][32] investigated the impact of eccentricity on gravitational wave detection, specifically the potential loss in the signal-tonoise ratio when "circular" waveform templates are applied to search for eccentric binaries.
Eccentric black hole binaries have also been studied with direct numerical simulations. Ref. [33] studied the variation of the signal to noise of the eccentric evolutions of intermediate mass binary black hole mergers as a function of mass and eccentricity. Ref. [34] presented binary black holes in zoom-whirl orbits where the waveforms are modulated by the harmonics of these zoom-whirls. In Ref. [35], the authors studied the transition from inspiral to plunge in general relativity by computing gravitational waveforms of eccentric nonspinning, equal mass black-hole binaries. They analyzed the radiation of energy and angular momentum in gravitational waves, the contribution of different multipolar components and the final spin of the remnant black hole. Ref. [36] presented results from numerical simulations of equal-mass, nonspinning binary black hole inspiral and merger for various eccentricities, and they measured the final mass and spin of the remnant black hole. Ref. [37] compared a numerical relativity simulation of an eccentric binary system with eccentricity 0.1 with corresponding post-Newtonian (PN) results. They found better agreement when the eccentric PN expressions are expanded in terms of the frequency-related parameter x ≡ (ΩM ) 2/3 , where Ω is orbital frequency and M is total mass of the binary, rather than the mean motion n = 2π/P , where P is the orbital period.
Beyond the Newtonian limit, the orbital eccentricity is not uniquely defined and a variety of definitions have appeared in the literature. Ref. [38] used a definition of the eccentricity for which a Newtonian orbit is momentarily tangent to the true orbit (the "osculating" eccentricity), while other authors [27][28][29]39] defined multiple "eccentricities" to encapsulate different aspects of noncircular orbits at PN order. Another useful definition for large eccentricity in numerical simulations is given in Refs. [40,41].
Similarly, numerical relativists [42][43][44][45][46] introduced several methods for defining and measuring the eccentricity using the residual oscillations in the orbital frequency, proper horizon separation and coordinate separation. These eccentricity definitions are necessary to compare the numerical waveforms with the waveforms produced by analytic techniques (i.e., PN methods). They behave differently depending on the magnitude of the eccentricity and details of the numerical simulation, like employed gauge conditions, or presence of numerical noise. This makes it important to specify the validity regimes of these definitions.
This paper deals with two related topics: First, we revisit many of the eccentricity definitions used so far in numerical work and compare them systematically. We find that for eccentricities of a few percent, most definitions work satisfactorily. However, for very small eccentricities, e ∼ 10 −4 , computation of the eccentricity based on the extracted gravitational waves is superior. In the second part of the paper, we measure decay of orbital eccentricity and periastron advance for inspiraling black hole binaries, and compare these measurements to post-Newtonian calculations.
Section II summarizes eccentricity definitions that are useful for measuring eccentricity in quasi-circular runs. In Section II, we compare these approaches, as well as some new ones, on the 15-orbit inspiral presented by Boyle et al. [14] and on the data of a new simulation of an eccentric (e = 0.05) nonspinning equal mass binary black hole. Next, by measuring the extrema in the eccentricity estimator, we estimate in Sec. III the decay of the eccentricity of these runs as well as the radial frequency. This allows us in Section IV to estimate the periastron advance for these runs from the ratio of the orbital frequency to the radial frequency as well as the periastron advance of a set of quasi-circular nonspinning binaries of mass ratios 2, 3, 4 and 6. The numerically estimated periastron advance is then compared to the 3PN formula of the periastron advance [29,39,47].

A. Definitions
For a non-precessing binary in an orbit with zero eccentricity, orbital variables and their time derivatives change monotonically as the holes inspiral to merger. In numerical simulations, however, a small eccentricity is introduced by imperfections of the initial data. As a result, small residual oscillations with amplitude proportional to the eccentricity are added to the monotonically changing orbital variables and their derivatives. To estimate the eccentricity, one needs to determine these residual oscillations.
Different methods to estimate the eccentricity [42,44,45] used the orbital frequency, separation between the holes (coordinate or proper separation), or some Newtonian formula containing both of these variables. Similarly, time derivatives of these variables could be used in these definitions of the eccentricity. Basically all approaches construct an eccentricity estimator e X (t) such that for Newtonian orbits e X (t) = e cos(Ω r t + φ), where e is the eccentricity 1 and Ω r is the frequency of radial oscillations in the quasi-circular orbit. The key property of e X (t) is that it is an oscillating function with amplitude equal to e. In order to define eccentricity for general relativistic inspirals, one computes a tentative eccentricity estimator e X (t), and checks its behavior. If it behaves as Eq. (1), one reads off the eccentricity e as the amplitude of the oscillations. The resulting eccentricity estimates are not local in time nor continuous functions of time but rather orbit-averaged quantities. Deviation from sinusoidal behavior indicates that particular eccentricity estimator is not reliable, and one must verify to what extent the eccentricity estimators behave as expected and to what extent they agree.
The estimated value of the eccentricity will differ slightly depending on the method used and the noise in the numerical data. In this paper, we compare typical eccentricity estimates using a Newtonian formula as in Ref. [42] or the orbital frequency and separation as in Ref. [45]. These eccentricities are also compared to new ones computed from the wave phase and frequency extracted at a given radius. Other definitions of the eccentricity could be used, but we restrict the study to these typical definitions.
To make this rather abstract discussion more concrete, consider the Newtonian formula for the radial distance d between the two objects with eccentricity e Newt Based on this formula, one can define the eccentricity estimator e d (t) where the average distanced equals d 0 in Newtonian gravity. For a general relativistic system, one obtains d(t) by a fit over several radial oscillation periods. If the residual d(t)−d(t) oscillates sinusoidally-which it indeed does for sufficiently large eccentricity-the amplitude of these oscillations defines an associated eccentricity e d . From the trajectory of the two objects, one can also use the orbital phase and frequency to define the corresponding eccentricity estimators using the following Newtonian relation [48]: where M is the mean anomaly and Φ is the orbital phase. Equivalent to Eq. (4) for numerical simulations is the relationship where Ω 0 is the average fitted orbital frequency and Φ 0 is some phase offset. Then the eccentricity estimator e Φ (t) is written as From the time derivative of Eq. (5) and the replacement Ω r → Ω 0 , we obtain an eccentricity estimator in terms of the orbital frequency (as in Ref. [45]) Notice that since the radial oscillation results from eccentricity, Ω r is different from Ω 0 , the average of the orbital frequency. The eccentricities of Eqs. (6) and (7) will differ by a factor Ω r /Ω 0 . For Newtonian orbits, Ω r /Ω 0 = 1, and this factor drops out. But for the binary black hole case, the factor is about 1.4, causing the difference between Figs. 4 and 5 below. This is easily seen by writing the eccentricity estimator from Eqs. (5) and (7) asΦ

B. Numerical data
Before introducing several further eccentricity estimators, let us briefly describe the numerical binary black hole simulations that we will analyze. All runs have been performed with the Spectral Einstein Code (SpEC) [49]. We will primarily analyze the 16 orbit long inspiral simulation of an equal mass, non-spinning black hole binary presented in Ref. [14] (specifically, the run labeled 30c-1). This run with eccentricity of about 6 × 10 −5 is used to compute the eccentricity data in Figs. 2, 3, 4 and 5. To compute eccentricity estimators, we use the orbital frequency Ω, the coordinate separation between the holes D, the proper horizon separation s (defined as the integrated distance between the holes along the coordinate axis, cf. Ref. [14]) as well as the gravitational wave phase φ and the gravitational wave frequency ω.
Furthermore, we utilize recent runs of quasi-circular nonspinning binaries [50] with mass ratios 2, 3, 4 (lasting 15 orbits) and mass ratio 6 (lasting 8 orbits). The eccentricity of these runs is also of the order of magnitude 10 −5 . The periastron advance and the resulting frequency modulation are estimated in Fig. 8.
As a separate check, another equal mass nonspinning binary with moderate eccentricity (e ∼ 0.05) is evolved to compare various eccentricity estimators and measure the periastron advance for a case that is not quasi-circular. Figure 1 shows the proper separation as well as the orbital frequency as a function of time for this eccentric binary.

C. A Newtonian definition
The first use of eccentricity estimators was by Buonnano, Cook & Pretorius [42], who consider the following relationship that holds for Newtonian orbits with eccentricity e Newt : Here Ω φ (t) and φ(t) denote orbital frequency and phase, respectively, and r is the separation of the masses. Motivated by Eq. (9), Buonnano, Cook & Pretorius define an eccentricity estimator where now Ω φ (t) and r(t) are extracted from the numerical simulation. To compute this eccentricity estimator e BCP , we fit the function Ω φ (t) 2 r(t) 3 /M to a polynomial in time, We found that a fifth order polynomial ensures a good fit. The polynomial order needs to be high enough to reliably capture the smooth inspiral trend in Ω(t) 2 r(t) 3 /M , but it should not capture the higher frequency oscillations due to eccentricity. When applying this procedure to a binary black hole inspiral, one has to decide how to generalize the Newtonian separation r(t) to curved space. We use two choices, the coordinate distance D(t) between the centers of the apparent horizons, and the proper separation s(t) between the apparent horizons, computed along a straight coordinate line connecting the centers of the apparent horizons. In Fig. 2, we plot the eccentricity estimator e BCP computed using the coordinate separation and proper horizon separation as described above. In the top panel, we plot e BCP (t) using the binary run with eccentricity e ∼ 0.05. Using the proper horizon separation s, the estimated initial eccentricity, 0.07, is larger by nearly a factor of 2 than in the case where the coordinate separation D is used (0.03). This is due to different numerical values for the distances, (s/D) 3 ∼ 1.8. Both eccentricity estimators are in phase during the whole time interval as expected. In both cases, the eccentricity magnitude decreases between t = 0 and t = 2500M . In this case, a clear decaying sinusoidal signal is obtained without any higher harmonics showing up at later times.
In the bottom panel, we examine the equal mass binary with eccentricity e ∼ 6 × 10 −5 . For this case, no clean sinusoidal signal is apparent. While e BCP computed from s(t) shows oscillations, they are faster than the orbital period, and can therefore not be attributed to orbital eccentricity. Because e BCP does not show the expected behavior, it is not meaningful to attribute a value of eccentricity to this analysis. For these small eccentricities, e BCP is dominated by other effects, possibly the coordinate dependence of the separation measurements.

D.
Eccentricity from orbital variables Husa et al [45] fitted directly the orbital frequency Ω(t) or the coordinate separation D(t) to a function of the form with fitting parameters t m , the coalescence time, and the coefficients a i . The eccentricity estimator is then defined as where X NR (t) is the numerical orbital variable and X fit (t) is the polynomial fit of X NR (t). We shall compute three eccentricity estimators using Eq. (12), which differ in the quantity being fitted: e s (t) and e D (t) are based on proper separation and coordinate separation between the black holes, with the value k = 1; e Ω (t) uses the orbital frequency, where k = 2. In the Newtonian limit, these estimators are identical to first order in eccentricity. Figure 3 shows these eccentricity estimators for a run with fairly large eccentricity and for a run with very small eccentricity. For large eccentricity e = 0.05, the various eccentricity estimators have a smooth decaying sinusoidal signal. This allows measuring a nearly identical value of the eccentricity for the three orbital variables from the amplitude of the residual oscillations. The phasing is also consistent between the different eccentricity estimates: The orbital frequency is a maximum when the separation is a minimum and vice-versa.
In the bottom panel of Fig. 3, we plot the eccentricity estimators applied to a simulation with much smaller ec- centricity e ∼ 6 × 10 −5 . The behavior of e D and e Ω is erratic. Higher-order harmonics are clearly visible, and the extrema are not monotonically decreasing, as one would expect from the circularizing effect of gravitational radiation. However, e s shows no increase in the eccentricity during the late stages of the inspiral, and no additional significant harmonics appears even at t = 3500M . The order of the polynomial fit depends on the time range of the fit. In this case, a fifth order polynomial was enough to capture the oscillatory behavior in the eccentricity estimator in the time range 500M < t < 3500M . Note that the orbital phase could also be used to measure the eccentricity estimator using Eq.(13) (but without division by X fit ).

E. Eccentricity from gravitational waves
All eccentricity estimators discussed so far utilize coordinate-dependent quantities like separation or orbital frequency. Therefore, one might suspect that the higher harmonics visible in Figs. 2 and 3 are caused by gauge effects. The gravitational radiation at future null infinity is expected to be gauge-invariant, removing the dependence on gauge-dependent quantities. These considerations motivate the use of the gravitational wave phase and frequency to define eccentricity.
Then the gravitational-wave frequency is defined as The waveforms extracted at finite radii are extrapolated to null infinity using the procedure in [51]. The wave phase φ and frequency ω are measured as a function of the retarded time t − r * , where r * is the tortoise-coordinate radius defined as where M ADM is the ADM mass of the initial data. At early times, the gravitational waveforms are contaminated by high frequency noise from imperfect initial data.
To measure the amplitudes and locations of the extrema in the eccentricity estimator more accurately, the residual functions are filtered using a low-pass Butterworth filter with the Matlab function "filtfilt" [51]. The filtered data can be used to measure the eccentricity for retarded time t − r * > ∼ 1000M . Based on the gravitational wave phase, we define the eccentricity estimator where an additional factor of 1/2 arises because the wave phase is approximately twice the orbital phase. In Fig. 4, we plot the eccentricity estimator computed from the gravitational wave phase of the (2,2) mode extracted at the radii r = 75M , r = 240M and extrapolated to infinity using terms up to 1/r 2 versus t− r * . The eccentricity estimate is independent of the radius value at which the wave is extracted, and various estimates agree to within 5% in both amplitude and phase for different radii of extraction.
Using the wave frequency we define the eccentricity estimator e ω (t) Computation of the gravitational wave frequency ω = dφ/dt requires a derivative of φ(t), which increases numerical noise. Given the small amplitude of the effect under consideration (the fractional change in ω is 2e = O(10 −4 )), the increased noise noticably affects e ω . It is usable only at finite extraction radius, and even there only for t − r * > ∼ 2000M . In Fig. 5, we compute the eccentricity estimator from the wave frequency extracted at r = 75M and r = 240M . The extrapolated data to infinity is not shown because of its sensitivity to noise. The two curves have a nearly sinusoidal behavior with the phase agreeing to within 10%. Eccentricity estimator e φ computed from the gravitational wave phase as a function of retarded time t−r * . In this plot, the eccentricity estimator is computed from the gravitational wave extracted at finite radii r = 75M and r = 240M and from data extrapolated to infinity. The three curves agree in amplitude and phase to within 5% in the retarded time interval 1000M < t − r * < 3000M .
However, the amplitude differs by 25% between the wave data measured at r = 75M and r = 240M . The reduced sensitivity to noise is an important advantage of e φ over e ω . For the binary with eccentricity 0.05, plots similar to Figs. 4 and 5 with smooth sinusoidal behavior could easily be obtained.
Computation of the eccentricity from gravitational radiation (e φ and e ω ) is better behaved than the methods using orbital variables. Only one harmonic mode appears in the data-even for the low-eccentricity run with e ∼ 6 × 10 −5 -and the eccentricity is decreasing as the binaries inspiral toward merger. We attribute this improvement to the disappearance of coordinate and gauge effects when the data are extracted further away from the holes.
The eccentricities extracted from e φ and e ω in Figs. 4 and 5 are inconsistent with each other; they differ by a factor Ω r /Ω φ as explained in Sec. II A.
One might also consider a definition of the eccentricity based on taking the time derivative of the wave frequency. From Eq. (4), the second time derivative of the orbital phase is given by: where the amplitude of the oscillatory part is 2e M2 +Ṁ 4 . Then, the eccentricity estimator computed from the time derivative of the wave frequency e dω is then defined as The main advantage of a such a definition is that it requires a lower order fitting polynomial. Unfortunately, the numerical derivatives necessary to computeφ amplify noise, and so this method becomes impractical for the numerical evolutions considered.
In Table I, we summarize the eccentricity definitions examined in this paper, the data range between t i /M and t f /M employed in the fits, and the order of the fitting polynomial n for the 15-orbits quasi-circular nonspinning binary. We also give an estimate of the eccentricity value at t/M = 1000, 2000 and 3000 and its estimated error δe/e for each method.

III. BEHAVIOR OF ECCENTRICITY DURING INSPIRAL
Radiation reaction reduces eccentricity during the inspiral of a binary compact object, as shown by the post-Newtonian calculation by Peters [16]. Using the quadrupole approximation, Peters derived the evolution of the orbital eccentricity during the inspiral caused by the emission of gravitational waves. In the limit of small eccentricity, the eccentricity is related to the semi-major axis a by e ∝ a 19/12 .   In Sec. II, we established that the eccentricity estimators e φ (wave phase) and e s (proper horizon separation) show the cleanest oscillatory behavior. Using these two eccentricity estimators, we compute as follows the eccentricity as a function of time for the much longer inspirals considered here. We first define the "average" eccentricity over one half of a radial oscillation as the difference between two consecutive extrema (from minimum to maximum, or vice versa) of the eccentricity estimator We further associate this eccentricity with the time halfway between the two extrema under consideration: At the time of this average eccentricity, the separation is measured numerically. In the case when gravitational wave data is used, the wave phase is approximated as a function of the separation by using the retarded time t − r * . The results are plotted in Fig. 6. Fitting a power-law log e = α + β log s (24) to the numerical data yields β ≈ 1.4. These decay estimates are in reasonable agreement with Peters' prediction (β = 19/12 ≈ 1.583), as can be seen by the indicated power-law in Fig. 6. The orbital eccentricity decays similarly in the two simulations with different eccentricity.

IV. PERIASTRON ADVANCE
The periastron advance is one of the new features for relativistic eccentric orbits that is not present in Newtonian gravity. It has been computed analytically in the post-Newtonian regime up to third order but-to our knowledge-it has never been estimated numerically in binary black hole simulations. Periastron advance will lead to a modulation of the gravitational wave signal for eccentric binaries and will impact gravitational wave detection strategies. Therefore, it is important to know what this frequency is and how it changes as a function of the mass ratio. The fractional periastron advance per orbit, K, is defined as where ∆Φ = Φ − 2π is the periastron advance per orbit. The dimensionless parameter K is related to the radial frequency Ω r and the orbital frequency Ω Φ through A. Numerical method for measuring the periastron advance From an eccentricity estimator e X , cf. Eq. (13), one can read off not only the eccentricity (via the amplitude of e X ), but also the frequency of the radial motion, Ω r (from the oscillation period). We shall define the period of the radial oscillation as twice the time interval between two consecutive extrema (from minimum to maximum, or vice versa) in the eccentricity estimator curve. We employ the following procedure to compute the periastron advance: 1. Choose a cleanly oscillating eccentricity estimator e X (t). We will use e φ , cf. Fig. 4.
3. Interpolate the orbital phase Φ to the times t k . Between neighboring data points, the orbital phase changes by Φ(t k+1 ) − Φ(t k−1 ), whereas the radial phase changes by 2π. Therefore, the ratio between orbital and radial phase increase is (Φ(t k+1 ) − Φ(t k−1 ))/2π, and so For the very low eccentricity simulation (e ∼ 5×10 −5 ), the periastron advance is very difficult to measure because the amplitude of e φ is so small. The uncertainty in the extracted Ω Φ /Ω r is about ±0.1 for 0.02 ≤ M Ω Φ ≤ 0.03. The error in the estimated periastron advance increases at higher frequencies as the binary evolves closer to merger. The eccentricity estimators depend on details of the polynomial fits, and it is more difficult to read off these small eccentricity residuals near the plunge. Therefore, Ω Φ /Ω r for the e ∼ 5 × 10 −5 run should not be trusted for M Ω Φ > ∼ 0.03. In the simulation with larger eccentricity e ∼ 0.05, by contrast, the periastron advance is easier to measure, because the amplitude of e Φ is proportional to the eccentricity. We obtain correspondingly smaller errors, about 3% at frequencies M Ω < ∼ 0.03. While we are able to extract Ω Φ /Ω r at higher frequencies for the simulation with e ∼ 0.05, recall that the numerical data are constructed from consecutive extrema of e φ . At late times (close to merger), there is an increasing amount of orbital evolution during such an interval, which renders ambiguous both the definition of Ω Φ /Ω r and its association with one orbital frequency.  Figure 7 shows the computed periastron advance for the two equal-mass simulations considered here. To facilitate comparison with analytical estimates (see next section), we plot Ω Φ /Ω r as a function of orbital frequency. The latter is approximated as half the gravitational wave frequency. (This is justified because the deviation from this value is much smaller than the error in estimating the eccentricity and the periastron advance.) We will discuss this figure in the next subsection.

B. Results
From Fig. 7 we see that Ω Φ /Ω r is positive (i.e. the fully general relativistic calculation produces indeed a periastron advance), and the periastron advance increases with increasing orbital frequency M Ω Φ , again consistent with expectations. The solid and the dashed lines in Fig. 7 indicate the periastron advance for a test-mass orbiting a Schwarzschild black hole, and for an equal-mass binary at 3rd post-Newtonian order (see Appendix for details), and we can now compare these calculations with the fully relativistic BBH simulations. The scatter in the numerical data Ω Φ /Ω r represents a measure of the uncertainty in the periastron advance of the numerical simulations. For the e ∼ 0.05 simulation, this scatter is much smaller than the difference from the 3PN calculation. Therefore, we have positively detected a difference between fully numeric simulations and 3PN calculations. M Ω Φ (Ω Φ /Ω r from the e ∼ 10 −5 simulation coincides with the data for the e ∼ 0.05 run, although with larger scatter, because of trying to extract much smaller variations in the numerical data.) The difference between the numerical periastron advance curve and the 3PN result is about 3% at Ω Φ = 0.02 and continually increases to about 5% at Ω Φ = 0.03. The fully NR periastron advance seems to follow more closely the test-mass calculation than the equal-mass 3-PN prediction. Note that comparing either of the two analytic results is imperfect: The 3-PN calculation is for equal masses, but because of the nature of post-Newtonian perturbation theory becomes increasingly less reliable for increasing frequency M Ω Φ . The test-mass limit, in contrast, is an exact calculation, but for a system different from an equal-mass binary. Unequal mass binaries with mass ratios very different from unity should result in better agreement with the testmass limit, and we will explore this case next.
Extracting the periastron advance from a series of nonspinning unequal mass simulations [50], we obtain the data plotted in Fig. 8. These simulations have very low eccentricity in order to accurately model circularized binaries for gravitational wave data-analysis, with eccentricities indicated in Fig. 8. The smallness of the eccentricity is unfortunate for our purposes, as this increases the errors in the extracted periastron advance. The periastron advance for q = 2 is very similar to the equal mass periastron advance data. For higher mass ratio, the numerically computed Ω Φ /Ω r seems to increase and approach the test-mass result; however, the large uncertainty in Ω Φ /Ω r for these runs prevents us from drawing strong conclusions.

C. Laplace-Runge-Lenz vector
The Laplace-Runge-Lenz vector points towards the periapsis of the orbit from the center of motion, and therefore it would seem that observing this vector during a simulation would result in an immediate measure of the periastron advance. This vector is defined in ADM coordinates in terms of the canonically conjugate position R and momentum P as [39]: where L = R × P and µ is the reduced mass. Unfortunately, the magnitude of this vector is proportional to e, i.e., it will typically be very small. Moreover, it is computed as the difference between two large terms that almost cancel each other, resulting in large numerical errors. Furthermore, relativistic effects, such as gauge effects, might affect the two terms in Eq. (28) differently, thus disproportionately affecting the small difference A.
Yet another obstacle is that the numerical data do not give the canonical position and momentum. For all these reasons, we found it impossible to measure the periastron advance from the Laplace-Runge-Lenz vector even for the binary run with e ∼ 0.05.

V. DISCUSSION
We have dealt with three aspects of eccentricity in binary black hole simulations: how to measure eccentricity, its decay during the inspiral, and periastron advance.
With regard to techniques to measure eccentricity, this paper provides a systematic comparison between several different estimators. The ones shown in Figs. 2, 3, 4 and 5 each displayed a different behavior, even though these definitions reduce precisely to the usual eccentricity e in the Newtonian limit. Differences appear mainly because the data corresponds to a binary in the last phase of the inspiral before merger when relativistic effects are significant-a regime in which the Newtonian relations between the orbital variables are no longer valid.
The eccentricity estimator e BCP (see Fig. 2) exhibits two very undesirable features: For the e ∼ 0.05 simulation, e BCP depends strongly on the choice of how separation between the black holes is measured (coordinate distance D vs. proper separation s). For small eccentricities e ∼ 5 × 10 −5 , no regular oscillatory behavior is apparent, rendering e BCP useless as an eccentricity estimator. This might be because it uses a definition where eccentricity comes in the next-to-leading term, and the leading order Newtonian expression is not satisfied. Also, the high power of the contribution of the orbital variable makes the eccentricity easily affected by high-order harmonic modes in the orbital variables. We have observed similar behavior when we explored alternative definitions of the eccentricity based on Newtonian formulas combining orbital variables.
Eccentricity measures based on orbital quantities (see Fig. 3) give the right amplitude (for t < 2500M in the case of e Ω and e D ), and the phasing is quite consistent between the different eccentricity estimators. For instance, the orbital frequency is maximal when the separation is minimal. However, for the low-eccentricity simulation (e ∼ 5 × 10 −5 ) higher-order harmonics are clearly visible as the binary approaches the merger, in particular for the coordinate separation e D and the orbital frequency e Ω . The eccentricity measured from the proper horizon separation(e s ) is affected least by these coordinate effects.
Eccentricity measures based on extracted gravitational waves (see Figs. 4 and 5) result in clean oscillatory behavior, even for eccentricities as small as considered here. No high-order harmonics are noticeable in the wave extrapolated to infinity during the time interval considered. The eccentricity is calculated from the maximum and minimum values in the oscillating function without concern for the coordinate location in the orbit. It is especially straightforward to calculate numerically the eccentricity from the wave phase extrapolated to infinity without resorting to any notions of "distance" between the holes. Computing eccentricity from the gravitational wave phase is therefore the preferred method. Unfortunately, the gravitational wave phase is not as easily accessible as orbital quantities: One needs to extract gravitational waves, the waveform is delayed by the light-travel time to the extraction radius, and, for best results, one may have to extrapolate to infinity. Therefore, in practice, eccentricity estimators based on orbital quantities may be useful for immediate diagnostics during a simulation, then confirmed and refined subsequently by eccentricity estimators based on gravitational wave properties.
Notice that the eccentricity measurement could be affected by noise sources such as the "junk radiation" early in the simulation or by poor boundary conditions causing radiation reflection at the outer boundary. These additional oscillations could easily be interpreted as eccentricity. In principle, however, one should be able to distinguish them from the eccentricity by the frequency of the oscillation.
The second part of this paper describes measuring the decay of orbital eccentricity during the inspiral of equal mass non-spinning black hole binaries, revisiting earlier work [44]. For both simulations considered, we find that eccentricity measured via proper separation (e s ) and via gravitational wave frequency decays with the same power of proper of separation, s β , with exponent β ≈ 1.4. This is somewhat smaller than the value predicted by post-Newtonian expansions, 19/12 ≈ 1.58. The earlier work, which was based on fewer data-points at closer separation, found a distinctively smaller exponent when computing eccentricity from proper separation rather than from the orbital frequency.
The third part of this paper presents a measurement of periastron advance for equal and unequal mass nonspinning black hole binaries. For eccentric binaries, periastron advance will result in a characteristic modulation of the observed GW signal, and hence it is important to quantify its frequency. We find that the numerically computed periastron advance Ω Φ /Ω r disagrees with both 3PN predictions for equal-mass binaries, as well as with the test-mass limit of geodesic motion in a Schwarzschild background. As shown in Fig. 7, the periastron advance for black hole binaries lies roughly halfway between these two analytic calculations. The unequal mass evolutions considered have very small eccentricities; this is unfortunate for our current purposes, as this made it impossible to measure periastron advance well enough to test reliably the approach to the test-mass limit with increasing mass-ratio. While the data appear to approach the test-mass limit as the mass ratio deviates from unity, cf. Fig. 8, detailed confirmation will have to await until this analysis is repeated with somewhat higher eccentricity runs in the future. Nevertheless, even the equal-mass case shows that periastron advance is yet another feature of fully numerical calculations that is not accurately predicted by post-Newtonian expansions. To achieve agreement, one may have to go to higher order post-Newtonian expansions, or one may have to incorporate finite-size effects. More pragmatically, for applications to gravitational wave data-analysis, one might also introduce fitting parameters into the post-Newtonian models, and choose these parameters to enhance agreement with the numerical waveforms.
symmetric mass ratio ν ≡ m 1 m 2 /(m 1 + m 2 ) 2 , where m 1 and m 2 are the masses of the two bodies, as