\section{Introduction}\label{sec:introduction}

The detections of gravitational wave signals \cite{PhysRevLett.116.061102, PhysRevLett.116.241103,Abbott:2017vtc,Abbott:2017gyy,Abbott:2017oio,TheLIGOScientific:2017qsa,LIGOScientific:2018mvr, Romero-Shaw:2019itr} have been found  to be consistent with models of the waveform emitted from the merger of compact objects under the assumption of quasicircularity of the binary's orbit prior to the merger. The assumption of quasicircularity motivated by the efficient circularization
of binaries as a consequence of the emission of gravitational waves \cite{PhysRev.131.435,PhysRev.136.B1224} simplifies significantly the complexity of the signal and has accelerated the development
 of inpiral-merger-ringdown (IMR) waveform models: several mature IMR models for quasicircular  coalescences, i.e.~neglecting eccentricity, are now publicly available \cite{Husa2015,PhysRevD.93.044007,phenomp,seobnrv4,london2018,PhysRevD.98.084028,Blackman:2017dfb,PhysRevResearch.1.033015, PhysRevD.99.064045,Khan:2018fmp,Khan:2019kot,phenX,phenXHM,phenT}, and are being used to search and infer the parameters of observed binary black hole systems \cite{LIGOScientific:2018mvr}.

Recently, population synthesis studies \cite{Belczynski:2016obo,Park:2017zgj,PhysRevD.97.103014,Samsing:2013kua} have shown that active galactic nuclei and globular clusters can host a population of moderate and highly eccentric binaries emitted in the frequency band of ground-based detectors. Therefore, the increase in sensitivity of the detectors will increase the likelihood of detecting binary systems with non-negligible eccentricities. The modeling of the gravitational waveforms from eccentric black hole binaries is complicated by the addition of a new timescale to the binary problem, the periastron precession \cite{1985AIHS...43..107D}. This new timescale induces oscillations in the waveforms due to the asymmetric emission of gravitational radiation between the apastron and periastron passages.

The orbits of eccentric black hole binaries are typically described using the quasi-Keplerian (QK) parametrization \cite{PhysRevD.70.064028}, which is currently known up to 3 post-Newtonain (PN) order \cite{PhysRevD.70.104011}. This parametrization has proved to be a key element in developing inspiral PN waveforms \cite{PhysRevD.93.064031,Moore:2018kvz,Moore:2019xkm,Loutrel:2018ydu,Moore_2019}. The generation of IMR eccentric models relies on the connection of an eccentric PN inspiral with a circular merger \cite{PhysRevD.97.024031,PhysRevD.98.044015}. Alternatively, one can substitute the PN waveform with one produced within the effective one body (EOB) formalism describing an eccentric inspiral \cite{PhysRevD.96.044028,PhysRevD.96.104048,Chiaramello:2020ehz}. Some eccentric IMR waveform models show good agreement with numerical waveforms up to $e \sim 0.2$ for nonspinning configurations \cite{PhysRevD.97.024031}. Recent work has shown possible extensions of the eccentric PN and EOB inspiral waveforms to include spin effects \cite{PhysRevD.98.104043,Chiaramello:2020ehz}.
A key step in the generation of IMR waveform models is the production of hybrid waveforms~\cite{PhysRevD.82.064016,Ohme_2012,PhysRevD.84.064029,MacDonald_2011,Ajith_2012,PhysRevD.84.064013,Bustillo:2015ova} between PN/EOB inspiral and numerical relativity (NR) waveforms. The hybridization procedure consists of smoothly attaching a PN/EOB inspiral waveform to a NR one in order to get the full description of the gravitational radiation of the binary system. The generation of datasets of hybrid waveforms has been used in the quasicircular case to calibrate and validate the accuracy of IMR waveform models \cite{phenX,phenXHM,PhysRevD.98.084028,Blackman:2017dfb,PhysRevResearch.1.033015}.

In this paper we present the input data and some key tools required  for the development of an IMR eccentric waveform
model calibrated to eccentric hybrid PN-NR waveforms.
In Sec. \ref{sec:NRcatalog} we first present our NR catalog of nonspinning and spinning eccentric binaries, computed with the private \texttt{BAM} code \cite{PhysRevD.77.024027} and the open-source EinsteinToolkit (\texttt{ET}) \cite{Loffler:2011ay,maria_babiuc_hamilton_2019_3522086}. This includes a discussion of our procedure to specify the initial parameters of the eccentric simulations in Sec.  \ref{sec:NRcatalogA}, a study of
the remnant quantities in Sec.~\ref{sec:NRcatalogB}, and  a new method for measuring the eccentricity of NR waveforms with arbitrarily high eccentricity  in Sec.~\ref{sec:NRcatalogC}. We find that the final spin and mass are consistent within the error estimates with the quasicircular case, which extends the study in \cite{PhysRevD.77.081502} to the eccentric spinning case.
We hybridize the dominant gravitational waveform $(l=2,m=2)$ mode between numerical relativity and post-Newtonian waveforms in Sec. \ref{sec:PNEcc}. This will provide the input data for future work on constructing waveform models that contain the inspiral, merger, and ringdown,  and it allows us to perform injections into detector noise which contain a long inspiral phase.
In Sec. \ref{sec:PE} we use such injections of hybrid waveforms, as well as those of pure numerical relativity waveforms, to study the parameter biases introduced when using quasicircular waveform models to estimate the parameters of highly eccentric spinning systems.
Unless explicitly noted, we are working in geometric units $G=c=1$. To simplify expressions we will also set the total mass of the system $M=1$ in Secs. \ref{sec:NRcatalog} and \ref{sec:PNEcc}. We define the mass ratio $q=m_1/m_2$ with the choice $m_1 > m_2$, so that $q>1$. We also introduce the symmetric mass ratio $\eta=q/(1+q)^2$, and we shall denote the black hole's dimensionless spin vectors by $\vec{\chi}_i= \vec{S}_i/m^2_i$ for $i=1,2$.

\section{Numerical Relativity dataset} \label{sec:NRcatalog}

\subsection{Overview} \label{sec:NRcatalog0}

We present a catalog of 60 eccentric NR simulations performed with the nonpublic \texttt{BAM} code \cite{PhysRevD.77.024027} and the open-source \texttt{ET} code \cite{Loffler:2011ay,maria_babiuc_hamilton_2019_3522086} with the multipatch \texttt{LLAMA} infrastructure \cite{Pollney:2009yz}. The numerical setup of both codes is the same as in \cite{PhysRevD.99.023003}. Most of the simulations are run with the \texttt{ET} code using the \texttt{LLAMA} module due to its ability to extract the waves at larger extraction radii.
The different simulations and their initial conditions are described in Table \ref{tab:tabNR3}. In Fig. \ref{fig:NRsims} we show our choices of mass ratio $q$, initial eccentricity $e_0$, and effective spin parameter, $\chi_{\text{eff}}= (m_1 \chi_{1,z}+m_2 \chi_{2,z})/(m_1+m_2)$. We have also added  20 public eccentric SXS simulations presented in \cite{PhysRevD.98.044015}.

\begin{figure}[ht!]
\centering
\includegraphics[scale=0.35]{plot3DBAMETSXSbo.pdf}
\caption{Initial eccentricity $e_0$, mass ratio $q$ and effective spin parameter $\chi_{\text{eff}}= (m_1 \chi_{1,z}+m_2 \chi_{2,z})/(m_1+m_2)$ for the numerical relativity simulations generated with the \texttt{BAM}, EinsteinToolkit and \texttt{SpEc} \cite{SpEC} codes. The thick black line represents the cases with $\chi_{\text{eff}}=0$.}
\label{fig:NRsims}
\end{figure}

\subsection{Initial parameters of eccentric NR simulations} \label{sec:NRcatalogA}

We use conformally flat Bowen-York initial data \cite{PhysRevD.21.2047} in the center-of-mass frame, where the free parameters are the spins and masses of the two black holes, the separation, and the momentum of one of the two black holes (the momentum of the second black hole is then equal in magnitude but opposite in direction). We first choose the masses and spins as displayed in Fig.~\ref{fig:NRsims}.
To be able to construct hybrid waveforms, the minimal separation, i.e. the separation at periastron, has to be large enough that
the PN approximation is still roughly valid. We then use a simple PN approximation as discussed below to compute the apastron separation required to achieve a chosen value of the eccentricity, and a further PN approximation to compute the appropriate value of the momentum corresponding to this value of the eccentricity. Owing to the simplicity, i.e., the low order, of the PN approximations used, neither the periastron separation nor the measured eccentricity will exactly coincide with the specified values.  In this study we choose our initial choice for the approximate periastron separation to be $r_{\text{min}} \sim 9M$, with slightly different values to account for mass ratio and spin effects which can significantly increase the computational cost of the simulations. We start our simulations at the apastron, where the PN approximation employed to specify the initial momentum and the agreement with the PN data that we use for hybridization, will be more accurate than during other points of the orbit.

Larger choices of eccentricity for the same configuration of masses and spins thus lead to a larger merger time and number of orbits, as one can see in Table \ref{tab:tabNR3}. For instance, focusing on simulations with identification numbers (IDs) 34, 35, and 36, one observes an increase in the merger time when increasing the initial eccentricity. This increase in merger time also implies an increase in the computational cost of the simulation.

Using the QK parametrization at Newtonian order, one can relate the initial minimum and maximum separations by

\begin{equation}
r_{\text{min}}= r_{\text{max}} \frac{1-e}{1+e}.
\label{eq:eq2}
\end{equation}
As stated above, for our simulations we choose $r_{\text{min}} \sim 9M$ such that the PN approximation is still roughly valid. Then  for $e_0 =  0.1,0.2,0.5$ Eq. \eqref{eq:eq2} implies that  $r_{\text{max}} = 11 M$, $13.5 M$, $27 M $, respectively. These values of  $r_{\text{max}}$ are rough estimates
based on a Newtonian order calculation; in practice, we slightly modify those values of initial separations to account for the increase of computational cost depending on the mass ratio and the spins of the simulations as observed in Table \ref{tab:tabNR3}. For instance, in the case of negative spin components the merger time is significantly reduced \cite{PhysRevD.97.084002}; thus, we increase $r_{\text{max}}$  for $e=0.1,0.2$ cases to produce longer NR waveforms which  are  easier to hybridize afterward.

