Eccentric binary black holes Comparing numerical relativity and small mass-ratio perturbation theory

,


I. INTRODUCTION
Binary black hole (BBH) mergers have dominated the gravitational wave (GW) observations of the LIGO and Virgo detectors [1,2] in the first, second and the third observing runs [3][4][5][6].One key parameter of these astrophysical systems is the mass ratio q = m 2 /m 1 ≤ 1 of the binaries' components.Current GW observations [5,[7][8][9] predominantly find mass ratios close to unity with a few observations showing support for low mass ratios [10,11].
With the increasing number of GW detections in the upcoming observing runs by ground-based detectors [7,12], and space-borne detectors, like the LISA observatory [13,14], it is likely that more binaries with mass asymmetries are found.In particular, LISA will be sensitive to binaries with mass ratios ranging from q ∼ 1, over intermediate mass-ratio systems (q ∼ 10 −3 ) to extreme mass ratio inspirals at q ∼ 10 −5 .Furthermore, third-generation ground-based detectors with improved low frequency sensitivity relative to today's ground-based detectors will be able to detect the capture of stellar mass black holes (BHs) by intermediate mass BHs with mass-ratios down to q ∼ 10 −3 [15].Thus, the modelling of GWs from BBHs at all mass ratios is of preeminent relevance for a correct detection and analysis of these sources.
Orbital eccentricity is another important parameter describing binary systems.While emission of gravitational waves generally does reduce eccentricity [23,24], orbital eccentricity is a key parameter to constrain the formation scenario of these binaries and thus the astrophysical origin of GW sources [25][26][27][28][29][30][31].For current ground-based detectors there are ongoing efforts to search for signatures of orbital eccentricity in the detected GW signals [32][33][34][35][36][37][38][39][40][41][42].For future GW detectors, especially for high-mass ratio binaries in LISA, it is expected that emission of GWs has not circularized most binaries yet.Therefore, the correct modelling of orbital eccentricity effects is fundamental to accurately describe such systems in future detectors, in particular for extreme mass-ratio inspirals, which are described by SMR perturbation theory.
Recently, Ref. [43] demonstrated that NR simulations at modest mass-ratios (q 0.1) can be used to gain insight into the accuracy of the SMR expansion, confirming the known leading-order term, and predicting next-to-leading order contributions.Ref. [43] considered non-eccentric (quasi-circular) binaries only, with both BHs non-spinning.Here, we begin to extend the analysis in [43] to eccentric BBHs, while still keeping both BHs non-spinning.The non-circularity of the binary's orbit introduces a new timescale to the two-body problem, the timescale of the periastron precession, which induces oscillations in the dynamical and GW quantities of the binary system complicating substantially the analysis relative to the quasi-circular case described in [43].
An additional difficulty arises from the fact that eccentricity is a gauge dependent parameter in general relativity, thus complicating the comparison between SMR evolutions and NR arXiv:2209.03390v1[gr-qc] 7 Sep 2022 simulations.In order to overcome this problem, we develop tools to extract gauge invariant quantities from both SMR and NR waveforms.Using PN theory, we derive relations among the eccentricity defined from the orbital and (2,2)-mode gravitational wave frequency, and the PN temporal eccentricity.We show that a commonly used definition of eccentricity based on the (2,2)-mode frequency -Eq.( 5) below-does not reduce to the Newtonian definition of eccentricity.We therefore adopt a new definition of eccentricity, e gw in Eq. ( 6) below.This new definition continues to be based on the frequency of the (2,2) GW-mode, but also satisfies the correct Newtonian limit.
NR simulations of BBHs have been routinely performed since more than a decade ago.Motivated by expectation of small eccentricities for most of the GW signals in the frequency band of ground-based detectors, most of the NR groups have focused on the production of simulations of quasi-circular BBH mergers [44][45][46][47][48][49][50][51][52][53][54][55], with the exception of a limited number of studies exploring BBH coalescences in eccentric orbits [55][56][57][58][59][60][61][62][63].The Spectral Einstein Code (SpEC) [64] is an accurate and efficient NR code that has been used to study quasi-circular inspirals in great depth [47,52].Eccentric inspirals at low eccentricity were studied to some extent [57][58][59][60].We expand SpEC's capabilities for accurate simulations of binaries with larger eccentricities, 0.2 e 0.8, which are characterized by large variations of GW frequency and amplitude between apastron and periastron passages.We have produced a set of 52 non-spinning eccentric simulations between mass ratios 1 : 10 and 1 : 1, and with initial eccentricities of up to 0.7.The number of orbits is typically 20, yielding a dataset with the longest evolutions and highest initial eccentricities up to date, which sets it also apart from the simulations of other groups [55,61,63].
The main purpose of this article is to compare NR and SMR calculations.This requires a map from the instantaneous state of a NR simulation to the geodesic on which the point-particle instantaneously moves in its motion around the central black hole.We characterize the instantaneous state of SMR and NR simulations by symmetric mass ratio, ν = m 1 m 2 /(m 1 + m 2 ) 2 , eccentricity, e gw , and orbit-averaged frequency of the 22mode, ω 22 .These quantities can be uniquely determined in SMR and NR configurations and they generate an unambiguous map between SMR and NR configurations, as described in Sec.V.
We find that the leading order prediction in the SMR expansion for the energy and angular momentum fluxes agree with the NR results to within 10%.The next-to-leading order SMR contributions to the fluxes can be estimated by rescaling the difference of the NR and leading order SMR contribution by a factor of the symmetric mass-ratio.The result has a very small dispersion in symmetric mass ratio, which implies that the next-to-next-to leading order SMR contribution is small, even for comparable masses.This is compatible with the findings of [65] in the quasi-circular case.Comparing the zeroeccentricity limit of our next-to-leading order estimate to the exact results of [65] we find the results to be comparable with an overall small shift likely due to the orbit-averaging procedure applied to extract quantities from the eccentric NR simu-lations.A similar analysis is done for the periastron advance.In this case the NR results are within 8% of the leading-order (geodesic) SMR result, the next-to-leading SMR contribution is compatible with previous exact calculations in the quasicircular limit [66,67], and the next-to-next-to-leading SMR contribution appears small in the comparable mass regime.
This article is organized as follows.In Sec.II we present a detailed description of the new dataset of eccentric nonspinning NR simulations produced for this work.We investigate the relations between different eccentricity definitions in Sec.III, and we provide a definition of eccentricity based on the (2,2)-mode frequency, which reduces to the Newtonian definition of eccentricity in the Newtonian limit.Section IV describes the SMR evolutions performed in this work, and in Sec.V discusses the mapping between SMR and NR configurations.In Sec.VI we compare the quantities extracted from the NR simulations to the SMR perturbation theory results, and provide constraints on the values of the next order terms in the SMR expansion for the GW energy and angular momentum fluxes, as well as the periastron advance.In Sec.VII we summarize our main conclusions and discuss future work.The appendices contain additional technical details: Appendix A describes our method to set the initial parameters in the NR simulations, in Appendix B we assess the quality of the NR waveforms, and in Appendix C we provide details of the derivation of the relations between different definitions of eccentricity using PN theory.

II. NUMERICAL RELATIVITY SIMULATIONS
The NR simulations produced in this work are performed with the SpEC code [64], utilizing numerical techniques summarized in [47,52].In particular, SpEC evolves a firstorder representation of the generalized harmonic evolution system [68] using a multi-domain spectral method [69][70][71][72].At the outer boundary constraint-preserving boundary conditions [68,73,74] are employed, whereas black hole excision is used inside the apparent horizons [70][71][72]75].The transition to ringdown is accomplished with the techniques described in Refs.[70,72].Initial data are constructed with the eXtended Conformal-Thin Sandwich (XCTS) approach [76][77][78], and we describe in Sec.II A and Appendix A how we achieve binaries with a desired value of orbital eccentricity.
For improved performance for eccentric systems we adopt part of the modifications developed in [79] to produce accurate simulations of hyperbolic encounters.Most notably, adaptive mesh refinement and GW output is triggered more frequently to adjust to periastron passages which happen on fast timescales, and which cause pulses of higher-frequency GWs that travel through the computational grid.

