Bayesian inference of binary black holes with inspiral-merger-ringdown waveforms using two eccentric parameters

Orbital eccentricity is a crucial physical effect to unveil the origin of compact-object binaries detected by ground- and spaced-based gravitational-wave (GW) observatories. Here, we perform for the first time a Bayesian inference study of inspiral-merger-ringdown eccentric waveforms for binary black holes with non-precessing spins using two (instead of one) eccentric parameters: eccentricity and relativistic anomaly. We employ for our study the multipolar effective-one-body (EOB) waveform model SEOBNRv4EHM, and use initial conditions such that the eccentric parameters are specified at an orbit-averaged frequency. We show that this new parametrization of the initial conditions leads to a more efficient sampling of the parameter space. We also assess the impact of the relativistic-anomaly parameter by performing mock-signal injections, and we show that neglecting such a parameter can lead to significant biases in several binary parameters. We validate our model with mock-signal injections based on numerical-relativity waveforms, and we demonstrate the ability of the model to accurately recover the injected parameters. Finally, using standard stochastic samplers employed by the LIGO-Virgo-KAGRA Collaboration, we analyze a set of real GW signals observed by the LIGO-Virgo detectors during the first and third runs. We do not find clear evidence of eccentricity in the signals analyzed, more specifically we measure $e^{\text{GW150914}}_{\text{gw, 10Hz}}= 0.08^{+0.09}_{-0.06}$, $e^{\text{GW151226}}_{\text{gw, 20Hz}}= {0.04}^{+0.05}_{-0.04} $, and $e^{\text{GW190521}}_{\text{gw, 5.5Hz}}= 0.15^{+0.12}_{-0.12}$.

In this paper, we perform a Bayesian inference study using the waveform model SEOBNRv4EHM [94,95], which describes eccentric effects using two eccentric parameters: the initial eccentricity and the relativistic anomaly.SEOBNRv4EHM is built upon the quasi-circular SEOBNRv4HM model [96] for binary black holes (BBHs) with non-precessing spins and includes eccentric corrections up to 2PN order [94] in the (l, |m|) = (2, 2), (2, 1), (3,3), (4,4), (5,5) multipoles.When restricting to the (l, |m|) = (2, 2) modes, we refer to the model as SEOBNRv4E.We modify the initial conditions of the original version of the SEOBNRv4EHM model, which was presented in Ref. [95], so that now they are specified at an orbit-averaged orbital frequency instead of an instantaneous orbital frequency.We show that this new parametrization simplifies the stochastic sampling across parameter space.Furthermore, in order to increase the efficiency in evolving the dynamics, we combine reduced tolerances of the Runge-Kutta integrator during the inspiral with the optimized Hamiltonian and integrator from Refs.[97,98], and assess the impact of these modifications in the accuracy of the waveform model across parameter space.We find that for most of the parameter space the loss of accuracy is < 1% in unfaithfulness, indicating that the more efficient version of the SEOBNRv4EHM is still highly accurate for our purposes.
With the more efficient model, we use a highly parallelizable nested sampler parallel Bilby [99] to assess firstly the impact of the different initial conditions, which available in the SEOBNRv4EHM model.We find that the initial condi-tions specified at an orbit-averaged frequency perform as accurately as the ones specified at an instantaneous orbital frequency, but with a more efficient sampling of the parameter space, which translates into shorter wall clock times for the parameter-estimation runs.Furthermore, we show, for the first time, the impact of the radial phase parameter (i.e., the relativistic anomaly), using an inspiral-merger-ringdown model and a full parameter-estimation code (see, e.g., Ref. [91] for a study of the argument of periapsis based on the overlaps).
We also validate the model SEOBNRv4EHM by performing zero-noise injections of public eccentric NR waveforms from the Simulating eXtreme Spacetimes (SXS) catalog [100,101].In particular we select a set of three equal-mass, nonspinning eccentric simulations with initial eccentricities measured from the orbital frequency at first periastron passage of 0.07, 0.13 and 0.25, respectively.In order to compare the eccentricity from NR waveforms and the SEOBNRv4EHM model, we choose a common definition of eccentricity based on the GW signal, e gw , as introduced in Ref. [102], and measure it from the waveform using the gw eccentricity package [103].We find that SEOBNRv4EHM is able to accurately recover the injected parameters of the NR signals for all the eccentricity values considered.
Finally, we analyze some GW events observed by the LVK collaboration with SEOBNRv4EHM.In particular, we analyze GW150914 [104], GW151226 [105] and GW190521 [106,107].The choice is based on the fact that the first GW signal is still one of the loudest GW events so far, and has been found as a non-eccentric binary by many studies in the literature [72,81,82,108].On the other hand, GW151226 and GW190521 have been found to have signatures of eccentricity by some studies in the literature [75,[78][79][80]82].We measure the following values of eccentricity for these GW events, e GW150914 gw, 10Hz = 0.08 +0.09 −0.06 , e GW151226 gw, 20Hz = 0.04 +0.05 −0.04 , and e GW190521 gw, 5.5Hz = 0.15 +0.12 −0.12 , and therefore, find no clear evidence of eccentricity in these signals.
This paper is organized as follows.In Sec.II, we provide an overview of the eccentric model SEOBNRv4EHM.In Sec.II B, we introduce new initial conditions specified at an orbit-averaged orbital frequency, and in Sec.II C, we introduce an optimized version of the model, SEOBNRv4EHM opt, and assess its performance.We present the methodology for parameter estimation in Sec.III A. In Sec.III B we demonstrate the accuracy of SEOBNRv4EHM opt in the quasi-circular limit, in Sec.III C we show the importance of the radial phase parameter with mock-signal injections, and in Sec.III D we asses the accuracy of the model against eccentric NR waveforms from the SXS catalog using zero noise.In Sec.III E, we analyze GW events detected by the LVK Collaboration.In Sec.IV, we summarize our main conclusions and discuss future work.Finally, in Appendix A we provide details about the derivation of the expressions used in the initial conditions specified at an orbit-averaged frequency.