To produce initial data for a desired eccentricity we then make use of Eq. (3.25) of \cite{PhysRevD.99.023003} to perturb the initial tangential momentum of the black holes by a factor $\lambda_t$ from its quasicircular value. The expression for $\lambda_t$ in terms of the eccentricity at 1PN order is
\begin{equation}
\lambda_t (r,e_0, \eta, \text{sign})= 1+ \frac{e_0 }{2}\times \text{sign}\times \left[ 1-\frac{ 1}{r}(\eta +2) \right],
\label{eq:eq1}
\end{equation}
where $\eta$ is the symmetric mass ratio, $r$ is the orbital separation and $\text{sign}=\pm 1$ depends on the initial phase of the eccentricity estimator \cite{PhysRevD.99.023003}.  We refer the reader to Sec. III D of \cite{PhysRevD.99.023003} for an explicit derivation of Eq. \eqref{eq:eq1}. Taking Eq. \eqref{eq:eq1}, we compute the correction factor applied to the momentum as the mean between the inverse of the expression with the plus sign plus the expression with the minus sign,

\begin{equation}
\begin{split}
\bar{\lambda}^0_t (r,e_0, \eta) & = \frac{1}{2} \left[ \lambda_t(r,e_0, \eta,+1)^{-1}+\lambda_t(r,e_0, \eta,-1) \right] \\
& = \frac{8 r^2-e_0^2 (\eta -r+2)^2}{4 r (e_0 (-\eta +r-2)+2 r)}.
\end{split}
\label{eq:eq3}
\end{equation}
We use the combination of factors in Eq. \eqref{eq:eq3} because we have experimentally tested it to see that it works more accurately than just specifying a value of $\lambda_t(r,e,\text{sign})$ with a given sign. In Table \ref{tab:tabNR3} one can compare the value of the desired initial eccentricity, $e_0$, specified in Eq. \eqref{eq:eq3}, and the actually measured initial eccentricity, $e_\omega $, from the orbital motion of the simulation. Both quantities are also displayed in Fig. \ref{fig:eccMeasEccDes}, where we have differentiated  among nonspinning and positive and negative spin simulations. The results point out that the use of Eq. \eqref{eq:eq3} produces differences of less than $10 \%$ between $e_\omega$ and $e_0$ in nonspinning cases at low eccentricities of the order of $0.1$. However, when spins are present or the eccentricities are higher, the inaccuracy of the formula becomes manifest, with differences of the order of $20\% -30 \%$, this is due to the fact that Eq. \eqref{eq:eq1} was derived assuming a nonspinning binary in the low eccentric limit. Additionally, one can check in Table \ref{tab:tabNR3} and Fig. \ref{fig:eccMeasEccDes} to see that the differences between $e_\omega$ and $e_0$ are smaller for the cases with positive spins than in cases with negative spins because in Eq. \eqref{eq:eq1} the radiation reaction effects, which are more significant for negative spins, are also not taken into account.

\subsection{Final state of spinning eccentric systems} \label{sec:NRcatalogB}

We compare the final state of the eccentric NR simulations with the predicted final mass and final spin of the quasicircular (QC) NR fits \cite{PhysRevD.95.064024} as an indicator of circularization of the coalescence process as the binary merges. This is basically an extension of \cite{PhysRevD.77.081502} to the eccentric spinning case with more moderate values of the eccentricity but longer NR evolutions.

The final mass and final spin of the simulations are computed using the apparent horizon (AH) of the remnant black hole and are shown in Table  \ref{tab:tabNR3}.
The magnitude $S$ of the angular momentum of the final black hole can be computed from the integral
\begin{equation}
S = \frac{1}{8 \pi} \oint_{\mbox{AH}} K_{ij} n^i \phi^j dA,
\label{eq:eq5}
\end{equation}
see the discussion in \cite{PhysRevD.75.064030, PhysRevD.67.024018}. Here for the \texttt{BAM} code \cite{PhysRevD.77.024027} the vector $\phi^j$ is a coordinate-based approximation to the (approximate) axial Killing vector of the black hole horizon as in \cite{PhysRevD.75.064030}, and for the EinsteinToolkit code  the  \texttt{QuasiLocalMeasures} module is used, which constructs
an approximate Killing vector with rotational symmetry around the spin axis as in \cite{PhysRevD.54.4899, Thornburg:2003sf}.
The vector $n^i$ is a spacelike unit normal to the horizon surface and $K_{ij}$ is the extrinsic curvature.
The final mass can be computed from the Christodoulou formula in terms of the black hole (BH) angular momentum and
AH area $A$ as
\begin{equation}
M_f = \sqrt{M^2_{irr}+ \frac{S^2}{4 M^2_{irr} }}, \quad M_{irr}=\sqrt{\frac{A}{16\pi}}.
\label{eq:eq4}
\end{equation}
where $M_{irr}$ is the irreducible mass. The dimensionless final spin can then be computed as $\chi_f = S/M_f^2$.

In Fig. \ref{fig:FinalStateErrors} we have computed the absolute and relative errors between the eccentric simulations and the quasicircular NR final mass and final spin fitting formulas \cite{PhysRevD.95.064024},
\begin{equation}
\Delta X = \left| \frac{X^{NR}}{X^{QC}}-1 \right| \times 100, \quad X=M_f \text{ or }\chi_f.
\label{eq:eq6}
\end{equation}
The results in Fig. \ref{fig:FinalStateErrors} show that the differences in the final spin are generally higher than in the final mass. However, the differences with respect to the quasicircular fitting values are as high as $\sim 1 \%$ which is entirely consistent with numerical errors and gauge artifacts on the apparent horizon surfaces and inaccuracies in the fits. Hence, we can conclude that within the current knowledge of systematic errors (compare with \cite{PhysRevD.95.064024}), the final state of the eccentric simulations, up to the values of eccentricity studied here, is consistent with the quasicircular values. Identifying small physical deviations between the quasicircular and eccentric final states will require numerical simulations with improved error estimates.

\begin{figure}[hbt!]
\centering
\begin{minipage}[hbt!]{\columnwidth}
\captionsetup{justification=centering}
\includegraphics[scale=0.42]{plotLogAbsErrorFinalStateBO.pdf}
\end{minipage}

\begin{minipage}[hbt!]{\linewidth}
\captionsetup{justification=centering}
\includegraphics[scale=0.42]{plotLogErrorFSbo.pdf}
\end{minipage}
\caption{ (Top panel) The absolute difference between the final mass and spin of the simulations and the QC NR fits as a function of the ID of the simulations in Table \ref{tab:tabNR3}. (Bottom panel) The absolute relative error for the phase and amplitude, $\Delta X = |X^{NR}/X^{QC}-1|\times 100$ for $X=M_f,\chi_f$, relative error of the final mass and final spin of the simulations against the QC NR fits as a function of the ID of the simulations in Table \ref{tab:tabNR3}.}
\label{fig:FinalStateErrors}
\end{figure}

\subsection{Measuring the eccentricity of highly eccentric systems} \label{sec:NRcatalogC}

This subsection aims to extend the discussion on the measurement of the eccentricity in NR presented in \cite{PhysRevD.99.023003} to highly eccentric systems. An eccentricity parameter is chosen to describe the noncircularity of orbits, such that for bound orbits its value ranges between 0 and 1, corresponding to circular and extremely elliptical configurations, respectively. Such an eccentricity can  be defined naturally only in Newtonian gravity, whereas in general relativity the eccentricity is a gauge dependent quantity. To measure the eccentricity in NR data one defines quantities known as eccentricity estimators, which estimate the eccentricity from the relative oscillations of a certain combination of dynamical quantities such as the orbital separation or orbital frequency, or wave quantities like the amplitude or frequency of the $(l,m)=(2,2)$ mode. All of these different estimators are usually defined such that they agree in the Newtonian limit and in the low eccentricity limit.

In \cite{PhysRevD.99.023003}, where we studied the reduction of residual eccentricity in initial datasets, we chose our eccentricity estimator based on the orbital frequency as
\begin{equation}
e_\omega(t)= \frac{\omega (t)-\omega(e=0)}{2 \omega(e=0)} \equiv  \frac{\omega (t)-\omega^{\text{fit}}(t)}{2 \omega^{\text{fit}}(t)},
\label{eq:eq8}
\end{equation}
 where $\omega (t)$ is the orbital frequency of the simulation and $\omega (e=0)$ is the orbital frequency in the quasicircular limit. We note that when dealing with numerical simulations, the quasicircular frequency in Eq. \eqref{eq:eq8} is typically replaced by a fit, $ \omega^{\text{fit}}(t)$, of the nonoscillatory part of the frequency~\cite{PhysRevD.99.023003}. This eccentricity estimator  is  used largely to measure the residual eccentricity of NR simulations of quasicircular black hole binaries. We remark that while in Eq. \eqref{eq:eq8} we decide to use the orbital frequency calculated from the BH motion, one can also  use the gravitational wave frequency extracted from the waves as discussed below.  Furthermore, gauge effects can impact the eccentricity measurement from the orbital frequency of NR codes as extensively discussed in the small eccentricity limit in \cite{Purrer:2012wy,PhysRevD.82.124016}. Here we follow the practice used in the literature \cite{Purrer:2012wy,PhysRevD.99.023003,PhysRevD.82.124016,PhysRevD.83.104034} to avoid contamination of the eccentricity measurement through the gauge quantities, like the choice of the value of the $\eta$-parameter in the gamma driver condition \cite{Alcubierre:2002kk}, which can  lead to residual oscillations in the orbital frequency complicating the determination of the eccentricity.

In \cite{PhysRevD.99.023003} we argued that the procedure shown there, based on Eq. \eqref{eq:eq8}, to measure the eccentricity is limited to values as high as $e\sim 0.1$ due to the lack of an accurate ansatz to fit the higher order contributions beyond the sinusoidal contribution. While the lack of an ansatz for high eccentricities is a clear limitation, the use of Eq. \eqref{eq:eq8} biases the eccentricity measurement due to its reliance on a noneccentric fit of the orbital frequency and due to the fact that Eq. \eqref{eq:eq8} for high eccentricities does not reduce to the common definition of eccentricity in the Newtonian limit.