A. Numerical relativity dataset
We have produced 52 new numerical relativity simulations of binary black holes on eccentric orbits.The simulations are summarized in Table I; for each of the mass-ratios q = Table I.Properties of the NR simulations used in this work.Columns 2-5 give the initial data parameters needed to reproduce each simulation (see main text), whereas columns 6-10 give some physical properties: the number of orbits, N orbits , the time to merger T merger /M, and the reference time T ref /M corresponding to a frequency of the (2,2)-mode ω ref 22 = 0.042, at which the eccentricity e ref gw and mean anomaly l ref /(2π) are extracted from the simulation.Here, M is the total mass of the binary after initial transients have settled down [52].These parameters and additional properties can be found at numerical precision in the metadata files accompanying each simulation.(*) Simulation performed with SKS initial data, differently from the rest of simulations, which used SHK initial data (see main text for details).SXS:BBH:2528, q = 1, egw = 0.45 SXS:BBH:2518, q = 1, egw = 0.03 SXS:BBH:2558, q = 1/6, egw = 0.46 SXS:BBH:2552, q = 1/6, egw = 0.02 SXS:BBH:2567, q = 1/10, egw = 0.34 SXS:BBH:2564, q = 1/10, egw = 0.01 0.025 0.000 0.025 t (s) Figure 1.Visualization of simulations at two eccentricities each for three different mass-ratios.Shown is h + at inclination angle ι = π/3 and coalescence phase φ = 0, for a binary of total mass of 60M at a distance of 430 Mpc.For ease of plotting, the waveforms are offset vertically.On each waveform, the location is marked where the orbit-averaged GW frequency ω 22 equals our reference value Mω ref 22 = 0.042; for M = 60M this corresponds to a GW frequency of 22.6 Hz, near the start of the frequency band of current GW detectors.The right panel enlarges the merger part of the signals.
1, 1/2, 1/3, 1/4, 1/6, 1/8 and 1/10, simulations with several different eccentricities e gw are computed.Within the XCTS formalism to construct initial data, the simulations in Table I were produced using superposed harmonic Kerr (SHK) initial data [80], except for the simulation SXS:BBH:2527, which used superposed Kerr-Schild (SKS) initial data [81] as this is the simulation with largest initial separation and eccentricity, and initial tests with SHK initial data were not successful 1 .
For each simulation, Table I reports on the parameters values necessary to reproduce the initial data with the techniques described in [82]: the inverse mass ratio 1/q = m 2 /m 1 ≥ 1, the orbital separation D 0 /M 0 , where M 0 is the initial ADM mass, the initial orbital frequency M 0 Ω 0 , and the initial radial velocity parameter a 0 [83,84].The procedure to determine the initial parameters of the simulations is described in Appendix A. The simulations are started at or very close to apastron due to limitations of the radial map used by the dual-frame method [85] employed to solve the Einstein equations in SpEC [86].Specifically, the radial mapping of Eq. ( 9) in [86] connecting the comoving and inertial frames does not allow the orbital separation to increase more than 1.5 times the initial separation.We note that this limitation has been recently overcome in SpEC by defining a new radial map, however, it is not applied for simulations in this publication, and we leave it to future work to report on this new feature.
To convey a sense of the physical properties of the BBHs studied, Table I also lists the number of orbits to merger, 1 After tuning some settings in the linear solvers of the SpEC initial data code, the SHK initial data was successfully computed, but in order to save computational resources the evolution with SHK initial data was not produced.
N orbits and the time to merger T merger /M, where M is the total mass.We also specify the time (before merger) T ref /M where the orbit averaged frequency of the (2,2)-mode reaches the value ω ref 22 = 0.042, as well as eccentricity e ref gw and mean anomaly l ref /(2π) at this reference time.These quantities are defined with the procedures outlined below in Sec.III.The reference frequency is chosen to be consistent with the length2 of the shortest simulation, which corresponds to SXS:BBH:2520 with 4963M of evolution and 18 orbits.Apart from this particular case, most of the simulations have typically a time to merger > 10 4 M.This makes our dataset of eccentric NR waveforms the one with the longest evolutions of eccentric binary black holes to date.
We extract the gravitational radiation from each simulation using the same techniques as in [52], and decompose Each mode h lm is further split into real amplitude and phase as with an associated GW mode frequency of A sample of the computed numerical waveforms are shown in Fig. 1.One can observe that the highly eccentric configurations develop a very complex structure in the waveform due to the eccentricity of the orbits followed by the BHs. Figure 1 also shows a zoom-in of the merger part on the waveforms, to highlight the similarity of the merger and ringdown parts of the waveform with different eccentricities 3 .Merger and ringdown of the high and low eccentricity inspirals agree well with each other, indicating that the circularization hypothesis is accurately fulfilled for our dataset, consistently with the findings in [56,61,63].We note that recently some unexpected dependence of the kick velocity on eccentricity was found in [87].A similar analysis of the kick velocity can be performed on our dataset, and we leave such study of the final velocity as well as other remnant properties for future work.

B. Eccentricity, azimuthal frequency & mean anomaly
We start with the eccentricity definition proposed by Mora & Will [88], -mode frequency at periastron and apastron are indicated with orange and black dots, respectively.These are used to compute the orbit-averaged frequency of the (2,2)-mode (solid red curve), and the eccentricity e gw (dashed green curve) through Eq. ( 6).Bottom panel: Time evolution of the mean anomaly (solid purple line) computed using Eq. ( 10) for the same simulation as in the top panel.The vertical dashed gray lines in both panels correspond to the times of the periastron passages.
where Ω p orb and Ω a orb are the values of the orbital frequency at consecutive periastron and apastron passages, i.e. maxima and minima of Ω orb (t).Equation ( 4) is easy to compute from orbital trajectories and reduces precisely to the normal eccentricity in the Newtonian limit [88].e Ω orb was for instance used in [58] to analyse generic precessing & eccentric BBH inspirals.To avoid the coordinate-dependence of Ω orb , recent papers (e.g.[63]) have applied Eq. (4) to frequencies directly defined from the gravitational radiation: where ω a , ω p refer to the (2,2)-mode frequency ω 22 at apastron and periastron, respectively.This procedure is illustrated in the top panel of Fig. 2: The time-dependent ω 22 (t) has maxima ω 22 p ,i and minima ω 22 a ,i indicated with the black and orange dots, where the integer i labels the extrema.The maxima and minima correspond to periastron and apastron passages, respectively, and occur at times t p i and t a i .We show below in Sec.III that e ω 22 disagrees with e Ω orb ; most notably, e ω 22 does not have the correct Newtonian limit.Therefore, we introduce a new eccentricity definition e gw measured from the frequency of the (2,2)-mode, which has the correct Newtonian limit: This new gravitational-wave frequency e gw is also plotted in the top panel of Fig. 2. The dashed curve for e gw is obtained by constructing interpolating functions through all maxima {ω 22 p i } and through all minima {ω 22 a i }, and then evaluating Eqs. ( 5) and (6) for these interpolating functions.
The average azimuthal frequency from the (2,2)-mode for the interval between the i-th and i + 1-th periastron passages is defined as We associate this frequency with the temporal midpoint and interpolate the discrete ( tp i , ω 22 i ) data to obtain a continuous ω 22 (t) curve.This curve is also included in Fig. 2.
The mean anomaly of the eccentric binary is defined as [89] where t p i and t p i+1 are the times of the periastron-passages immediately before and after the time t of interest, and is plotted in the lower panel of Fig. 2.
The NR quantities introduced so far are used in Fig. 3 to illustrate the entire NR dataset produced in this work.Figure 3 shows the tracks of each simulation in the parameter space spanned by the orbit-averaged (2,2)-mode frequency, ω 22 and the eccentricity, e gw .Each simulation is color-coded by its mass ratio.We also indicate the value of the mean anomaly at the reference frequency used to perform the analysis.One can observe that the mean anomaly at the reference frequency is randomly distributed.We assess the accuracy of the simulations by computing the unfaithfulness between waveforms at different resolution in Appendix B, and we obtain that our dataset of simulations has a median maximum mismatch between different resolutions of < 10 −3 , indicating a convergent behavior of the waveforms with increasing resolution.

C. Quantities for comparisons with small mass-ratio theory
In our comparisons with small mass-ratio perturbation theory, we will also utilize several more quantities extracted from the NR simulations.We define an orbit-averaged radial frequency based on the periastron passages as which is interpolated to a continuous Ω r 22 (t) curve.From this, we compute periastron advance K as the ratio between the azimuthal and radial frequencies [90], The instantaneous energy and angular momentum fluxes are computed from the GW modes, h lm , using the expressions [91] Ėgw = 1 16π where ḣ = dh/dt, the indicates the imaginary part and ḣ * lm denotes the complex conjugate of ḣlm .In the case of nonspinning binaries only the z-component of the angular momentum flux is non-zero.Analogous to Eq. ( 7) we define the orbit average of either of these fluxes as, We associate these discrete averages over each radial oscillation period with the mid-time tp i , and interpolate to obtain continuous functions Ėgw (t) and Jgw z (t).A first estimate of the peaks is computed using an envelope subtraction method as in [58].Each estimate of the peak is used to set a window of ∼ 30M on which a polynomial fit is performed.Finally, this polynomial fit is used to compute the value of the peak.