NOTATION
In this paper, we use geometric units, setting G = c = 1 unless otherwise specified.
We consider a binary with masses m 1,2 , with m 1 ≥ m 2 , and spins S S S 1,2 .We define the following combination of masses A relevant combination of masses for GW data analysis is the chirp mass defined as [109] M = ν 3/5 M. ( For non-precessing configurations the spin components are aligned/anti-aligned with the orbital angular momentum, L L L, with direction L L L. In this case the dimensionless spin components can be defined as where i = 1, 2. It is convenient to define the effective-spin parameter χ eff [86,110,111], and the spin combinations,

II. ECCENTRIC EFFECTIVE-ONE-BODY WAVEFORM MODEL WITH NON-PRECESSING SPINS
In this section we describe the waveform model, SEOBNRv4EHM, used throughout the paper.We provide an overview of the model, develop a new procedure to specify the initial conditions at an orbit-averaged frequency, and introduce modifications to increase the computational efficiency of the model.

A. Overview
The SEOBNRv4EHM model was presented in Ref. [95], and we refer therein for further details.SEOBNRv4EHM is constructed within the EOB formalism [83][84][85][86][87], and thus it is constituted by three main building blocks: • EOB Hamiltonian: The EOB conservative dynamics is determined by the EOB Hamiltonian constructed from the effective Hamiltonian, H eff , as described in Refs.[112,113], augmented with the parameters (K, d SO , d SS , ∆t 22 peak ) calibrated to NR waveforms from Ref. [114], through the energy map [83] For spins anti-aligned/aligned with the orbital angular momentum the motion is restricted to a plane.As a consequence, the dynamical variables in the Hamiltonian are the (dimensionless) radial separation r ≡ R/M, the orbital phase ϕ, and their (dimensionless) conjugate momenta p r ≡ P r /µ and p ϕ ≡ P ϕ /µ.
• Radiation-reaction force: The dissipative effects in the EOB dynamics are described by a radiation-reaction (RR) force F , which enters the Hamilton equations of motion, as [115,116] where the dot represents the time derivative d/d t, with respect to the dimensionless time t ≡ T/M, ĤEOB ≡ H EOB /µ, and Fϕ ≡ F ϕ /M.The equations are expressed in terms of p r * ≡ p r ξ(r), which is the conjugate momentum to the tortoisecoordinate r * , and ξ(r) ≡ dr/dr * can be expressed in terms of the potentials of the effective Hamiltonian [115].The components of the RR force are computed using the following relations [84,87] where ω = φ is the (dimensionless) orbital frequency, and Φ E is the energy flux for quasi-circular orbits written as a sum over waveform modes using [117,118] where d L is the luminosity distance between the binary system and the observer.
where t ℓm match is defined from the peak of the amplitude of the (2, 2)-mode (see Ref. [96] for details).The mergerringdown modes of SEOBNRv4EHM are the same as in the underlying SEOBNRv4HM model [96,114], and thus we assume that the effects of eccentricity at mergerringdown are negligible.The inspiral-plunge multipoles are constructed upon the NR-calibrated quasicircular SEOBNRv4HM multipoles [96] by incorporating the 2PN eccentric corrections, including spin-orbit and spin-spin effects, as derived in Ref. [94], and are augmented by an orbit-averaged procedure to compute the nonquasi-circular (NQC) terms (see Sec. II B from Ref. [95] for details).

B. New eccentric initial conditions
The initial conditions of SEOBNRv4EHM are expressed in terms of the eccentricity, e, and the relativistic anomaly ζ defined in the Keplerian parametrization of the orbit [95] where u p is the inverse semilatus rectum.Given the masses, spins, the initial instantaneous orbital frequency ω 0 , initial eccentricity e 0 and relativistic anomaly ζ 0 , the initial conditions for r 0 and p ϕ 0 , in absence of radiation reaction, can be obtained by solving the equations [95] ∂ ĤEOB with p r (p ϕ , e, ζ) and ṗr (p ϕ , e, ζ) given by the 2PN-order expressions in Eqs.(C3) and (C4) of the Appendix C in Ref. [95].
The initial condition for p r 0 can be computed using the solution for r 0 and p ϕ 0 , and numerically solving where ṙ(0) is the 2PN-order expression for ṙ at zeroth order in the RR effects (see Eq. (C5) in Ref. [95]), while ṙ(1) is the first-order term in the RR part of ṙ, for which we use the quasi-circular expression derived in Ref. [87] ṙ being Φ qc E the quasi-circular energy flux given in Eq. ( 9).Finally, the initial value p r 0 is converted into the tortoisecoordinate conjugate momentum p r * 0 , using the relations in Sec.II A of Ref. [95], so that together with r 0 and p ϕ 0 , it can be introduced in Eqs.(7) to evolve the EOB equations of motion.
The specification of eccentricity and relativistic anomaly at an instantaneous orbital frequency, which enters the righthand side (RHS) of Eq. ( 12), implies significant variations of the length of the EOB dynamics (and thus of the waveform length), as one of the eccentric parameters is modified.This effect is displayed in Fig. 1 for the orbital frequency evolution for a configuration with mass ratio q = 3, dimensionless spins χ 1 = 0.3, χ 2 = 0.5, and total mass 70M ⊙ at a starting frequency of 20Hz.In particular, we show in the upper panel of Fig. 1, the instantaneous orbital frequency at a fixed initial eccentricity, e 0 = 0.2, for three values of the initial relativistic anomaly, ζ 0 = [0, π/3, π] (dashed lines).Imposing that the initial eccentricity and relativistic anomaly are specified at an instantaneous orbital frequency causes necessarily that the length of the evolution is substantially different when the initial relativistic anomaly is at periastron (ζ 0 = 0) or apastron (ζ 0 = π).In the bottom panel of Fig. 1, we display the instantaneous orbital frequency for three values of initial eccentricity, e 0 = [0.01,0.1, 0.2], at a fixed value of the initial relativistic anomaly (ζ 0 = π/3).The evolutions (dashed lines) show that the higher the eccentricity the longer the evolution due to the larger oscillations in the instantaneous orbital frequency, and therefore, the chosen value of ω 0 is crossed at earlier times in the inspiral.
We note that specifying the initial eccentricity and rela-tivistic anomaly at an instantaneous orbital frequency is a particular parametrization of elliptical orbits.Other possible parameterizations exist in the literature.For instance, in post-Newtonian (PN) models based on the quasi-Keplerian parametrization [119][120][121][122][123][124][125][126][127][128][129][130][131][132][133][134][135], the initial conditions for the evolution of a binary in elliptical orbits are specified at an initial orbit-averaged orbital frequency.In these evolutions the initial parameters are the orbit-averaged orbital frequency, eccentricity and radial phase (typically the mean anomaly).One of the consequences of this parametrization is that at a fixed value of eccentricity, changes in the radial phase correspond to different positions in the same orbit.While an increase (decrease) of the value of eccentricity at fixed value of the radial phase causes a decrease (increase) of the length of the evolution.An additional motivation to explore a new parameterization of the initial conditions for SEOBNRv4EHM is the application of the model to Bayesian inference studies.Particularly, when performing parameter estimation with SEOBNRv4EHM, the usage of the initial conditions based on the instantaneous orbital frequency produce an increase in the structure of the log-likelihood surface as can be observed in Fig. 2. There, we show the log-likelihood computed for a zero-noise injection of a SEOBNRv4E waveform with initial eccentricity e 0 = 0.1, total mass 65M ⊙ , and dimensionless spins χ 1 = 0.3, χ 2 = 0 at a starting frequency of 20Hz, and recovering with SEOBNRv4E with all the parameters fixed to injected values except for the initial eccentricity and relativistic anomaly.We consider 5000 random points in the parameter space ζ 0 ∈ [0, 2π] and e 0 ∈ [0, 0.3], and use the parameter estimation code Bilby [136,137] to compute the likelihood.In the uppermost panel we fix the value of the initial relativistic anomaly to periastron (ζ 0 = 0) and sample only in eccentricity, while in the mid panel we fix ζ 0 = 1 and sample both in initial eccentricity and relativistic anomaly.In both panels the specification of the initial conditions at an instantaneous orbital frequency leads to a complex structure in the log-likelihood values, which can pose a challenge for stochastic samplers as shown in Sec.III C.
Therefore, we introduce a new parametrization of the EOB initial conditions where the initial eccentricity and relativistic anomaly are specified at an orbit-averaged orbital frequency, ω 0 .The orbit-averaged orbital frequency can be defined as where T r is the radial period.The integral in Eq. ( 15) can be computed using the 2PN expressions in the Keplerian parametrization from Ref. [94], and a detailed derivation can be found in Appendix A.
The expression for ω 0 can be inverted so that the instantaneous orbital frequency is expressed in terms of the orbitaveraged frequency, eccentricity, relativistic anomaly, mass ratio and spins.Therefore, the instantaneous initial orbital frequency, ω 0 , entering the RHS of Eq. (12), can be expressed as 1 − e Given masses, spins, initial eccentricity and relativistic anomaly at a particular orbit-averaged frequency ω 0 , we use Eq. ( 16) to compute the corresponding instantaneous orbital frequency ω 0 , which enters the RHS of Eq. (12).The rest of the procedure to compute the EOB initial conditions is not modified.The impact of this new parameterization in the EOB dynamics, can be observed in Fig. 1, where the solid lines correspond to the orbital frequency evolution of SEOBNRv4EHM using initial conditions specified at an orbit-averaged frequency of 20Hz.The main effect of this new parameterization is the almost constant merger time when varying the initial relativistic anomaly at a fixed value of initial eccentricity (top panel), and the reduction of the length of the evolution with increasing e 0 at a fixed value of ζ 0 (bottom panel) due to the larger emission of radiation at periastron passages.This behavior resembles the one in PN models based on the quasi-Keplerian parametrization [126,128,131,134,135], except for the fact that the EOB initial conditions are not adiabatic as in PN, and thus causing small variations of the merger time with different values of ζ 0 at a fixed value of e 0 (see top panel of Fig. 1).
Finally, sampling in e 0 and ζ 0 at ω 0 , produces notably simpler log-likelihood structures as shown in Fig. 2. In the top panel, where the sampling is performed only in e 0 , the number of oscillations in the log-likelihood curves is reduced with respect to the initial conditions based on the instantaneous orbital frequency.The results of the orbit-average initial conditions create a pattern with an easy to identify maximum in log-likelihood (green curve) at the injected value of e 0 = 0.1.While in the bottom panel, where the sampling in e 0 and ζ 0 is performed, it emerges a clearly defined pattern, and the values of e 0 and ζ 0 at the maximum log-likelihood point for the orbit-averaged initial conditions are e orb-avg log L max = 0.099, ζ orb-avg log L max = 0.77, while for the instantaneous initial conditions are e orb-avg log L max = 0.086, ζ orb-avg log L max = 4.829.An accurate recovery of the injected values requires more than 5000 points in the 2D parameter, however, the results already indicate that even with a low number of points the initial conditions based on the orbit-averaged orbital frequency are closer to the injected value, and may need less points than the recovery using the initial conditions based on the instantaneous orbital frequency.Thus, the initial conditions at an orbit-averaged frequency may be more adequate for data analysis applications such as parameter estimation, and in Sec.III C we further explore the consequences of these different initial conditions with stochastic sampling techniques.

C. Computational performance
One of the main applications of waveform models is the inference of the source parameters using Bayesian inference methods.These methods typically require of the order of ∼ 10 7 − 10 8 or more waveform evaluations.The SEOBNRv4EHM model is built upon the SEOBNRv4HM model [96], and thus it inherits the low computational efficiency of the previous generation of SEOBNR models 1 .Several techniques exist to increase the computational efficiency of EOB waveforms, such as reduced-order or surrogate models [143][144][145][146][147][148][149][150][151], the postadiabatic approximation [152,153], as well as methods targeting specifically parameter estimation, such as reduced order quadratures [154][155][156][157] or relative binning [158].
Here, we decide to increase the efficiency of SEOBNRv4EHM by reducing the absolute and relative tolerances of the 4thorder Runge-Kutta integrator from 10 −10 and 10 −9 to 10 −8 and 10 −8 , respectively 2 .Furthermore, as in the SEOBNRv4EHM model the Hamiltonian and the radiation-reaction force are the same as in the SEOBNRv4HM model, we also use the optimized Hamiltonian and integrator from Refs.[97,98].In order to ease notation we refer to the model with reduced tolerances and optimizations as SEOBNRv4EHM opt, and SEOBNRv4E opt when referring to the model containing only the (l, |m|) = (2, 2) multipoles.
The reduction of the ODE tolerances implies an increase in the efficiency of the model with waveform evaluation timescales of the order of O(100 ms), while decreasing the accuracy of the model.In order to quantify the latter across parameter space, we compute the unfaithfulness between the SEOBNRv4EHM and SEOBNRv4EHM opt models for 4500 points in the parameter space q ∈ [1, 50], χ 1,2 ∈ [−0.9.0.9], e 0 ∈ [0, 0.5], ζ 0 ∈ [0.2π] for a dimensionless starting frequency of Mω 0 = 0.023.We define the inner product between two waveforms, h A and h B [159,160] where a tilde indicates Fourier transform, a star complex conjugation and S n ( f ) the power spectral density (PSD) of the detector noise.In this work, we employ for the PSD the LIGO's "zero-detuned high-power" design sensitivity curve [161].Similarly, as in Ref. [141], we use f in = 10Hz and f max = 2048Hz.
To assess the agreement between two waveforms -for instance, the signal, h s , and the template, h t , observed by a detector, we define the faithfulness function [96,162], 1 The computational efficiency of the SEOBNR models has been recently significantly improved by the new generation of SEOBNRv5 models [138][139][140][141][142]. 2 The reduction of the tolerances is a similar approach to the one followed in Ref. [78] to improve the efficiency of the TEOBResumS-Dali model [93].
where λ λ λ = {m 1,2 , χ 1,2 , e 0 , ζ 0 } denotes the set of intrinsic parameters of the binary.When comparing waveforms, we choose the same inclination angle for the signal and template ι s = ι t = π/3, while the coalescence time and azimuthal angles of the template, (t 0 t , φ 0 t ), are adjusted to maximize the faithfulness of the template.The maximizations over the coalescence time t c and coalescence phase φ 0t are performed numerically.Similarly as in Sec.IV of Ref. [141] we set a grid of 8 points in the coalescence phase of the signal φ 0 s ∈ [0, 2π], and average over it to compute F .Finally, we introduce the unfaithfulness or mismatch as In Fig. 3 we show the distribution of median unfaithfulness over the total mass range [20 − 300]M ⊙ between the models containing only the (l, |m|) = (2, 2) modes, SEOBNRv4E and SEOBNRv4E opt models, as well as for the models including also the subdominant harmonics, i.e., (l, |m|) = (2, 2), (2, 1), (3,3), (4, 4), (5,5) multipoles.The results demonstrate a remarkable good agreement between the optimized (SEOBNRv4EHM opt and SEOBNRv4E opt) and the original models (SEOBNRv4EHM and SEOBNRv4E), with a median of unfaithfulness of 7.7 × 10 −6 for the dominant mode models, and 2.1 × 10 −5 for the models including higher order modes.As expected the differences between models with higher order modes are larger than for the dominant-mode models, due to the fact that small changes in the termination of the dynamics caused by the modifications in SEOBNRv4EHM opt impact more significantly the higher multipoles.In particular, the reduced tolerances in the integration can lead to small differences in the non-quasicircular coefficients computed from the input values (see Ref. [114] for details), which affect more the higher modes due to their low power.Furthermore, we also observe a tail of cases with unfaithfulness larger than 1% between the original and optimized models.These cases correspond to the more challenging parts of the parameter space considered with e 0 > 0.3 − 0.5 and high spins (χ 1,2 > 0.8 − 0.9), where the models have known limitations, such as the orbit-averaged procedure (see Appendix B of Ref. [95]), and the use of non-quasicircular corrections calibrated to quasi-circular binaries, which may increase the differences between the models.Finally, we assess the computational efficiency of the SEOBNRv4EHM opt model by timing the waveform generation and comparing it to the original SEOBNRv4EHM model, as well as the SEOBNRv4HM model in the quasi-circular limit.We consider binary's configurations with mass ratio q = 3, dimensionless spins χ 1 = −0.5, χ 2 = 0.3, initial relativistic anomaly ζ 0 = 1, total mass range M ∈ [20,200]M ⊙ , start-ing frequency f start = 10Hz and two different initial eccentricities e 0 = 0, 0.2.The results of the walltimes to generate the waveforms are shown in Fig. 4, where we are including all the modes up to l = 4, and a fixed sampling rate of 8192Hz for all the total masses considered 3 .The outcome of the benchmark demonstrates the significant increase in speed of the SEOBNRv4EHM opt model with respect to the SEOBNRv4EHM and SEOBNRv4HM models.For the configurations considered, we observe approximately a factor 2-5 improvement in speed.In the quasi-circular limit SEOBNRv4EHM is on average, over the total mass considered, slightly (∼ 1.5 times) faster than SEOBNRv4HM due to a more efficient implementation of some operations involved in the waveform calculation of SEOBNRv4EHM, while SEOBNRv4EHM opt is a factor ∼ 2.7 faster than SEOBNRv4EHM.Regarding models with only the (l, |m|) = (2, 2) modes, the hierarchy of curves is similar with the SEOBNRv4E opt being faster than SEOBNRv4E by a factor ∼ 3.1.This increase in speed compared to models with higher-order modes is due to the lack of common operations performed for models with higher-order modes, which limit the speed of the latter.When considering an initial eccentricity e 0 = 0.2 in the bottom panel of Fig. 4, we find an average (over total masses) increase of speed of ∼ 7.7 and ∼ 3.8 for the SEOBNRv4E opt and SEOBNRv4EHM opt models.The speed increase is more significant at low total masses as the main cost of waveform generation comes from the evaluation of the EOB dynamics, which is computed less frequently in the optimized model due to the reduced integration tolerances.
In summary, the SEOBNRv4E opt and SEOBNRv4EHM opt models imply a significant acceleration in waveform evaluation with respect to their original counterparts with a minor reduction in accuracy, in the region of parameter space of interest for our study.Furthermore, the sampling rate used for the benchmarks here is rather high4 (8192Hz), and typical applications for data analysis may use lower ones, which implies that the walltimes for waveform evaluation may be further reduced.As a consequence the optimized models reach speeds competitive for parameter estimation, and we show in Sec.III that they can be used to perform parameter-estimation runs in the order of hours and days.

III. BAYESIAN INFERENCE STUDY
The main application of the SEOBNRv4EHM opt waveform model is the Bayesian inference of source parameters of GWs emitted by BBHs.Thus, we introduce the methods and parameter-estimation codes used to infer the binary parameters in Sec.III A, we show the accuracy of the model in the quasi-circular limit in Sec.III B, we assess in Sec.III C the impact of the different initial conditions discussed in Sec.II B, as well as the importance of the relativistic-anomaly parameter.In Sec.III D we further investigate the accuracy of the model by performing a series of synthetic NR signal injections into zero detector noise.Finally, in Sec.III E, we analyze three real GW events detected by the LVK collaboration (GW150914, GW151226 and GW190521), and compare with results from the literature.

A. Methodology for parameter estimation
For the parameter-estimation study we employ parallel Bilby5 [99], a highly parallelized version of the Bayesian inference Python package Bilby [136,137], incorporating the nested sampler dynesty [163].Based on previous experience with parallel Bilby [141], we use a number of auto-correlation times nact = 30, number of live points nlive = 2048, and the remaining sampling parameters with their default values, unless otherwise specified.Furthermore, the runs are performed using distance marginalization as implemented in Bilby, and the phase marginalization is activated when using the (l, |m|) = (2, 2)-mode models to further reduce the computational cost.
For the choice of priors, we follow broadly Refs.[1,141,164].We choose a prior in inverse mass ratio, 1/q, and chirp mass, M, such that it is uniform in component masses.The priors in initial eccentricity, e 0 , and relativistic anomaly, ζ 0 ∈ [0, 2π], are chosen to be uniform.In order to facilitate the comparison with precessing-spin results, the priors on the spin-components, χ 1,2 , are chosen such that they correspond to the projections of a uniform and isotropic spin distribution along the ẑ-direction [164].The luminosity distance prior is chosen to be proportional to ∝ d 2 L , unless otherwise specified.The rest of the priors are set according to Appendix C of Ref. [1].The specific values of the prior boundaries for the different parameters vary depending on the application, and we specify them in the subsequent sections.

B. Quasi-circular limit
Eccentricity is a parameter, which defines the ellipticity of an orbit between two limits: the parabolic and the circular case.In this section we consider the latter, and demonstrate that the SEOBNRv4EHM opt model is able to accurately describe GWs from quasi-circular BBHs.In Sec.III of Ref. [95] it is shown that SEOBNRv4EHM has a comparable accuracy to SEOBNRv4HM by computing the unfaithfulness against quasi-circular NR waveforms.
Here, we assess the accuracy of SEOBNRv4EHM opt in the zero-eccentricity limit by computing the unfaithfulness, as defined in Sec.II C, against the accurate quasi-circular SEOBNRv4HM model for 4500 random points distributed in the following parameter space: q ∈ [1, 50], χ 1,2 ∈ [−0.9.0.9], with an inclination angle of ι = π/3, for a total mass range [20 − 300]M ⊙ and starting frequency of Mω 0 = 0.023 in geometric units.In Fig. 5 we show the distribution of the median unfaithfulness over the total mass range considered when comparing the (l, |m|) = (2, 2)-mode models, SEOBNRv4 and SEOBNRv4E opt, as well as the corresponding models including higher multipoles, SEOBNRv4HM and SEOBNRv4EHM opt.For SEOBNRv4E opt the median value of unfaithfulness is 8.1 × 10 −6 , while when including higher order modes it degrades to 3.8 × 10 −5 .In both cases there are no configurations with unfaithfulness larger than 1%.Therefore, both the SEOBNRv4E opt and SEOBNRv4EHM opt models are faithful across parameter space to the SEOBNRv4 and SEOBNRv4HM models, respectively, considering that SEOBNRv4 was calibrated to NR with unfaithfulness below 1%.)-modes models, SEOBNRv4E opt and SEOBNRv4 (blue), as well as between the higher-order mode models, SEOBNRv4EHM opt and SEOBNRv4HM (green) for 4500 configurations in the parameter space q ∈ [1, 50], χ 1,2 ∈ [−0.9.0.9] for a dimensionless starting frequency of Mω 0 = 0.023.The vertical dashed lines indicate the median values of the distribution.
We further investigate the implications of these unfaithfulness results in parameter estimation by performing a mocksignal injection into zero-detector noise.With zero noise, and flat priors, the likelihood will peak at the true parameters when using the same model for injection and recovery.This makes it easier to see biases that are arising from model differences.We use the SEOBNRv4 model as a signal, and recover the injected parameters with the reduced order model (ROM) SEOBNRv4 ROM [114] and the SEOBNRv4E opt model.For the latter we also sample in initial eccentricity and relativistic anomaly.We consider a configuration with mass ratio q = 4, total mass M = 90.08M⊙ and BH's dimensionless spins χ 1 = 0.5 and χ 2 = −0.1 defined at 20Hz.
For this injection we choose the inclination with respect to the line of sight of the BBH to be ι = 0.1 rad, coalescence   I. and polarization phases are ϕ = 0.6 rad and ψ = 0.33 rad, respectively.The luminosity distance to the source is chosen to be 600 Mpc, which produces a three-detector (LIGO Hanford, LIGO Livingston and Virgo) network-SNR of ρ N mf = 67.9 when using the LIGO and Virgo PSD at design sensitivity [161].
We choose a uniform prior in inverse mass ratio and chirp mass, with ranges 1/q ∈ [0.05, 1] and M ∈ [5, 100]M ⊙ .The priors on the magnitudes of the dimensionless z-components of the spins are a i ∈ [0, 0.99].For SEOBNRv4E opt we take a uniform prior in the initial eccentricity e 0 ∈ [0, 0.3], and uniform in the initial relativistic anomaly The resulting posteriors for the inverse mass ratio, chirp mass, effective-spin parameter, luminosity distance, initial eccentricity and relativistic anomaly are shown in Fig. 6.We find remarkable agreement between the posteriors of SEOBNRv4E opt and SEOBNRv4E ROM for all the parameters.The injected values and recovered parameters are displayed in Table I, where additional parameters are shown.Even for a relatively high SNR injection (ρ N mf ∼ 68) both models are able to accurately recover the injected parameters for the inverse mass ratio, chirp mass and effective-spin parameter within the 90% credible intervals, while the luminosity distance presents a small bias due to the limited mode content of both models (only the (l, |m|) = (2, 2) multipoles), which creates a degeneracy between the inclination angle and the luminosity distance and thus, complicates the measurement of both quantities [96,[165][166][167].This degeneracy and, thus the bias, can be removed by including higher multipoles in the waveform, but we do not use models with higher multipoles, as the focus of this injection study is the assessment of the agreement between the SEOBNRv4E opt and SEOBNRv4E ROM models.
Regarding initial eccentricity and relativistic anomaly, SEOBNRv4E opt measures an initial eccentricity consistent with zero, e 0 = 0.01 +0.02 −0.01 , while the initial relativistic anomaly becomes an uninformative parameter as for quasi-circular orbits the radial phase provides no information about the position of the binary components.
In conclusion, the SEOBNRv4EHM opt model has an accurate zero eccentricity limit comparable to the underlying quasi-circular SEOBNRv4HM model.This is extremely important as eccentricity is a gauge dependent parameter in General Relativity between two well defined limits (circular and parabolic), and being able to unambiguously recover one of them ensures that the values of eccentricity can be quoted with respect to a uniquely defined physical configuration 6 .

C. Eccentric case
Two different initial conditions (ICs) for the SEOBNRv4EHM model were presented in Sec.II B: 1) one based on the specification of the initial eccentricity and relativistic anomaly, (e 0 , ζ 0 ) at an instantaneous orbital frequency ω 0 , hereafter referred as instantaneous ICs, and 2) one based on the specification of (e 0 , ζ 0 ) at an orbit-averaged orbital frequency ω 0 , hereafter called orbit-averaged ICs.
Here, we study the implications of these ICs by performing two mock-signal injections in zero detector noise using the SEOBNRv4E opt model as a signal with two different initial eccentricities e 0 = [0.1,0.2], and recovering with the SEOBNRv4E opt model.We consider a configuration with mass ratio q = 3, initial relativistic anomaly ζ 0 = 1.2, total mass M = 76.4M⊙ and BH's dimensionless spins χ 1 = 0.5 and χ 2 = −0.1 defined at a starting frequency of 20Hz 7The priors on the sampling parameters are chosen as in Sec.III B. The injection and recovery of the parameters are performed using the SEOBNRv4E opt model with orbit-averaged and instantaneous ICs.Furthermore, we also assess the importance of the relativistic anomaly parameter, ζ 0 , by recovering the parameters by setting the initial starting point of the orbit at periastron (ζ 0 = 0), which is different from the injected value ζ 0 = 1.2.For completeness we also measure the binary parameters in the quasi-circular limit (e 0 = 0).The resulting posterior distributions of the different cases are shown in Fig. 7 for the inverse mass-ratio, initial eccentricity and initial relativistic anomaly.The injected values, the median values and 90% credible intervals of the posterior distributions of some parameters are provided in Table II.
The left side of the violin plots in Fig. 7 shows the comparison of the orbit-averaged and instantaneous ICs as well as the zero-eccentricity case.We find that both ICs are able to accurately recover the injected parameters for all the parameters of the two injected eccentricities, e 0 = 0.1 and e 0 = 0.2.
In the case of the instantaneous ICs the ζ 0 posterior is sharply peaked at the injected value due to the fact that for these ICs a change in ζ 0 at a fixed value of e 0 provides a completely different evolution, while for the orbit-averaged ICs variations of ζ 0 at a fixed value e 0 describe different values of the radial phase in the same orbit as shown in Fig. 1.These different descriptions of the orbits have also implications in the efficiency of the sampler and the cost of the parameter estimation runs.For instance, the e 0 = 0.2 injections were performed with the same sampler settings as described in Sec.III A and using 6 nodes of 32 cores each8 with an averaged wall clock time of 11.8 and 35.7 hours for the orbit-averaged and instantaneous ICs, respectively.This demonstrates that the orbit-averaged ICs lead to a more computationally efficient sampling of the parameter space than the instantaneous ICs without a loss of accuracy.
We focus now on the right side of the violin plots in Fig. 7, where the results of setting the initial relativistic-anomaly parameter to a different value (ζ 0 = 0) from the injected one (ζ 0 = 1.2) for both the instantaneous and the orbit-averaged ICs are displayed.For the orbit-averaged ICs, neglecting the relativistic anomaly parameter can lead to biases in the quasi-circular parameters like the inverse mass ratio, while the eccentricity parameter is accurately recovered (see right-side distributions in Fig. 7).This can be explained by the fact that during the parameter-estimation run the template waveform compensates the lack of another eccentric degree of freedom by modifying the rest of the parameters in order to describe more accurately the signal.This points out the relevance of taking into account the radial-phase parameter when using orbit-averaged ICs, and is in agreement with the findings of Ref. [170] during the construction of an eccentric surrogatewaveform model.Regarding the instantaneous ICs, the quasicircular parameters are recovered with no biases, but a multimodal posterior in the eccentricity parameter is found in the middle panels of Fig. 7.These multimodalities can be explained by a degeneracy in the waveforms with fixed ζ 0 at a relatively high total mass (M ∼ 76M ⊙ ), and the fact that the model assumes circularization at merger-ringdown.As shown in Sec.III D these multimodalities can sometimes be an artifact of the specific parameterization chosen, and be removed by going to a definition of eccentricity based on the waveform.However, this is not the case for the injections shown in this section, and the multimodalities remain even after a redefinition of the eccentric parameters, indicating that for these particular injections these are true multi-modalities in the posterior distribution.Therefore, the instantaneous ICs when neglecting ζ 0 can substantially complicate the sampling as well as the measurement of the eccentricity parameter due to the multimodalities.
Finally, we also show in Fig. 7 and .The median values also report the 90% credible intervals.The binary parameters correspond to the total mass M, chirp mass M, inverse mass ratio 1/q, effective spin parameter χ eff , initial eccentricity e 0 , initial relativistic anomaly ζ 0 , angle between the total angular momentum and the line of sight θ JN , luminosity distance d L , coalescence phase ϕ ref and the network matched-filtered SNR for LIGO-Hanford/Livingston and Virgo detectors ρ N mf .For each injection the recovery is done with the SEOBNRv4E opt model using orbitaveraged ICs (ω 0 ICs), instantaneous ICs (ω 0 ICs), and setting ζ 0 = 0 for both the orbit-averaged and instantaneous ICs, ω 0 ICs (ζ 0 = 0) and ω 0 ICs (ζ 0 = 0), respectively.Additionally, the values recovered setting the initial eccentricity to zero, e 0 = 0, are also shown.injection with small eccentricity at the injected SNR (∼ 46) the quasi-circular model is still able to recover the parameters accurately within the 90% credible intervals.However, for the high eccentricity injection (e 0 = 0.2) this is no longer the case and the biases in parameters like the chirp mass can be as large as 8% with respect to the injected value.
In summary, the orbit-averaged ICs provide a more efficient sampling of the eccentric parameter space than the instantaneous ones without a loss of accuracy.As a consequence we adopt the orbit-averaged ICs hereafter for the analysis using the SEOBNRv4EHM opt model in this paper.Furthermore, we have shown that neglecting the radial-phase parameter, as currently done in the TEOBResumS-Dali and the SEOBNREHM [90] models can lead to biases in the recovered parameters for both instantaneous and orbit-averaged ICs, unless one varies additional parameters.For instance, in Ref. [81] to avoid biases in the posteriors with TEOBResumS-Dali, which only employs the eccentricity parameter to describe elliptical orbits, the starting frequency of the waveform is also sampled during the parameter estimation runs.However, we find more natural to keep the starting frequency of the waveform fixed as done in the LVK analysis of the Gravitational Wave Transient Catalogs [2][3][4], and thus, use two eccentric parameters, eccentricity and relativistic anomaly, which can vary freely during the parameter-estimation run.