Therefore, we decide to change to another estimator \cite{PhysRevD.66.101501}, constructed also from the orbital frequency,
\begin{equation}
e_{\omega}(t)= \frac{\omega^{1/2}_p - \omega^{1/2}_a}{\omega^{1/2}_p + \omega^{1/2}_a},
\label{eq:eq9}
\end{equation}
where $\omega_a,\omega_p$ are the orbital frequencies at apastron and periastron, respectively. The eccentricity estimator in Eq.  \eqref{eq:eq9} does not depend on any noneccentric fit of the orbital frequency. Furthermore, as shown in Appendix \ref{sec:AppendixB} the eccentricity estimator from Eq. \eqref{eq:eq8} in the Newtonian limit at high eccentricities does not reduce to the eccentricity parameter and it is not normalized, while the eccentricity estimator from Eq.  \eqref{eq:eq9} fulfills all these conditions.

\begin{figure}[!]
\centering
\captionsetup{justification=centering}
\includegraphics[scale=0.45]{EccMinMaxbo.pdf}
\caption{ Time evolution of the orbital frequency $M \omega_{orb}$, the orbital frequency at apastron $M \omega_a$, the orbital frequency at periastron $M \omega_p$ and the eccentricity estimator $e_\omega$ defined in Eq. \eqref{eq:eq9}.}
\label{fig:EccMinMaxbo}
\end{figure}

 We measure the eccentricity from the maxima and minima of the orbital frequency corresponding to the periastron and apastron passages, respectively. Additionally, we produce an interpolated function from the maxima $\omega_p$  and the minima $\omega_a$ and substitute them into Eq. \eqref{eq:eq9}, so that one can estimate the evolution of the eccentricity from those points. The interpolation is calculated using the Hermite method implemented in the function \texttt{Interpolation} in \textit{Mathematica} \cite{Mathematica}.  The new procedure to measure the eccentricity is shown in Fig. \ref{fig:EccMinMaxbo}, where the time evolution of the orbital frequency, the interpolated functions of the maxima and minima of the orbital frequency and the eccentricity are shown for the configuration with ID 60 from Table \ref{tab:tabNR3}. As expected the eccentricity is a monotonically decaying function, whose value at $t=200M$, after the burst of junk radiation, is $e_\omega=0.415 \pm 0.005$. The error in the eccentricity, $\delta e_\omega$, is computed using error propagation: from Eq. \eqref{eq:eq9} we obtain
\begin{equation}
\delta e_{\omega}= \frac{  \delta \omega}{\left( \omega_a^{1/2} + \omega_p^{1/2} \right)^2} \left[ \frac{\omega_a^{1/2}}{\omega_p^{1/2}} + \frac{\omega_p^{1/2}}{\omega_a^{1/2}} \right],
\label{eq:eq91}
\end{equation}
where we have assumed $\delta \omega_a =\delta \omega_p= \delta \omega$. Motivated by the results of the error in the convergence analysis of the orbital frequency in \cite{NRUIB} we have taken as a conservative estimate   $\delta  \omega =0.0001$.  The error estimate of Eq. \eqref{eq:eq91} is the statistical error associated with the eccentricity measurement taking into account the error of the orbital frequency from different resolutions of the NR simulations. We remark that this error does not take into account systematics coming from the use of a different eccentricity estimator, nor contributions from the interpolation error when the number of minima and maxima is small due to the short length of the simulations. Because of the difficulties in quantifying the systematics associated with the choice of eccentricity estimator estimator and the fact that the interpolation error is a subdominant effect for most of the simulations here, we restrict for simplicity our eccentricity error calculation to Eq. \eqref{eq:eq91}.

The main drawback of this method is that when the simulations are so short that there is only one minimum and one maximum it becomes inefficient and inaccurate. Furthermore, one could choose the frequency of the $(l,m)=(2,2)$ mode and compute the orbital frequency as $\omega_{orb} \approx \omega_{22}/2$, employing the same method discussed in this section. Nevertheless, as pointed out in \cite{Purrer:2012wy} the usage of the orbital frequency from the $(2,2)$ mode requires additional postprocessing of the data due to the presence of high frequency noise when taking a time derivative of the phase of the $(2,2)$ mode.  As a conclusion, if one has long enough highly eccentric simulations, the method introduced in this section allows one to  measure the eccentricity as a monotonically decaying function for the whole inspiral, which is a key tool to be used to construct a time domain eccentric waveform model.

\section{Hybridization of eccentric waveforms} \label{sec:PNEcc}
In the eccentric case the hybrization of the PN-NR waveforms is a challenging problem. The higher the eccentricity the stronger the interaction between the binary components at each periastron passage, which can break the post-Newtonian, weak-field and low velocity, approximation and generate a secular dephasing between both waveforms. Moreover, the lack of a general description in PN theory of eccentric black hole binary systems poses the main difficulty.
 Therefore, we briefly review the status of the PN theory for eccentric systems in Sec. \ref{sec:ReviewPN}. In Sec. \ref{sec:Exhyb} we show an example of our procedure to hybridize eccentric PN-NR waveforms.

\subsection{Review of eccentric post-Newtonian theory} \label{sec:ReviewPN}

As far as the authors know by the time of writing this communication, the orbital averaged gravitational wave energy flux for eccentric binaries is known up to 3PN order \cite{PhysRevD.80.124018} using the 3PN QK parametrization \cite{PhysRevD.70.104011}. Our strategy consists in evolving the $3.5$PN Hamilton's equations of motion in Arnowitt, Deser, and Misner transverse-traceless (ADMTT) gauge\cite{PhysRev.116.1322,Deser:1960zzc,PhysRev.117.1595} for a point particle binary,
\begin{equation}
\frac{d \bm{X}}{dt}=\frac{\partial \mathcal{H}}{\partial \bm{P}}, \quad \frac{d \bm{P}}{dt}=-\frac{\partial \mathcal{H}}{\partial \bm{X}}+\bm{F}, \quad \frac{d \bm{S}_i}{dt}=\frac{\partial \mathcal{H}}{\partial \bm{S}_i}\times \bm{S}_i, \quad i=1,2,
\label{eq:eq10}
\end{equation}
with $\bm{X}$,$\bm{P}$, and $\bm{S}_i$ being the position, momentum and spin vectors in the center-of-mass frame, $\mathcal{H}$ the Hamiltonian described in Sec. II of \cite{PhysRevD.99.023003} and $\bm{F}$ the radiation reaction force described in \cite{Buonanno:2005xu} enhanced with the eccentric contribution to the energy flux from \cite{PhysRevD.80.124018}. The eccentric term in the flux is expressed in the QK parametrization and depends only on the orbital frequency $\omega$, which is computed while evolving the system, and the eccentricity $e_t$, for which we use its 3PN expression in terms of the orbital energy and the angular momentum of the system, which are variables computed at each time step.

The solution of the PN point particle equations, Eq. \eqref{eq:eq10}, can be used to compute the gravitational radiation emitted by the system. Here, the lack of general PN expressions for the waveforms of point particles evolving on quasielliptical orbits sets a strong limitation. The instantaneous terms of the waveform multipoles are  known up to 3PN order for general non-spinning systems with arbitrary eccentricity \cite{PhysRevD.91.084040}. Recently, the complete description of the 3PN nonspinning multipoles was computed including tail, tail-of-tail, and memory terms within the QK parametrization for low eccentricities \cite{Boetzel:2019nfw, PhysRevD.100.084043}. At this point using the 3PN instantaneous terms  only \cite{PhysRevD.91.084040} introduces more errors than the quadrupole order due to the missing tail and tail-of-tail terms that enter at 1.5PN, 2.5PN and 3PN orders, respectively. Additionally, the translation of the generic solution we obtain from solving Eq. \eqref{eq:eq10} to the QK form of the waveform modes in \cite{Boetzel:2019nfw, PhysRevD.100.084043} is more involved due to the fact that they split the dynamical variables into adiabatic and postadiabatic contributions. Therefore, we will restrict here to the quadrupole formula to generate the $(l,m)=(2,2)$ mode and leave for future work the generation of full 3PN waveforms, which will additionally allow us to construct multimode eccentric hybrids.

\subsection{Hybridization example} \label{sec:Exhyb}

The hybridization of PN and NR waveforms consists of determining the time shift and phase offset which minimizes the difference between both waveforms in a certain time window. This hybridization procedure is well established in the quasicircular case \cite{Bustillo:2015ova,PhysRevD.82.064016,MacDonald:2011ne,PhysRevD.87.024009,PhysRevD.84.064013}.  The time shift is usually computed by  minimizing  a suitable quantity
that measures disagreement of the two waveforms, such as an overlap integral \cite{PhysRevD.87.024009,PhysRevD.84.064013}, or the deviation between phase or frequency of the $(2,2)$ mode \cite{Bustillo:2015ova}. However, in the eccentric case the calculation of the time shift requires alignment of the peaks due to eccentricity of both waveforms in the hybridization window.  This alignment is complicated to obtain with the phase because the peaks corresponding to each periastron passage are not very pronounced and they are difficult to estimate. One could use the frequency of the $(2,2)$ mode. However, it is a quantity obtained from a time derivative of the phase, which for NR waveforms tends to be noisy.  As a consequence, for simplicity we use the amplitude of the $(2,2)$ mode to determine the time shift of the waveform because it is a clean quantity with clearly defined peaks. We remark that aligning the oscillations of the amplitude of PN and NR waveforms in a certain hybridization window is equivalent to minimizing their difference as the maximum agreement between both quantities is obtained when they are aligned at the eccentric peaks.

As an example, we take the NR simulation with ID 60 of Table \ref{tab:tabNR3}, which is a mass ratio $q=4$ nonspinning configuration with an initial eccentricity of $e^0_\omega=0.415 \pm 0.005$ and initial orbital separation at apastron $D_0=27.5M$. We take the initial conditions of the NR simulation defined by the initial position vector, momenta (velocities in the case of \texttt{SpEc} waveforms), and dimensionless spin vectors: $\{\bm{X},\bm{P}/\bm{V},\bm{S}_{1},\bm{S}_{2}\}_{t=0}$. The fact that PN and NR coordinates for the initial data agree up to 1.5PN order \cite{Tichy:2002ec,Yunes:2006iw,PhysRevD.74.104011} makes this identification a good approximation. However, we have checked to see that the differences between the PN and NR initial conditions can produce discrepancies between the NR and PN waveforms of the order of $10 \%$.