III. DISCUSSION ABOUT ECCENTRICITY DEFINITIONS
There is a large variety of measures of eccentricity in use in general relativity [92].Many of these measures derive from the trajectories of the binaries, and are therefore coordinate dependent.This makes them generally unsuitable for comparisons between different modelling approaches, which may be computed in different gauges or where there may be no well-defined notion of trajectory at all.However, one gauge invariant observable common to all approaches to modelling gravitational waves from compact binaries is the waveform itself.In this sense, it may seem more reasonable to define eccentricity in terms of gravitational wave quantities rather than quantities dependent on the trajectories of the black holes.
A gravitational wave mode (see Sec. II), has an instantaneous frequency ω lm = φlm , which can be related in the inspiral regime to the instantaneous orbital Ω orb = φorb by the approximation [16] However, as eccentricity increases the approximation of Eq. ( 16) is no longer valid as can be observed in the top panel of Fig. 4, where the left and right-hand sides of Eq. ( 16) in the case of the (l, m) = (2, 2) multipole are displayed.In the top plot of Fig. 4, the upper and bottom panels correspond to a q = 1/6 configuration with two different initial eccentricities e 0 ω 22 = 0.03, 0.63, respectively.The relation between the orbital and the (2,2)-mode frequency is no longer the simple factor 2, as in the quasi-circular case.In order to derive the relation between both frequencies in the more generic eccentric case, we use PN theory.Specifically, we compute ω 22 at 1PN order using the instantaneous gravitational modes from [93].We obtain a 1PN-accurate expression for ω 22 in harmonic coordinates of the form, where r denotes two time derivatives on r.The explicit expression for F is given in Eq. (C6) in Appendix C 1 together with details of the derivation.Because of φ = Ω orb , Eq. ( 17) is a relation between ω 22 and Ω orb .
The top panel of Fig. 4 shows that the use of Eq. ( 17) with NR coordinates (ω 1PN, c NR 22 in the figure) agrees notably better with the (2,2)-mode NR frequency than 2Ω orb .The deviations in Eq. ( 16) increase with eccentricity.The relative error can be larger than 10%, whereas Eq. ( 17) leads to differences smaller than 1%.
It is important to note that the scaling relation in Eq. ( 16) between orbital and gravitational wave frequencies is still satisfied in an orbit-averaged sense.This is shown in the upper panel of Fig. 4 for the orbit-averaged frequencies ω 22 (solid red lines) and Ω orb (blue dots).
Let us now turn to eccentricity defined from the extrema of a frequency.Equation ( 5) can be evaluated from ω 22 (as written), or from the orbital frequency Ω orb .Because ω 22 (t) and Ω orb (t) have modulations of different amplitude (as seen in the top panels of Fig. 4), the corresponding eccentricities e ω 22 and e Ω orb are also different, as visible in the lower panels   I) with two different initial eccentricities.For each simulation twice the orbital frequency 2Ω orb (blue solid lines), the frequency of the (2,2)-mode (red solid lines), ω 22 , and the 1PN expression for the frequency of the (2,2)-mode from Eq. ( 17) evaluated using the NR coordinates, ω 1PN, c NR 22 , (black dashed lines) are shown.Additionally, the orbit-averaged values of the frequency of the (2,2)-mode (red dashdotted lines), ω 22 , and twice the orbit-averaged orbital frequency (blue dots), 2 Ω orb , are displayed for each configuration.Bottom panel: Eccentricity evolution computed from the orbital and (2,2)-mode frequencies using Eqs.( 5) and ( 6), and the 1PN expression for the eccentricity of the (2,2)-mode computed from (18) using NR coordinates, e 1PN, c NR ω 22 (black dashed lines). of Fig. 4. 4 Given the remarkable agreement of the PN approximation to ω 22 with respect to the NR values, one can insert Eq. ( 17) into the right hand side of Eq. ( 5), expand the corresponding expressions up to 1PN, and obtain an approximation 4 The eccentricity curves in the lower panels show a spurious bump close to merger arising from the interpolation of the maxima and minima close to the plunge.Our analysis focuses on the inspiral regime and is not affected by this feature.We leave to future work the improvement of the eccentricity measurement in the transition from inspiral to plunge.
for e ω 22 in terms of the coordinates as, where the expression for G is given by (C14) in Appendix C 1, and the subscripts/superscripts a, p refer to the apastron and periastron, respectively.The bottom panel of Fig. 4 shows that Eq. ( 18) successfully reproduces e ω 22 .Given the overall agreement, we do not pursue to explore higher PN orders, or possible resummations of this PN expression to improve its behavior in the strong field regime, and we leave possible extensions of these expressions, like the inclusion of spin effects, for future work.
The relations in Eqs. ( 17) and ( 18) allow one to obtain an estimate of the eccentricity measured from the (2,2)-mode frequency from the coordinates of the system.This can be useful, for instance, to set an eccentricity reduction or eccentricity control procedure based on the eccentricity measured from the waveforms instead of the trajectories without having to evolve the system such that the gravitational waves reach the extraction radii, and thus, saving computational time.
A detailed derivation of the relation e Ω orb − e t up to 3PN order for non-spinning binaries can be found in Appendix C 2. We focus here on the relation e ω 22 − e t , which is derived up to 1PN order in Appendix C 3, providing where x = Ω 2/3 orb , and γ = 1/c 2 is a bookkeeping parameter identifying the 1-PN corrections.At Newtonian order, Eq. ( 19) reduces to While Eq. ( 20) achieves the right limits for circular and parabolic orbits -e 0PN ω 22 (e t = 0) = 0 and e 0PN ω 22 (e t = 1) = 1it disagrees otherwise.This can be easily seen by expanding Eq. ( 20) for small eccentricities, which explicitly demonstrates that for small eccentricities in the Newtonian limit, e ω 22 does not reduce to e t , but rather to 3/4e t . in the Newtonian limit, but its relation is e ω 22 ∼ 3e t /4.An expansion of Eq. (20) in the large eccentricity limit 1 − e t 1 yields which also exhibits a wrong slope ( √ 3) for e t near 1.Equations (21) and (22) show that the definition of eccentricity based on the (2, 2)-mode frequency will be different from the Newtonian definition of eccentricity in the two limits of the bound case.Additional PN orders will introduce higher frequency corrections to the Newtonian behaviour, whose impact in the leading Newtonian correction factors between the eccentricity will depend on the region of the parameter space considered.
The relation e t (e ω 22 ) at Newtonian order can be obtained by inverting Eq. ( 20), Applying Eq. ( 23) to e ω 22 will yield an eccentricity-definition that reduces to the Newtonian definition of eccentricity.
As a consequence of the previous analysis we propose a new definition of eccentricity measured from the frequency of the (2,2)-mode, which corrects the naive result e ω 22 obtained from the extrema of ω 22 by Eq. ( 23), By construction, e gw reduces to the Newtonian definition of eccentricity in the Newtonian limit.In the bottom panel of Fig. 4, e gw is shown to be closer to e Ω orb than e ω 22 .Both e gw and e Ω orb have the correct Newtonian limit, and the differences may be explained due to coordinate effects affecting e Ω orb , and higher PN terms, as e gw is obtained from Eq. ( 20).This new definition of eccentricity is adopted throughout the rest of the paper, and its applications to data analysis are further investigated in upcoming work [103].

IV. SMR THEORY AND DATA
In the small mass-ratio (SMR) limit, the dynamics of a black hole binary can be described through the gravitational self-force formalism.For the inspiral part of the waveform, this formalism leads to a systematic expansion of the waveform in integer powers of the symmetric mass-ratio ν.This expansion is known as the post-adiabatic (PA) expansion.In this section, we introduce the necessary parts of this formalism to produce SMR eccentric inspirals for comparison to our NR data.For a more in depth review of the formalism see e.g.[21,104].