D. Numerical-relativity injections
In Ref. [95], the SEOBNRv4EHM model was shown to be accurate to with an unfaithfulness below 1% for a dataset of public eccentric NR waveforms from the SXS catalog [100,171].In this section we further investigate the accuracy of the SEOBNRv4EHM model against NR waveforms by performing zero-noise injections, and recovering the parameters with the SEOBNRv4E opt model.
We consider a set of 3 public eccentric simulations SXS:BBH:1355, SXS:BBH:1359 and SXS:BBH:1363, which correspond to equal-mass, nonspinning configurations with initial eccentricities measured from the orbital frequency at first periastron passage of 0.07, 0.13 and 0.25, respectively (see Table I of Ref. [95] for details).For these injections we choose a total mass M = 70M ⊙ , inclination with respect to the line of sight of the BBH ι = 0 rad, coalescence phase ϕ ref = 0 rad, and luminosity distance d L =2307 Mpc, which produces a three-detector (LIGO Hanford, LIGO Livingston and Virgo) network-SNR of ρ N mf = 20 when using the LIGO and Virgo PSD at design sensitivity.The priors are the same as in Sec. on the right side of the x-axis to ease the visualization of the results.The dashed vertical lines indicate the 90% credible intervals.The solid vertical lines correspond to the injected parameters, which are also shown in Table I.Bottom row: Same as in the upper row but for an injected signal with initial eccentricity e 0 = 0.2 at 20Hz.III C, with the only exception that for the run corresponding to the SXS:BBH:1363 NR waveform we set a larger upper bound for the eccentricity prior of e 0 ∈ [0, 0.5] in order to avoid railing of the posterior against the upper bound of the prior.Before analysing the results of the injections, we consider the problem of mapping the eccentric parameters, eccentricity and radial phase, from the NR waveforms to the SEOBNRv4E opt model.This problem stems from the gaugedependency of the eccentricity parameter in General Relativity, and can be avoided by adopting a common definition of the parameters defining elliptical orbits in the NR simulation and the SEOBNR model.In particular, we adopt a definition of eccentricity, e gw , measured from the frequency of the (2,2)mode with the correct Newtonian limit [102] with where ω 22,a , ω 22,p refer to the values of the (2,2)-mode frequency at apastron and periastron, respectively.As the radial-phase parameter describing a binary in an elliptic orbit, we use the mean anomaly, with the following definition [172] where t p i is the time of the i-th periastron passage measured from the (2,2)-mode frequency.
These definitions of eccentricity and mean anomaly can be applied to the posterior distributions from the parameterestimation runs as a post-processing step employing its highly efficient implementation in the gw eccentricity Python package9 [103].The procedure consists in evaluating the waveform for each sample of the posterior distributions, and applying the gw eccentricity package to measure the eccentricity and mean anomaly at a desired point in the evolution.The process can be parallelized to measure the eccentricity and mean anomaly in a much smaller timescale than an actual parameter-estimation run 10 .In the SEOBNRv4EHM model the eccentricity is specified at an orbit-averaged orbital frequency, however, the eccentricity definition introduced in Eq. ( 20) is based on the (2,2)-mode frequency.As shown in Ref. [102], for eccentric binaries the instantaneous orbital and (2,2)-mode frequencies are not related by a simple factor of 2, as in the quasi-circular case, and this can cause that for some samples in the posterior distribution, the instantaneous (2,2)-mode frequency does not reach the starting frequency specified in the initial conditions of the SEOBNRv4EHM model.In order to avoid that this situation prevents the eccentricity measurement for some samples we have also implemented an option to integrate the EOB dynamics backwards in time from a specific starting frequency.For the rest of the calculations involving gw eccentricity in this paper, we integrate 2000M backwards in time.
In Fig. 8, we summarize the parameter-estimation results of the injections.We report the marginalized 1D and 2D posteriors for the chirp mass M and the effective-spin parameter χ eff , the initial eccentricity and relativistic anomaly, and the GW eccentricity and mean anomaly measured at 20Hz.In Table III we provide the values of the injected parameters and the median of the inferred posterior distribution with the 90% confidence intervals for both models.The results show that SEOBNRv4E opt is able to recover M and χ eff , as well as the mass ratio and total mass for all the NR injections within the 90% confidence intervals.The NR waveforms contain all the multipoles up to l ≤ 8, while the SEOBNRv4E opt contains only the (l, |m|) = (2, 2) modes.This difference in mode content explains why there are some small biases in the luminosity distance and θ JN parameter.
Regarding the eccentric parameters, the middle panel of Fig. 8 shows the initial eccentricity and relativistic anomaly posterior distributions at 20Hz for the three different injections.For the highest eccentricity injection (SXS:BBH:1363), the eccentricity posterior is bimodal with two modes centered at e 0 ∼ 0.25 and e 0 ∼ 0.35.However, when moving to the definition of eccentricity and mean anomaly based on the (2,2)mode frequency introduced in Eq. ( 20), we find a unimodal posterior in e gw and l gw , which indicates that the bimodality in (e 0 , l 0 ) is simply a consequence of the parametrization of the EOB initial conditions used in the SEOBNRv4E opt model.Furthermore, the e gw and l gw parameters at 20Hz are measured from the NR waveforms, and shown as vertical lines in the right panel of Fig. 8.The posterior distributions for e gw and l gw are consistent within the 90% credible intervals with the injected values for the three NR waveforms.
The results of the injections demonstrate that SEOBNRv4E opt is able to accurately recover the eccentric and non-eccentric parameters of the injected NR waveforms, and they are consistent with the low unfaithfulness values of SEOBNRv4E against NR waveforms reported in Ref. [95].Further studies of the accuracy of the model will require larger datasets of eccentric NR waveforms including spins and higher eccentricities, and we leave for future work investigating the waveform systematics of the SEOBNRv4EHM opt model and its biases against NR waveforms.