To leverage these differences we decide to modify the initial condition vector of the PN evolution by modifying the initial separation by a $\delta r$ such that the difference in the amplitude of the Newman-Penrose scalar, $\psi_4$, for the $(2,2)$ mode between PN and NR is minimal. In our example we obtained $\delta r = 0.08$. The outcome of such a calculation can be observed in the top panel of Fig. \ref{fig:HybridErrors}, where the time domain amplitudes of the PN and NR waveforms are shown. We do not show the full time domain range of the hybrid waveform in the top panel of  Fig. \ref{fig:HybridErrors} to better display the matching PN/NR region. The procedure is also applied to eccentric aligned-spin configurations. We find that initial highly eccentric configurations require larger $\delta r$ than low eccentric ones, and that the hybridization errors for high negative spins, where radiation reaction plays a dominant role, are 1 order of magnitude higher than for nonspinning or low spins due to the lack of expressions for PN spinning eccentric waveforms.

The procedure to construct the hybrid waveform is similar to the one presented in \cite{Bustillo:2015ova}.  We first choose the matching region to be after the junk radiation burst; in our particular case we take $t/M \in (275,375)$, which corresponds to less than one gravitational wave cycle as shown in the top panel of Fig.       \ref{fig:HybridErrors}.  To understand the choice of this short hybridization window for eccentric waveforms, we first explain the criteria for hybridizing quasicircular ones, following \cite{QCHybrids}. Quasicircular waveforms are  hybridized over several cycles as the low frequency approximant, typically EOB, is very accurate and resembles faithfully the NR behavior during the late inspiral. Furthermore, hybridization over several cycles is required to accurately compute the time alignment between waveforms by averaging out residual oscillations due to eccentricity and high frequency numerical noise coming from NR.
In the eccentric case, the time alignment is much easier to compute as the peaks in the gravitational wave (GW) frequency ease such an alignment, so there is no need to use several cycles. Moreover, the inaccuracy of the current low frequency eccentric approximants also sets a clear limitation to faithfully reproduce the NR waveforms along several cycles. Hence, we have chosen a small hybridization window to ensure small errors in the GW amplitude and frequency between PN and NR. Choices of hybridization window including several cycles make that error increase to $10\%$ or larger depending on the case, due to the inaccuracy of the PN approximant. We have also checked to see that the election of different peaks for hybridization in the GW amplitude (in the inspiral regime)  does not significantly  change the errors, maintaining them below $1 \%$  as quoted in the lower panel of Fig. 4.

