Eccentric binary black hole surrogate models for the gravitational waveform and remnant properties: comparable mass, nonspinning case

We develop new strategies to build numerical relativity surrogate models for eccentric binary black hole systems, which are expected to play an increasingly important role in current and future gravitational-wave detectors. We introduce a new surrogate waveform model, \texttt{NRSur2dq1Ecc}, using 47 nonspinning, equal-mass waveforms with eccentricities up to $0.2$ when measured at a reference time of $5500M$ before merger. This is the first waveform model that is directly trained on eccentric numerical relativity simulations and does not require that the binary circularizes before merger. The model includes the $(2,2)$, $(3,2)$, and $(4,4)$ spin-weighted spherical harmonic modes. We also build a final black hole model, \texttt{NRSur2dq1EccRemnant}, which models the mass, and spin of the remnant black hole. We show that our waveform model can accurately predict numerical relativity waveforms with mismatches $\approx 10^{-3}$, while the remnant model can recover the final mass and dimensionless spin with errors smaller than $\approx 5 \times 10^{-4}M$ and $\approx 2 \times10^{-3}$ respectively. We demonstrate that the waveform model can also recover subtle effects like mode-mixing in the ringdown signal without any special ad-hoc modeling steps. Finally, we show that despite being trained only on equal-mass binaries, \texttt{NRSur2dq1Ecc} can be reasonably extended up to mass ratio $q\approx3$ with mismatches $\simeq 10^{-2}$ for eccentricities smaller than $\sim 0.05$ as measured at a reference time of $2000M$ before merger. The methods developed here should prove useful in the building of future eccentric surrogate models over larger regions of the parameter space.


I. INTRODUCTION
Detection of gravitational waves (GWs) [1,2] by the LIGO [3] and Virgo [4] detectors has opened a new window in astrophysics to probe binary compact objects -binary black holes (BBHs) being the most abundant source for these detectors. Both detection and extraction of source properties from the GW signal relies on the availability of accurate inspiral-merger-ringdown (IMR) waveform models for BBHs. While numerical relativity (NR) provides the most accurate gravitational waveforms for BBHs, they are computationally expensive, taking weeks to months to generate a single waveform. Data-driven surrogate modeling strategies [5][6][7][8][9][10][11][12][13] have been shown to be capable of producing waveforms that are nearly indistinguishable from NR with evaluation times of less than 0.1 seconds. While NR surrogate waveform models for nonspinning [7], aligned-spin [10], and precessing BBHs [9,13] are well developed, NR surrogate modeling of eccentric systems is completely unexplored. lead to systematic biases if the actual signal corresponds to an eccentric system [43]. Such biases can also lead to eccentric systems being misidentified as a violation of general relativity (GR). Even if all binaries are found to be circular, eccentric models are necessary to place bounds on the eccentricity. Therefore, including eccentricity in our GW models is important, especially as the detectors become more sensitive.
Notably, all of these models rely on the assumption that the binary circularizes by the merger time. While this is approximately true for many expected sources [52,60], this necessarily places a limit on the range of validity of these models. In addition, none of these models are calibrated on eccentric NR simulations, even though their accuracy is tested by comparing against eccentric simulations.
Apart from the waveform prediction, BBH remnant modeling from eccentric sources is also of crucial astrophysical importance [61][62][63][64]. For example, recoils from eccentric mergers can be up to 25% higher than the circular case [63,64], which result in a higher likelihood of ejections from astrophysical hosts like star clusters and galaxies.
It is, therefore, timely to invest in building faithful eccentric BBHs waveform and remnant models that address some of these limitations. In this paper, we develop a detailed framework for constructing a surrogate model with eccentric NR data. We then build a two-dimensional surrogate model, NRSur2dq1Ecc, over parameters that describe eccentricity for equal-mass, nonspinning systems to demonstrate the efficacy of the proposed methods. This is the first eccentric waveform that is directly trained on eccentric NR simulations and does not need to assume that the binary circularizes before merger. The model can produce waveforms that are of comparable accuracy to the NR simulations used to train it. Furthermore, despite being trained only on equal-mass eccentric BBHs, we find that the model can be reasonably evaluated beyond its training range upto mass ratio q ≈ 3 provided the eccentricities are small.
In addition to the waveform model, we build a surrogate model for the remnant mass and spin, NRSur2dq1EccRemnant, which can provide accurate predictions for the final state of eccentric binary mergers. This work paves the way forward for building future eccentric surrogate models: we expect that the methods developed here can be applied straightforwardly to aligned-spin eccentric BBHs, while the precessing case requires significantly more work.
The rest of the paper is organized as follows. Sec. II describes the NR simulations. Sec. III describes data decomposition, parameterization and construction of the surrogate model. In Sec. IV, we test the surrogate model by comparing against NR waveforms. We end with some concluding remarks in Sec. V.