E. Analysis of GW events
In this section, we analyze 3 GW events recorded by the LIGO and Virgo detectors [1,3,4] during the first and third observing runs: GW150914, GW151226 and GW190521.We employ strain data from the Gravitational Wave Open Source Catalog (GWOSC) [173], and the released PSD and calibration envelopes included in the Gravitational Wave Transient Catalogs GWTC-2.1 [3], and their respective parameter- estimation samples releases 11 .
We analyze GW150914, GW151226 and GW190521 using SEOBNRv4EHM opt and SEOBNRv4E opt with parallel Bilby and the settings described in Sec.III A. For SEOBNRv4EHM opt, we restrict to a mode content l ≤ 4 in order to avoid the increase of computational cost due to the high sampling rates necessary to resolve the (5, 5)-mode, as at the current SNRs the impact of this mode on accuracy is limited.

GW150914
The first observation of GWs from a BBH coalescence, GW150914, had one of the highest SNRs (∼ 23.7) of the GW 11 For GW190521 we employ the samples of the SEOBNRv4PHM model from Refs.[106,107], which were produced using parallel Bilby.
events observed during the first three observing runs [3,104] of the LVK.The binary parameters are consistent with a nonspinning binary with comparable masses [108].We choose priors in inverse mass ratio, 1/q ∈ [0.05, 1], and chirp mass, M ∈ [20,50]M ⊙ , such that the induced priors in component masses are uniform.Uniform priors were also used for initial eccentricity, e 0 ∈ [0, 0.3], and the initial relativistic anomaly, ζ 0 ∈ [0., 2π].The other priors are chosen as in Sec.III D. For the analysis we use a starting frequency of 10Hz at which the waveform is generated and, thus at which the initial eccentricity and relativistic anomaly are defined.This choice of starting frequency ensures that the higher order modes in SEOBNRv4EHM opt are in band at the minimum frequency of 20Hz at which the likelihood calculation starts.We also remark the importance that both SEOBNRv4EHM opt and SEOBNRv4E opt are generated at the same starting frequency, so that both have the same priors in initial eccentricity.
The posterior distributions of chirp mass, effective-spin parameter, GW eccentricity and GW mean anomaly are displayed in the top row of Fig. 9.In Table IV we also report the median values and the 90% credible intervals of the posterior distributions for other binary parameters.For comparisons we include in our analysis the samples of SEOBNRv4PHM [162] from the GWTC-2.1 catalog [3].We find that binary parameters like chirp mass, effective-spin parameter and mass ratio measured by SEOBNRv4E opt and SEOBNRv4EHM opt are consistent with the ones measured with SEOBNRv4PHM.This is expected as GW150914 is consistent with a nonspinning binary, and thus the effects of spin-precession which are accurately described by SEOBNRv4PHM are negligible.Regarding the eccentric parameters, although both SEOBNRv4E opt and SEOBNRv4EHM opt have median values of eccentricity distinct from zero, e 10Hz gw = 0.08 +0.1 −0.07 and e 10Hz gw = 0.08 +0.09 −0.06 , respectively, the posterior distributions have a strong support in the zero eccentricity region, which is in agreement with other analyses of GW150914 with eccentric waveforms [72,81,82,108].Additionally, we have produced a run with SEOBNRv4E opt setting e 0 = 0, and obtained a signalto-noise natural log Bayes factor of 284.6 +0.1 −0.1 which is slightly larger than 283.9 +0.1 −0.1 obtained for SEOBNRv4E opt when sampling in (e 0 , ζ 0 ) (see Table IV).Unfortunately, this information is not available for the SEOBNRv4PHM model in the public GWOSC results, but GWTC2-1 does include a log Bayes factor for the similar IMRPhenomXPHM model, of log BF = 303.45+0.14 −0.14 .These results indicate that the non-precessing eccentric hypothesis is disfavoured against the precessing-spin quasi-circular one, and that GW150914 is more consistent with a quasi-circular BBH merger.