A. Equations of motion
In the SMR limit an eccentric inspiral of non-spinning black holes can be described as a series of evolving (perturbed) eccentric orbits in a Schwarzschild background.Eccentric orbits in Schwarzschild are often identified by their semi-latus rectum p and geodesic eccentricity e g , which in turn are defined through the periastron and apastron positions, r p and r a , p = 2r a r p r a + r p , and The position along the eccentric orbit is tracked by a phase q r conjugate to the radial action, defined such that q r = 0 mod 2π corresponds to the orbit being at periastron.The equations of motion for the evolution of the inspiral can be described as an expansion in the symmetric mass ratio ν (keeping the total mass M fixed), where t is retarded time at future null infinity, Ω r geo and Ω geo are the geodesic radial and azimuthal frequencies (w.r.t t), and the F's and f are the first order (gravitational self-force) corrections to the equations of motion.
By applying a near-identity (averaging) transformation Eqs.(27) can be put in an orbit averaged form (without loss of generality) [105].The leading terms give rise to the adiabatic (or 0-post-adiabatic, 0PA) approximation to the inspiral equations of motion, The next order in ν in the approximation -the 1-postadiabatic or 1PA order -requires knowledge of the average parts of the second order F p and F e g , i.e. the second order gravitational self-force.Despite major progress in calculating the second order self-force and corresponding 1PA corrections for non-spinning quasi-circular inspirals [65,106,107], there are no second-order self-force results yet for eccentric inspirals.Without the input of the second order self-force, any 1PA corrections based purely on the conservative part of the firstorder self-force are not gauge invariant [108], and not suitable for comparison with NR.Consequently, for this work we will focus on comparisons with the adiabatic (0PA) SMR results.

B. Gravitational wave strain
The gravitational wave strain produced by a test particle orbiting a Schwarzschild black hole can be found by solving the Teukolsky equation for ψ 4 .We write ψ 4 at future null infinity as where mn = mΩ geo + nΩ r geo , and Z lmn are the mode amplitudes.The strain-modes at infinity, Eq. ( 1), are then given as where in the last step we have written the strain explicitly in terms of the variable evolved by Eq. ( 28).To obtain the strain produced by an adiabatic (0PA) inspiral, one simply elevates the geodesics variables (p, e g , q r , q φ ) in Eq. ( 31) to their inspiral (evolving) counterparts in Eq. ( 28).

C. SMR data and interpolation
To produce SMR 0PA waveforms 5 we need the various quantities appearing on the right-hand sides of Eqs.(28) and (31).The 0th order "frequencies" Ω are known analytically [109], while the F p , F e g , and A lmn need to be calculated numerically.All three may be obtained by solving the Teukolsky equation sourced by a test mass following an eccentric geodesic to obtain the Z lmn 's in Eq. ( 29), which we do using the arbitrary precision frequency domain code developed in [110][111][112].
Specifically, we calculate F p , F e g , and Z lmn on a grid of Chebyshev nodes in x = (MΩ geo ) 2/3 (18 nodes between 0.001 and 0.130) and e g (12 nodes between 0 and 0.5), and interpolate the results using Chebyshev polynomials.The resulting interpolant has a typical relative interpolation error of about 10 −5 .
Note that the SMR 0PA inspiral waveforms generated here could in principle have been generated with Fast EMRI Waveforms (FEW) framework [113][114][115].We chose a different approach because FEW was not yet publicly available when this project started and to retain a better control over numerical errors in the model.In particular, the FEW model was not designed to faithfully reproduce the minima and maxima of the waveform frequency ω 22 .
From a (0PA) SMR inspiral we have two distinct ways of obtaining the average orbital and radial frequencies.We can apply the procedure of Secs.II B and II C to extract the average orbital ω 22 and radial frequencies Ω r 22 from the SMR 0PA waveform.We will denote these frequencies ω 0PA and Ω r 0PA .Alternatively, we have the instantaneous geodesic frequencies Ω geo and Ω r geo as they appeared in Eq. ( 28).In the ν → 0 limit, i.e. when there is no inspiral, Eq. ( 31) gives the following expression of the waveform frequency ω lm , = n (m dφ dt + n dq r dt )A lmn (p, e g )e −i(mφ+nq r ) n A lmn (p, e g )e −i(mφ+nq r ) (33) From this we note that the waveform frequency is exactly 2π periodic in q r , and consequently the radial period is exactly 2π/Ω r geo .A less obvious observation is that the average of the second term in (34) vanishes after averaging over a radial period.A sufficient condition for this to be true is since this guarantees that n A lmn (p, e g )e −inq r is confined to a half of the complex plane and must return to the same complex argument after one period.The condition ( 35) is clearly satisfied for low eccentricity orbits since A lmn = O(e g n ).However, condition (35) is easily violated by high eccentricity zoomwhirl orbits.Nonetheless, we observe empirically that the average of the second term (34) vanishes in all geodesic waveforms used in this work.
We thus find that in the ν → 0 limit we have exactly, This, of course, does not come as a surprise, since this is precisely what the frequency recovery procedure of Secs.II B and II C was designed to achieve.However, using the SMR 0PA inspiral waveforms we can now investigate what happens for finite values of ν when the system is evolving.Figure 5 shows both the frequencies, Ω r 0PA and ω 0PA /2, recovered from a SMR 0PA waveform at equal mass (ν = 1/4) and the geodesic frequencies, Ω r geo and Ω geo , inferred from the underlying inspiral dynamics.Even at equal mass there is hardly any perceivable difference between the two sets of frequencies.
To compare the frequencies obtained through the two procedures more closely we pick three frequencies along the adiabatic inspiral depicted in Fig. 5.For each of these frequencies we generate a series of adiabatic inspirals with symmetric mass-ratios varying between ν = 10 −3 and ν = 1/4 going through that point (and randomized initial values of q r ).For each of these inspirals we extract the azimuthal and radial frequency from the waveform using the procedure of Secs.II B and II C. Figure 6 shows the difference between these frequencies and the corresponding values obtained directly from the underlying geodesic.We observe a small, but measurable, difference between the two sets of frequencies, which appears to grow linearly with ν and is larger for higher frequencies.Since the SMR 0PA waveform contains no higher order frequency corrections, this difference arises purely from unintended side effects of the frequency recovery procedure.Some contributing factors are the averaging over a radial period while the inspiral is evolving, and limitations in establishing a radial period in the first place.

E. Eccentricity
To calculate the gauge invariant eccentricty e gw for a SMR 0PA waveform we again have two options.First, we can follow the procedure of Secs.II B and II C to determine the minima and maxima of ω 22 of the SMR 0PA waveform, and compute e gw using Eqs.( 5) and ( 6).We will refer to this as e 0PA gw .Alternatively, we want to obtain e gw directly from the dynamical variables p and e g .Unfortunately, there is no analytic closed form expression for e gw in terms of p and e g .Instead we start from the (numerical) "snapshot" waveform generated by a test particle going around a geodesic with fixed p and e g .The snapshot waveform h lm is a biperoidic function of the radial and azimuthal phases q r and φ as described by Eq. (31).Using the expression for ω 22 in Eq. ( 34), we find the minima and maxima of the frequency with respect to q r and calculate the corresponding e ω 22 , which can be input to (6) to provide e gw .We will refer to this quantity as e Figure 7 explores the difference between e geo gw and e 0PA gw for adiabatic inspirals.As expected, the difference between these two approaches for obtaining e gw vanishes in the ν → 0 limit.For ν 0, this difference grows again proportional to ν, similarly to Fig. 6.

F. Geodesic snapshot vs. inspiral waveform quantities
The preceding subsections have explored the difference between extracting the frequencies Ω φ/r , and eccentricity e gw from evolving adiabatic (0PA) waveforms and extracting the same information from geodesic "snapshots" that are not evolving at all.The relative difference between the two methods is found to be O(ν).For the comparisons in the rest of this work we choose to work with the geodesic snapshot SMR quantities, since these can in general be obtained more efficiently and reliably.For the leading order comparisons this will not make a difference.However, for any higher order corrections that we infer, we must be aware that these also contain a next-to-leading order correction due to comparing NR quantities from an evolving waveform with SMR geodesic snapshot quantities.