II. NUMERICAL RELATIVITY DATA
NR simulations for this work are performed using the Spectral Einstein Code (SpEC) [65] developed by the Simulating eXterme Spacetimes (SXS) collaboration [66]. We follow the procedure outlined in Ref. [67] to construct initial orbital parameters that result in a desired eccentricity. The constraint equations are solved employing the extended conformal thin sandwich formalism [68,69] with superposed harmonic Kerr free data [70]. The evolution equations are solved employing the generalized harmonic formulation [71,72]. The time steps during the simulations are chosen nonuniformly using an adaptive time-stepper [73]. Further details can be found in Ref. [73] and references within. We perform 47 new eccentric NR simulations that have been assigned the identifiers SXS:BBH:2266 -SXS:BBH:2312, and will be made available through the SXS public catalog [74].
The component BH masses, m 1 and m 2 , and dimensionless spins, χ 1 and χ 2 , are measured on the apparent horizons [73] of the BHs, where index 1 (2) corresponds to the heavier (lighter) BH. The component masses at the relaxation time [73] are used to define the mass ratio q = m 1 /m 2 ≥ 1 and total mass M = m 1 + m 2 . Unless otherwise specified, all masses in this paper are given in units of the total mass. When training the surrogate model, we restrict ourselves to q = 1, χ 1 , χ 2 = 0 in this work.
The waveform is extracted at several extraction spheres at varying finite radii from the origin and then extrapolated to future null infinity [73,75]. These extrapolated waveforms are then corrected to account for the initial drift of the center of mass [76,77]. The spin-weighted spherical harmonic modes at future null infinity, scaled to unit mass and unit distance, are denoted as h m (t) in this paper.
The complex strain h = h + − ih × is given by: where h + (h × ) is the plus (cross) polarization of the waveform, −2 Y m are the spin = −2 weighted spherical harmonics, and ι and ϕ 0 are the polar and azimuthal angles on the sky in the source frame. We model modes with ( , m) = (2, 2), (3,2), (4,4). Because of the symmetries of equal-mass, nonspinning BBHs, all odd-m modes are identically zero, and the m < 0 modes can be obtained from the m > 0 modes. Therefore, we model all non-zero ≤ 3 and (4, ±4) modes, except the m = 0 modes. We exclude m = 0 memory modes because (non-oscillatory) Christodoulou memory is not accumulated sufficiently in our NR simulations [78]; this defect was recently addressed in both Cauchy characteristic extraction (CCE) [79][80][81] and extrapolation [82] approaches. The (4, 2) mode, on the other hand, was found to have significant numerical error in the extrapolation procedure [73,75]. We expect this issue to be resolved with CCE as well. Therefore, in future models, we should be able to include the m = 0 modes as well as modes like the (4,2) mode.
The remnant mass m f and spin χ f are determined from the common apparent horizon long after ringdown [73]. For nonprecessing systems like the ones considered here, the final spin is directed along the direction of the orbital angular momentum. Unlike previous surrogate models [13,83,84], we do not model the recoil kick in this work, as the symmetries of equal-mass, nonspinning BBHs restrict the kick to be zero.

III. SURROGATE METHODOLOGY FOR ECCENTRIC WAVEFORMS
In this section, we describe our new framework to build NR surrogate models for eccentric BBHs. We begin by applying the following post processing steps that simplify the modeling procedure.