GW151226
GW151226 is one of the GW events with lowest total mass observed in the first observing run, and it was identified in the GWTC-1 catalog [1] to exclude support for χ eff = 0 at > 90% probability.Furthermore, Ref. [78] analyzed this event with the TEOBResumS-Dali model [169] and parallel Bilby, with a uniform prior in eccentricity, and constrained the ini- .Effective-spin parameter, chirp mass, GW eccentricity and GW mean anomaly parameters inferred for the real GW events analysed with SEOBNRv4E opt and SEOBNRv4EHM opt.Comparisons are presented with the SEOBNRv4PHM (when available) and the IMRPhenomXPHM models from the GWTC-2.1 catalog [3], except for GW190521 where samples are from Refs.[106,107].The waveform starting frequency for GW150914 is f start = 10Hz, for the low total mass event GW151226 is f start = 20Hz, in order to reduce the computational cost, while for GW190521 f start = 5.5Hz in order to have the mode content in band at the likelihood minimum frequency ( f min = 11Hz).
tial eccentricity to be e 0 < 0.15 at 90% at a starting frequency of 10Hz.Moreover, Ref. [72] using the SEOBNRE model [89] and the reweighing technique with a log-uniform prior in eccentricity found a much tighter constraint e 0 < 0.04 at 10Hz.
For our analysis we use a uniform prior in initial eccentricity e 0 ∈ [0, 0.3] and relativistic anomaly ζ 0 ∈ [0, 2π], and priors in inverse mass ratio, 1/q ∈ [0.125, 1], and chirp mass, M ∈ [5, 100]M ⊙ , such that they are uniform in component masses.The rest of the priors are chosen as in the analysis of GW150914.We use a starting frequency of 20Hz at which the waveform is generated and, at which e 0 and ζ 0 are defined.Due to the low total mass of the event, we restrict to a higher starting frequency than in the case of GW150914, and we only use the SEOBNRv4E opt model in order to reduce the  IV.Median and 90% credible intervals for the total mass, chirp mass, inverse mass ratio, effective-spin parameter, GW eccentricity, GW mean anomaly, θ JN , luminosity distance and signal-versus-noise (natural) log Bayes factor, log BF , measured with the SEOBNRv4E opt and SEOBNRv4EHM opt models for GW150914, GW151226 and GW190521.Additionally, we include the results of SEOBNRv4PHM for GW150914 and GW190521 and of IMRPhenomXPHM for GW151226 from the GWTC-2.1 catalog [3], except for GW190521 where samples are from Refs.[106,107].The waveform starting frequency for GW150914 is f start = 10Hz, for the low total mass event GW151226 is f start = 20Hz, in order to reduce the computational cost, while for GW190521 is f start = 5.5Hz in order to have the mode content in band at the likelihood minimum frequency ( f min = 11Hz).)) or SEOBNRv4EHM opt ( HMs ≡ higher modes) to perform the run.Sampling rate (srate) and data segment duration (seglen) are specified in the data settings, while the number of autocorrelaction times (nact) and number of live points (nlive) are specified in the sampler settings.The time reported is walltime, while the total computational cost in CPU hours can be obtained multiplying this time by the reported number of CPU cores employed.