V. CHOOSING "INDEPENDENT" VARIABLES
In this section we study several options for the variables describing the state of a binary inspiral.We present some of the choices of variables made in the literature when comparing SMR and NR results, discuss their applicability in the eccentric case and, finally, describe the choice of variables which better adapt to our study.
The instantaneous state of a non-spinning eccentric binary is captured by four dynamical variables.For example, in the SMR setup these are the (p, e g , q r , φ) that appear in Eq. (27).In this work we are interested in comparing quantities that are observable for a distant observer.Since the instantaneous value of φ is completely degenerate with the position of this observer, it carries no useful information about the state of the binary.Moreover, we are presently interested in observables that are integrated over a radial cycle, eliminating q r .Thus we can identify the instantaneous state of the binary with two variables, like (p, e g ).Of course, (p, e g ) are not gauge invariant, and therefore not useful to find an NR simulation in the same instantaneous state.In order to compare the SMR and NR results, we need a set of two variables that can fix the instantaneous state of the binary and be unambiguously computed both in NR and the SMR formalism.
One pair of variables extensively used in the literature [58,66] are the azimuthal and radial frequencies Ω orb and Ω r orb defined in analogy to Eqs. ( 7) and ( 11) from the orbital frequency Ω orb .These two frequencies can be calculated analytically at geodesic order in the SMR formalism [109], and they can be extracted from NR data [58].However, since they are derived from the coordinate trajectories in NR, they are not fully gauge invariant (e.g.Fig. 17 of [116]).
A second possibility are frequencies computed from the gravitational radiation instead, e.g.ω 22 and Ω r 22 , which are manifestly gauge invariant.Then, as shown in Secs.III and IV, the orbit-averaged azimuthal frequencies from the waveform and the trajectories can be related by a factor 2, while the radial frequency stays the same.These frequencies are plotted in the top panel of Fig.   I, whereas the grey shaded area indicates the region covered by Schwarzschild geodesics, bounded by the diagonal black curve representing quasi-circular geodesics.Bottom panel: Eccentricity, e gw , computed using Eqs.( 5) and ( 6), as a function of ω 22 , for the same NR simulations as in the upper panel.While geodesics exist at all eccentricities, we have only generated SMR configurations in the grey shaded area.In both panels each NR simulation has been color-coded according to its symmetric mass ratio ν.
the corresponding circular orbit at the same ω 22 , i.e. the NR frequencies fall outside the range spanned by geodesics.It might be possible to rectify this situation by applying a linear mass-ratio gravitational self-force correction to the NR frequencies.This would, however, result in a very convoluted analysis requiring SMR inputs on the NR side of the comparison.Thus, to avoid such a complication we discard the radial and azimuthal frequencies as independent variables to describe both NR and SMR eccentric inspirals.
A third possibility as a pair of independent variables are the binding energy, E b , and the dimensionless angular momentum, j.This pair of variables has been extensively used for comparisons between NR simulations and effective-onebody (EOB) evolutions [117][118][119].Both quantities can be analytically calculated at geodesic order in the SMR formalism [109], and they can also be extracted from the NR simulations.Nonetheless, the computation of the reduced angular momentum and the binding energy from NR simulations requires the application of some unknown offsets to both quantities.This is due to the fact that E b and j are reconstructed by integrating the fluxes to infinity and using the initial (ADM) or final mass and angular momentum.However, the fluxes at the start of the simulations, due to junk radiation, and at the end, due to the exponential power decay during ringdown, are not very well resolved.Consequently, the obtained E bj curves are generally off by a shift in the E bj plane [118].Even after that shift is applied, it is possible that the NR data exist in the region of the E bj plane that is inaccessible by geodesics.Thus, in order to avoid the introduction of systematics from the determination of the offset in E b and j, here we do not consider them as independent variables for the mapping between SMR and NR configurations.
Finally, we present a combination of variables such that Schwarzschild geodesics and NR simulations lie in the same region of parameter space.These variables are the eccentricity measure e gw , defined in Eqs. ( 5) and ( 6), and the orbitaveraged azimuthal frequency computed from the (2,2)-mode, ω 22 , defined in Eq. ( 7).The lower panel of Fig. 8 displays e gw as a function of ω 22 , for the NR simulations in Table I and for Schwarzschild geodesics.The e gw -ω 22 plane is naturally overlapping for both NR and Schwarzschild geodesics, without the need of any shifts or rescalings.This eccentricity definition by construction spans the range from 0 (circular) to 1 (parabolic orbit).This remains true at any level of the SMR approximation.Additionally, we note that given a (2,2)-mode waveform e gw is uniquely determined, and thus, it is a gauge invariant observable.
However, we note that the leading order contribution to e gw cannot be computed analytically in the SMR formalism, but requires solving the first order field equations numerically.Similarly, the next-to-leading order in mass-ratio correction to e gw requires the second order metric perturbation, which has not yet been calculated for eccentric orbits.Consequently, calculating the next-to-leading order contribution to the expansion of any observable at fixed e gw and ω 22 requires obtaining the second order metric perturbation.(The only exception to this are quantities at fixed e gw = 0 or e gw = 1, since the higher order corrections to these values are zero by construction.) In the limit e gw → 0, fixing e gw and ω 22 reduces to the usual comparisons done for quasi-circular inspirals.Hence, we consider these two variables, e gw and ω 22 , as our independent variables for the comparison of NR and SMR inspirals.

VI. RESULTS
In this section we compare the energy and angular momentum fluxes, as well as the periastron advance, obtained from NR and SMR adiabatic evolutions, and provide constraints on the magnitude of the next order term in the SMR expansion for the mass ratios considered here.We consider the orbitaveraged energy and angular momentum fluxes from the NR simulations in Table I, computed using Eqs.( 13)-( 15), as well as the periastron advance K, computed using Eq. ( 12).
The orbit-averaged fluxes extracted from the NR simulations are illustrated in the top two rows of Fig. 9.In order to reduce the dynamical range of the fluxes, we rescale them with the Newtonian (0PN) quasi-circular values for these quantities [16] ĖQC,0PN gw