After choosing the hybridization region, we have to compute the time shift  $\tau$  and phase offset $\varphi_0$, which reduce the difference between the PN and NR waveforms in the matching window,
\begin{equation}
h^{PN}(t)= e^{i \varphi_0} h^{NR}(t+\tau).
\label{eq:eq101}
\end{equation}
To align the waveforms in time we choose $\tau$ such that it minimizes the amplitude difference along the matching window. For the phase offset  we decide to align the phases at the beginning of the window, $\varphi_0=\phi^{NR}(t_0 - \tau)-\phi^{PN}(t_0)$, where $t_0$ is the initial time of the window. Once $\tau$ and $\varphi_0$ are calculated the hybrid waveform is constructed as a piecewise function
\begin{equation}
h^{hyb}(t)=
\begin{cases
     e^{i \varphi_0} h^{PN}(t+\tau)  \hspace{2.95cm}   \text{ if } t< t_1 \\
     w^{-}(t)e^{i \varphi_0} h^{PN}(t+\tau)+  w^{+}(t) h^{NR}(t) \quad \text{if } t_1 < t < t_2  \\
 h^{NR}(t)    \hspace{4.12cm} \text{if }   t  > t_2
\end{cases}
\label{eq:eq11}
\end{equation}
where $t_1=275M$ and $t_2=375M$. The functions $w^{\pm}(t)$ denote the blending functions defined in the interval $\left[t_1,t_2 \right]$ that monotonically go from 0 to 1 and 1 to 0, respectively,
\begin{equation}
w^+(t)_{[t_1,t_2]}= \frac{t-t_1}{t_2-t_1}, \quad w^-(t)_{[t_1,t_2]}= 1- w^+(t).
\label{eq:eq12}
\end{equation}
 The result of the application of such a hybridization procedure can be observed in the bottom panel of Fig. \ref{fig:HybridErrors}, where the absolute value of the relative errors between the hybrid and NR amplitudes and frequencies are shown. The quantity $\Delta X$ is defined as $\Delta X= |X^{hyb}/X^{NR}-1|\times 100$ for $X=A_{22},\omega_{22}$. The errors in the gravitational wave frequency and amplitude are both below  $1 \%$, with those for the amplitude being slightly smaller due to the choice of the amplitude as the quantity with which to minimize the agreement between PN and NR waveforms.

\begin{figure}[hbt!]
\centering
\begin{minipage}[hbt!]{\columnwidth}
\captionsetup{justification=centering}
\includegraphics[scale=0.45]{plotHybridNRPNboDefzoom.pdf}
\end{minipage}

\begin{minipage}[hbt!]{\linewidth}
\captionsetup{justification=centering}
\includegraphics[scale=0.45]{plotHybridErrorsboOmega.pdf}
\end{minipage}
\caption{ (Top panel) The time domain amplitude of the $| r \psi^{2,2}_{4} |$ mode. The blue thick and the red dashed curves represent the PN and NR waveforms, respectively, and the vertical black lines highlight the hybridization window. (Bottom panel) The absolute value of the relative error for the gravitational wave frequency and amplitude, $\Delta X = |X^{hyb}/X^{NR}-1| \times 100$ for $X=\omega_{22},A_{22}$, of the hybrid against the NR waveform in the matching region is displayed.}
\label{fig:HybridErrors}
\end{figure}

Finally, note that the PN waveform used to produce the hybrid is evolved backward in time from $D_0/M=27.5$ to $D_f/M=60$.  This makes the initial eccentricity increase with respect to the NR waveform.
 Next, we explicitly show the systematics affecting the measurement of the initial eccentricity of the hybrid. We display in Fig. \ref{fig:OmegaComparisons} the time evolution of the orbital frequency for the same hybrid waveform of Fig. \ref{fig:HybridErrors},  $\omega_{orb} \approx \dot{\phi}_{22}/2$, computed from the phase of the $(2,2)$ mode of the Newman-Penrose scalar and  the strain computed using the fixed-frequency integration algorithm \cite{Reisswig:2010di}. We also compute the orbital frequency from the PN dynamics as
\begin{equation}
\omega=\left| \frac{\bm{v}\times \bm{r}}{r^2}\right|
\label{eq:eq100}
\end{equation}
where $r=|\bm{r}|$, and $\bm{v},\bm{r}$ are the velocity and the position vectors in the center-of-mass frame. The curves from Fig. \ref{fig:OmegaComparisons} indicate that the orbital frequencies computed from $\psi_4$ and $h$ overestimate and underestimate, respectively, the values of eccentricity with respect to the ones from the dynamics. This is confirmed from the values for the initial eccentricity one obtains from the orbital frequency of the strain, $\psi_4$ and the dynamics, $e_0^h=0.55\pm 0.01$, $e_0^{\psi_4}=0.84\pm 0.03$, and  $e_0^{\text{dyn}}=0.65\pm 0.01$, respectively. These three values of eccentricity are measured at the same initial time,  $t=600M$. These results lead to the conclusion that the eccentricity measured from the frequency of the $(2,2)$ mode  is higher for $\psi_4$ than for $h$. This can be understood using the fact that $h \approx \int \int \psi_4 dt' dt $; therefore, $h$ is a smoother function than $\psi_4$. As shown in Fig. \ref{fig:OmegaComparisons}, this is not a particular result of our procedure to measure the eccentricity but a general fact which can be reproduced by any method to measure the eccentricity based on the oscillations of the frequency of the $(2,2)$ mode.  We have  decided to show the orbital frequency from the PN dynamics, as it contains more cycles and eases the visualization of the effect, but the same effect can be obtained with the orbital frequency from the BH motion of a NR simulation. Moreover, we remark that these differences were also noted in  \cite{Purrer:2012wy}, where Puerrer et al. explicitly computed the factor between the eccentricity estimator calculated from the gravitational wave frequency of $h$ and $\psi_4$ in the low eccentric limit. Thus, one expects to see these discrepancies even augmented as the eccentricity increases, as is the case for the waveforms studied in this article.
We also note that we choose not to integrate backward too long in the past of the binary due to the inaccuracy of the eccentric PN fluxes, which makes the solutions inaccurate for extremely high eccentricities and the inaccuracy of the PN expressions for the waveform, which also become more and more inaccurate for high eccentricities.

\begin{figure}[!]
\centering
\captionsetup{justification=centering}
\includegraphics[scale=0.45]{omegaComparisonsBO.pdf}
\caption{ Time evolution of the orbital frequency $M \omega_{orb}$ computed from the phase of $r \psi^{2,2}_4$, the orbital frequency computed from the dynamics, $\omega=|\frac{\bm{v}\times \bm{r}}{r^2}|$, and from the phase of the strain $h^{2,2}$.}
\label{fig:OmegaComparisons}
\end{figure}

\section{Parameter estimation with eccentric signals} \label{sec:PE}

In this section we employ the waveforms introduced in Secs. \ref{sec:NRcatalog} and \ref{sec:PNEcc} for data analysis studies. First, we analyze the impact of the eccentricity when computing overlaps against quasicircular models. Second, we perform parameter estimation studies injecting eccentric NR and hybrid waveforms into detector noise and compute parameter biases   using three different IMR  quasicircular models available in the LIGO libraries, \texttt{LALSUITE} \cite{lalsuite}.

\subsection{Match calculation}
A generic black hole binary evolving in a quasielliptical orbit is described by 17 parameters.  The intrinsic parameters are the individual masses of the binary $m_1$, $m_2$, the six components of the two spin vectors $\vec{S}_1$ and $\vec{S}_2$, the orbital eccentricity $e$ and the argument of the periapsis $\Omega$. The extrinsic parameters describing the sky position of the binary with respect to the detector are the distance from the detector to the source $r$, the coalescence time $t_c$, the inclination $\iota$, the azimuthal angle $\varphi$, the right ascension ($\phi$), the declination ($\theta$), and the polarization angle  ($\psi$). These parameters together describe the strain induced in a detector from a passing gravitational wave \cite{PhysRevD.47.2198};
\begin{equation}
\begin{split}
h(t,\zeta,\Theta)& =  \left[ F_+ (\theta,\phi,\psi)h_+(t-t_c;\iota,\varphi, \zeta)  \right. \\
& \left.  +  F_\times (\theta,\phi,\psi)h_\times (t-t_c;\iota,\varphi, \zeta)  \right].
\end{split}
\label{eq:eq13}
\end{equation}
 Where $\Theta=\{ t_c, r, \theta, \phi,\iota, \varphi, \psi  \}$ is the set of extrinsic parameters, $\zeta=\{m_1, m_2, \vec{S}_1, \vec{S}_2, e, \Omega \}$ are the intrinsic parameters and $F_+$, $F_\times$ are the antenna patterns functions defined in \cite{PhysRevD.47.2198} . The detector response is written in terms of the waveform polarizations $(h_+,h_\times)$ which combine to define the complex waveform strain
\begin{equation}
h(t)=h_+ - i h_\times = \sum_{l=2}^{\infty}\sum_{m=-l}^{l} Y^{-2}_{lm}(\iota, \varphi ) h_{lm} (t-t_c;\zeta),
\label{eq:eq14}
\end{equation}
 where $Y^{-2}_{lm}(\iota, \varphi )$ are spin-weighted -2 spherical harmonics, and where $h_{lm}$ refers to the $(l,m)$ waveform mode.
The comparison between two waveforms is usually quantified by an overlap integral, which is a noise-weighted inner product between  signals \cite{Jaranowski2012}, and which can be maximized over subsets or all of the parameters of the signal. Given a real-valued detector response, the inner product between the signal $h_{resp}^S(t)$ and the model $h_{resp}^M(t)$ is defined as
\begin{equation}
	\braket{h^S_{resp}|h^M_{resp}}= 4 \text{Re} \int^{+ \infty}_{0} \frac{\tilde{h}^S_{resp}(f) \tilde{h}^{M*}_{resp}(f)}{S_n(|f|)}df,
\label{eq:eq15}
\end{equation}
where $\tilde{h}$ denotes the Fourier transform of $h$, $h^*$ denotes the complex conjugate of $h$, and $S_n(|f|)$ is the one sided noise power spectral density (PSD) of the detector.

 The normalized match optimized over a relative time shift and the initial orbital phase can be written as
\begin{equation}
\mathcal{M}(\iota_S,\varphi_{0_S}) = \max_{t_c, \varphi_{0_S}} \left[ 	\frac{\braket{h^S_{resp}|h^M_{resp}}}{\sqrt{\braket{h^S_{resp}|h^S_{resp}}\braket{h^M_{resp}|h^M_{resp}}}}\right] .
\label{eq:eq16}
\end{equation}
The match is close to 1 when the model is able to faithfully reproduce the signal, while values of the match close to 0 indicate large disagreement between the two waveforms.
 In Eq. \eqref{eq:eq16} the match is computed for given values of the angles $(\iota_S,\varphi_{0_S})$ of the signal and maximizing over phase and time shifts. We will take only the $h_{22}$ mode of the eccentric hybrids and a quasicircular (QC) waveform model and compute single mode mismatches maximized over a time shift $t_0$ and a phase offset $\phi_0$  as
\begin{equation}
\mathcal{M}\mathcal{M} = \max_{t_0, \phi_{0}} \left[ 	\frac{\braket{h^{\text{hyb}}_{22}|h^{\text{QC}}_{22}}}{\sqrt{\braket{h^{\text{hyb}}_{22}|h^{\text{hyb}}_{22}}\braket{h^{\text{QC}}_{22}|h^{\text{QC}}_{22}}}}\right] .
\label{eq:eq116}
\end{equation}
To simplify the comparisons, we introduce the mismatch,	$1- \mathcal{M}\mathcal{M}$. Values of the mismatch close to zero indicate good agreement between the signal and the model, while the higher the mismatch the larger the difference between both waveforms, indicating that the model is not able to accurately represent the signal.

Having set the notation for the calculation of the mismatch, we compute the mismatch between the eccentric $(2,2)$ mode hybrids computed in Sec. \ref{sec:PNEcc} and the quasicircular model PhenomX \cite{phenX, phenXHM}, which is an upgrade of the aligned-spin PhenomD model \cite{PhysRevD.93.044006,PhysRevD.93.044007}, with calibration to a larger NR dataset and also to extreme-mass ratio waveforms. We employ the Advanced LIGO' `zero detuned high power' PSD \cite{LIGOPSD} to compute the overlap in Eq. \eqref{eq:eq15}. The integral of Eq. \eqref{eq:eq15} is evaluated within a frequency range of $20$-$2000$ Hz.  The nonmonotonic behavior of the GW frequency of eccentric systems complicates the determination of the frequency range of a signal in the detector band. The ideal case would be one in which the initial apastron and periastron frequencies are below $20$Hz. This would mean that the whole waveform starts before the cutoff frequencies of the detectors and one observes the complete eccentric inspiral of the binary. Another possibility is that both frequencies are above $20$ Hz, then the signal is very short and much of the inspiral waveform is lost. Finally, it is also possible that during some part of the waveform the periastron frequencies are above $20$ Hz and apastron frequencies are below $20$ Hz. The latter is typically the case for our hybrid waveforms.

\begin{figure}[hbt!]
\centering
\captionsetup{justification=centering}
\includegraphics[scale=0.5]{mmMatch22.pdf}
\caption{ Mismatches for the $(l,m)=(2,2)$ mode between the eccentric hybrid waveforms corresponding to the cases presented in Table \ref{tab:tabNR3} and the quasicircular PhenomX waveform model as a function of the total mass of the system. The green, blue, black, and red lines correspond to eccentric PN-NR hybrid waveforms with initial eccentricities $e \leq 0.3$, $0.3 \leq e \leq 0.5$, and $0.5 \leq e \leq 0.84$, respectively. }
\label{fig:mmMatch22}
\end{figure}

In Fig. \ref{fig:mmMatch22} we show the single mode mismatches between the eccentric hybrids and PhenomX for a range of total mass for a system where $M_T \in\{20,200\} M_\odot$. As expected, for larger total masses of the system most of the waveform in the frequency band of the detector is in the merger and ringdown parts and the mismatches are even below the $3 \%$ threshold. This is consistent with the results obtained in Sec. \ref{sec:NRcatalogB}, which show the agreement in the final state between the eccentric simulations and the quasicircular fits. However, the lower the total mass the higher the mismatch, this is due to the fact that at low frequencies there is more inspiral part of the waveform in the frequency band, and therefore the inability of the quasicircular model to resemble the eccentric inspiral becomes notorious. One can also appreciate in Fig. \ref{fig:mmMatch22} that, generally, the higher the initial eccentricity the higher the mismatch for the whole mass range.

\subsection{Eccentric injections into detector noise}
In this section we show some applications of the eccentric waveform dataset to parameter estimation. We inject eccentric hybrids into Gaussian noise realization recolored to match the spectral density of the Advanced LIGO and Virgo detectors at design sensitivity. All simulations with the same injected signal are performed with the same noise realization.  We do not address the challenge of detecting eccentric signals and instead assume that the signal has been detected by standard CBC search pipelines \cite{PhysRevD.81.024007}. We perform parameter estimation using the python-based \texttt{BILBY} code \cite{Ashton:2018jfp}. Of the numerous stochastic samplers implemented in \texttt{BILBY}, we employ the nested sampler \texttt{CPNEST} \cite{JohnVeitchCpnest} and use waveform approximants implemented in \texttt{LALSUITE} \cite{lalsuite} as the model templates. The eccentric hybrids make use of the numerical relativity injection infrastructure \cite{Galley:2016mvy,Schmidt:2017btt}.

State-of-the-art non-spinning eccentric IMR models \cite{PhysRevD.97.024031,PhysRevD.98.044015,PhysRevD.96.044028,PhysRevD.96.104048} have not yet been implemented in \texttt{LALSUITE}. The only eccentric waveform models in \texttt{LALSUITE} are inspiral nonspinning frequency domain approximants \cite{PhysRevLett.112.101101,PhysRevD.90.084016,PhysRevD.93.064031}. We have decided not to use such inspiral waveform models to avoid bias induced by the sharp cutoff at the end of the waveform \cite{Mandel:2014tca}. For a study of the eccentricity measurement using such inspiral approximants, see \cite{PhysRevD.98.083028,Moore:2019vjj}. We restrict to IMR quasicircular approximants and perform parameter estimation analysis on the injected eccentric signals by sampling in the 15 parameters of a quasicircular black hole binary.

\begin{table}[h!]
\begin{center}
 \def1.1 {1.3 }
\begin{tabular}{c c c c c c c c c c c c c c c c }
\hline
\hline
ID &  & Simulation &  &    q  & &${\chi}_{1,z}$ & & $\chi_{2,z}$ &    &$D/M$ & &  $e_\omega \pm \delta e_\omega$  \\
\hline
61&   & \texttt{SXS:BBH:1355}&  &     1& & $0$ & & $0$ & &   12.97      & &   $0.090 \pm 0.003$ \\
62&   & \texttt{SXS:BBH:1359}&  &    1&  &$0$ &  &$0$  & &    15.73 	  & &   $0.146 \pm 0.003$ \\
63&   & \texttt{SXS:BBH:1361}&  &   1&  &$0$ & & $0$   & &    16.69 	  & &   $0.209 \pm 0.003$ \\
\hline
\hline
\end{tabular}
\end{center}
\caption{Summary of the injected  NR simulations. The first column denotes the identifier of the simulation, and the second column indicates the name of the simulation as presented in \cite{PhysRevD.98.044015}. The next columns show the mass ratio, the z component of the dimensionless spin vectors, the initial orbital separation and the initial orbital eccentricity measured using the procedure detailed in Sec. \ref{sec:NRcatalogC}.}
\label{tab:tabNRPE}
\end{table}

We inject three  NR equal mass nonspinning simulations described in Table \ref{tab:tabNRPE} into a network of gravitational wave detectors composed of the LIGO-Hanford, LIGO-Livingston \cite{TheLIGOScientific:2014jea}, and Virgo interferometers \cite{TheVirgo:2014hva}, each operating at design sensitivity. We set a reference frequency of $f_{\text{ref}}=20 \text{ Hz}$, where the waveforms start. Some injected parameters are displayed in Table \ref{tab:tabNRbias}, while the declination is $\delta=-1.21$ rad, the right ascension $\alpha=1.37$ rad, and the coalescence phase $\phi=0$ rad. From these simulations, the $\{(l,m\}=\{(2, \pm 2), (3,\pm 2), (4, \pm 4), (5, \pm 4), (6, \pm 6)\}$ modes are used. We do not inject odd m modes because they are zero by symmetry. For the injected signal we choose the luminosity distance $D_L=430$ Mpc, which is similar to the first detection of a gravitational wave signal, GW150914 \cite{LIGOScientific:2018mvr}, which produces a high network signal-to-noise ratio (SNR) as shown in Tables \ref{tab:tabNRbias} and \ref{tab:tabHybridbias}.

We employ a uniform-in-volume prior on the luminosity distance, $p(D_L|H)\propto D_L^2 $,  between 50 and 1500 Mpc.  The inclination and polarization angles  both  have uniform priors between $(0,\pi)$. We use the standard priors for the extrinsic variables, as in Table I of \cite{Ashton:2018jfp}. Instead of sampling in the component masses, we sample in   mass ratio q and   chirp mass $\mathcal{M}_c$, with ranges $(0.05,1)$ and $(15, 60) M_\odot$, respectively. The spin priors are set differently according to the approximant. If the approximant is nonprecessing, we set the option of \texttt{aligned\_spin=True} in the \texttt{BBHPriorDict} function of \texttt{BILBY}, which samples in the dimensionless spin z components between -0.8 and +0.8. For precessing approximants, we sample in the tilt angles $(\theta_1,\theta_2)$, the angle between the spin vectors $\phi_{12}$, the angle between \textbf{J} and \textbf{L} $\phi_{JL}$, and the dimensionless spin magnitudes $(a_1,a_2)$. The priors for  $a_1,a_2,\theta_1,\theta_2,\phi_{JL}$, and $\phi_{12}$ are the same as in Table I of \cite{Ashton:2018jfp}.
We also define a uniform prior for the coalescence time of 2s centered at the injection time.

We take three quasicircular models as approximants: 1) IMRPhenomD \cite{PhysRevD.93.044006,PhysRevD.93.044007}, a nonprecessing model with only the $(2,\pm 2)$ modes, 2)  IMRPhenomHM \cite{london2018},  nonprecessing model including higher order modes, and 3) IMRPhenomPv2 \cite{phenomp},  an effective precessing model.