GW event sampler
computational cost.
The results are shown in the middle row of Fig. 9, while in Table IV the median values and the 90% credible intervals of the posterior distributions are reported.For comparisons we include in our analysis the samples of IMRPhenomXPHM [174] from the GWTC-2.1 catalog [3].The quasi-circular binary parameters like chirp mass, effective-spin parameter and mass ratio measured by SEOBNRv4E opt show differences with respect to the ones measured by IMRPhenomXPHM.
This can be explained because of the different physical content of each model, as IMRPhenomXPHM includes higherorder modes and describes precessing-spin quasi-circular binaries, while SEOBNRv4E opt includes only the (l, |m|) = (2, 2) modes and describes non-precessing eccentric binaries.This translates into a preference of IMRPhenomXPHM for unequal masses and larger effective spin values than SEOBNRv4E opt.Apart from these differences, the posterior distributions of SEOBNRv4E opt and IMRPhenomXPHM have large overlapping regions as can be observed in the left plot in the middle row of Fig. 9.
The value of eccentricity measured by SEOBNRv4E opt is e 20Hz gw = 0.04 +0.05 −0.04 , and the posterior distribution of GW eccentricity (middle and right plots in Fig. 9) have strong support at zero eccentricity.This combined with an uninformative posterior distribution for the GW mean anomaly indicate that GW151226 is more consistent with a quasi-circular binary than a non-precessing eccentric binary with e 0 ≤ 0.3.The comparison between the value of eccentricity in Table IV with the ones in Refs.[72,78] is challenging due to the following issues: the value of eccentricity that we report is computed at a different starting frequency (20Hz), than the one (10Hz) used in Refs.[72,78], and the definition of eccentricity that we are using is based on the waveforms generated from the samples, while Refs.[72,78] report the initial eccentricities considered in the definitions of the initial conditions in the TEOBResumS-Dali and SEOBNRE models.In order to estimate the eccentricity value at some other frequency we can apply the Newtonian relation between the eccentricity and frequency [23] e = e ref ( f / f ref ) −19/18 to our eccentricity measurement to obtain e 10Hz gw ∼ 0.1, which is closer to the value reported in Ref. [78].This fact can be explained probably due to the fact that Ref. [72] uses a log-uniform prior in eccen-tricity which puts more weight on the low eccentricity region than a uniform prior.
Finally, we have also produced a run with SEOBNRv4E opt setting e 0 = 0, and obtained a signal-to-noise natural log Bayes factor of 39.8 +0.1 −0.1 which is slightly smaller than 40.4 +0.1 −0.1 obtained for SEOBNRv4E opt when sampling in (e 0 , ζ 0 ) (see Table IV).The difference in log Bayes factors is ∼ 0.6 which indicates a minor preference for the eccentric hypothesis.When comparing to the quasi-circular precessing-spin results of IMRPhenomXPHM from the GWTC-2.1 catalog we find log BF = 47.59 +0.14 −0.14 , which indicates that the eccentric non-precessing spin hypothesis is less preferred than the precessing-spin quasi-circular one with a log Bayes factors of ∼ 7.2 in favour of the latter.

GW190521
GW190521 is particularly intriguing, with only 4 cycles in the band of the detectors, thus being consistent with a mergerringdown dominated signal.It has been attributed to a variety of physical systems including a head-on collision of exotic compact objects [175], a nonspinning hyperbolic capture [80] and an eccentric binary [74,75], although some other recent studies do not find clear evidence for eccentricity [82].
We analyze GW190521 using the SEOBNRv4EHM opt and SEOBNRv4E opt models with a prior uniform in initial eccentricity, e 0 ∈ [0, 0.3], and relativistic anomaly, ζ 0 ∈ [0, 2π].We employ priors in inverse mass ratio 1/q ∈ [0.05, 1], and chirp mass M ∈ [60,200]M ⊙ such that the induced priors are uniform in component masses.The rest of the priors are as in the analysis of GW150914, except for the luminosity distance prior which is chosen to be uniform in comoving volume instead of ∝ d 2 L , in order to match the settings of Ref. [141].The starting frequency of waveform generation is 5.5Hz such that the higher-order modes are in band at the minimum frequency of the likelihood evaluation (11Hz).As discussed in Ref. [176], the analysis of GW190521 using different sampling methods can lead to systematics in recovering some modes in the posterior distribution.In order to avoid that, we use a higher number of live points of 8192 for the parallel Bilby runs of this event for SEOBNRv4E opt, while for SEOBNRv4EHM opt we use a lower number of live points 2048 to reduce the computational cost 12 .
The results of our analysis are shown in the bottom row of Fig. 9, while in Table IV the median values and the 90% credible intervals of the posterior distributions are reported.For comparisons we include in our analysis the samples of SEOBNRv4PHM from Refs.[106,107].The quasi-circular parameters like chirp mass or total mass show the largest differences with respect to the SEOBNRv4PHM, while other pa-rameters like the effective-spin parameter and the mass ratio are quite consistent in the median.However, as in the case of non-precessing quasi-circular models (see Ref. [176]), the SEOBNRv4EHM opt and SEOBNRv4E opt models are not able to reproduce the secondary mode in the inverse mass ratio posterior.These differences between SEOBNRv4PHM and SEOBNRv4EHM opt are expected due to the GW190521 signal being merger-ringdown dominated, and these waveform models include distinct physical effects.
Focusing on the eccentric parameters, we find median values of the GW eccentricity of e 5.5Hz gw = 0.15 +0.12 −0.12 for both SEOBNRv4EHM opt and SEOBNRv4E opt.The large median values of eccentricity contrast with the uniformative posterior distribution of eccentricity and the large uncertainty in the 90% credible intervals, which combined with an uninformative posterior distribution of the GW mean anomaly (see bottom panels of Fig. 9), indicates that the eccentricity parameter is poorly constrained in GW190521.This situation can be explained by the fact that the signal is extremely short, and thus merger-ringdown dominated.However, in SEOBNRv4EHM opt the eccentricity effects are included only through the inspiralmerger EOB modes, while at merger-ringdown the binary is assumed to have circularized, and the merger input values and the ringdown model are the same as in SEOBNRv4HM (see Sec. II for details).For moderate eccentricities, as the ones considered here, and non-precessing spins it has been shown in Refs.[177][178][179] that the effects of eccentricity in the final mass and spin of NR waveforms are subdominant.However, non-precessing large eccentricity, as well as precessing-spin eccentric cases are fairly unexplored.Therefore, in order to clearly measure eccentricity in high-mass systems, the binary should have large enough eccentricity at merger.
Apart from this limitation of the SEOBNRv4EHM opt model 13 , we have also attempted to estimate the validity of the non-precessing eccentric hypothesis versus the non-precessing quasi-circular one by producing a run with SEOBNRv4E opt setting e 0 = 0, and comparing the signalto-noise Bayes factors.For the zero-eccentricity run with SEOBNRv4E opt we find log BF = 77.7 +0.1 −0.1 , which is consistent within the error bars with the values obtained for the eccentric run with SEOBNRv4E opt, log BF = 77.8+0.1 −0.1 .This points out that the non-precessing eccentric hypothesis is equally favoured as the non-precessing quasi-circular one.However, when comparing to the quasi-circular precessingspin results from the discovery paper of GW190521 14 [106] we observe that for the NRSur7dq4 model [180] the log Bayes factor is log BF = 84.49.This produces a log Bayes factor between the quasi-circular precessing-spin and nonprecessing eccentric hypothesis of ∼ 6.7, indicating that the quasi-circular precessing-spin hypothesis is preferred over the non-precessing eccentric one with a prior in e 0 ∈ [0, 0.3].
The comparison of Bayes factors from different waveform models can be complicated as the result may be affected not only by the different physical effects included in the models, but also by the waveform systematics between the different waveform approximants.Different waveform models including the same physical effects can lead to different log Bayes factors (see for instance Table II of Ref. [176]).As a consequence, we consider the comparison of the eccentric nonprecessing against quasi-circular precessing-spin hypothesis using log Bayes factors from other waveform families as an approximate estimate, and leave for future work a more detailed study [181].
Besides calibration to eccentric NR waveforms one of the main limitations in the analysis of GW190521 with SEOBNRv4EHM opt is the lack of inclusion of spin-precession effects, which impact significantly the morphology of the templates at merger-ringdown, and may substantially modify the measured parameters.This points out the necessity to produce waveform models, which include both eccentricity and spin-precession effects, and there is ongoing work [182] to include both effects in the new generation of SEOBNR models [140,141] built within the new pySEOBNR python infrastructure [183].
Finally, the SEOBNRv4EHM opt results have been obtained on the order of a few days or a week using parallel Bilby (see Table V).This makes SEOBNRv4EHM opt a standard tool that can be used with a highly parallelizable nested sampler like parallel Bilby, and we plan to extend the Bayesian inference study presented here, using the machine-learning code DINGO [184][185][186], to all the GW events observed during the third-observing run [181].