JQC,0PN
z,gw Here, Ω orb denotes the orbital frequency for which we substitute ω 22 /2, see Sec.III for details.The rescaling by Eqs.(37) produces a smooth dependence of the fluxes in parameter space, with practically no curves crossing each other.This is because most of the mass ratio dependence is already accounted for by the rescaling factors.In the right-hand panels, the data are plotted as function of e ω 22 only.This projection highlights how well the normalization accounts for the νand ω 22 -dependence, with only the eccentricity-dependence remaining.The eccentricity dependence qualitatively resembles the expected analytical behaviour for the energy flux for eccentric binaries with corrections of the form ∼ 1 (1−e 2 ) x (1 + ae 2 + be 4 + ...), where x = 7/2 or 2, for the energy and angular momentum fluxes respectively, while a and b are coefficients which can be found in [120,121].We do not introduce eccentric corrections to the rescaling factors as our eccentricity definition, e gw , only reduces to the temporal eccentricity, e t , at Newtonian order, while higher PN order corrections may be important to reproduce the eccentricity dependence of the NR fluxes, especially at the end of the inspiral regime.Hence, we leave the exploration of the eccentricity dependence of the NR fluxes for future work.
The periastron advance is not rescaled, since the Newtonian value is simply 1 by Kepler's first law.As a consequence, a larger dependence of this quantity on mass ratio is observed when projected into the K NR − e gw plane.While the values corresponding to a fixed reference of ω ref 22 = 0.042 (red circles) as a function of eccentricity show a similar behavior as the fluxes.Overall, the inspection of the NR curves in Fig. 9 indicates that most of the mass ratio dependence may be already captured by the leading order mass ratio contribution.
Moving to the comparison of NR against SMR results, the NR and SMR fluxes rescaled by the leading order symmetric mass ratio squared as well as the periastron advance, are shown in Figs. 10, 11 and 12 for three different reference frequencies representative of the full inspiral, Mω ref 22 = 0.034, 0.042, 0.063.The SMR fluxes are determined numerically from geodesic snapshots at the quantities selected in Sec.V, (ν, e gw , ω ref 22 ), as explained in Sec.IV.The SMR values for the periastron advance correspond to the analytic geodesic result for the periastron advance, which can be readily obtained from expressions available in the Black Hole Perturbation Toolkit [122], where K is the complete elliptical integral of the first kind, and p and e g have been evaluated at the corresponding values of ω ref 22 and e gw .
In the case of the fluxes (Figs. 10 and 11), both NR and SMR show qualitatively good agreement, which is maintained with increasing eccentricity.This indicates that the effect of eccentricity is well captured by the SMR calculations.Ad-ditionally, the dependence on mass ratio when rescaling by the leading order symmetric mass ratio contribution is small.The qualitative agreement between the NR and SMR degrades with increasing reference frequency, as expected because higher order mass ratio corrections are larger in the strong field.The periastron advance, shown in Fig. 12, has  stronger dependence on mass ratio than the fluxes, especially at high frequencies.As in the case of the fluxes, with increasing eccentricity the agreement of the periastron advance between NR and SMR does not substantially degrade, indicating that eccentric effects are accurately described within SMR theory using adiabatic evolutions.Overall, for both fluxes and periastron advance the SMR curves overestimate the NR results for all frequencies, mass ratios and eccentricities.
Before proceeding to a more quantitative comparison of the difference between NR and SMR, we assess the accuracy of the NR values shown in Figs.10-12 by comparing NR data obtained with different numerical settings.The data in Figs.10-12 was obtained from the highest numerical reso- lution (Lev3) with applied center-of-mass (CoM) correction 6 , extrapolation order 4, and all the spin-weighted spherical harmonic modes up to l ≤ 8.For the particular reference frequency of ω ref 22 = 0.042, we show in Fig. 13 the absolute difference between the quantities computed from this reference waveform against the same quantities calculated from a waveform, where one of the previous conditions is modified at a time.Precisely, the differences are computed against a waveform without CoM correction; using extrapolation order 3; a lower resolution (Lev2), and also in the case in which only l ≤ 4 modes are included.The largest differences for the three quantities typically occur when comparing against the lower resolution (Lev2).Furthermore, the individual errors are summed in quadrature for an overall error estimate for the subsequent analysis.
We now perform a more quantitative comparison of SMR and NR results for the particular reference frequency, ω ref 22 = 0.042.The difference between the SMR and NR fluxes rescaled by the leading order power of symmetric mass ratio as a function of eccentricity are shown in the top and mid left panels of Fig. 14, while in the bottom panels the differences for periastron advance are displayed.Each data-point carries the error bar determined through the analysis in Fig. 13.We see that at 0PA order there is already good agreement between NR and SMR results, with relative differences typically of the order 10%, with the largest discrepancies occurring at equal masses, as expected from a small mass ratio expansion.
tions.The right panels of Fig. 14 show that this scaling collapses the three quantities into 1-dimensional curves.These 1-dimensional curves represent the next-to-leading order 1PA contribution to the respective quantity, as a function of eccen-tricity.The small residual spread in mass ratio in these curves represents unknown yet higher order terms.The fact that the right panels of Fig. 14 collapse to quasi 1-dimensional curves indicates that such ≥ 2PA contributions are small compared to the 1PA contribution.
Additionally, we have added to the right top and mid panels of Fig. 14 the second order self-force results for the quasicircular fluxes from [65].The agreement is good, but a small shift is noticeable between the quasi-circular results with respect to our eccentric results.This feature may be a consequence due to the fact that the results from [65] are based on a two-timescale expansion computed from the self-force dynamics, while our fluxes are averaged over a radial period of an evolving SMR waveform.However, a more detailed study is required to determine the source of this small discrepancy, which is within the error bars.
In the case of the periastron advance, when rescaling by an additional power of symmetric mass ratio in the bottom right panel of Fig. 14, we include also the quasi-circular 1PA SMR result [66,67] and the two quasi-circular non-spinning NR simulations (q = 1, 1/8) from [90].The eccentric results have comparatively large error bars at small eccentricities because the amplitude of the oscillations in ω 22 (from which all quantities are derived) becomes small and more difficult to resolve.A similar shift as in the case of the fluxes is present in the bottom right panel of Fig. 14 between the quasi-circular results from [90] and the low-eccentricity data-points of our new analysis.We leave for future work the precise determination of such small differences between our results for the fluxes and the periastron advance, and the existing quasi-circular results from the literature.
Finally, we remark that the dependence of the fluxes and periastron advance with eccentricity resembles a functional form as expected from the PN results [120,121,127], where for instance the eccentric corrections to the fluxes are of the form ∼ 1 (1−e 2 ) x (1 + ae 2 + ...), where x and a are coefficients to be determined.This suggests that fitting such results as a function of eccentricity and mass ratio could provide some phenomenological expressions for the unknown 1PA SMR terms as a function of eccentricity, mass ratio and frequency, as a similar eccentricity dependence is observed for other frequencies in the inspiral.We leave such a task for future work, as well as the production of new eccentric NR simulations at smaller mass ratios, which may help assess the contributions of the unknown higher order terms in the SMR perturbation theory for eccentric non-spinning binaries.

VII. CONCLUSIONS
We have presented a new set of BBH NR simulations produced with the SpEC code with the objective of exploring the accuracy of the small mass-ratio expansion for eccentric nonspinning binary black holes.In particular, our study aims to extend recent work [43] on assessing the accuracy of the SMR theory for non-spinning quasi-circular BBH to non-spinning eccentric BBHs.
The simulations produced in this work cover mass ratios, The energy and angular momentum fluxes (top and mid panels) have been rescaled by the leading order symmetric mass ratio power, ν −2 .Right panels: Same quantity as in the corresponding left plot rescaled by an additional power in symmetric mass ratio.In all panels each point is color-coded by symmetric mass ratio and carries the error bar computed in Fig. 13.The red dots in right top and mid panels corresponds to the quasi-circular second order self-force results from [65].In the right bottom plot the gray dot refers to the SMR prediction for quasi-circular binaries from [90], and the dots circled by magenta disks correspond to the periastron advance values for the q = 1, 1/8 quasi-circular NR simulations computed in [90].
q ∈ [0.1, 1], initial eccentricities, e 0 gw ∈ [0.01, 0.7], and initial mean anomalies close to apastron, l 0 ∼ π.Each simulation is performed at three different resolutions, and most of them have 20 orbits, which makes our dataset the one with the longest eccentric BBH simulations to date.These simulations are compared to waveforms produced using the gravitational self-force formalism.Using an existing frequency domain Teukolsky code [110][111][112], we have generated eccentric inspirals in a Schwarzschild background that are accurate to leading order in the SMR expansion.
As a first step towards comparing the NR and SMR results, we adapted the orbit-average method from [58] to extract the radial and azimuthal frequencies, the energy and angular momentum fluxes, and measure the eccentricity from waveforms.We have validated this procedure to extract orbit-averaged frequencies by using the 0PA inspirals, where the geodesic azimuthal and radial frequencies are provided as an outcome of performing such evolutions.We find that the procedure of ex-tracting the frequencies, and eccentricities produces relative differences of 10 −5 in the early inspiral, while the discrepancies increase up to ∼ 10 −2 close to merger due to a combination of the boundary effects and the rapid increase of the frequencies, which is a clear limitation of the procedure.Thus, we restrict this study to the inspiral part of the waveform, and leave for future work an improvement of the extraction procedure to describe more faithfully the transition from inspiral to plunge of the signal.
We investigated eccentricity e Ω orb defined from the orbital frequency and eccentricity e ω 22 defined from the gravitational wave (2,2) mode, and found them to systematically differ.Using PN theory we have derived relations between different definitions of eccentricity.The instantaneous orbital and (2,2)-mode frequency are not related by the simple factor 2 for eccentric binaries, as is the case of the orbit-averaged frequencies, and thus, we have provided PN-accurate expressions relating both, which produce relative differences of ∼ 10 −2 when tested on NR simulations.Furthermore, we have provided PN-accurate expressions relating e Ω orb , e ω 22 and the temporal eccentricity, e t .We show that in the Newtonian limit e ω 22 ∼ 3e t /4, so that e ω 22 does not have the correct Newtonian limit.In Eq. ( 6), we propose a new eccentricity definition e gw based on the (2,2)-mode frequency, which reduces to e t in the Newtonian limit.
Comparisons between NR and SMR require a map which associates a SMR inspiral with the instantaneous state of an eccentric NR inspiral.We investigated several proposals in the literature for variables that identify the same inspiral in the NR simulations and SMR evolutions.We find that some choices used in the literature lead to the NR simulations lying outside the range spanned by the geodesic results, hampering comparisons.We propose to use as variables the orbit-averaged azimuthal frequency, ω 22 , and eccentricity e gw , measured both from the instantaneous frequency of the (2,2)-mode, which do not suffer from this limitation.
Moving to the comparison between NR and SMR results, we have focused on the energy and angular momentum fluxes, as well as the periastron advance.Overall, we find good agreement between the NR and SMR values, with relative differences typically 10%, and no particular degradation with increasing eccentricity.
We assess the contributions coming from the unknown higher order term in the PA expansion (1PA) by considering the difference between the NR and SMR fluxes and periastron advance.After rescaling by the symmetric mass-ratio cubed, we find that the differences collapse to one dimensional curves as a function of eccentricity with very small spread in mass ratio, see Fig. 14.This behavior indicates that the next order term in the SMR expansion (2PA) has a very small contribution compared to the 1PA term.Furthermore, we compare these differences for the fluxes and the periastron against available results in the literature for quasi-circular binaries from [65,90], and find that the results are consistent with our findings, except for small shifts which are within the error bars.We leave for future work the precise determination of the origin of this small feature.
The eccentricity dependence of the fluxes and periastron advance rescaled by symmetric mass ratio also suggests a functional form similar to the one predicted by the known PN results [120,121,127].An interesting extension of the work presented here would be the modelling of these differences between the adiabatic SMR inspirals and the NR simulations, by fitting them as a function of mass ratio, eccentricity and orbit-averaged frequency {q, e gw , ω 22 }, and provide some phenomenological expressions which can be used to compute the unknown 1PA term for the fluxes and periastron advance.Another possible future direction is to focus on comparing the phasing between NR and SMR, and extend previous studies for quasi-circular binaries [43] to the eccentric case.Future work will also include extending our set of simulations to higher mass ratios, and to gradually incorporate spins.Other applications of the simulations will include the calculation of the redshift factor [128], extending current studies on quasi-circular binaries [129] to the eccentric case.Finally, these simulations will also be of paramount relevance to assess the accuracy of the currently existing inspiral-mergerringdown eccentric waveform models [59,[130][131][132][133].
1.5.We will avoid this problem by starting NR simulation near apastron. 7Second, there are previous results on how the tangential momentum in an eccentric binary varies with the eccentricity.Specifically, we will utilize the correction factor of the tangential momentum [51] λ 0 t (r, e, ν, sign) where r is the orbital separation and sign = ±1 is the sign of the correction [51].While this correction has been derived in the low eccentricity limit, it has been shown [63] to be useful to determine the initial parameters in eccentric moving puncture simulations.We average the correction with both signs to arrive at λ0 t (r, e, ν) = as in Eq. ( 2.3) of [63].The third consideration concerns the choice of coordinates: our SpEC simulations start from superposed harmonic Kerr (SHK) data [80], whereas Eq. (A1) was derived in Arnowitt-Desner-Misner transverse-traceless (ADMTT) coordinates.Therefore, we will also employ a coordinate transformation from ADMTT coordinates to harmonic coordinates.Overall, we proceed as follows: 1. Choose mass-ratio q ≤ 1, and a desired eccentricity e 0 .Set spins χ χ χ i = 0, masses m 1 = 1/(1 + q), m 2 = q/(1 + q) (so that M 0 = m 1 + m 2 = 1), and ν = m 1 m 2 /M 2 0 = q/(1 + q) 2 .2. Choose a tentative initial separation as the apastron distance of a Newtonian binary with periastron distance of r p = 9M 0 , i.e.D0 = (1 + e 0 )(1 − e 0 ) −1 r p .If a PN evolution with the same parameters indicates that the inspiral may be too short, increase D0 .3. Compute 3.5PN quasi-circular estimates for the tangential and radial momenta, p 0 t , p 0 r in ADMTT coordinates using Eqs.(A2) and (2.16) in [63].4. Calculate the correction factor λt using Eq.(A2). 5. Construct the ADMTT position and momentum vectors in Cartesian coordinates, x x x ADMTT = (0, D0 , 0), (A3) Here, we placed the black holes on the y-axis.Note that λt ≤ 1, so that the tangential momentum is reduced, consistent with our goal to start at apastron.6. Apply the transformation from ADM to harmonic coordinates [134] to obtain the position and velocity vector in harmonic coordinates, A6) 7 Very recently a new radial map has been developed and implemented in SpEC, which avoids these restrictions.where Y Y Y and V V V are operators mapping the ADM coordinates to harmonic coordinates expanded up to 3PN order [134].(Note that the expressions in [134] are restricted to non-spinning binaries.)7. Read of the initial data parameters from the position and velocity vectors in harmonic coordinates, ) where Euclidean vector operations are used.
Figure 15 compares the target eccentricity e 0 with the actual eccentricity e 0 gw achieved near the start of each simulation.There is an offset between these eccentricities.We note that specially for high eccentricities the use of the correction factor is not accurate due to the fact that it is an expression derived in the low eccentricity limit.Furthermore, we attribute the larger differences between our target and measured initial eccentricities as compared to other studies like [63] due to the assumptions on the identification we made between harmonic and superposed harmonic coordinates, and inaccuracies in the PN expressions for the eccentric corrections being amplified due to the transformation from the ADM to the harmonic gauge.
The calculation of the initial parameters presented in this section is useful for placing points in the eccentric parameter space with a limited accuracy.In the future we plan to adopt an iterative procedure to specify the desired initial eccentricity and mean anomaly as done in [60,135], to accurately and efficiently populate the eccentric parameter space.