We plot the posterior probability distribution for the chirp mass, mass ratio, effective spin parameter and luminosity distance for the PhenomD approximant in Fig. \ref{fig:posteriorsNRPhD} with $90 \%$ credible intervals specified by the dashed vertical lines and the injected values by the  thick vertical magenta lines. As a control case, we also show in Fig. \ref{fig:posteriorsNRPhD} the posterior distribution of an equal mass nonspinning zero-eccentricity injection performed using the hybridized surrogate model NRHybSur3dq8 \cite{PhysRevD.99.064045} with the same injected parameters  as in Table \ref{tab:tabNRbias} and recovered with the PhenomD model.  The NRHybSur3dq8 injected waveform contains all higher order modes up to $l=4$, which in this case seems to cause the small bias one observes in the luminosity distance when recovering with the PhenomD model, which contains only the $(2,|2|)$ modes. For the rest of the parameters, like the mass ratio, the chirp mass and the effective spin parameter, we obtain results consistent with the accuracy of the PhenomD model for parameter estimation of injected signals as shown in \cite{phenX}.

The posterior distributions for the rest of the approximants are shown in Fig. \ref{fig:posteriorsNR}.  The same information is summarized in Fig. \ref{fig:errorbarsPE}, where the median and the error bars corresponding to the $90 \%$ credible intervals of the posterior distribution are shown as a function of the initial eccentricity. Note that the bars corresponding to the same initial eccentricity but different approximants have been separated by a small amount to ease the visualization of the results. For the lowest initial eccentricity, $e_0=0.09$, the results for the four quantities are pretty different. The chirp mass and the effective spin parameter produce similar distributions for the three approximants, while for the mass ratio and the luminosity distance, PhenomHM distributions are closer to the injected values than PhenomD and PhenomPv2.

Furthermore, for $e_0=0.14$ and $e_0=0.2$ we observe increasingly poor agreement with the injected values, except for the mass ratio where the lowest initial eccentricity signal produces wider distributions than those with higher initial eccentricities. This can also be checked in Table \ref{tab:tabNRbias}, where the recovered parameters, median values and $90 \%$ credible intervals, are compared with the injected values. Regarding the effective spin parameter and the chirp mass, the increase of initial eccentricity in the injected signal shifts the posteriors for the three quasicircular models, while for the mass ratio the increase of initial eccentricity reduces the bias on the measurement of the mass ratio, probably as a consequence of the shift in the chirp mass distribution, as displayed in Fig. \ref{fig:q1PhDcontourPlot}  for the all injections recovered with the PhenomD model.

One observes also that PhenomHM recovers the injected parameters more effectively than PhenomD and PhenomPv2. For the luminosity distance the probability densities tend to flatten and be closer to the prior distributions for high initial eccentricities. One notes again that PhenomHM has less parameter bias than PhenomD and PhenomPv2. Injected values of the sky position like the right ascension $\alpha=1.375$ rad and  $\delta=-1.21$ rad are well recovered for all nine runs, probably due to the expensive parameter estimation (PE) settings described in Appendix \ref{sec:AppendixC} : $\alpha=1.37^{+0.01}_{-0.01}$ rad and $\delta=-1.21^{+0.01}_{-0.01}$ rad.

Furthermore, we have computed the recovered matched-filter SNR for the detector network, $\rho_{\text{Match}}$, for each simulation. This quantity,  $\rho_{\text{Match}}$, is computed calculating the matched filter between the detector data with the eccentric signal injected and the waveform of the approximant waveform model with the parameters corresponding to the highest log-likelihood value of the posterior distribution. The results of such a calculation are shown in Table \ref{tab:tabNRPE}. As expected, we observe that the zero-eccentricity injection recovers much more SNR than the eccentric injections, with decreasing values of the recovered SNR with increasing eccentricity.

Additionally, we display the values of the log Bayes factor for each simulation. The Bayes factor is computed here as the ratio between the signal and null evidence [see Eq. (13) of \cite{Thrane_2019}]. One can observe that both the recovered matched-filter SNR and the log Bayes factor decrease more the higher the initial eccentricity of the injected signal is. The matched-filter SNR produces similar values between models for simulations with the same initial eccentricity. However, the log Bayes factor tends to be slightly higher for the aligned-spin waveform models, PhenomD and PhenomHM, for the lowest initial eccentric injected signal, while for higher initial eccentricities the precessing model IMRPhenomPv2 shows slightly greater log Bayes factors than the aligned-spin ones. The highest log Bayes factor is obtained for the zero-eccentricity injection.

\begin{table*}
\begin{center}
 \def1.1 {1.4}
\begin{tabular}{  c  c  c   c   c  c  c  c   c  c c c c }
\hline
\hline
$e_0$& Model &  $m_1/M_\odot$&  $m_2/M_\odot$& $\mathcal{M}_c/M_\odot$&  $q$ &     $D_L/$Mpc&  $\chi_{\text{eff}}$ &$\psi$ (rad) &$\iota$ (rad)& $\rho_{\text{Match}}$  & $\log \mathcal{B}$ \\
\hline
\multirow{3}{*}{$0.09$} & PhenomD   &  $35.06^{+2.55}_{-1.92}$&  $31.29^{+1.10}_{-1.33}$ &   $28.40^{+0.17}_{-0.17}$&  $0.87_{-0.12}^{+0.10}$&   $384^{+49}_{-82}$ &  $0.00^{+0.02}_{-0.02}$  & $1.60^{+1.28}_{-1.34}$  &$0.54^{+0.32}_{-0.32}$    & 89.40  & 3463.79 \\
   & PhenomHM   &  $34.05^{+2.14}_{-1.16}$&  $31.79^{+0.69}_{-1.17}$ &   $28.38^{+0.16}_{-0.16}$&  $0.92_{-0.11}^{+0.07}$&   $429^{+16}_{-33}$ &  $-0.01^{+0.02}_{-0.02}$ & $2.01^{+0.97}_{-1.82}$ & $0.28^{+0.21}_{-0.17}$    & 89.28 & 3463.78 \\
   & PhenomPv2   &  $35.26^{+2.97}_{-2.06}$&  $31.28^{+1.18}_{-1.53}$ &   $28.44^{+0.21}_{-0.18}$&  $0.86^{+0.11}_{-0.13}$&   $412^{+24}_{-66}$ &  $0.00^{+0.02}_{-0.02}$   & $1.65^{+1.22}_{-1.32}$& $0.39^{+0.32}_{-0.22}$   & 89.19  & 3459.54\\
 \hline

\multirow{3}{*}{$0.14$} & PhenomD   &  $34.03^{+1.34}_{-0.72}$&  $32.63^{+0.44}_{-0.73}$ &   $28.86^{+0.15}_{-0.15}$&  $0.95^{+0.04}_{-0.07}$&   $407^{+53}_{-84}$ &  $0.02^{+0.02}_{-0.02}$ & $1.58^{+1.21}_{-1.23}$ & $0.54^{+0.32}_{-0.32}$   & 84.87   & 3288.25 \\
   & PhenomHM   &  $33.76^{+0.96}_{-0.54}$&  $32.73^{+0.35}_{-0.56}$ &   $28.82^{+0.16}_{-0.14}$&  $0.96^{+0.03}_{-0.05}$&   $408^{+46}_{-52}$ &  $0.02^{+0.02}_{-0.02}$  & $1.91^{+0.46}_{-0.58}$ & $0.54^{+0.19}_{-0.25}$     & 84.74  & 3283.61 \\
   & PhenomPv2   &  $34.22^{+1.48}_{-0.89}$&  $32.54^{+0.54}_{-0.82}$ &   $28.87^{+0.19}_{-0.21}$&  $0.94^{+0.05}_{-0.08}$&   $389^{+33}_{-60}$ &  $0.01^{+0.02}_{-0.03}$  & $1.70^{+1.09}_{-1.07}$  &$0.64^{+0.25}_{-0.18}$ & 85.08  & 3302.37 \\

 \hline