A. Processing the training data
In order to construct parametric fits (cf. Sec. III D) for the surrogate model, it is necessary to align all the waveforms such that their peaks occur at the same time. We define the peak of each waveform, τ peak , to be the time when the quadrature sum, reaches its maximum. Here the summation is taken over all the modes being modeled. We then choose a new time coordinate, such that A tot (t) for each waveform peaks at t = 0. Next, we use cubic splines to interpolate the real and imaginary parts of the waveform modes onto a common time grid of [−5500M , 75M ] with a uniform time spacing of dt = 0.1M . The initial time of −5500M is chosen so that we can safely eliminate spurious initial transients in the waveform, also known as junk radiation [73], for each waveform in our dataset.
Once all the waveforms are interpolated onto a common time grid, we perform a frame rotation of the waveform modes about the z-axis such that the orbital phase is zero at t = −5500M . The orbital phase is obtained from the (2, 2) mode [cf. Eq. (14)]. Because of the symmetry of the equal-mass, equal-spin systems considered here, the oddm modes are identically zero and so we need not worry about remaining φ orb → φ orb + π rotational freedom as was necessary in Refs. [7][8][9][10]13]. This preprocessing of time and phase ensures that the waveform varies smoothly across the parameter space, which in turn makes modeling easier.

B. Measuring eccentricity and mean anomaly
Departure of NR orbits from circularity is measured by a time-dependent eccentricity and mean anomaly. Eccentricity takes values between [0, 1] where the boundary values correspond to a quasicircular binary and an unbound orbit [85], respectively. Mean anomaly, on the other hand, is bounded by [0, 2π). While it may seem most natural to estimate orbital parameters from the BH trajectories, this task is complicated by the fact that any such measurement will be impacted by the gauge conditions chosen by the NR simulation. We instead choose to estimate eccentricity and anomaly parameters directly from the waveform data at future null infinity.