IV. CONCLUSIONS
In this paper we have improved and validated the multipolar non-precessing eccentric SEOBNRv4EHM model presented in Ref. [95], and shown its applicability to Bayesian inference studies.
Within the SEOBNRv4EHM model, elliptical orbits are described by two parameters, initial eccentricity, e 0 , and initial relativistic anomaly, ζ 0 , specified at an instantaneous orbital frequency ω 0 .Here, we present a new parametrization of the initial conditions, where e 0 and ζ 0 are specified at an orbitaveraged orbital frequency ω 0 .This new orbit-averaged initial conditions lead to smoother variations of e 0 and ζ 0 across parameter space, and as a consequence more efficient sampling of the parameter space.
The improvement in sampling efficiency due to the new initial conditions has also being accompanied with the de-velopment of SEOBNRv4EHM opt, a faster version of the SEOBNRv4EHM model.The SEOBNRv4EHM opt model combines a reduction of the absolute and relative tolerances of the Runge Kutta integrator from 10 −10 and 10 −9 , to 10 −8 and 10 −8 , with the use of the optimized Hamiltonian and integrator from Refs.[97,98].These modifications lead to a factor of ∼ 3 − 7 speed-up depending on the binary parameters.The reduction of the tolerances implies a reduction in accuracy of the SEOBNRv4EHM opt model in the corners of parameter space (i.e., high eccentricities and high spins), where the model usage is limited, because it is very sensitive to the attachment point of the inspiral and merger-ringdown EOB modes.The trade-off between accuracy and efficiency of the new SEOBNRv4EHM opt model, with a waveform evaluation time of O(100)ms, makes it a competitive model for use in parameter-estimation studies.
Given the accuracy and computational efficiency of SEOBNRv4EHM opt, we have performed a Bayesian inference study on mock signals and real GW events detected by the LVK collaboration.We have first investigated the quasicircular limit of the SEOBNRv4EHM opt model by computing the unfaithfulness against the quasi-circular SEOBNRv4HM model [96] for 4500 random configurations in the parameter space q ∈ [1, 50], spins χ 1,2 ∈ [−0.9, 0.9] and total masses M ∈ [20, 300]M ⊙ , at a dimensionless starting frequency of Mω = 0.023.The results show that SEOBNRv4EHM opt has a median unfaithfulness of 3.8 × 10 −5 and no cases with unfaithfulness > 1%, indicating that the quasi-circular limit is accurately recovered.Furthermore, we have performed a mocksignal injection into zero noise using SEOBNRv4 [114] as a signal, and SEOBNRv4E opt and SEOBNRv4 ROM [114] as templates.We have considered a configuration with mass ratio q = 4, total mass M = 90.08M⊙ and BH's dimensionless spins χ 1 = 0.5 and χ 2 = −0.1 defined at 20Hz.The recovery of quasi-circular parameters, like mass ratio, chirp mass or the effective-spin parameter by SEOBNRv4E opt agrees remarkably well with the ones from SEOBNRv4 ROM.While the initial eccentricity and relativistic anomaly measured by SEOBNRv4E opt are e 0 = 0.01 +0.02 −0.01 and ζ 0 = 3.09 +2.57−2.48 , which indicate that the signal is compatible with zero eccentricity.Therefore, SEOBNRv4EHM opt is able to correctly reproduce the zero-eccentricity limit, with an accuracy comparable to the underlying quasi-circular model SEOBNRv4HM.
Moving to the eccentric sector, we have studied with zeronoise injections, using the SEOBNRv4E opt as a signal and template, the impact of the initial conditions based on (e 0 , ζ 0 ) specified at ω 0 , and the new prescription, where (e 0 , ζ 0 ) are specified at ω 0 .For this study we have chosen a configuration with two different initial eccentricities e 0 = [0.1,0.2], mass ratio q = 3, initial relativistic anomaly ζ 0 = 1.2, total mass M = 76.4M⊙ and BH's dimensionless spins χ 1 = 0.5 and χ 2 = −0.1 defined at 20Hz.The results show that both the orbit-averaged and the instantaneous initial conditions are able to accurately recover the corresponding injected signal, however, the averaged wall clock-time of the runs using the instantaneous initial conditions is a factor 3 slower than the orbit-averaged initial conditions, due to the complicated structure of the likelihood across the (e 0 , ζ 0 ) parameter space in the case of the instantaneous initial conditions (see Fig. 2).As a consequence we adopt the orbit-averaged initial conditions as the default ones in the SEOBNRv4EHM and SEOBNRv4EHM opt models.
Moreover, we have investigated the impact of neglecting the radial-phase parameter, relativistic anomaly, in the previous model injections by starting the orbits of the templates at periastron (i.e., ζ 0 = 0).The posterior distributions of the orbit-averaged ICs show larger biases than the instantaneous ICs in the recovery of the quasi-circular parameters, for instance, 8% bias in the case of the chirp mass, while in terms of the eccentricity parameter the instantaneous ICs may develop multimodalities in the posteriors, as is the case of the injections considered here, indicating that the parametrization cannot adequately reproduce the injected signal.For the orbitaveraged ICs the eccentricity parameter is compatible with the injected value within the 90% credible intervals.This indicates that neglecting the radial-phase parameter when performing parameter estimation of eccentric signals can induce biases not only in the measurement of the eccentricity, but also in the estimation of other quasi-circular parameters like mass ratio or the spins due to the strong correlation of eccentricity with these parameters.
The accuracy of the SEOBNRv4EHM model in the eccentric case was investigated in Ref. [95] by computing the unfaithfulness of the model against a set of public eccentric NR waveforms from the SXS catalog [100,171].Here we have also validated the accuracy of the SEOBNRv4EHM opt model by performing a set of injections of synthetic NR signals into a network of LIGO-Virgo detectors at design sensitivity.We have injected into zero-detector noise three eccentric NR waveforms, SXS:BBH:1355, SXS:BBH:1359 and SXS:BBH:1363, corresponding to equal-mass, nonspinning configurations with initial eccentricities measured from the orbital frequency at first periastron passage of 0.07, 0.13 and 0.25, respectively.For these injections we choose a total mass M = 70M ⊙ , inclination ι = 0 and SNR= 20.The results are summarized in Fig. 8 and Table III.In order to compare the eccentricity from the NR waveforms and the SEOBNRv4EHM opt model, we have adopted a common definition of eccentricity and radial phase, based on the frequency of the (2,2)-mode [102], and used its efficient implementation in the open-source gw eccentricity Python package [103] to post-process the parameter estimation runs.We have found that the recovery of the parameters with SEOBNRv4E opt does not produce significant biases, and that the measurement of the GW eccentricity, e gw , and GW mean anomaly, l gw , is consistent with the injected values for the three injections considered.A more comprehensive Bayesian inference study will be required to assess the modeling inaccuracies and how they translate into biases in both eccentric and quasi-circular parameters.Here, new methods of inference such as machine learning techniques, like DINGO [184][185][186], may offer an alternative method to efficiently perform large-scale injections campaigns with moderate computational cost [181].
Besides injection studies, we have demonstrated that SEOBNRv4EHM opt can be used as a standard tool in Bayesian inference studies of real GW events.We have analyzed three GW events (GW150914, GW151226 and GW190521) detected by the LVK Collaboration in the first and third observing runs.The eccentricity measured for the three events is e GW150914 gw, 10Hz = 0.08 +0.09 −0.06 , e GW151226 gw, 20Hz = 0.04 +0.05 −0.04 , and e GW190521 gw, 5.5Hz = 0.15 +0.12 −0.12 .As a consequence, we do not find clear evidence of orbital eccentricity in any of the GW events considered, when using a non-precessing eccentric model with initial eccentricities e 0 ∈ [0, 0.3].For the GW150914 and GW151226 we have compared with the quasi-circular results from the GWTC-2.1 catalog [3], while for GW190521 with the results from Refs.[106,107] (the precessing-spin SEOBNRv4PHM [162] and IMRPhenomXPHM [174] models).For GW150914 we find good agreement with SEOBNRv4PHM due to the fact that the event is consistent with a nonspinning binary, while for GW151226 and GW190521, we find some discrepancies, which are likely due to the inclusion of spin-precession effects in the quasi-circular models, which are not included in the SEOBNRv4EHM opt model.This is a clear limitation of the SEOBNRv4EHM opt model as well as all the current existing inspiral-merger-ringdown eccentric waveform models.There is ongoing work [182] to include such effects in the new generation of SEOBNR models [140,141] developed within the new pySEOBNR infrastructure [183].
Regarding the analysis of real GW events, we plan in the future to extend the Bayesian inference study presented here, using the machine-learning code DINGO, to all the GW events detected during the third-observing run [181] in order to set constraints on the eccentricity of the observed population of BBHs.