\multirow{3}{*}{$0.2$} & PhenomD   &  $35.65^{+1.52}_{-0.85}$&  $34.01^{+0.51}_{-0.82}$ &   $30.13^{+0.16}_{-0.16}$&  $0.94^{+0.05}_{-0.07}$&   $420^{+72}_{-109}$ &  $0.07^{+0.02}_{-0.02}$ & $1.57^{+1.36}_{-1.18}$  & $0.61^{+0.41}_{-0.37}$   &   81.88  & 3102.70 \\
   & PhenomHM    &  $35.47^{+1.36}_{-0.78}$&  $33.97^{+0.46}_{-0.72}$ &   $30.06^{+0.16}_{-0.15}$&  $0.95^{+0.04}_{-0.07}$&   $438^{+43}_{-47}$ &  $0.06^{+0.02}_{-0.02}$ & $0.42^{+0.90}_{-0.29}$ & $0.54^{+0.16}_{-0.20}$    & 81.97  &  3101.79 \\
   & PhenomPv2   &  $37.13^{+2.11}_{-1.76}$&  $33.20^{+1.00}_{-1.12}$ &   $30.12^{+0.21}_{-0.22}$&  $0.87^{0.09}_{-0.09}$&   $414^{+41}_{-69}$ &  $0.06^{+0.02}_{-0.02}$  & $1.62^{+0.98}_{-1.29}$ & $0.66^{+0.25}_{-0.18}$ &  82.05  & 3112.97  \\

 \hline
0 & PhenomD &  $34.07^{+2.05}_{-1.29}$&  $30.91^{+1.24}_{-1.77}$ &   $28.24^{+0.16}_{-0.16}$&  $0.91^{+0.07}_{-0.10}$&   $375^{+48}_{-75}$ &  $0.0^{+0.02}_{-0.02}$ & $1.58^{+1.28}_{-1.24}$  & $0.53^{+0.33}_{-0.32}$   &  173.16   &  3632.19 \\
 \hline
& Injected &  $32.5$&  $32.5$& $28.29$&  $1$ &      430 &  $0$ &$0.33$  & $0.3$ & \\
\hline
\hline
 \end{tabular}

\end{center}
\caption{Black hole binary recovered parameters for the three NR simulations from Table \ref{tab:tabNRPE}. The last row corresponds to the injected parameters. In the penultimate row we show the recovered parameters of the zero-eccentricity injection performed with the NRHybSur3dq8 model. The first column describes the initial eccentricity of the injected signal. Then we specify the approximant, the component masses, the chirp mass, the mass ratio, the luminosity distance, the effective spin parameter, the polarization angle, the inclination, the recovered matched-filter SNR for the detector network, and the log of the Bayes factor.}
\label{tab:tabNRbias}
\end{table*}

\begin{figure*}[!]

\noindent\begin{minipage}{\textwidth}

\noindent\begin{minipage}[h]{.45\textwidth}
\includegraphics[scale=0.42]{q1NSMcPhDe0.pdf}
\end{minipage}
\hspace{1cm}
\begin{minipage}[h]{.45\textwidth}
\includegraphics[scale=0.42]{q1NSmassratioPhDe0.pdf}
\end{minipage}

\end{minipage}

\noindent\begin{minipage}{\textwidth}

\noindent\begin{minipage}[h]{.45\textwidth}
\includegraphics[scale=0.42]{q1NSDLPhDe0.pdf}
\end{minipage}
\hspace{1cm}
\begin{minipage}[h]{.45\textwidth}
\includegraphics[scale=0.42]{q1NSchieffPhDe0.pdf}
\end{minipage}

\end{minipage}

\caption{Posterior probability distributions for the injected NR simulations of Table \ref{tab:tabNRPE} and a zero-eccentricity injection using the NRHybSur3dq8 model. The vertical dashed lines correspond to $90 \%$ credible regions. The magenta thick vertical line represents the injected value. The green, black, blue and red curves represent distributions sampled using the IMRPhenomD approximant with injected initial eccentricities $e_0= 0.0 ,0.09,0.14,0.2$, respectively.  }

\label{fig:posteriorsNRPhD}
\end{figure*}

\begin{figure*}[!]

\noindent\begin{minipage}{\textwidth}

\noindent\begin{minipage}[h]{.45\textwidth}
\includegraphics[scale=0.42]{q1NSerrbarMc.pdf}
\end{minipage}
\hspace{1cm}
\begin{minipage}[h]{.45\textwidth}
\includegraphics[scale=0.42]{q1NSerrbarq.pdf}
\end{minipage}

\end{minipage}

\noindent\begin{minipage}{\textwidth}

\noindent\begin{minipage}[h]{.45\textwidth}
\includegraphics[scale=0.42]{q1NSerrbarDL.pdf}
\end{minipage}
\hspace{1cm}
\begin{minipage}[h]{.45\textwidth}
\includegraphics[scale=0.42]{q1NSerrbarchieff.pdf}
\end{minipage}

\end{minipage}

\caption{Median values and error bars corresponding to $90 \%$ credible regions  of the posterior probabibility distributions for the injected NR simulations of Table \ref{tab:tabNRPE}. The vertical magenta line represents the injected value. The black, blue and red segments represent the median values and errors bars of the distributions sampled using the IMRPhenomD, IMRPhenomHM and IMRPhenomPv2 approximants, respectively. The cases are represented for three initial eccentricities of the injected signal, $e_0=0.09,0.14,0.2$. To ease the visualization of the horizontal bars, cases with the same initial eccentricity and run with different approximants have been separated a $\Delta e =0.003$.}

\label{fig:errorbarsPE}
\end{figure*}

\begin{figure*}[!]

\noindent\begin{minipage}{\textwidth}

\noindent\begin{minipage}[h]{.45\textwidth}
\includegraphics[scale=0.42]{hybridMcbo.pdf}
\end{minipage}
\hspace{1cm}
\begin{minipage}[h]{.45\textwidth}
\includegraphics[scale=0.42]{hybridmassratiobo.pdf}
\end{minipage}

\end{minipage}

\noindent\begin{minipage}{\textwidth}

\noindent\begin{minipage}[h]{.45\textwidth}
\includegraphics[scale=0.42]{hybridDLbo.pdf}
\end{minipage}
\hspace{1cm}
\begin{minipage}[h]{.45\textwidth}
\includegraphics[scale=0.42]{hybridchieffbo.pdf}
\end{minipage}

\end{minipage}

\caption{ Posterior probability distributions for the injected spinning eccentric hybrid waveform, with initial eccentricity   $e_0=0.420\pm 0.006$. The vertical dashed lines correspond to $90 \%$ credible regions. The thick vertical magenta line represents the injected value. The black, blue and red curves represent distributions sampled using the IMRPhenomD, IMRPhenomHM and IMRPhenomPv2 approximants, respectively.}
\label{fig:posteriorsHybrid}
\end{figure*}

We repeat the same procedure while injecting a hybrid waveform, including only the $(l,m)=(2,\pm 2)$ modes, of an eccentric spinning waveform with ID 8 of Table \ref{tab:tabNR3}. This is an equal mass with the z component of the dimensionless spin vectors $\chi_{1z}=\chi_{2z}=-0.25$ and initial eccentricity $e_0=0.420\pm 0.006$. The posterior distribution for the chirp mass, mass ratio, luminosity distance and $\chi_{\text{eff}}$ are shown in Fig. \ref{fig:posteriorsHybrid} for IMRPhenomD, IMRPhenomHM,  and IMRPhenomPv2 as waveform models. In this case, the parameter biases are much higher than in the previous injection study mainly due to the fact that the injected signal has a much higher initial eccentricity.

The values of the recovered parameters as well as the injected values are shown in Table \ref{tab:tabHybridbias}. The injected values of the sky position like the right ascension $\alpha=1.375$ rad and  $\delta=-1.21$ rad are again well recovered parameters for the three runs $\alpha=1.37^{+0.01}_{-0.01}$ rad and $\delta=-1.21^{+0.01}_{-0.01}$ rad. The bias in the chirp mass is $\sim 4 M_\odot$ for the three models. Here one again observes  the correlation between chirp mass and mass ratio. The shift in chirp mass posteriors with respect to the injected value translates into a better determined mass ratio distribution, which is clearly the case for PhenomHM, which performs unexpectedly well in recovering the mass ratio parameter, while PhenomD and PhenomPv2 show much wider distributions and much larger credible intervals.

The posteriors of the luminosity distance show also large error bars for the three models, where PhenomHM again reduces the bias with respect to PhenomD and PhenomPv2. The recovered effective spin parameter is completely off with respect to the injected value for the three approximants. The recovered $\chi_{\text{eff}}$ is positive, while the injected one is negative. The bias in the effective spin parameter is  approximately $-0.3$ for the three models, indicating the inability of the quasicircular models to estimate the spin parameter of highly eccentric spinning binaries with quasicircular models. Regarding the  recovered matched-filter SNR and the log Bayes factor displayed in Table \ref{tab:tabHybridbias}, one can observe that while the SNR provides comparable values among models, the values of the log Bayes factor indicate that PhenomPv2 fits  the data scarcely better than PhenomHM and PhenomD.

\begin{table*}
\begin{center}
 \def1.1 {1.4 }