Appendix B: Numerical relativity waveform quality
In this appendix we assess the accuracy of the NR waveforms listed in Table I.For each simulation SpEC employs multiple subdomains.The shape, size and number of subdomains is dynamically varied during the simulations according to the spectral adaptive mesh refinement (AMR) procedure [81,136].The accuracy of the simulations is controlled by a tolerance parameter which determines when AMR should add or remove grid points within a given subdomain, and when a subdomain should be split into two, or when two neighboring subdomains should be combined into one.As a consequence, it is difficult to obtain strict convergence as a function of the AMR tolerance parameter.Convergence may fail, for instance, due to two identical simulation having different AMR tolerances in a particular subdomain modifying the number of grid points in it, or different subdomain boundaries in a particular time.Notwithstanding these issues, most simulations in the SXS catalog show convergence with the AMR tolerance [52].
In this work we have run each simulation at three different AMR tolerances, henceforth called different resolutions for brevity.This appendix extends the error analysis of our main results with calculations of the mismatch between waveforms obtained at the highest and second highest resolutions.Following Ossokine et al 2020 [137], we compute the SNR-weighted mismatch between waveforms computed from the highest and the next-to-highest resolutions.The mismatches are computed for binary masses 20M ≤ M ≤ 200M , and using as a Power Spectral Density (PSD), the Advanced LIGO's zero-detuned high-power design sensitivity curve [138].When both waveforms are in band, we use f min = 10Hz and f max = 2048Hz, as the lower and upper bounds of the integral.For waveforms where this is not the case, we set f min = 1.05 f start , where f start is the starting frequency of the NR waveform.To represent dependence on M, we compute the mean and the maximum over M. The results of the mismatch calculation are shown in Fig. 16.The vertical dashed lines denote the median values of each distribution.We note that the median value of the maximum mismatch is below 10 −3 , while for the mean mismatch is ∼ 10 −4 .Three simulations (SXS:BBH:2517, SXS:BBH:2525, SXS:BBH:2564) have maximum mismatches above 1%, MSNR max = 0.011, 0.065, 0.041, respectively.The highest mismatch occurs for SXS:BBH:2525 which is both the shortest evolution in our dataset (making it more prone to systematics due to the ringdown transition in SpEC [70,72]), and which was also the first simulation produced in our dataset, so it does not take into account some improvements in SpEC, which have been introduced during this project (see Sec. II for details).Overall, the mismatches are comparable to the ones obtained in the SXS catalog for quasi-circular binaries [52] (see Fig. 9 there, but note that [52] uses a flat PSD).This indicates that SpEC is capable to perform simulation of eccentric BBH with a numerical error comparable to the quasi-circular case.In this Appendix we use PN theory to investigate the relation among e ω 22 , e Ω orb and the post-Newtonian e t [139].In the following, we set the total mass, M = 1, to ease the notation.

Relation e ω 22 − e Ω orb
Section III showed that the differences between e ω 22 and e Ω orb can be explained within PN theory.This appendix derives the relations used there at 1PN using harmonic coordinates.As a first step, we calculate e ω 22 from h 22 at the 1PN order [93], where γ = 1 c 2 is the PN order bookkeeping parameter, i is the imaginary unit, r is the radial separation, φ is the orbital phase and the overdot represents a time derivative.
Taking the complex argument of Eq. (C1) and expanding to 1PN order yields The where ) This result is used in the main text in Eq. ( 17).At the turning points apastron and periastron, ṙ = 0 and φ = 0, and Eq.(C6) simplifies to This result is used in the main text in Eq. ( 18).The expressions derived above can be useful to estimate ω 22 , or the eccentricity e ω 22 , for NR simulations, where the trajectories are output in Cartesian or polar coordinates 8 , and are typically cleaner quantities than the frequencies of the extracted waveform modes, especially for finite difference codes [46].Another application of Eq. (C10) is for eccentricity reduction/control purposes, where short evolutions are done to iteratively converge to the desired value of eccentricity.In these methods [46,51,83,84] one typically chooses a trajectory-based definition of eccentricity instead of a waveform-based one due to the extra computational cost, which involves the evolution of the gravitational radiation reaching the extraction radii.However, with the expressions provided in Eq. (C10), one can obtain an approximation of e ω 22 from the coordinates.