ACKNOWLEDGMENTS
It is a pleasure to thank Mohammed Khalil for providing the expressions of the orbit-averaged orbital frequency used throughout the paper.We also would like to thank Sergei Ossokine and Harald Pfeiffer for helpful discussions about the implementation of the model and the initial conditions, as well as Marta Colleoni for useful comments to improve the manuscript.The computational work for this manuscript was carried out on the computer cluster Hypatia at the Max Planck Institute for Gravitational Physics in Potsdam.
This research has made use of data or software obtained from the Gravitational Wave Open Science Center (gwosc.org), a service of LIGO Laboratory, the LIGO Scientific Collaboration, the Virgo Collaboration, and KAGRA.LIGO Laboratory and Advanced LIGO are funded by the United States National Science Foundation (NSF) as well as the Science and Technology Facilities Council (STFC) of the United Kingdom, the Max-Planck-Society (MPS), and the State of Niedersachsen/Germany for support of the construction of Advanced LIGO and construction and operation of the GEO600 detector.Additional support for Advanced LIGO was provided by the Australian Research Council.This material is based upon work supported by NSF's LIGO Laboratory which is a major facility fully funded by the National Science Foundation.Virgo is funded, through the European Gravitational Observatory (EGO), by the French Centre National de Recherche Scientifique (CNRS), the Italian Istituto Nazionale di Fisica Nucleare (INFN) and the Dutch Nikhef, with contributions by institutions from Belgium, Germany, Greece, Hungary, Ireland, Japan, Monaco, Poland, Portugal, Spain.KA-GRA is supported by Ministry of Education, Culture, Sports, Science and Technology (MEXT), Japan Society for the Promotion of Science (JSPS) in Japan; National Research Foundation (NRF) and Ministry of Science and ICT (MSIT) in Korea; Academia Sinica (AS) and National Science and Technology Council (NSTC) in Taiwan.The initial conditions for eccentric orbits were derived in Sec.II.C of Ref. [95], given an initial (instantaneous) orbital frequency ω, eccentricity e, and relativistic anomaly ζ.In order to perform orbital evolutions starting from different points on (roughly) the same eccentric orbit, we derive an expression for the instantaneous frequency in terms of the orbit-averaged frequency.
We use the Keplerian parametrization Thus, we start with a given initial orbit-averaged frequency ω, eccentricity e, and relativistic anomaly ζ.Then, use Eq.(A7) to compute the instantaneous frequency ω(ω, e, ζ), and follow the same procedure as in Ref. [95] to obtain the initial conditions for the dynamical variables.

Figure 1 .
Figure 1.Top panel: Orbital frequency evolution of SEOBNRv4EHM for a configuration with mass ratio q = 3, dimensionless spins χ 1 = 0.3, χ 2 = 0.5, e 0 = 0.2, total mass 70M ⊙ at a starting frequency of 20Hz, and three different values of initial relativistic anomaly ζ 0 = [0, π/3, π].The dashed curves correspond to the specification of (e 0 , ζ 0 ) at an instantaneous orbital frequency, ω 0 , while the solid lines correspond to the initial conditions specified at an orbit-averaged orbital frequency ω 0 .Bottom panel: Same configuration as in the upper panel, but fixing ζ 0 = π/3, and using three different values of initial eccentricity, e 0 = [0.01,0.1, 0.2], for both types of initial conditions.The horizontal black solid line in both panels indicates a dimensionless frequency of Mω 0 = 0.0216, corresponding to 20Hz at 70M ⊙ .

Figure 2 .
Figure 2. Zero-noise injection of a SEOBNRv4E waveform with initial eccentricity e 0 = 0.1, total mass 65M ⊙ and spins χ 1 = 0.3, χ 2 = 0 at a starting frequency of 20Hz, and recovering with SEOBNRv4E with all the parameters fixed to injected values except for the eccentric ones (e 0 , ζ 0 ).In all the plots we consider 5000 points randomly distributed in the parameter space ζ 0 ∈ [0, 2π] and e 0 ∈ [0, 0.3].In the top panel, the relativistic anomaly is fixed to ζ 0 = 0, and only e 0 is sampled.The blue curve corresponds to using initial conditions based on the instantaneous orbital frequency, ω 0 , while the green curve corresponds to using the orbit-averaged orbital frequency, ω 0 .The mid and bottom plots correspond to an injected value of ζ 0 = 1, and the use of ω 0 and ω 0 in the initial conditions, respectively.In both mid and bottom panels, e 0 and ζ 0 are sampled, and each point in parameter space is color-coded by its log-likelihood value.

Figure 6 .
Figure 6.Inverse mass ratio, chirp mass, effective spin parameter, luminosity distance, initial eccentricity and relativistic anomaly posterior distributions for a synthetic quasi-circular BBH signal injection using the SEOBNRv4 model, and recovering the parameters with SEOBNRv4E opt (blue) and SEOBNRv4 ROM (orange).The dashed vertical lines indicate the 90% credible intervals.The solid vertical lines correspond to the injected parameters, which are also shown in TableI.

2 Figure 7 .
Figure 7. Top row: Violin plots for the posterior distributions of the inverse mass ratio, initial eccentricity and relativistic anomaly for a synthetic BBH signal injection using the SEOBNRv4E opt model with e 0 = 0.1 and ζ 0 = 1.2 at 20Hz.The posteriors are computed using the SEOBNRv4E opt model with orbit-averaged ICs (blue), instantaneous ICs (orange), setting ζ 0 = 0 for the orbit-averaged (brown) and instantaneous (pink) ICs, and for the zero-eccentricity case e 0 = 0 (gray).The cases where ζ is fixed during sampling (ζ 0 = 0) have been placed on the right side of the x-axis to ease the visualization of the results.The dashed vertical lines indicate the 90% credible intervals.The solid vertical lines correspond to the injected parameters, which are also shown in TableI.Bottom row: Same as in the upper row but for an injected signal with initial eccentricity e 0 = 0.2 at 20Hz.

Figure 8 .
Figure 8. 2D and posterior distributions for some relevant parameters from the equal mass nonspinning synthetic BBH signals with total mass 70M ⊙ .The signal waveforms correspond to the NR waveforms from the public SXS catalog SXS:BBH:1355, SXS:BBH:1359 and SXS:BBH:1363 with GW eccentricities e gw = 0.05, 0.1, 0.25 and GW mean anomalies e gw = 0.05, 0.1, 0.25 defined at 20Hz, respectively.The other parameters are specified in Table III.In the 2D posteriors the solid contours represent the 90% credible intervals and black dots show the values of the parameters of the injected signal.In the 1D posteriors they are represented by dashed and solid vertical lines, respectively.The parameter estimation is performed with the SEOBNRv4E opt model.Left: Chirp mass and effective-spin parameter.Middle: Initial eccentricity and relativistic anomaly at 20Hz.Right: GW eccentricity and GW mean anomaly at 20Hz.

Figure 9
Figure 9. Effective-spin parameter, chirp mass, GW eccentricity and GW mean anomaly parameters inferred for the real GW events analysed with SEOBNRv4E opt and SEOBNRv4EHM opt.Comparisons are presented with the SEOBNRv4PHM (when available) and the IMRPhenomXPHM models from the GWTC-2.1 catalog[3], except for GW190521 where samples are from Refs.[106,107].The waveform starting frequency for GW150914 is f start = 10Hz, for the low total mass event GW151226 is f start = 20Hz, in order to reduce the computational cost, while for GW190521 f start = 5.5Hz in order to have the mode content in band at the likelihood minimum frequency ( f min = 11Hz).
Appendix A: Derivation of the orbit-averaged orbital frequency

•
Waveform multipoles: The GW multipoles are composed by two main parts: the inspiral-plunge multipoles h

Table I .
Injected and median values of the posterior distributions for the synthetic injection with the SEOBNRv4 model, recovered with SEOBNRv4E opt and SEOBNRv4 ROM.The median values also report the 90% credible intervals.The binary parameters correspond to the total mass M, chirp mass M, inverse mass ratio 1/q, effectivespin parameter χ eff , initial eccentricity e 0 , initial relativistic anomaly ζ 0 , angle between the total angular momentum and the line of sight θ JN , luminosity distance d L , coalescence phase ϕ ref and the network matched-filtered SNR for LIGO-Hanford/Livingston and Virgo detectors ρ N mf .
Table II the results of sampling with the initial eccentricity set to zero.For the Table II.Injected and median values of the posterior distributions for two synthetic injections with the SEOBNRv4E opt model with initial eccentricities e 0 = [0.1,0.2]

Table III .
Injected and median values of the posterior distributions for three synthetic NR injections with the same quasi-circular parameters and different initial eccentricities, and recovered with SEOBNRv4E opt.The median values also report the 90% credible intervals.The binary parameters correspond to the total mass M, chirp mass M, inverse mass ratio 1/q, effective-spin parameter χ eff , initial eccentricity e 0 , initial relativistic anomaly ζ 0 , angle between the total angular momentum and the line of sight θ JN , luminosity distance d L , coalescence phase ϕ ref and the network matched-filtered SNR for LIGO-Hanford/Livingston and Virgo detectors ρ N mf .At the bottom of the table the injected and measured GW eccentricity, e gw , and GW mean anomaly, l gw , are also reported.

Table V .
Settings and evaluation time for the different parameter estimation runs on real GW events.The mode content indicates the use of SEOBNRv4E opt ((2, |2|