\begin{tabular}{  c   c   c   c   c  c  c  c   c  c c c c }
\hline
\hline
$e_0$& Model &  $m_1/M_\odot$&  $m_2/M_\odot$& $\mathcal{M}_c/M_\odot$&  $q$ &     $D_L/$Mpc&  $\chi_{\text{eff}}$ &$\psi$ (rad) &$\iota$ (rad) & $\rho_{\text{Match}}$ & $\log \mathcal{B}$\\
\hline
\multirow{3}{*}{$0.42$} & PhenomD   &  $37.52^{+1.30}_{-0.76}$&  $36.04^{+0.49}_{-0.73}$ &   $31.86^{+0.19}_{-0.2}$&  $0.95^{+0.04}_{-0.06}$&   $474^{+62}_{-101}$ &  $0.06^{+0.02}_{-0.02}$ &$2.60^{+0.31}_{-0.33}$  & $1.54^{+1.22}_{-1.19}$   &  82.68  & 2895.91 \\
   & PhenomHM   &  $37.23^{+0.95}_{-0.37}$&  $36.62^{+0.30}_{-0.81}$ &   $32.07^{+0.18}_{-0.23}$&  $0.98^{+0.02}_{-0.06}$&   $384^{+54}_{-45}$ &  $0.08^{+0.02}_{-0.02}$ & $2.28^{+0.18}_{-0.16}$  & $1.04^{+1.12}_{-0.26}$  &  82.54 & 2894.17 \\
   & PhenomPv2   &  $39.15^{+2.08}_{-1.62}$&  $35.20^{+0.84}_{-1.06}$ &   $31.87^{+0.23}_{-0.26}$&  $0.88^{+0.07}_{-0.08}$&   $413^{+77}_{-110}$ &  $0.05^{+0.03}_{-0.03}$ &$2.33^{+0.33}_{-0.44}$  & $1.46^{+1.36}_{-0.45}$ &  82.62  & 2910.28 \\
 \hline
& Injected &  $32.5$&  $32.5$& $28.29$&  $1$ &  $430$ &  $-0.25$ &$0.33$  & $0.3$ & & \\
\hline
\hline
 \end{tabular}

\end{center}
\caption{Black hole binary recovered parameters for the spinning hybrid waveform from Fig. 	\ref{fig:posteriorsHybrid}. The last row corresponds to the injected parameters. The first column describes the initial eccentricity of the injected signal. Then we specify the approximant, the component masses, the chirp mass, the mass ratio,  the luminosity distance, the effective spin parameter, the polarization angle, the inclination, the recovered matched-filter SNR for the detector network, and the log of the Bayes factor.}
\label{tab:tabHybridbias}
\end{table*}

This section shows examples of the kind of study that one is able to perform with the current eccentric waveform dataset.  We have shown the limitations of the current IMR quasicircular to estimate the parameters of moderately eccentric waveforms including a moderately spinning case. For the cases studied in this section, we have found that although the use of quasicircular models to estimate parameters of eccentric signals leads to inevitable biases, aligned-spin quasicircular models with higher order modes leverage the impact of these biases for the mass ratio and the luminosity distance when compared to aligned-spin models with only the $(2,\pm 2)$ modes or precessing models. Because of the computational cost of the PE runs and the amount of eccentric waveforms available, we leave for future work a detailed study of the whole dataset using not only quasicircular models but also eccentric waveform approximants.

\section{Summary and conclusions} \label{sec:summary}
In this paper we have presented the first parameter study of numerical relativity simulations of eccentric spinning black hole binaries.
We have presented a simple procedure to set up the initial parameters of eccentric simulations. The higher the initial eccentricity of the simulation, the longer the initial separation has to be  to avoid the immediate plunge of the binary due to the strong interactions at the periastron. This increases the computational cost of the simulations of Table \ref{tab:tabNR3} with $e_0 \sim 0.4$, which is roughly double the one with $e_0 \sim 0.2$, as can be observed in their merger times. Additionally, longer initial separations produce long enough waveforms, which allows one to avoid the breakdown of the post-Newtonian approximation and ease the posterior construction of PN-NR hybrid waveforms. As part of the postprocessing step, we have computed the final mass and final spin of the 60 new simulations presented in Table \ref{tab:tabNR3}. We have compared the final mass and final spin of those simulations with quasicircular NR fits   \cite{PhysRevD.95.064024} and found that relative differences are as high as $1 \%$, which is completely consistent with the inaccuracies of the fitting formulas and gauge transient in the apparent horizon quantities. Therefore, we have extended previous work \cite{PhysRevD.77.081502} on the circularization of eccentric nonspinning numerical relativity simulations to the eccentric spinning case. Note that the eccentricities of the simulations presented in this communication have more moderate values than the ones presented in \cite{PhysRevD.77.081502}, although ours are much longer and include spins.

A crucial part of this  work has been to extend the low eccentric procedure to measure the eccentricity in NR \cite{PhysRevD.99.023003} to the arbitrary high eccentric limit. We have shown that the eccentricity estimator used in \cite{PhysRevD.99.023003} cannot be used for high eccentricities because it does not reduce to the Newtonian definition of the eccentricity. Additionally, its reliance on a noneccentric fit makes it numerically inaccurate, and it can produce eccentricity values higher than 1. As a consequence, we have decided to use another eccentricity estimator \cite{PhysRevD.66.101501} that is also constructed upon the orbital frequency and which does not rely on any noneccentric fit. This eccentricity estimator reduces to the Newtonian definition of eccentricity for arbitrarily high eccentricities. We have shown that with this eccentricity estimator we are able to robustly measure the eccentricity for the whole evolution, which will be a key result for generating a future eccentric waveform model.

We have then taken the NR waveforms and hybridized the $(2,2)$ mode with PN waveforms. The production of the eccentric PN waveforms has required to solve the point particle 3.5PN equations of motion in ADMTT coordinates  \cite{PhysRevD.99.023003} enhanced with the eccentric contribution to the energy flux from \cite{PhysRevD.80.124018}. The absence of complete generic PN expressions for the waveform modes has caused the inaccuracy of the PN waveforms to dominate the error in the hybridization procedure. The use of the instantaneous terms at 3PN order \cite{PhysRevD.91.084040} produces inaccurate waveforms due to the lack of low order tail terms, while the full 3PN expressions in  \cite{Boetzel:2019nfw} are restricted to the QK parametrization and rely on a certain decomposition of the dynamical variables which complicates their combination with the generic numerical solution of the equations of motion \eqref{eq:eq10}. Therefore, we have restricted to the use of the quadrupole formula with a correction procedure for the initial orbital separation. We have developed a procedure which corrects the initial orbital separation of the PN evolution code for a certain $\delta r$, such that it minimizes the difference in amplitude between the PN and NR $(2,2)$  waveform modes. We have shown that with that procedure we are left with relative errors in the amplitude and phase below $1 \%$ in the hybridization region. These errors in amplitude and phase are high compared to the quasicircular ones \cite{QCHybrids}, where the PN knowledge is wider. Therefore, we expect that in the future an improvement in  knowledge about post-Newtonian waveforms will allow us to construct more accurate hybrid waveforms, not only for the $(2,2)$ mode  but also for the higher order modes.

We have also compared the hybrid waveforms with quasicircular IMR waveform models. This has been done by first computing the mismatch of the eccentric hybrid dataset against the quasicircular nonprecessing  PhenomX model \cite{phenX,phenXHM}. We find that the mismatches become much higher than $3 \%$ for binaries with a total mass lower than $100 M_\odot$, while for total masses higher than $150 M_\odot$ the mismatch lowers below $3 \%$ due to the fact that most of the eccentric waveforms in the frequency band of the detector are in the merger-ringdown  parts, which as shown in Sec. \ref{sec:NRcatalogB}, due to circularization agrees really well with the quasicircular model.

Additionally, we have made a set of injections into Gaussian detector noise colored to match the LIGO and Virgo design detector sensitivities. We have studied the parameter biases on recovered parameters when using quasicircular models as approximants. We have used three different quasicircular models to recover the parameters and have shown that, although the use of quasicircular models leads to inevitable biases in parameters like the effective spin parameter or the chirp mass, where the biases are similar among the three models, others like the mass ratio and the luminosity distance present lower biases when using quasicircular aligned-spin models including higher order modes. Another important feature is the correlation between chirp mass and mass ratio, the better the measurement of the chirp mass the worse the determination of the mass ratio, and vice versa. This can be clearly observed in Figs. \ref{fig:posteriorsHybrid}  and \ref{fig:posteriorsNR}  where for initial eccentricities $0.09$ the chirp mass is well measured for the three models  but the mass ratio distributions are not. As the initial eccentricity increases, so does the shift in the chirp mass distribution, and  the mass ratio is generally better determined. In the case of the spinning eccentric hybrid, the high initial eccentricity produces clear biases in all quantities and, unexpectedly, PhenomHM recovers the injected value of the mass ratio well and performs the best for the luminosity distance. The study of this phenomenology for the different cases that we have available is ongoing and we leave for a future communication the extension of these results to the whole parameter space.

The work presented in this communication is a natural extension of \cite{PhysRevD.99.023003}. We have set up the current infrastructure  of our group for quasicircular waveform modeling to the eccentric case. As shown in this paper, we have developed new methods to produce a set of spinning eccentric hybrid waveforms which can actually be used for data analysis purposes. The next natural step is to use this hybrid dataset to produce a calibrated eccentric IMR waveform, which can be used for the detection and parameter estimation of eccentric blackhole binaries.

\section{Acknowledgements}  \label{sec:Acknowledgements}

We thank Frank Ohme for the useful comments about manuscript and Maria Haney for the valuable discussions. This work was supported by European Union FEDER funds, the Ministry of Science, Innovation and Universities, and  Spanish Agencia Estatal de Investigación Grants No. FPA2016-76821-P, FPA2017-90687-REDC, FPA2017-90566-REDC. FIS2016-81770-REDT, FPA2015-68783-REDT,  Spanish Ministry of Education, Culture and Sport grants  FPU15/03344, FPU15/01319, Vicepresid`encia i Conselleria d’Innovació, Recerca i Turisme, Conselleria d’Educació, i Universitats del Govern de les Illes Balears i Fons Social Europeu, Generalitat Valenciana (PROMETEO/2019/071),  EU COST Actions CA18108 , CA17137, CA16214, and CA16104, H2020-MSCA-IF-2016. Marie Skłodowska-Curie individual fellowships Proposal No. 751492. The authors thankfully acknowledge the computer resources at MareNostrum and the technical support provided by Barcelona Supercomputing Center (BSC) through Grants No. AECT-2019-2-0010, AECT-2019-1-0022, AECT-2018-3-0017, AECT-2018-2-0022, AECT-2018-1-0009, AECT-2017-3-0013, AECT-2017-2-0017, AECT-2017-1-0017, AECT-2016-3-0014, AECT2016-2-0009,  from the Red Española de Supercomputación (RES) and PRACE (Grant No. 2015133131). \texttt{BAM} and \texttt{ET} simulations were carried out on the BSC MareNostrum computer under PRACE and RES  allocations and on the FONER computer at the University of the Balearic Islands.