Measuring eccentricity
Various methods to extract the eccentricity from NR simulations have been proposed in the literature [60,[86][87][88]. As the eccentricity evolves during the binary's orbit [21], these methods use dynamical quantities such as some combination of the (2, 2) mode's amplitude, phase, or frequency. All of these methods reduce to the eccentricity parameter in the Newtonian limit. The estimated value of the eccentricity may differ slightly depending on the method used and the noise in the numerical data. However, as long as they provide a consistent measurement of eccentricity that decays monotonically with time, one can use any of the eccentricity estimators for constructing a surrogate waveform model. For this work, we use the following definition of eccentricity based on orbital frequency [89]: where ω a and ω p are the orbital frequencies at apocenter (i.e. point of furthest approach) and pericenter (i.e. point of closest approach), respectively. Unlike several other eccentricity estimators proposed in literature [60,[86][87][88], the one defined in Eq. (4) is normalized and reduces to the eccentricity parameter in the Newtonian limit at both low and high eccentricities [43].
ω a (t) We first compute the orbital frequency, where φ orb is the orbital phase inferred from the (2,2) mode (cf. Eq. (14)), and the derivative is approximated using second-order finite differences. We then find the times where ω orb passes through a local maxima (minima) and associate those to pericenter (apocenter) passages, to obtain ω p (ω a ). We find that using the local maxima/minima of the amplitude of the (2, 2) mode to identify the pericenter/apocenter times leads to a consistent value for the eccentricity. We then interpolate ω p and ω a onto the full time grid using cubic splines. This gives us ω p (t) and ω a (t), which are used in Eq. (4). Figure 1 shows an example of the measured eccentricity for the NR simulation SXS:BBH:2304. We see that our method provides a smooth, monotonically decreasing e(t). The estimate become unreliable near merger where finding local maxima/minima in ω orb becomes problematic as the orbit transitions from inspiral to plunge. The estimate also becomes problematic whenever the eccentricity is extremely small, thereby preventing the appearance of an identifiable local maxima/minima. This does not affect our modeling, however, as we only require an eccentricity value at a reference time while the binary is still in the inspiral phase. We select a reference time of t ref = −5500M and parameterize our waveform model by While estimating e ref , we include the data segment slightly before t ref as this allows us to interpolate, rather than extrapolate, when constructing e(t) in Eq. (4).

Measuring mean anomaly
In the Newtonian context, the mean anomaly l of an eccentric orbit is defined as where t 0 is a time corresponding to the previous pericenter passage and P is the radial period, which is defined to be the time between two successive pericenter passages. In the Newtonian case P is a constant, but in GR it changes as the binary inspirals. However, one can continue to use Eq. (7) as a meaningful measurement of the radial oscillation's phase for the purpose of constructing a waveform model [51].
For each NR waveform, we compute the times for all pericenter passages using the same procedure as in Sec. III B 1. We divide the time array into different orbital windows defined as [t peri is the time for i th pericenter passage. The orbital period in each window is given by P i = t peri i+1 − t peri i , and the mean anomaly by Note that each l i (t) grows linearly with time over [0, 2π) for the window [t peri i , t peri i+1 ). To obtain the full l(t), we simply join each l i (t) for consecutive orbits. Finally, the value for mean anomaly parameterizing our waveform model is then simply the evaluation of the mean anomaly (9) Figure 2 shows an example application of our method to estimate the mean anomaly of the NR simulation SXS:BBH:2304.

Targeted parameter space
In Fig. 3, we show the measured values for eccentricity and mean anomaly at t ref for all 47 NR waveforms, which leads to the following 2d parameter space for our model: 3 shows a large gap in the parameter space, which reflects an inherent limitation in our current approach to achieve target eccentricity parameters from the initial data. The method we use to construct initial orbital parameters [67] seeks to achieve target values of (e ref , l ref ) at a time 500M after the start of the simulation. The initial orbital frequency is chosen such that time to merger is 6000M , as predicted by a leading-order PN calculation. Unfortunately, this is only approximate, leading to different merger times for different simulations. Consequently, when we estimate the eccentricity parameters at t ref = −5500M , this is no longer a fixed time from the start of the simulation. The eccentricity parameters evolve differently for different simulations during this time, leading to the clustering in Fig. 3. In the future, we plan to resolve this using a higher order PN expression, or an eccentric waveform model [52][53][54] to predict the time to merger.

C. Waveform data decomposition
Building a surrogate model becomes more challenging for oscillatory and complicated waveform data. One solution is to transform or decompose the waveform data into several simpler "waveform data pieces" that also vary smoothly over the parameter space. These simpler data pieces can then be modeled more easily and recombined to get back the original waveform. Successful decomposition strategies have been developed for quasi-circular NR surrogates [7][8][9][10]13]. In order to develop similar The colors indicate the maximum mismatch, which systematically increases near the high eccentricity boundary where few training data points are available.
strategies for eccentric waveform data, we have pursued a variety of options. We now summarize the most successful decomposition technique we have tried, while relegating some alternatives to Appendix A.

Decomposing the quadrupolar mode h22
The complex (2, 2) waveform mode, can be decomposed into an amplitude, A 22 , and phase, φ 22 . For non-precessing systems in quasicircular orbit, A 22 and φ 22 are slowly varying functions of time, and have therefore been used as waveform data pieces for many modeling efforts. For eccentric waveforms, however, both amplitude and phase show highly oscillatory modulations on the orbital time scale (cf. Figs. 1 and 2 for the frequency, which is a time-derivative of the phase). This demands further decomposition of the waveforms into even simpler data pieces. One natural solution could have been to build interpolated functions of the local maxima and minima of A 22 and φ 22 . The secular trend of these functions can then be subtracted out from the original amplitude and phase. The resulting residual amplitude and phase data may be easier to model. Unfortunately, as mentioned in Sec. III B 1, finding the local maxima/minima becomes problematic near the merger. We instead follow a simpler approach whereby the amplitude and phase of a quasicircular q = 1, nonspinning NR waveform (SXS:BBH:1155) is used as a proxy for the secular trend of the amplitude and phase. We then compute the residual amplitude and phase, where  22 . Note that noneccentric waveform data is plentiful [73] and accurate surrogate models have been built for noneccentric NR waveforms [10,13]. So extending the residual amplitude and phase computation to spinning, unequal-mass systems is straightforward. For instance the surrogate model of Ref. [10] can be used to generate A 0 22 and φ 0 22 for generic aligned-spin systems.

Decomposing the higher order modes
In this paper, we model the quadrupolar mode and the higher-order modes differently. For h 22 , we model data pieces closely associated with the amplitude and phase as described above. On the other hand, for higher order modes, we first transform the waveform into a co-orbital frame in which the waveform is described by a much simpler and slowly varying function. This is done by applying a time-dependent rotation given by the instantaneous orbital phase: where φ 22 is the phase of the (2, 2) mode (cf. Eq. (10)), φ orb is the orbital phase, and h C m represents the complex modes in the co-orbital frame.
We use the real and imaginary parts of h C m as our waveform data pieces for the nonquadrupole modes. As shown in Fig. 5, the h C m data have less structure, making them easier to model. We find that using quasicircular h C m to subtract off the secular trend does not provide any modeling advantage. We, therefore, model the real and imaginary parts of h C m without any further data decomposition.

Summary of waveform data pieces
To summarize, the full set of waveform data pieces we model is as follows: ∆A 22 , ∆φ 22 for the (2,2) mode, and real and imaginary parts of h C m for the (3,2) and (4,4) modes.

D. Building the waveform model
We decompose the inertial frame waveform data into many waveform data pieces as summarized in Sec. III C 3. For each of these data pieces, we now build a surrogate model using reduced basis, empirical interpolation, and parametric fits across the parameter space. The detailed procedure is outlined in Refs. [5,9], which we only briefly describe here.
For each waveform data piece, we employ a greedy algorithm to construct a reduced basis [90] such that the projection errors (cf. Eq. (5) of Ref. [9]) for the entire data set onto this basis are below a given tolerance. We use a basis tolerance of 10 −2 radians for ∆φ 22 , 1.5 × 10 −3 for ∆A 22 and 2 × 10 −5 for the real part of h C 32 . For all other data pieces, basis tolerance is set to 5 × 10 −5 .
These choices are made so that we include sufficient number of basis functions for each data piece [9 for ∆A 22 , 12 for ∆φ 22 , 7 (5) for the real (imaginary) part of h C 32 and 10 (6) for the real (imaginary) part of h C 44 ] to capture the underlying physical features in the simulations while avoiding over fitting. We perform additional visual inspection of the basis functions to ensure that they are not noisy in which case modeling accuracy can become comprised (cf. Appendix B of Ref. [9]). The next step is to construct an empirical interpolant in time using a greedy algorithm which picks the most representative time nodes [5,[91][92][93]. The number of the time nodes for each data piece is equal to the number of basis functions used. The final surrogate-building step is to construct parametric fits for each data piece at each of the empirical time nodes across the two-dimensional parameter space {e ref , l ref }. We do this using the Gaussian process regression (GPR) fitting method as described in Refs. [83,94].

E. Evaluating the waveform surrogate
To evaluate the NRSur2dq1Ecc surrogate model, we provide the eccentricity e ref and mean anomaly l ref as inputs. We then evaluate the parametric fits for each waveform data pieces at each time node. Next, the empirical interpolant is used to reconstruct the full waveform data pieces (cf. Sec. III C 3).
We compute the amplitude and phase of the (2, 2) mode,  The final mass and spin fits are also constructed using the GPR fitting method as described in Refs. [83,94].

IV. RESULTS
In this section we demonstrate the accuracy of NRSur2dq1Ecc and NRSur2dq1EccRemnant by comparing against the eccentric NR simulations described in Sec. II. We do this by performing a leave-one-out cross-validation study. In this study, we hold out one NR waveform from the training set and build a trial surrogate from the remaining 46 eccentric NR waveforms. We then evaluate the trial surrogate at the parameter value corresponding to the held out data, and compare its prediction with the highest-resolution NR waveform. We refer to the errors obtained by comparing against the left-out NR waveforms as cross-validation errors. These represent conservative error estimates for the surrogate models against NR. Since we have 47 eccentric NR waveforms, we build 47 trial surrogates for each error study. We compare these errors to the NR resolution error, estimated by comparing the two highest available NR simulations.  (4,4) in the co-orbital frame are shown in the middle and lower panels respectively. The waveform is aligned such that the peak of the amplitude occurs at t = 0 and the orbital phase is zero at t ref = −5500M .

Time domain error without time/phase optimization
In order to quantify the accuracy of NRSur2dq1Ecc, we first compute the normalized L 2 -norm between the NR data and surrogate approximation where h(t) andh(t) correspond to the complex strain for NR and NRSur2dq1Ecc waveforms, respectively. Here, t 1 and t 2 denote the start and end of the waveform data. As the NR waveforms are already aligned in time and phase, the surrogate reproduces this alignment. Therefore, we compute the time-domain error E without any further time/phase shifts. In Fig. 6, we report both the full waveform and individual mode errors for NRSur2dq1Ecc. For comparison, we also show the NR resolution errors. When computing the full waveform error we use all modes included in the surrogate model ( , m) = (2, 2), (3, 2), (4, 4) in  Eq. (17). To compute errors for individual modes, we restrict the sum in Eq. (17) to only the mode of interest. The NRSur2dq1Ecc errors are comparable to the NR errors in Fig. 6. However, we find that the surrogate errors have an extended tail around two orders of magnitude larger than the largest NR mismatch. While this could imply overfitting, we find that highest mismatches correspond to the parameter space adjacent to the higher eccentricity e ref boundary where only few (to none) training waveforms are used. As will be discussed in Sec. IV A 2, the sparsely  We further note in Fig. 6 that the highest error in each mode corresponds to the same point in the parameter space indicating consistency in our modeling. Furthermore, as we only deal with mass ratio q = 1 waveforms, the contribution of the higher modes are expected to be negligible compared to the dominant (2, 2) mode (see for example, Ref. [95]). Therefore, even though the (3,3) and (4, 4) modes have larger relative errors compared to the (2, 2) mode, their contribution to the total error is much smaller. This can be verified by comparing the full waveform errors to the (2,2) mode errors in Fig. 6.

Frequency domain mismatch with time/phase optimization
In this section, we estimate leave-one-out crossvalidation errors by computing mismatches between the NR waveform and the trial surrogate waveform in the frequency domain. The frequency domain mismatch between two waveforms, h 1 and h 2 is defined as: whereh(f ) indicates the Fourier transform of the complex strain h(t), * indicates a complex conjugation, Re indicates the real part, and S n (f ) is the one-sided power spectral density of a GW detector. Before transforming the time domain waveform to the frequency domain, we first taper the time domain waveform using a Planck window [96], and then zero-pad to the nearest power of two. The tapering at the start of the waveform is done over 1.5 cycles of the (2, 2) mode. The tapering at the end is done over the last 20M . Once we obtain the frequency domain waveforms, we compute mismatches following the procedure described in Appendix D of Ref. [9]. The mismatches are optimized over shifts in time, polarization angle, and initial orbital phase. We compute the mismatches at 37 points uniformly distributed on the sky of the source frame, and use all available modes for the surrogate model.
We consider a flat noise curve S n (f ) = 1 as well as the Advanced-LIGO design sensitivity Zero-Detuned-HighP noise curve from Ref. [97]. We take f min to be the frequency of the (2, 2) mode at the end of the initial tapering window while f max is set at 4f peak 22 , where f peak 22 is the frequency of the (2, 2) mode at its peak. This ensures that the peak frequencies of all modes considered in our model are captured well, and we have confirmed that our mismatch values do not change for larger values of f max . Note that when computing mismatches using Advanced LIGO noise curve, for masses below ∼ 70M , f min is greater than 20Hz, meaning that the signal starts within the detector sensitivity band.
The mismatches computed using the flat noise curve are shown in the left panel of Fig. 7. The histograms include mismatches for all 47 NR waveforms and source-frame sky locations. We find that the typical surrogate mismatches are 10 −5 − 10 −3 , which are comparable to but larger than the NR errors. As an example, Fig. 8 shows the surrogate and NR waveforms for the case that leads to the largest mismatch in the left panel of Fig. 7.
In Fig. 3 Fig. 7 shows the mismatches computed using advanced LIGO design sensitivity noise curve [97] for different total masses M of the binary. For each M , we compute the mismatches for all 47 NR waveforms and source-frame sky locations and show the distribution of mismatches using vertical histograms known as violin plots. Over the mass range 20 − 180M , the surrogate mismatches are at the level of ∼ 10 −4 − 10 −3 but with an extending tail as before. However, we note that these errors are typically smaller than the mismatches for other eccentric waveform models [52][53][54].

B. Mode mixing
NR waveforms are extracted as spin-weighted spherical harmonic modes [98,99]. However, during the ringdown, the system can be considered a single Kerr black hole perturbed by quasinormal modes; perturbation theory tells us that the angular eigenfunctions for these modes are the spin-weighted spheroidal harmonics [100,101]. A different decay rates. This more complicated behavior is referred to as mode mixing, since power flows between different spherical-harmonic modes [102]. This mixing is particularly evident in the (3, 2) mode as significant power of the dominant (2, 2) spherical-harmonic mode can leak into the (3, 2) spherical-harmonic mode. As the surrogate accurately reproduces the spherical harmonic modes from the NR simulations, it is also expected to capture the effect of mode mixing without any additional effort [10]. We demonstrate this for an example case in Fig. 9 where we plot the amplitude of individual modes of the waveform during the ringdown. We show that the mode mixing in the (3, 2) mode is effectively recovered by the surrogate model.

C. NRSur2dq1EccRemnant errors
In addition to the waveform surrogate, we also build a remnant surrogate model, NRSur2dq1EccRemnant, that predicts the mass and spin of the final BH left behind after the merger. This is the first such model for eccentric BBHs (but see e.g. Refs. [63,64]). Figure 10 shows the cross-validation errors of NRSur2dq1EccRemnant in predicting the remnant mass and spin. We find that NRSur2dq1EccRemnant can predict the final mass and spin with an accuracy of 5 × 10 −4 M and 2 × 10 −3 respectively. We further compute the errors for a noneccentric remnant model, NRSur3dq8Remnant [83], when compared against the same eccentric NR simulations, finding that errors in NRSur3dq8Remnant are comparable with NRSur2dq1EccRemnant errors. This suggests that noneccentric remnant models may be sufficient for equalmass nonspinning binaries with eccentricities e ref ≤ 0.2. However, we expect such models to disagree with eccentric simulations in the more general case of unequal-mass, spinning binaries (see for e.g. Ref. [43]).
Even though these amplitude and phase differences are computed at q = 1, we treat them as a proxy for the modulations due to eccentricity at any q. We then add these modulations to the amplitude and phase of h 0 m (q), the noneccentric surrogate model evaluated at the given q, to get the full amplitude and phase: The final surrogate prediction, which we view as a new, simple model NRSur2dq1Ecc+, is then: To assess the accuracy of NRSur2dq1Ecc+ we compare against eight publicly available eccentric NR simulations with q = 2 and q = 3 [51,73]. These NR waveforms are shorter in length than the ones used to train our surrogate model. To ensure fair comparison between surrogate predictions and NR waveforms, we build a test surrogate 2 which is parameterized by e ref and l ref at In Fig. 11, we show mismatches computed using the advanced LIGO design sensitivity noise curve, between the NRSur2dq1Ecc+ model and eccentric NR data at q = 2, 3. We include all modes available in the model while computing the mismatch. For simplicity, we only consider a single point in the source-frame sky, with an inclination angle of π/3. These differences cannot be accounted for by a time or phase shift, therefore, mean anomaly is an important parameter to include for waveform modeling and data analysis of eccentric binaries.
waveforms for small eccentricities. However, we advise caution with extrapolation-type procedures in general.

E. Importance of mean anomaly for data analysis
Many existing waveform models [52][53][54][55] for eccentric binaries parameterize eccentric characteristics of the waveform by only one parameter e ref while keeping l ref fixed. We, however, use both e ref and l ref as parameters in our model. We find that not allowing l ref as an independent parameter results in large modeling error, indicating that the mean anomaly is important to consider when modeling the GW signal from eccentric binaries.
To demonstrate the importance of mean anomaly also in data analysis, we present a simple study. We generate As we already account for allowed time and frame shifts when computing the mismatch, ignoring this difference can lead to modeling errors or biased parameter estimation. In the right panel of Fig. 13, we show the waveform amplitude for the cases with l ref = 0 and l ref = π. The clear differences in the amplitude reinforce our assertion that this mismatch cannot be accounted for by a time or frame shift.

V. CONCLUSION
We present NRSur2dq1Ecc, the first eccentric NR surrogate waveform model. This model is trained on 47 NR waveforms of equal-mass nonspinning BBH systems with eccentricity e ref ≤ 0.2, defined at a reference time t ref = −5500M before the waveform peak. The model includes the (2, 2), (3,2) and (4, 4) spin-weighted spherical harmonic modes. Due to the symmetries of the equal-mass, nonspinning systems considered here, this is equivalent to including all ≤ 3 and (4, ±4) modes, except the m = 0 modes. This is the first eccentric BBH model that is directly trained on eccentric NR simulations and does not require that the binary circularizes before merger. We also present NRSur2dq1EccRemnant, the first NR surrogate model for the final BH properties of eccentric BBH mergers. This model is also trained on the same set of simulations. We use Gaussian process regression to construct the parametric fits for both models. Both NRSur2dq1Ecc and NRSur2dq1EccRemnant will be made publicly available in the near future.
Through a leave-one-out cross-validation study, we show that NRSur2dq1Ecc accurately reproduces NR waveforms with a typical mismatch of ∼ 10 −3 . We further demonstrate that our remnant model, NRSur2dq1EccRemnant, can accurately predict the final mass and spin of the merger remnant with errors 5 × 10 −4 M and 2 × 10 −3 respectively. We showed that despite being trained on equal-mass binaries, NRSur2dq1Ecc can be reasonably extended up to mass ratio q ≈ 3 with mismatches 10 −2 for eccentricities e ref 0.05 at t ref = −2000M . Finally, we demonstrate that the mean anomaly, which is often ignored in waveform modeling and parameter estimation of eccentric binaries, is an important parameter to include. Exclusion of mean anomaly can result in poor modeling accuracy and/or biased parameter inference.
The NR simulations used for this work were performed using the Spectral Einstein Code (SpEC) [65]. SpEC's development efforts have been primarily focused on evolutions of binary black hole systems in quasi-circular orbits [73]. To efficiently generate accurate training data for high eccentricity systems, it may be necessary to improve certain algorithmic subroutines. For example, as noted in Sec. III B 3, we found it difficult to achieve target values of (e ref , l ref ) at a reference time before merger. We also noticed that the waveform's numerical error was noticeably larger near pericenters, suggesting better adaptive mesh refinement algorithms [103] may be necessary for highly eccentric simulations. We have also explored several data decomposition techniques and parametrizations for building eccentric NR surrogate models, which can guide strategies for future models. Our final framework for building eccentric NR surrogates is quite general, and we expect that it can be applied straightforwardly to higher dimensional parameter spaces including unequal masses and aligned-spins. We leave these explorations to future work.

APPENDIX
In this Appendix, we describe various alternate modeling strategies we pursued before deciding on the formalism presented in the main text.

Appendix A: Choice of data decomposition
In this work, we have modeled the amplitude (A 22 ) and phase (φ 22 ) of the (2, 2) mode by modeling the residual (∆A 22 , ∆φ 22 ) of these quantities with respect to a quasicircular NR waveform (cf. Sec. III C). Alternatively, one could instead model the amplitude and frequency (or their residuals), and then integrate the frequency to obtain the phase. The frequency of the (2,2) mode is given by where φ 22 is defined in Eq. (10). The corresponding residual is given by: where ω 0 22 is the frequency of (2, 2) mode for the quasicircular NR waveform.
We, therefore, explore four different data decomposition strategies for the (2, 2) mode, summarized below: In order to explore the effectiveness of these strategies, we build a separate surrogate model using each strategy. When building the frequency surrogates (ω 22 or ∆ω 22 ) we use a basis tolerance of 10 −3 rad/M . For A 22 (φ 22 ) we use the same tolerance as used for ∆A 22 (∆φ 22 ) in Section III D. We compute the normalized L 2 -norm between the NR data and each surrogate approximation using Eq. (17). In Fig. 14, we show the surrogate errors for all four different strategies. We find that modeling the frequency ω 22 or residual frequency ∆ω 22 yields at least two-to-three orders of magnitude larger E than when we model the phase φ 22 or ∆φ 22 . Furthermore, modeling  17)) for the four different decomposition strategies we consider. We find that modeling the residual amplitude ∆A22 and residual phase ∆φ22 yields the least errors.
the residual amplitude ∆A 22 proves to be slightly more accurate than the case where we model the amplitude A 22 directly. Therefore, in the main text, we build surrogate models of the residual {∆A 22 , ∆φ 22 } (cf. Sec. III C).

Appendix B: Choice of fit parameterization
When building the surrogate models in main text, fits across parameter space are required for the waveform model as well as the remnant model (cf. Sec. III D). which can be useful if the eccentricity varies over several orders of magnitude in the NR dataset.
Similarly to the previous section, to explore the effectiveness of these strategies, we build a separate surrogate model using each strategy. Here, however, we consider all modes included ( , m) = (2, 2), (3,2), (4,4) and evaluate E errors [cf. Eq. (17)]. In Fig. 15, we show E errors for each parameterization strategy. We find that while the alternative strategies using either log 10