Relation e Ω orb − e t
Most eccentric waveform models for compact binaries use PN theory to describe the inspiral regime and/or their initial parameters [94-102, 131, 133, 140].A commonly used description of eccentric orbits is the quasi-Keplerian parameterization [94] where three different eccentricity parameters e t , e r and e φ describe the orbit [139].These three eccentricities are not independent from each other, and they are all related at a given PN order.Eccentric PN waveform models typically use the temporal eccentricity, e t , as the eccentricity parameter.In the following, PN-accurate expressions between the eccentricity defined from the orbital frequency, e Ω orb , and e t , are computed.
In order to perform this calculation we use the 3PN expression for Ω orb , which can be found in Appendix A of [56].The calculation of e Ω orb requires the values of the orbital frequency at periastron and apastron, which correspond to values of the eccentric anomaly of u = 0 and u = π, respectively.Thus, at the turning points, the orbital frequency can be expressed as where γ is the PN bookkeeping parameter.Using the abbreviation ε ≡ 1 − e 2 t , the contributions at different PN order can be written as ) where the upper sign corresponds to apastron and the lower sign corresponds to periastron.To derive Eq. (C20) we assumed that the value of x is the same at apastron and periastron as it corresponds to an orbit-averaged frequency, which is evolved using the radiation reaction equations in an adiabatic evolution.This approximation may not be accurately fulfilled when post-adiabatic effects become more relevant as in the case of the binary close to merger.Substituting Eq. (C20) into Eq.( 4) and PN-expanding the result to 3PN order, one obtains We note that in the derivation of Eq. (C21) only the instantaneous contributions to the orbital frequency up to 3PN order have been used, while tail contributions or spin terms, which would appear beyond 1PN order, have been neglected.We leave for future work including spin effects, as well as contributions from the tail terms.

Relation e ω 22 − e t
This derivation proceeds similarly to appendix C 1, but starting from the quasi-Keplerian parametrization.We start with 1PN expressions for the (2,2)-mode waveform, h 22 , in the quasi-Keplerian parameterization [93], (C33) The time derivatives ẋ, ėt and φ can be found in [56,141], while for the eccentric anomaly, u, we use the Kepler equation at Newtonian order to write9 u = l + ėt sin u where l is the mean anomaly, and an expression for l in the quasi-Keplerian parametrization can be found in [56].We note that the 3PN Kepler equation can be found in [56,141], however, we restrict to low PN order for simplicity of the calculations, and to avoid the introduction of the true anomaly, which substantially complicates the higher order calculations [142].
At 1PN order, one can write the following expression for the frequency of the (2,2)-mode ω

Figure 2 .
Figure 2. Top panel: Time evolution of the frequency of the (2,2)mode (solid blue line) for the simulation SXS:BBH:2558.The values of the (2,2)-mode frequency at periastron and apastron are indicated with orange and black dots, respectively.These are used to compute the orbit-averaged frequency of the (2,2)-mode (solid red curve), and the eccentricity e gw (dashed green curve) through Eq. (6).Bottom panel: Time evolution of the mean anomaly (solid purple line) computed using Eq.(10) for the same simulation as in the top panel.The vertical dashed gray lines in both panels correspond to the times of the periastron passages.

Figure 5 .
Figure 5. Frequency extraction procedure applied to an SMR waveform at 0PA at equal mass.The solid curves arise directly from the SMR inspiral and its dynamics.The filled circles are the result of applying our frequency extraction procedure to the maxima of the instantaneous frequency ω 22 (t).Even at equal mass where the inspiral is fastest, the recovered orbital averaged azimuthal and radial frequencies agree well with the geodesic frequencies of the underlying SMR dynamics.This figure is analogous to Fig.4.

MFigure 6 .Figure 7 .
Figure 6.Absolute relative difference between the frequencies from the waveform (azimuthal ω 0PA /2 and radial Ω r 0PA ), and the frequencies from the geodesic inspiral (azimuthal Ω geo and radial Ω r geo ), as a function of symmetric mass ratio ν at three selected points along the inspiral from Fig.5.The frequencies extracted from the waveforms, ω 0PA and Ω r 0PA , have been computed using the orbit-average procedure of Secs.II B and II C employing the periastron passages.The gray line indicates a linearly increasing ν-dependence.

Figure 8 .
Figure 8. Top panel: Ratio of the orbit-averaged azimuthal and radial frequencies computed from the (2,2)-mode , ω 22 /(2 Ω r 22 ), as a function of the orbit-averaged frequency, ω 22 .The colored curves represent the NR simulations in TableI, whereas the grey shaded area indicates the region covered by Schwarzschild geodesics, bounded by the diagonal black curve representing quasi-circular geodesics.Bottom panel: Eccentricity, e gw , computed using Eqs.(5) and (6), as a function of ω 22 , for the same NR simulations as in the upper panel.While geodesics exist at all eccentricities, we have only generated SMR configurations in the grey shaded area.In both panels each NR simulation has been color-coded according to its symmetric mass ratio ν.

Figure 9 .
Figure 9. Left column: From top to bottom energy flux, angular momentum flux and periastron advance extracted from the NR simulations in Table I as a function of eccentricity, e gw , and orbit-averaged azimuthal frequency, ω 22 .Each curve corresponds to a NR simulation in Table I, and is color-coded by symmetric mass ratio ν.The energy and angular momentum fluxes are rescaled by the quasi-circular Newtonian expressions in Eqs.(37).The red planes indicate the reference frequency Mω ref 22 = 0.042.Right column: Projection of the left plots, in the Z-e ω 22 plane, where Z indicates the quantity in the z-axis (fluxes or periastron advance).The red-white circles indicate the points where each NR simulation passes the reference frequency Mω ref 22 = 0.042.

Figure 10 .Figure 11 .
Figure 10.Orbit-averaged energy flux rescaled by the leading order symmetric mass ratio dependence (ν −2 ) as a function of eccentricity at three different reference frequencies, Mω ref 22 = 0.034, 0.042, 0.063.Each marker corresponds to a NR simulation at the specified reference frequency, and it is color coded by mass ratio.The solid lines are the leading order SMR energy flux.

Figure 12 .
Figure 12.Periastron advance as a function of eccentricity at three different frequencies, Mω ref 22 = 0.034, 0.042, 0.063.Each marker corresponds to a NR simulation at the specified reference frequency, and it is color coded by mass ratio.The solid lines correspond to joining the values of the geodesic periastron precession at the same (ν, e gw , ω ref 22 ) values as the NR configurations.

Figure 13 .
Figure 13.Error estimates for the energy and angular momentum fluxes, and the periastron advance computed from the NR simulations in Table I at a reference frequency of ω ref 22 = 0.042.Taking as a reference data computed from the waveform computed with highest numerical resolution (Lev3), with extrapolation order 4, CoM correction, and all the modes up to l ≤ 8, we compute the absolute difference that arises when each one of these conditions is changed, i.e., comparing to the values computed from the waveform with extrapolation order 3; without CoM correction; against a lower resolution; and against the waveform with incomplete modes only up to l ≤ 4. In the case of periastron advance the impact of higher order modes is not assessed as this quantity is computed from the (2,2)-mode.The orange circles represent the quadrature sum of the individual error contributions.

Figure 14 .
Figure 14.Left panels: Difference between the SMR and NR fluxes and periastron advance as a function of eccentricity at a reference frequency of ω ref 22 = 0.042.The energy and angular momentum fluxes (top and mid panels) have been rescaled by the leading order symmetric mass ratio power, ν −2 .Right panels: Same quantity as in the corresponding left plot rescaled by an additional power in symmetric mass ratio.In all panels each point is color-coded by symmetric mass ratio and carries the error bar computed in Fig.13.The red dots in right top and mid panels corresponds to the quasi-circular second order self-force results from[65].In the right bottom plot the gray dot refers to the SMR prediction for quasi-circular binaries from[90], and the dots circled by magenta disks correspond to the periastron advance values for the q = 1, 1/8 quasi-circular NR simulations computed in[90].

Figure 15 .
Figure 15.Initial eccentricity (dots), e 0gw , as defined in Eq. (6), measured from the NR simulations in TableI, and initial eccentricity (crosses), e 0 , specified in Eq. (A2), as a function of the merger time of the simulations.Each simulation is color-coded according to its inverse mass ratio.

Figure 16 .
Figure16. of the SNR-weighted mismatch between the two highest resolutions for each simulation in TableI.The orange and green distributions correspond to the mean and maximum mismatch over the total mass range considered M =[20, 200]M .The vertical dashed lines correspond to the median values of the distributions.