# A higher-multipole gravitational waveform model for an eccentric binary black holes based on the effective-one-body-numerical-relativity formalism

aa r X i v : . [ g r- q c ] F e b A higher-multipole gravitational waveform model for an eccentric binary black holesbased on the eﬀective-one-body-numerical-relativity formalism

Xiaolin Liu, Zhoujian Cao ∗ ,

1, 2, † and Zong-Hong Zhu Department of Astronomy, Beijing Normal University, Beijing 100875, China School of Fundamental Physics and Mathematical Sciences,Hangzhou Institute for Advanced Study, UCAS, Hangzhou 310024, China

We construct a new factorized waveform including ( l, | m | ) = (2 , , (2 , , (3 , , (4 ,

4) modes basedon eﬀective-one-body (EOB) formalism, which is valid for spinning binary black holes (BBH) ingeneral equatorial orbit. When combined with the dynamics of

SEOBNRv4 , the ( l, | m | ) = (2 ,

2) modewaveform generated by this new waveform can ﬁt the original

SEOBNRv4 waveform very well in thecase of a quasi-circular orbit. We have calibrated our new waveform model to the Simulating eXtremeSpacetimes (SXS) catalog. The comparison is done for BBH with total mass in (20 , M ⊙ usingAdvanced LIGO designed sensitivity. For the quasi-circular cases we have compared our (2 ,

2) modewaveforms to the 281 numerical relativity (NR) simulations of BBH along quasi-circular orbits. All ofthe matching factors are bigger than 98%. For the elliptical cases, 24 numerical relativity simulationsof BBH along an elliptic orbit are used. For each elliptical BBH system, we compare our modeledgravitational polarizations against the NR results for diﬀerent combinations of the inclination angle,the initial orbit phase and the source localization in the sky. We use the the minimal matching factorrespect to the inclination angle, the initial orbit phase and the source localization to quantify theperformance of the higher modes waveform. We found that after introducing the high modes, theminimum of the minimal matching factor among the 24 tested elliptical BBHs increases from 90% to98%. Following our previous

SEOBNRE waveform model, we call our new waveform model

SEOBNREHM .Our

SEOBNREHM waveform model can match all tested 305 SXS waveforms better than 98% includinghighly spinning ( χ = 0 .

99) BBH, highly eccentric ( e ≈ .

15) BBH and large mass ratio ( q = 10)BBH. I. INTRODUCTION

LIGO [1] and VIRGO [2] have achieved the detection ofgravitational waves (GW) [3, 4]. The gravitational wavesfrom compact binary systems are the only detected GWsources, and they will be also the most likely gravitationalwave events to be detected in the near future. For thereported more than 50 events in the O1/O2/O3a data,the accurate waveform model of the gravitational wavesignal has played an important role. Based on the mostaccurate numerical relativistic (NR) simulation of BBHs,several inspiral-merger-ringdown (IMR) waveform mod-els have been constructed. Eﬀective-one-body (EOB) ap-proach [5] is a widely used method. The

EOBNR model isbased on EOB method and resums the results of postNewtonian (PN) approximation. Through the calibra-tion to NR waveforms EOBNR obtains an accurate modelof the time domain BBH waveform. Currently, the

EOBNR model [6] has gradually developed to

SEOBNRv4 model[7–10], which can accurately describe the BBH gravita-tional wave waveform of general spin-aligned BBH mov-ing along an equatorial quasi-circular orbit.In the source frame the gravitational wave can be de-composed by a set of spin weighted -2 spherical harmonicsrespect to diﬀerent directions. For a two-body system,source frame with z direction pointing to the orbital an- ∗ corresponding author † Zhoujian Cao: [email protected] gular momentum is the most preferred frame choice. Insuch a frame the quadrupole components contribute mostof the gravitational wave energy. In another word, the( l, | m | ) = (2 ,

2) spin weighted -2 spherical harmonic com-ponents are the dominant modes. Consequently othermodes than ( l, | m | ) = (2 ,

2) are called higher modes.When the diﬀerence between the two components masses,the spin precession and/or the orbit eccentricity becomeslarger, the impacts of higher-order modes of gravitationalwaves appears.For the quasi-circular orbit case, there are alreadymany models describing the waveforms of higher-ordermodes, which include

IMRPhenomHM , SEOBNRv4HM [11, 12]and others [13]. So far, most works about the detectionof gravitational wave signal and parameter estimation donot consider the inﬂuence of the orbit eccentricity (butsee [14–17]). This is because people suspect that the orbitof a BBH system has become circular before it enter theLIGO frequency band. Recently the event GW190521[18] makes people rethink the orbit eccentricity problemfor ground based detector [19].There are many works investigating the binary sys-tem moving along an eccentric orbit. The properties ofenergy and angular momentum diﬀusion of eccentric bi-nary system was seminally studied by Peters [20]. Thepost-Newton (PN) waveforms for gravitational radiationhave also been widely studied [21–23]. In Ref. [24], theauthors assumed a small eccentricity condition and gota post-circular (PC) waveform model based on the lowPN order waveform in frequency domain. Later the PCmodel was improved by a phenomenological method anddeveloped into an enhanced post-circular (EPC) model[25], which recovers the TaylorF2 model in quasi-circularcases. The x-model [26] is a PN waveform with variable x ≡ ( M ω ) / inspired by numerical relativity. In [27] Ec-centric, Nonspinning, Inspiral-Gaussian-process MergerApproximant ( ENIGMA ) model was proposed, which com-bines analytical and NR results using machine learn-ing algorithms. There are also some studies of ellipticalBBH systems based on EOB formalism. In Ref. [28], anEOB eccentric waveform model was constructed basedon an adiabatic approximation rather than NR calibra-tion. We proposed

SEOBNRE [29] waveform model basedon

SEOBNRv1 through adding elliptical correction termsto the EOB factorized waveform.

SEOBNRE model hasbeen validated to NR catalog in [30]. The authors of[31] modiﬁed

TEOBiResumS SM , a highly NR-faithful EOBmultipolar waveform model for quasi-circular orbits, toincorporate orbit eccentricity. They found that such awaveform model can also work for hyperbolic encounterBBHs [32, 33].The existing

SEOBNRE waveform model admits only(2 ,

2) mode. In this paper, we resum the higher modes ofgravitational waveforms for an equatorial orbit and com-bine them with the existing higher PN order EOB fac-torized waveforms to obtain a waveform model of highermodes for

SEOBNRE . At the mean time some improve-ments on

SEOBNRE model have been implemented. Forconvenient reference we call the new waveform modelproposed in the current work

SEOBNREHM to distinguishthe previous

SEOBNRE model. In Sec. II, we will showthe construction details of our new waveform. Then inSec. III we combine this new waveform with the dynam-ical system of

SEOBNRv4 to get

SEOBNREHM model. Thenwe verify that the ( l, m ) = (2 ,

2) mode of our new modelwaveform can match the waveform of

SEOBNRv4 well inthe case of quasi-circular orbits. Meanwhile when com-paring ( l, m ) = (2 ,

2) against NR catalogs, for both quasi-circular orbits cases and elliptical orbits cases, the ﬁttingfactors are bigger than 98%. In Sec. IV we test the per-formance of the higher modes in both quasi-circular andelliptical cases, we ﬁnd that the matching factor betweenour model waveform and the NR waveform is bigger than98%. Based on our new waveform model we analyze thepower fraction of higher modes in the Sec. V. A summaryis given in the last section. We adopt the geometric units c = G = 1 throughout this paper. II. EFFECTIVE-ONE-BODY MULTIPOLARWAVEFORM MODEL FOR ECCENTRICBINARY BLACK HOLEA. Eﬀective-one-body dynamics

For a binary system with component masses m and m , the EOB conservative dynamics for the binary can be described by the Hamiltonian [5] H EOB = M s ν (cid:18) H eﬀ µ − (cid:19) , (1)where µ = m m /M is the reduced mass, ν = µ/M is the symmetric mass ratio and M = m + m is thetotal mass. For a spin-aligned BBH we are concerned,we denote the spin angular momentums by S i = χ i m i ˆ L ,where χ i is dimensionless spin parameters and ˆ L = L / | L | is the direction of the orbital angular momentum. Here i = 1 , m ≥ m . The explicit expres-sion of the eﬀective Hamiltonian H eﬀ has been derived in[34, 35]. In the EOB framework, the spin-aligned orbitalevolution is described by orbit phase φ , radial distance r , and the corresponding canonical momentums p φ and p r . The other four parameters ( K, d SO , d SS , δ ) canbe determined through calibration against numerical rel-ativity simulations [10]. Here we adopt the conservativedynamics, i.e., the Hamiltonian from SEOBNRv4 model.The gravitational waves would take away energy andangular momentum. Such eﬀect can be described by theradiation-reaction force [10]. For elliptical orbits, we setthe initial conditions of the orbital evolution as follows ∂H∂ ˆ p φ = πf , (2) ∂H∂r = − e r , (3)where ˆ p = p/µ , and e ∈ (0 ,

1) is the initial eccentricity.When e = 0, the orbit is circular. When e > B. Eﬀective-one-body gravitational waveformmodes

The direction of the observer respect to the gravita-tional wave source frame corresponds to the angles ( ι, ϕ c )where the inclination angle ι represent the angle betweenthe line of sight and the direction of the angular momen-tum of the binary system, and ϕ c denotes coalescenceorbital phase. Consequently the gravitational wave canbe decomposed by spin-weighted -2 spherical harmonicsas h ( ι, ϕ c ; t ) = h + ( ι, ϕ c ; t ) − ih × ( ι, ϕ c ; t )= ∞ X l =2 l X m = − l − Y lm ( ι, ϕ c ) h lm ( t ) (4)Following [37], we express the full inspiral-merger-ringdown EOB waveform modes as h IMR lm ( t ) = ( h insp lm ( t ) , t < t lm match h RD lm ( t ) , t > t lm match (5)The deﬁnition of t lm match , and the ringdown waveformmodes h RD lm can be found in [12].The quasi-circular inspiral-plunge modes h insp lm can bewritten in a factorized form [38, 39] h insp lm = h F lm N lm , (6)where h F lm is a resummation of post-Newtonian waveformfor quasi-circular orbit [40–42]. The non-quasi-circular(NQC) term N lm is used to adjust EOB waveform modesin the late inspiral and merger phase. The Ref. [10] hasdetermined N as N = " (cid:18) ˆ p r ∗ r Ω (cid:19) a h + a h r + a h r / ! × exp (cid:26) i ˆ p r ∗ r Ω (cid:16) b h + b h ˆ p r ∗ (cid:17)(cid:27) (7)where Ω = ˙ φ and ˆ p r ∗ is the tortoise radial angular mo-mentum. The coeﬃcients a h and b h will be deter-mined by ﬁtting amplitude, frequency and their deriva-tives of h insp lm at t lm match to the values predicted by NRsimulations.In a quasi-circular case, the prefactor ˆ p r ∗ / ( r Ω) is smallduring early inspiral. Consequently N is about 1 andthe NQC does not take eﬀect. For an elliptical orbitˆ p r ∗ / ( r Ω) may be quite large near the perihelion point,even for early inspiral. In order to make sure the NQCterm takes eﬀect only near merger, we introduce a win-dow function to limit the scope of this term as N lm ( t ) = " (cid:18) W ( t ) ˆ p r ∗ r Ω (cid:19) a h lm + a h lm r + a h lm r / ! × exp (cid:26) W ( t ) i ˆ p r ∗ r Ω (cid:16) b h lm + b h lm ˆ p r ∗ (cid:17)(cid:27) , (8)where the window function W ( t ) reads W ( t ) = 11 + exp { [ − ¯ ω ( t − t w )] } . (9)The t w denotes the location of the window. Since p r = 0at perihelion and aphelion, p r can be used as an indica-tor where the test particle locates respect to the orbit.We choose t w equals to the time after which p r < t lm match to the time where p r = 0 and set it as our t w . When e = 0 which corresponds to the quasi-circularcase, the p r < t w = −∞ under these circumstances, which means W ( t ) ≡

1, and our NQC recovers original NQC termin

SEOBNRv4 model. The another parameter ¯ ω controlsthe width of the window function. In principle this ad-justable parameter should be determined by the calibra-tion to the numerical relativity simulations. In the cur-rent work we did not do such extensive investigation. Weinstead choose ¯ ω = ln 23 M based on the comparison to the24 elliptical BBH simulations of SXS. Diﬀerent choicesof ¯ ω may result in a diﬀerent matching factor. But thematching factor diﬀerence is less than 1%. C. The factorized gravitational waveform modesfor an eccentric orbit

The factorized resummation of waveform modes pro-posed in [38, 39] are composed of ν = 0 circular wave-form terms till 3PN and ν = 0 circular waveform termstill 5PN. In this work, we replace the 2PN order cir-cular parts with the 2PN general equatorial form [29].Compare to the original SEOBNRE waveform model wheredirect PN waveform for 2PN terms is used, this resumedwaveform can hopefully improve the performance. Andas we will show in the following part of the current paper,this resummation form improves quite much.The factorized waveform modes read [38, 39] h F lm = h ( N,ǫ ) lm ˆ S ( ǫ )eﬀ T lm e iδ lm f lm (10)where ǫ = 0 when l + m is even, otherwise ǫ = 1. Theamplitude correction factor f lm is resummed by f lm =( ρ lm ) l when l is even and f lm = ( ρ NS lm ) l + f S lm when l isodd. The explicit expression of these terms can be foundin [38, 39].Quasi-circular Newtonian waveform modes h ( N,ǫ ) lm arecontrolled by the parameter x ≡ ( M Ω) / . We generalizethem to eccentric orbit cases [21, 41, 42, 44] h Nlm = M νR n ( ǫ ) lm c l + ǫ ( ν ) Y l, − m ( π , φ h )ˆ h lm (11)ˆ h = y − y y iy y (12)ˆ h = y y (13)ˆ h = − iy − y y + 49 iy y + 23 iy y + 2 y y y (14)ˆ h = 7 y − y y + 5164 y y + 2732 iy y y + 3 y

32+ 3 y

32 + 38 iy y − y y − iy y (15)The parameters y , y and y are deﬁned as y = r r h , y = r h · v h r h , y = r h Ω h , (16)where the subscript h means harmonic coordinates, usedto distinguish it from the coordinates in EOB dynam-ics. Corresponding to the above harmonic coordinates,we can also use EOB coordinates to deﬁne three newvariables x = r r , x = r · v r , x = r Ω . (17)Following [45], we apply a 2PN gauge transformationfrom the harmonic coordinates ( y , y , y , φ h ) to the EOB

20 40 60 80 100 120 140 160 180 200M tot [M ⊙ ]10 −4 −3 −2 −2 −2 −1 - FF ID,(q,χ ,χ )

20 40 60 80 100 120 140 160 180 200M tot [M ⊙ ]10 −4 −3 −2 −2 −2 −1 - FF ID,(q,χ ,χ ) FIG. 1: The matching factor FF of the ( l, m ) = (2 ,

2) mode waveform for the spin-aligned BBHs between the new waveformmodel

SEOBNREHM proposed in the current work and the SXS catalog [43]. The Advanced LIGO designed PSD is used in thisplot. We highlight the ones whose matching factor is less than 99%. Left panel: 281 quasi-circular cases. Right panel: 24eccentric cases. χ A χ e ff −2.00−1.95−1.90−1.85−1.80 l o g ( − FF m i n ) FIG. 2: Matching factors of the ( l, m ) = (2 ,

2) mode waveform for the eccentric spin-aligned BBH between the new waveformmodel

SEOBNREHM proposed in the current work and the SXS catalog [43]. The color means the minimal matching factor respectto the total mass in the range (20 , M ⊙ . This ﬁgure corresponds to the right panel of the Fig. 1. coordinates ( x , x , x , φ ) (see the Appendix. A). We re-place the Newtonian waveform (11) with h ( N,ǫ,e ) lm = M νR n ( ǫ ) lm c l + ǫ ( ν ) Y l − ǫ, − m (cid:16) π , φ (cid:17) x l + ǫ . (18)We keep original factors ˆ S ( ǫ )eﬀ , T lm and e iδ lm unchanged,while change f lm in the following way. We introduce aTaylor expansion form ζ lm = c h lm + c h lm x + c h lm x + ... (19)We list the explicit expression of c h lm , , in the Ap-pendix. B. The above Taylor expansion has poor con-vergence near merger. So we resum it by (0 ,

4) Padeapproximation f elm = ( P [ ζ lm ] , l is even P (cid:2) ζ NS lm (cid:3) + ζ S lm , l is odd (20) Here NS (Non-Spin) means dropping out the spin terms,while S (Spin) means solely spin terms. Explicitly P means P [ A + A x + A x + A x + A x ]= A B x + B x + B x + B x (21) B = − A A (22) B = A − A A A (23) B = − A + 2 A A A − A A A (24) B = A − A A A + A A + 2 A A A − A A A . (25)We replace the 2PN resummed parts in the SEOBNRv4HM waveform with the new 2PN form. The full factorizedwaveform modes are h F lm = h (F ,e ) lm + h (F ,c ) lm (26) h (F ,e ) lm = h ( N,ǫ,e ) lm ˆ S ( ǫ )eﬀ T lm e iδ lm f elm (27) h (F ,c ) lm = h ( N,ǫ ) lm ˆ S ( ǫ )eﬀ T lm e iδ lm (cid:0) f lm − f lm (cid:1) , (28)where f lm is f lm = (cid:0) ρ lm (cid:1) l , l is even (cid:16) ρ NS lm (cid:17) l + f S lm , l is odd (29)The index 2PN represents 2PN parts. III. THE (2,2) MODE PERFORMANCE OF

SEOBNREHM

Given two waveforms h ( t ) and h ( t ), we deﬁne match-ing factor (FF) by inner product h h | h i with respect toa detector noise h h | h i = 4 R Z f max f min ˜ h ( f )˜ h ∗ ( f ) S n ( f ) d f FF ≡ max t c ,ϕ c h h | h i p h h | h i h h | h i (30)where S n ( f ) is the one-sided power spectral density(PSD) of the detector noise. The detector frequencyband ( f min , f max ) = (10Hz , A. Comparison with numerical relativitywaveforms

In this subsection we compare our new (2 ,

2) modewaveforms to the numerical relativity results catalog [43].Here we use the BBH numerical waveform from Simu-lating eXtreme Spacetime (SXS) generated by SpectralEinstein code (

SpEC ) [47, 48]. Since the numerical rel-ativity results can be scaled to any total mass, we plotthe comparison results respect to diﬀerent total massesin the left panel of the Fig. 1 for the 281 quasi-circularBBHs listed in the Appendix. C. Each line corresponds toa parameter combination ν , χ z and χ z . Most match-ing factors are greater than 99%. Only 10 parametercombinations among the 281 combinations in all admitmatching factor in the range (98%,99%) for some totalmasses. These 10 cases have been highlighted out in theleft panel of the Fig. 1. We ﬁnd out the minimal match-ing factor FF respect to the total mass in the range M ∈ [20 , M ⊙ for each line and list them respect tothe parameter combinations in the Appendix. C. -1.0-0.50.00.51.0 χ A ν -1.0-0.50.00.51.0 χ e ff −2.0−1.8−1.6−1.4−1.2−1.0−0.8 l o g ( − FF m i n ) FIG. 3: Matching factors of the ( l, m ) = (2 ,

2) mode wave-form for the quasi-circular spin-aligned BBH between the newwaveform model proposed in the current work and

SEOBNRv4 .Uniformly distributed 10000 points in the parameter space ν , χ z and χ z are chosen. The plot convention is the same tothat of the Fig. 2. For elliptical cases, we ﬁrstly estimate the initial eccen-tricity based on the ﬁrst orbit of NR simulation through e NR = r i ¨ r | r = r i (31)where subscript i means the initial perihelion. Regardingto our SEOBNREHM model we set the reference frequencyto the frequency when the eccentricity is estimated bythe Eq. (31) and search the initial eccentricity whichmaximizes the matching factor between the our model(2,2) mode waveform and the numerical relativity wave-form. We should notice that the eccentricity deﬁned bythe Eq. (31) equivalents to the eccentricity deﬁned in theEq. (3) only at Newtonian order. Similar analysis hasbeen done before [29, 30]. This fact explains part of thedeviations of the

SEOBNREHM estimated eccentricity e EOB from the NR estimated one e NR . We have listed the re-sults in the Table. I. We plot the (2,2) mode comparisonresults in the right panel of the Fig. 1. At the mean timewe ﬁnd out the minimal matching factor for each lineand plot the results in the Fig. 2. The minimal matchingfactors are also listed in the Table. I. TABLE I: SXS Numerical relativity waveforms of eccentric BBH used for the waveform validation in the current work. IDcorresponds to the id number in the SXS catalog. f o is the orbital frequency at the ﬁrst perihelion from the NR simulation. e NR is the estimated initial eccentricity by the Eq. (31). e EOB is the initial eccentricity at reference frequency Mf = 2 Mf o ofour SEOBNREHM model which maximizes the matching factor between the

SEOBNREHM waveform and the NR waveform for (2,2)mode. Matching factors between the

SEOBNEHM waveform model and the SXS Numerical relativity waveforms are also listedhere. FF means the minimal matching factor for the (2,2) mode respect to the total mass. FF and FF min are respectivelythe averaged matching factor and the minimal matching factor respect to the three angles ( κ, ι, ϕ ) (check the Sec. IV for detailexplanation).ID q χ χ Mf o e NR e EOB FF FF FF min . × − .

76% 99 .

85% 99 . . × − .

51% 99 .

63% 99 . . × − .

46% 99 .

61% 99 . . × − .

58% 99 .

60% 99 . . × − .

53% 99 .

56% 99 . . × − .

25% 99 .

15% 99 . . × − .

14% 99 .

15% 99 . . × − .

43% 98 .

30% 98 . . × − .

24% 98 .

21% 98 . . × − .

65% 99 .

59% 99 . . × − .

64% 99 .

64% 99 . . × − .

67% 99 .

56% 99 . . × − .

44% 99 .

53% 99 . . × − .

54% 99 .

49% 99 . . × − .

87% 98 .

69% 98 . . × − .

48% 98 .

54% 98 . . × − .

74% 99 .

65% 99 . . × − .

69% 99 .

78% 99 . . × − .

65% 99 .

60% 99 . . × − .

23% 99 .

08% 98 . . × − .

77% 99 .

23% 98 . . × − .

62% 99 .

21% 98 . . × − .

69% 99 .

33% 98 . . × − .

29% 98 .

44% 98 . B. Comparison with

SEOBNRv4

SEOBNRv4 waveform model is one of the most widelyused waveform model by LIGO data analysis. In thissubsection, we compare our new (2 ,

2) mode waveformsto the corresponding ones given by

SEOBNRv4 [10] forquasi-circular spin-aligned BBHs. We uniformly chose10000 samples at random from the parameter space ofsymmetric mass ratio ν ∈ [0 , . χ z ∈ [ − . , . χ z ∈ [ − . , . ν , χ z and χ z we scan the total mass in therange M ∈ [20 , M ⊙ to search the minimal match-ing factor FF min respect to the total mass. We plot theFF min respect to diﬀerent mass ratio and black hole spinin the Fig. 3. For convenient comparison to the ﬁguresin [10] we have used variables χ A = ( χ − χ ) / χ eﬀ = ( m χ + m χ ) /M in the plot.The above comparison results indicate that our newwaveform model can recover SEOBNRv4 very well inthe case of quasi-circular spin-aligned BBH. Especiallythe ﬁtting factor between our new waveform model

SEOBNREHM and the

SEOBNRv4 waveform model is biggerthan 99.9% for spinless BBHs. At the mean time we ﬁndthat the worst ﬁtting factor respect to the black holes’spin depends on the mass ratio. In the Fig. 4 we plot suchdependence between min(FF min ) ≡ min χ z ,χ z (FF min )and the mass ratio ν ( q ≡ m m ≥ ν & . q .

20) the ﬁtting factor is bigger than 99%; when ν & .

01 ( q .

60) the ﬁtting factor is bigger than 90%.

C. Comparison with

SEOBNRE

In [29, 30] we have constructed

SEOBNRE waveformmodel which extends the waveform model

SEOBNRv1 forspin-aligned circular BBH to spin-aligned eccentric cases.

SEOBNRE waveform model works only for (2,2) modewaveform. Our current

SEOBNREHM waveform model notonly extends

SEOBNRE to high modes but also improvesthe accuracy for (2,2) mode due to the resummation skilland the calibration to

SEOBNRv4 waveforms. In this sub-section we compare our new (2,2) mode waveform to ourprevious

SEOBNRE model.

TABLE II: Matching factors between the models waveform and the SXS numerical relativity waveform for the 5 mostlymismatched cases between

SEOBNREHM and SXS waveforms among the 281 quasi-circular simulations listed in the Appendix. C.ID corresponds to the id number in the SXS catalog. FF and FF min are respectively the averaged matching factor and theminimal matching factor respect to the three angles ( κ, ι, ϕ ). ‘E’ means model

SEOBNREHM and ‘C’ means model

SEOBNRv4HM .ID q χ χ FF( C ) FF min ( C ) FF( E ) FF min ( E )1422 7.95 -0.80 -0.46 98.87% 96.95% 98.75% 96.97%1424 6.46 -0.66 -0.80 99.03% 97.94% 98.87% 96.98%1425 6.12 -0.80 +0.67 99.22% 97.12% 98.98% 97.06%1427 7.41 -0.61 -0.73 99.39% 98.32% 99.27% 98.26%1428 5.52 -0.80 -0.70 98.95% 97.18% 98.69% 96.91% −3 −2 −1 − ( FF m i n ) m i n

100 10 4 2 1q

FIG. 4: The dependence of the minimal FF min respect toblack holes’ spin on the mass ratio, where the subscript minin FF min means the minimum respect to the total mass rangedin (20,200)M ⊙ . The minimal operation respect to black holes’spin is done based on the data shown in the Fig. 3. In order to understand the diﬀerence between thewaveform model

SEOBNRE and

SEOBNREHM , we use the nu-merical relativity waveform as a reference. Accordinglywe calculate the waveforms for the parameters combi-nation shown in the Fig. 2. We show the comparisonresults in the Fig. 5. Firstly we note that some pointsshown in the Fig. 2 are missing in the Fig. 5. This isbecause

SEOBNRE model breaks down for such high spincases. In addition the matching behavior to the NR wave-form of

SEOBNREHM improves in general than

SEOBNRE . Inmost cases the improvement is less than 1%. For sev-eral cases including SXS:BBH:0292, SXS:BBH:1439 andSXS:BBH:1432, the improvement is more than 10%.

IV. THE MULTIPOLAR WAVEFORMPERFORMANCE OF

SEOBNREHM

SEOBNRv4HM waveform model is one of the most widelyused higher modes waveform models used in LIGO dataanalysis. Here we calibrate our new

SEOBNREHM wave- form model to

SEOBNRv4HM for quasi-circular BBHs. Theeccentric BBH waveform model

TEOBiResumS SM also in-volves high modes. But unfortunately we can not accesssuch waveform data.Since the waveform of higher modes is very weak, thenumerical error of NR results is large for higher modes.In order to treat this issue, previous works on highermodes used combined waveform comparison to numericalrelativity simulations. In the current work we adopt thesame strategy.

A. Comparison with numerical relativitywaveforms

The gravitational wave strain detected by a detectorcan be described as [12], h ( ι, ϕ c , t c , θ , α, δ, ψ ; t ) = F + ( α, δ, ψ ) h + ( ι, ϕ c , t c , θ ; t ) + F × ( α, δ, ψ ) h × ( ι, ϕ c , t c , θ ; t )(32) F + ( α, δ, ψ ) ≡

12 (1 + cos α ) cos 2 δ cos 2 ψ − cos α sin 2 δ sin 2 ψ, (33) F × ( α, δ, ψ ) ≡ + 12 (1 + cos α ) cos 2 δ sin 2 ψ + cos α sin 2 δ cos 2 ψ, (34)where θ = ( m , m , χ , χ , e ) is the intrinsic parame-ters of a BBH. ( α, δ ) denotes the sky location of theGW source, and ψ is the polarization angle. Thesethree angles determine the antenna pattern functions F + , F × [49, 50]. ( ι, ϕ c ) corresponds to the direction ofthe observation line respect to the gravitational sourceframe (check the Eq. (4)) and t c is the coalescence time.Following [12], we deﬁne the angle κ ∈ [0 , π ) (eﬀectivepolarization [51]) and the amplitude A as e iκ ( α,δ,ψ ) = F + ( α, δ, ψ ) + iF × ( α, δ, ψ ) A ( α, δ ) (35) A ( α, δ ) = q F ( α, δ, ψ ) + F × ( α, δ, ψ ) . (36) χ A χ e ff −2.00−1.95−1.90−1.85−1.80 l o g ( − FF m i n ) χ A χ e ff −2.00−1.95−1.90−1.85−1.80 l o g ( − FF m i n ) FIG. 5: Matching factors of the ( l, m ) = (2 ,

2) mode waveform for the eccentric spin-aligned BBH for the

SEOBNRE waveformmodel. Top row is for the matching factors between the

SEOBNRE waveform model and the SXS catalog [43]. The bottom rowis for the matching factors between the

SEOBNRE waveform model and the

SEOBNREHM waveform model. The plot convention isthe same to the Fig. 2.

Firstly we show one comparison example for eccen-tric BBH between our

SEOBNREHM waveform and the NRwaveform. The Fig. 6 shows the waveform comparisonfor the most eccentric BBHs in the SXS catalog whichcorresponds to SXS:BBH:1374. We can see the consis-tence of the modes amplitude and of the real part ofeach mode between the

SEOBNREHM waveforms and theNR waveforms. Regarding to the combined gravitationalwave strain, the higher modes than (2,2) improves thematching behavior near merger.When comparing a waveform with higher modes, wedeﬁne the matching factor as [12]FF( κ, ϕ, ι, θ ) = max t c ,κ ′ ,ϕ ′ h h | h i p h h | h i h h | h i (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) κ κ ′ κ κϕc ϕ ′ ϕc ϕι ι ι θ θ θ . (37)In this work, the sky location of the source is set to bethe same for EOB/NR waveforms. For a given intrinsic parameter θ the matching factor is a function of angles( κ, ϕ, ι ). As an example we show the matching factorbehavior corresponding to SXS:BBH:1374 respect to theangles ( ϕ, ι ) in the Fig. 7. Along with cos ι decreasesfrom 1 to 0, the orientation of the orbital plane is fromface-on to edge-on with respect to the observer. For thespin-aligned BBHs, the cos ι < ι > ι = 0always results in better matching than ι = π/ SEOBNREHM model to NR waveform. This feature is truefor both with and without higher modes. We havealso checked that this feature is common for other spin-aligned BBHs. This feature has also been reported in[12]. We understand this fact in the following way. Fol-lowing (2,2) mode, the next strongest mode is (3,3). Not-ing that the spin weighted spherical harmonic functions −3 −2 −1 D L | h l m |/ M NREOB (l,m)=(2,2)EOB (l,m)=(2,1)EOB (l,m)=(3,3)EOB (l,m)=(4,4) -0.250.000.25-0.0250.0000.025-0.050.000.05-0.0150.0000.015 0 600 1200 1800 2400 3000 3600−0.10.00.1 D L h + / M t/M 3650 3700 3750 3800 NR (l ≤ 4, m ≠ 0)EOB (l,|m|)=(2,2),(2,1)(3,3),(4,4)EOB (l,|m|)=(2,2)

FIG. 6: Waveform comparison between our

SEOBNREHM waveform model and SXS:BBH:1374, corresponding to spinless BBHwith mass ratio q = 3. The top panel shows the comparison of the amplitudes for the ( l, m ) = (2 , , (2 , , (3 , , (4 ,

4) modes.The second to ﬁfth row corresponds to the comparison of the real part of the ( l, m ) = (2 , , (2 , , (3 , , (4 ,

4) modes respectivelybetween the

SEOBNREHM waveform and the NR waveform. The bottom row shows h + ( ι, ϕ, κ ) with choice ι = π/ , ϕ = π, κ = 0.The red line means only (2,2) mode is used to combine the gravitational wave strain. The green line means the (2,2) mode andthe higher modes (2,1), (3,3) and (4,4) are used to combine the gravitational wave strain. read − Y ± ( ι, ϕ c ) = 18 e − iϕ c r π (1 ∓ cos ι ) , (38) − Y ± ( ι, ϕ c ) = ± e − iϕ c r π (1 ∓ cos ι ) sin ι. (39)So we have (cid:12)(cid:12)(cid:12) − Y ± − Y ± (cid:12)(cid:12)(cid:12) ∝ ι which means the contributionof (3,3) increases when ι increases from 0 to π/

2. Unlessthe waveform model for (3,3) mode is as accurate as (2,2)mode, the matching factor may decrease when ι increasesfrom 0 to π/ SEOBNREHM model waveform to the NR waveform we havedeﬁned an averaged and a minimal matching factor over ( κ, ϕ, ι ) respectively as [12]FF( θ ) ≡ π Z π d κ Z π d ϕ Z − FF( κ, ϕ, ι, θ ) d(cos ι ) , (40)FF min ( θ ) ≡ min κ,ϕ,ι FF( κ, ϕ, ι, θ ) . (41)Firstly we check the quasi-circular BBHs whose pa-rameters have been listed in the Appendix. C. In allwe have checked 281 NR simulations for quasi-circularBBH. The NR multipolar waveform was constructed byusing l ≤ , m = 0 modes. Respect to diﬀerent totalmass, we plot the resulted matching factor in the Fig. 8.Since the higher modes are stronger when the mass ratioia larger, we have divided the eccentric BBH cases into0 c o s ι . . . . . . . . . c o s ι . . . . . . . . . . . FIG. 7: The dependence on the angles ( ϕ, ι ) of the matching factor between the

SEOBNREHM model waveform and the SXS NR( l ≤ , m = 0) waveform. The contour lines together with the contour values of the matching factor FF deﬁned in the Eq. (37)are shown in the ﬁgure. The intrinsic parameters except the total mass corresponds to the case SXS:BBH:1374 (correspondingto the waveform shown in the Fig. 6). In addition κ = 0 and the Advanced LIGO designed PSD are used in this ﬁgure. Thetop row corresponds to the total mass M = 40 M ⊙ and the bottom row corresponds to the total mass M = 180 M ⊙ . In theleft column, only ( l, | m | ) = (2 ,

2) mode of the

SEOBNREHM model is used to combine the gravitational wave strain. In the rightcolumn, modes ( l, | m | ) = (2 , , (2 , , (3 , , (4 ,

4) of the

SEOBNREHM model are used. three groups, including q <

2, 2 ≤ q ≤ q > min are bigger than 98% for all casesexcept SXS:BBH:1422, SXS:BBH:1424, SXS:BBH:1425and SXS:BBH:1428. Together with SXS:BBH:1427 welist these 5 smallest FF min cases among the 281 quasi-circular BBHs in the Table. II. These 5 simulationswere not available when the paper [12] was published.In the Table. II we also list the matching results for SEOBNRv4HM waveform model. We ﬁnd that the behaviorof

SEOBNRv4HM is similar to that of

SEOBNREHM for these5 cases.We are more interested in the performance of highermodes of the new waveform in elliptical cases. We alsocalculate the averaged matching factor FF and the min- imal matching factor FF min between the new highermodes waveform model

SEOBNREHM and the NR wave-forms. The comparison results are plotted in the Fig. 9.Here we again divide the BBHs into three groups basedon mass ratio, including q < q = 2 and q = 3. Asexpected, when q ≥ q < e EOB < . B. Comparison with

SEOBNRv4HM

In this subsection, we compare our new (2,2), (2,1),(3,3) and (4,4) modes waveforms to the corresponding1

20 40 60 80 100 120 140 160 180 200M tot [M ⊙ ]10 −⊙ −2 −2 ⊙ × 10 −2 −1 − FF

20 40 60 80 100 120 140 160 180 200M tot [M ⊙ ]10 −⊙ −2 −2 ⊙ × 10 −2 −1 q < 22 ≤ q ≤ 5q > 5

20 40 60 80 100 120 140 160 180 200M tot [M ⊙ ]10 −⊙ −2 −2 ⊙ × 10 −2 −1 − FF m i n

20 40 60 80 100 120 140 160 180 200M tot [M ⊙ ]10 −⊙ −2 −2 ⊙ × 10 −2 −1 FIG. 8: The matching factors between the

SEOBNREHM waveform model and the SXS NR waveforms for the quasi-circularBBHs listed in the Appendix. C. The NR waveform modes ( l ≤ , m = 0) and the Advanced LIGO design zero-detunedhigh-power noise PSD are used in this ﬁgure. The top row is for the averaged matching factor over angles κ, ϕ, ι whichis deﬁned in the Eq. (40). The bottom row is for the minimal matching factor which is deﬁned in the Eq. (41). In theleft column only ( l, | m | ) = (2 , SEOBNREHM waveform mode is used. In the right column the

SEOBNREHM waveform modes( l, | m | ) = (2 , , (2 , , (3 , , (4 ,

4) are used.TABLE III: Matching factors between the models waveform and the SXS numerical relativity waveform for the 10 mostlymismatched cases between

SEOBNREHM and

SEOBNRv4HM . ID corresponds to the id number in the SXS catalog. FF and FF min are respectively the averaged matching factor and the minimal matching factor respect to the three angles ( κ, ι, ϕ ). ‘E’ meansmodel

SEOBNREHM and ‘C’ means model

SEOBNRv4HM .ID q χ χ FF( C ) FF min ( C ) FF( E ) FF min ( E )0025 1.50 +0.50 -0.50 99.77% 99.48% 99.53% 99.04%0254 2.00 +0.60 -0.60 99.82% 99.55% 99.69% 99.41%0305 1.22 +0.33 -0.44 99.58% 99.54% 99.35% 98.83%1453 2.35 +0.80 -0.78 99.74% 99.40% 99.71% 99.44%1474 1.28 +0.72 -0.80 99.04% 98.53% 98.10% 96.42%1481 1.00 +0.73 +0.79 99.66% 99.53% 99.63% 99.48%1495 1.00 +0.78 +0.53 99.52% 99.37% 99.53% 99.41%1496 1.16 +0.80 +0.03 98.95% 98.75% 98.78% 98.21%1497 1.00 +0.68 +0.67 99.68% 99.58% 99.63% 99.54%1498 1.03 +0.22 -0.78 99.78% 99.49% 99.40% 98.69% ones given by SEOBNRv4HM [12] for quasi-circular spin-aligned BBHs. Similar to the above (2,2) mode com-parison to the

SEOBNRv4 model, we use the matchingfactor FF deﬁned in the Eq. (30) for each mode to dothe comparison. In order to use the NR simulations as a reference we check the BBH parameters corresponding tothe SXS quasi-circular BBH simulations listed in the Ap-pendix. C. For each given combination of ν , χ z and χ z we scan the total mass in the range M ∈ [20 , M ⊙ foreach mode to search the minimal matching factor FF min

20 40 60 80 100 120 140 160 180 200M tot [M ⊙ ]10 −⊙ −2 −2 ⊙ × 10 −2 −1 − FF

20 40 60 80 100 120 140 160 180 200M tot [M ⊙ ]10 −⊙ −2 −2 ⊙ × 10 −2 −1 q<2q=⊙q=2

20 40 60 80 100 120 140 160 180 200M tot [M ⊙ ]10 −2 −2 ⊙ × 10 −2 −1 − FF m i n

20 40 60 80 100 120 140 160 180 200M tot [M ⊙ ]10 −2 −2 ⊙ × 10 −2 −1 q<2q=⊙q=2 FIG. 9: The matching factors between the

SEOBNREHM waveform model and the SXS NR waveforms for the eccentric BBHslisted in the Table. I. The NR waveform modes ( l ≤ , m = 0) and the Advanced LIGO design zero-detuned high-power noisePSD are used in this ﬁgure. The top row is for the averaged matching factor over angles κ, ϕ, ι which is deﬁned in the Eq. (40).The bottom row is for the minimal matching factor which is deﬁned in the Eq. (41). In the left column only ( l, | m | ) = (2 , SEOBNREHM waveform mode is used. In the right column the

SEOBNREHM waveform modes ( l, | m | ) = (2 , , (2 , , (3 , , (4 ,

4) areused. When including higher modes, all of the matching factors are bigger than 98%. Here we have used the initial eccentricity e EOB listed in the Table. I to generate the

SEOBNREHM waveform. respect to the total mass.We plot the FF min respect to diﬀerent mass ratio andblack hole spin in the Fig. 10. Following the Fig. 3 wehave used variables χ A = ( χ − χ ) / χ eﬀ = ( m χ + m χ ) /M in the plot.For spinless cases, all matching factors for all (2,2),(2,1), (3,3) and (4,4) modes are bigger than 99.9%. Re-garding to (2,2) and (4,4) modes, all matching factors arebigger than 99% even for highly spinning BBHs.For (3,3) mode, there are 4 cases, includingSXS:BBH:0305, SXS:BBH:1474, SXS:BBH:1496 andSXS:BBH:1498, admit minimal matching factor less than99%. All other minimal matching factors are biggerthan 99%. Among these 4 cases, the matching fac-tor of SXS:BBH:1474 is smallest, which is 96.6%. Weplot the waveform comparison for SXS:BBH:1474 andSXS:BBH:1498 respectively in the Figs. 11 and 12. Re-garding to the waveform alignment for diﬀerent modes,we would like to mention one point. All modes wave-forms share the unique shifted time which is determinedby (2,2) mode. In another word, we did not maximize thematching factors through adjusting parameter t c except (2,2) mode.For (2,1) mode, all minimal matching factors are big-ger than 90% except 6 cases including SXS:BBH:0025,SXS:BBH:0254, SXS:BBH:1453, SXS:BBH:1481,SXS:BBH:1495 and SXS:BBH:1497. Among them, thematching factor of SXS:BBH:1481 is 53.7% which isthe smallest one. We plot the waveform comparison forSXS:BBH:1481 and SXS:BBH:0025 in the Figs. 13 and14. At the mean time we note that the corresponding(2,1) mode and (3,3) mode for the above mentioned 10cases are quite weak. This fact makes it hard to extractinformation from NR results and accordingly constructand tune model waveform. On the other hand this factmakes the mismatch mentioned above weakly aﬀect thewaveform matching to NR results directly. Speciﬁcally,the resulted averaged matching factor FF and minimalmatching factor FF min respect to the angles κ, ϕ, ι arequite similar between SEOBNREHM and

SEOBNRv4HM . Forconvenient reference we list these matching factors inthe Table. III.3 -1.0-0.50.00.51.0 χ A ν -1.0-0.50.00.51.0 χ e ff -1.0-0.50.00.51.0 χ A ν -1.0-0.50.00.51.0 χ e ff -1.0-0.50.00.51.0 χ A ν -1.0-0.50.00.51.0 χ e ff -1.0-0.50.00.51.0 χ A ν -1.0-0.50.00.51.0 χ e ff −2.0−1.8−1.6−1.4−1.2−1.0−0.8 l o g ( − FF m i n ) FIG. 10: Matching factors of the ( l, m ) = (2 , , (2 , , (3 , , (4 ,

4) modes waveform for the quasi-circular spin-aligned BBHbetween the new waveform model

SEOBNREHM proposed in the current work and

SEOBNRv4HM for the BBH parameters listed inthe Appendix. C. From left to right, the FF min respect to diﬀerent mass ratio for the ( l, m ) = (2 , , (2 , , (3 , , (4 ,

4) modesare shown respectively. The plot convention is the same to the Fig. 2.

V. POWER FRACTION OF HIGHER MODESFOR SPIN-ALIGNED BBHS

As we have seen in the previous sections that when thecontribution of higher modes increases the higher modesbecome more and more important to gravitational wave-form template for data analysis. Roughly we expect thecontribution of higher modes increases along with themass ratio q and orbital eccentricity increase. In thissection we would like to investigate the power fraction ofhigher modes.We deﬁne the power of each mode as [29] P lm = 116 π (cid:12)(cid:12)(cid:12) ˙ h lm (cid:12)(cid:12)(cid:12) . (42)At about the merger time P lm admits a peak. But usuallydiﬀerent mode approaches such peak at diﬀerent time.As an estimation of power fraction of higher modes, wecalculate the ratio between the peak P lm and the peak P η lm ≡ max t P lm max t P , (43)where the maximum is taken respect to time correspond-ing to the peak value.For the quasi-circular BBHs, we plot η , η and η in the top row of the Fig. 15 which corresponds to the281 SXS simulation results listed in the Appendix. C.In the second row we plot the results of the SEOBNREHM waveform model with e = 0 and BBH parameters cor-responding to the top row. We can see that the ratiogotten by the SEOBNREHM waveform model is very closeto the one gotten by the NR simulations. As expectedthe ratio for each mode increases with mass ratio q . Re-garding to the eﬀect of BH spin, the ratio for each mode increase with the magnitude of spin on the one hand. Onthe other hand we interestingly ﬁnd that positive spinrespect to the direction of orbital angular momentum re-sults in much bigger ratio for (3,3) and (4,4) modes. In-stead negative spin results in much bigger ratio for (2,1)mode.In order to check the eﬀect of eccentricity on the powerratio we have calculated diﬀerent initial eccentricity butkeep other BBH parameters unchanged. We found theratio η lm roughly does not change. As an example weplot the result with e = 0 . M f = 0 .

002 in the bottom row of the Fig. 15. Wesuspect this fact implies that the gravitational radiationstrongly circularizes the orbit near the merger and thequasi-circular state is approached before the ﬁnal BHforms.Since the eccentricity does not change the power ratio η lm , it is interesting to ask how about the following powerratio as a function of time ξ lm ( t ) ≡ P lm ( t ) P ( t ) . (44)We pick the cases with largest η lm shown in the Fig. 15as examples. For (2,1) mode, SXS:BBH:1433 is the onewith largest η ≈ .

24. SXS:BBH:1441 is the case withlargest η lm for both (3,3) ( η ≈ .

44) and (4,4) mode( η ≈ . ξ ( t ) in theFig. 16. Compared to quasi-circular case, the eccentric-ity just introduces an oscillation for ξ ( t ). There is noampliﬁcation of power ratio ξ introduced by the eccen-tricity.In the top panel of the Fig. 16, there is a minor butinteresting feature. Near the periastron point, ξ showstwo peaks. This is because eccentricity introduce a cor-rection to the phase of the waveform and consequently4 −3 −2 −1 D L | h l m |/ M NREOB(E) (l,m)=(2,2)EOB(E) (l,m)=(2,1)EOB(E) (l,m)=(3,3)EOB(E) (l,m)=(4,4)EOB(C) (l,m)=(2,2)EOB(C) (l,m)=(2,1)EOB(C) (l,m)=(3,3)EOB(C) (l,m)=(4,4) −0.250.000.25−0.0250.0000.025−0.0250.0000.025-0.0150.0000.015 0 1000 2000 3000 4000−0.10.00.1 D L h + / M t/M 4800 4900 5000 NR (l ≤ 4, m ≠ 0)EOB(E) (l,|m|)=(2,2),(2,1)(3,3),(4,4)EOB(C) (l,|m|)=(2,2),(2,1)(3,3),(4,4)

FIG. 11: Waveform comparison between model waveform and NR waveform for case SXS:BBH:1474. EOB(C) means

SEOBNRv4HM and EOB(E) means

SEOBNREHM . a correction to the wave frequency of each mode. Suchcorrection makes the frequency of ( l, m ) mode a little dif-ferent to the m times of the orbital frequency. When theeccentricity is strong enough, such correction results inthe mentioned two-peak behavior of ξ .So we conclude that eccentricity does not enhance thehigher modes than (2,2) of the gravitational wave. In-stead the eccentricity excite higher harmonic tones thanthe harmonic tone with frequency ω ≈ | m | ω orb for mode( l, m ), where ω orb means the orbital mean frequency[52, 53]. The Fig. 4 of [52] and the Fig. 5 of [53] haveshown the higher harmonics excitation behavior by ec-centricity based on post-Newtonian approximation. Herewe correspondingly show the full general relativity re-sults on such higher harmonics excitation in the Figs. 17and 18. For the eccentric cases we can see the dot dashbehavior of the power spectrum which is resulted fromthe burst near the periastron point from time to time.We once again see that the eccentricity does not am- plify higher spherical-harmonic mode, instead excites thehigher harmonics of each spherical-harmonic mode. Thespherical-harmonic modes of quasci-circular BBHs ad-mit solely the fundamental harmonic. In contrast, thespherical-harmonic modes of eccentric BBHs admit boththe fundamental harmonic and other higher harmonics[54]. VI. SUMMARY

The recent found gravitational wave event GW190521[18, 55, 56] inspired more interesting on the eccentric bi-nary black holes [17, 57] for ground based GW detec-tion. In order to properly treat the eccentricity duringthe gravitational wave data analysis, accurate waveformmodel for eccentric binary black hole is needed. Due tothe GW190521 like events, the demand of waveform tem-plate for eccentric binary black holes increases. In addi-5 −3 −2 −1 D L | h l m |/ M NREOB(E) (l,m)=(2,2)EOB(E) (l,m)=(2,1)EOB(E) (l,m)=(3,3)EOB(E) (l,m)=(4,4)EOB(C) (l,m)=(2,2)EOB(C) (l,m)=(2,1)EOB(C) (l,m)=(3,3)EOB(C) (l,m)=(4,4) −0.250.000.25−0.0250.0000.025-0.0050.0000.005-0.0150.0000.015 0 1000 2000 3000 4000−0.10.00.1 D L h + / M t/M 4800 4900 5000 NR (l ≤ 4, m ≠ 0)EOB(E) (l,|m|)=(2,2),(2,1)(3,3),(4,4)EOB(C) (l,|m|)=(2,2),(2,1)(3,3),(4,4)

FIG. 12: Waveform comparison between model waveform and NR waveform for case SXS:BBH:1498. EOB(C) means

SEOBNRv4HM and EOB(E) means

SEOBNREHM . tion, for the future gravitational wave detection projectsin space, the residual eccentricity of BBH events in theirfrequency bands is still relatively large. For such GWsources, using gravitational wave template of circular or-bit to detect gravitational wave signal will lead to theloss of signal-to-noise ratio [30].In the last few years, a handful of eccentric inspi-ral only [25, 52, 58–60] and Inspiral-Merger-Ringdown(IMR) models [27–29, 31, 61, 62] have been con-structed. Among these waveform models, only the mod-iﬁed TEOBiResumS SM [31] is valid for higher modes than(2,2).In this paper, we extended our previous

SEOBNRE wave-form model [29] to

SEOBNREHM to cover higher modes in-cluding (2,1), (3,3) and (4,4). At the mean time, theEOB factorized waveform is extended from quasi-circularcases to general equatorial cases, inspired by the work[38]. As the orbit of the binary system gradually circu-larizes along with the evolution, we just replace the low order post-Newtonian terms with the ones corrected byeccentricity.Our new

SEOBNREHM waveform model can recover the

SEOBNRv4HM waveforms of quasi-circular BBHs very well.In the eccentric cases of BBH, the matching factor be-tween the

SEOBNREHM waveform and the numerical rela-tivistic waveform is greater than 98%.The Hamiltonian related to the

SEOBNREHM dynamicsand the remnant black hole properties used in the ring-down part waveform are borrowed from the

SEOBNRv4 model. After more numerical relativity simulations abouteccentric BBHs are available [57], it is straight for-ward to ﬁnetune the related parameters involved in theHamiltonian and the remnant properties to improve the

SEOBNREHM ’s performance. Even without such ﬁnetun-ing, we have shown that

SEOBNREHM can match all SXSwaveforms including highly eccentric ones as good as the

SEOBNRv4 waveforms match the circular ones. Hopefullyour

SEOBNREHM can facilitate the gravitational wave data6 −4 −2 D L | h l m |/ M NREOB(E) (l,m)=(2,2)EOB(E) (l,m)=(2,1)EOB(E) (l,m)=(3,3)EOB(E) (l,m)=(4,4)EOB(C) (l,m)=(2,2)EOB(C) (l,m)=(2,1)EOB(C) (l,m)=(3,3)EOB(C) (l,m)=(4,4) −0.250.000.250.00000.0025−0.00050.00000.0005-0.0150.0000.015 0 1000 2000 3000 4000−0.10.00.1 D L h + / M t/M 4800 5000 NR (l ≤ 4, m ≠ 0)EOB(E) (l,|m|)=(2,2),(2,1)(3,3),(4,4)EOB(C) (l,|m|)=(2,2),(2,1)(3,3),(4,4)

FIG. 13: Waveform comparison between model waveform and NR waveform for case SXS:BBH:1481. EOB(C) means

SEOBNRv4HM and EOB(E) means

SEOBNREHM . analysis about eccentric BBHs [15, 17]. Acknowledgments

We would like to thank Alessandra Buonanno, Li-jing Shao and Antoni Ramos Buades for many helpfuldiscussions. This work was supported by the NSFC (No. 11690023, No. 11633001, No. 11920101003 andNo. 12021003) and by the Collaborative research pro-gram of the Institute for Cosmic Ray Research (ICRR),the University of Tokyo. Z. Cao and Z.-H. Zhu weresupported by “the Interdiscipline Research Funds of Bei-jing Normal University” and the Strategic Priority Re-search Program of the Chinese Academy of Sciences(No. XDB23000000 and No. XDB23040100).

Appendix A: Gauge transformation between the EOB coordinates and the harmonic coordinates

Ref. [45] has given the coordinate transformation among the harmonic coordinates, the ADM coordinates and theEOB coordinates, and the transformation among diﬀerent canonical velocities and canonical momentums. Accordingto such transformations we have y = x (cid:20) − νx − ν (3 ν + 14) x + 132 (cid:0) − ν − ν (cid:1) x (A1)7 −4 −3 −2 −1 D L | h l m |/ M NREOB(E) (l,m)=(2,2)EOB(E) (l,m)=(2,1)EOB(E) (l,m)=(3,3)EOB(E) (l,m)=(4,4)EOB(C) (l,m)=(2,2)EOB(C) (l,m)=(2,1)EOB(C) (l,m)=(3,3)EOB(C) (l,m)=(4,4) −0.250.000.25−0.010.000.01−0.0250.0000.025-0.0150.0000.015 0 1000 2000 3000 4000 5000 6000−0.10.00.1 D L h + / M t/M 6000 6200 6400 NR (l ≤ 4, m ≠ 0)EOB(E) (l,|m|)=(2,2),(2,1)(3,3),(4,4)EOB(C) (l,|m|)=(2,2),(2,1)(3,3),(4,4)

FIG. 14: Waveform comparison between model waveform and NR waveform for case SXS:BBH:0025. EOB(C) means

SEOBNRv4HM and EOB(E) means

SEOBNREHM . + x (cid:18) − ν − ν (13 ν + 10) x (cid:19) + x (cid:18)

132 (8 ν + 16) + 12 νx χ S + 132 (cid:0) ν − ν (cid:1) x + 132 (cid:0) ν − ν (cid:1) x (cid:19) + 132 x (cid:0) ν (cid:0) √ − νχ A χ S + χ A + χ S − (cid:1) − (cid:0) √ − νχ A χ S + 2 χ A + 2 χ S − (cid:1) + ν (cid:0) − χ S (cid:1)(cid:1)(cid:21) y = x (cid:20) x (cid:0) ν (cid:0) √ − νχ A χ S + χ A + χ S − (cid:1) − (cid:0) √ − νχ A χ S + 2 χ A + 2 χ S − (cid:1) + ν (cid:0) − χ S (cid:1)(cid:1) (A2)+ x (cid:18) − ν − ν (9 ν − x − ν (29 ν + 17) x (cid:19) − x (cid:0) − ν (cid:0) √ − νχ A χ S + 8 χ A + 8 χ S − (cid:1) + 2 (cid:0) √ − νχ A χ S + χ A + χ S (cid:1) + ν (cid:0) χ S − (cid:1)(cid:1)(cid:21) y = x (cid:20) (cid:0) ν + 7 ν (cid:1) x + 38 ν ( ν + 1) x + x (cid:18) ν ν (3 ν + 5) x (cid:19) + 3 νx x (cid:18)

18 ( − ν −

8) + 18 (cid:0) ν − ν (cid:1) x + 18 (cid:0) − ν − ν (cid:1) x (cid:19) -1.0-0.50.00.51.0 χ A -1.0-0.50.00.51.0 χ e ff -1.0-0.50.00.51.0 χ A -1.0-0.50.00.51.0 χ e ff -1.0-0.50.00.51.0 χ A -1.0-0.50.00.51.0 χ e ff -1.0-0.50.00.51.0 χ A -1.0-0.50.00.51.0 χ e ff -1.0-0.50.00.51.0 χ A -1.0-0.50.00.51.0 χ e ff -1.0-0.50.00.51.0 χ A -1.0-0.50.00.51.0 χ e ff -1.0-0.50.00.51.0 χ A -1.0-0.50.00.51.0 χ e ff χ A -1.0-0.50.00.51.0 χ e ff χ A -1.0-0.50.00.51.0 χ e ff FIG. 15: Power ratio η lm for modes (2,1), (3,3) and (4,4) respectively from the left to the right. The top row corresponds to thequasi-circular SXS simulations listed in the Appendix. C. The middle row corresponds to the SEOBNREHM waveform model resultswith e = 0 and BBH parameters corresponding to the top row. The bottom row corresponds to the SEOBNREHM waveformmodel results with e = 0 . Mf = 0 .

002 and other BBH parameters the same to the top row. −1 t/10 M10 −3 −2 −1 | ̇ h | /| ̇ h | e0=0̇0e0=0̇1e0=0̇2 e0=0̇3e0=0̇4e0=0̇510 −1 t/10 M10 −2 −1 | ̇ h | /| ̇ h | e0=0̇0e0=0̇1e0=0̇2 e0=0̇3e0=0̇4e0=0̇510 −1 t/10 M10 −3 −2 −1 | ̇ h | /| ̇ h | e0=0̇0e0=0̇1e0=0̇2 e0=0̇3e0=0̇4e0=0̇5 FIG. 16: The top panel: power ratio ξ BBH mass ratio and spin parameters corresponding to SXS:BBH:1433 but diﬀerentinitial eccentricity shown in the legend at reference frequency Mf = 0 . ξ and SXS:BBH:1441. The bottom panel: similar to the top panel but for ξ and SXS:BBH:1441. + 12 x (cid:0) − ν (cid:0) √ − νχ A χ S + 4 χ A + 4 χ S − (cid:1) + 2 √ − νχ A χ S + χ A + ν (cid:0) χ S − (cid:1) + χ S (cid:1)(cid:21) − νx χ S e − imφ h = e − imφ (cid:20) − m ν x x − imνx x − im (1 − ν ) νx x x (A4)+ x (cid:18) x (cid:0) − imν − imν (cid:1) − imνx (cid:19)(cid:21) where χ S = ( χ + χ ) / , χ A = ( χ − χ ) / FIG. 17: The normalized time-frequency behavior through constant-Q-transformation for the gravitational wave modes. TheBBH mass ratio and spin parameters corresponding to SXS:BBH:1433. The top row corresponds to the quasi-circular case.The bottom row corresponds to the eccentric case with initial eccentricity e = 0 . Mf = 0 . Appendix B: 2PN factorized waveform

Here we introduce two auxiliary parameters A = x x , A = x x . (B1)We list the ζ elm in the following. Firstly, (2 ,

2) mode: ζ e = c h + c h x + c h x + c h x (B2) c h = − A iA A + A c h = −

32 + ν A (cid:18) ν − (cid:19) + A A (cid:18) i − iν (cid:19) + A (cid:18) ν

28 + 114 (cid:19) (B4)1+ A (cid:18) ν − (cid:19) + A (cid:18) − ν − (cid:19) + A A (cid:18) iν

14 + i (cid:19) + A A (cid:18) iν

14 + i (cid:19) c h = A (cid:18) −√ − νχ A + 7 νχ S − χ S (cid:19) − iA (cid:0) √ − νχ A + (3 − ν ) χ S (cid:1) (B5) c h = 205 ν − ν

36 + 65252 + A (cid:18) ν − ν − (cid:19) + A (cid:18) − ν − ν − (cid:19) (B6)+ A A (cid:18) − iν − iν − i (cid:19) + A A (cid:18) − iν

63 + 935 iν

252 + 71 i (cid:19) + A A (cid:18) iν

24 + 179 iν

168 + i (cid:19) + A A (cid:18) − ν − ν

56 + 2384 (cid:19) + A A (cid:18) − ν − ν − (cid:19) + A A (cid:18) − iν

21 + 29 iν − i (cid:19) + A A (cid:18) iν

12 + 179 iν

84 + 2 i (cid:19) + A (cid:18) − ν − ν (cid:19) + A A (cid:18) ν

48 + 179 ν

336 + 142 (cid:19) + A (cid:18) ν

48 + 179 ν

336 + 142 (cid:19) + A A (cid:18) iν

24 + 179 iν

168 + i (cid:19) + A (cid:18) − ν

63 + 1355 ν − (cid:19) + A (cid:18) ν − ν − (cid:19) + A + A (cid:0) √ − νχ A + (1 − ν ) χ S (cid:1) + 12 h(cid:0) √ − νχ A + (1 + ν ) χ S (cid:1) − ν χ S i . (2 ,

1) mode: ζ e = c h + c h x + c h x (B7) A c h = A (B8) A c h = − χ A √ − ν − χ S A c h = A (cid:18) − ν (cid:19) + A A (cid:18) i − iν (cid:19) + A A (cid:18) − ν − (cid:19) + A (cid:18) ν − (cid:19) (B10) A c h = 114 (cid:18) (84 − ν ) χ A √ − ν + (84 − ν ) χ S (cid:19) + i A A (cid:18) (104 ν − √ − ν χ A + (8 ν − χ S (cid:19) (B11)+ 128 A (cid:18) (95 ν + 126) χ A √ − ν + 3(29 ν + 42) χ S (cid:19) + 128 A (cid:18) (327 ν − χ A √ − ν + (107 ν − χ S (cid:19) . (3 ,

3) mode: ζ e = c h + c h x + c h x (B12) c h = − iA − A A + 23 iA A + 4 iA A A c h = A (cid:18) − iν − i (cid:19) + A (cid:18) iν − i (cid:19) + A (cid:18) − iν − i (cid:19) (B14)+ A A (cid:18) − ν − (cid:19) + A A (cid:18) ν − (cid:19) + A (cid:18) ν − (cid:19) + A A (cid:18) iν

27 + 4 i (cid:19) + A A (cid:18) i − iν (cid:19) + A (cid:18) ν − (cid:19) + A (cid:18) ν

27 + 227 (cid:19) + A A (cid:18) iν i (cid:19) + A A (cid:18) − ν − (cid:19) c h = 19 (cid:18) ν − χ A √ − ν + (3 ν − χ S (cid:19) + 19 A (cid:18) − ν ) χ S − (25 ν − χ A √ − ν (cid:19) (B15)+ iA A (cid:18) (77 ν − √ − ν χ A + (31 ν − χ S (cid:19) + 118 A (cid:18) (119 ν − χ A √ − ν + (35 ν − χ S (cid:19) . (4 ,

4) mode: ζ e = c h + c h x (B16)2 c h = 3 A − iA A − A A − A

32 + 38 iA A + 2732 iA A + 3 A

32 + 51 A

64 + 764 (B17) c h = 11 − ν " (cid:0) − ν + 1358 ν − (cid:1) + 3 A (cid:0) A (cid:0) − ν + 30 ν + 60 (cid:1) − iA (cid:0) ν − ν + 2348 (cid:1)(cid:1) − A (cid:0) A (cid:0) A (cid:0) ν − ν + 859 (cid:1) + 12910 ν − ν + 17971 (cid:1) − iA (cid:0) ν − ν − (cid:1)(cid:1) (cid:0) A (cid:0) − ν + 2 ν + 4 (cid:1) + A (cid:0) A (cid:0) ν − ν − (cid:1) + 109 ν − ν + 46 (cid:1)(cid:1) + A (cid:0) A (cid:0) ν − ν − (cid:1) − A (cid:0) ν − ν + 5426 (cid:1) − (cid:0) ν − ν + 579 (cid:1)(cid:1) − iA A (cid:0) A (cid:0) ν − ν − (cid:1) − A (cid:0) ν − ν + 577 (cid:1) − ν − ν + 2832 (cid:1) . Appendix C: Numerical-relativity waveforms forquasi-circular BBH used in the current work

We have used 281 numerical relativity simulations inthe current paper for quasi-circular BBH correspondingdiﬀerent parameters combination of symmetric mass ra-tio ν and the black holes’ spin parameters χ z and χ z . Inthe following we list the parameters for these numericalrelativity simulations, the corresponding minimal match-ing factor FF for (2,2) mode (check the Sec. IIIfor detail explanation) and the averaged matching fac-tor FF( θ ), the minimal matching factor FF min ( θ ) respectto the three angles ( κ, ι, ϕ ) (check the Sec. IV for detailexplanation). ID q χ χ FF FF( θ ) FF min ( θ )0001 1 .

00 0 .

00 0 .

00 99 .

58% 99 .

58% 99 . .

00 0 .

00 0 .

00 99 .

58% 99 .

58% 99 . . − .

50 0 .

00 99 .

77% 99 .

68% 99 . .

00 +0 .

50 0 .

00 99 .

88% 99 .

72% 99 . .

50 0 .

00 0 .

00 99 .

49% 99 .

49% 99 . .

50 0 .

00 0 .

00 99 .

57% 99 .

57% 99 . .

50 +0 .

50 0 .

00 99 .

71% 99 .

62% 99 . .

50 +0 .

50 0 .

00 99 .

86% 99 .

63% 99 . . − .

50 0 .

00 99 .

57% 99 .

31% 98 . . − .

50 0 .

00 99 .

56% 99 .

38% 99 . . − .

50 +0 .

50 99 .

70% 99 .

59% 99 . .

50 +0 . − .

50 99 .

47% 99 .

53% 99 . .

00 0 .

00 0 .

00 99 .

66% 99 .

72% 99 . .

00 +0 .

50 0 .

00 99 .

52% 99 .

56% 99 . . − .

50 0 .

00 99 .

95% 99 .

87% 99 . . − .

50 0 .

00 99 .

69% 99 .

79% 99 . . − .

50 0 .

00 99 .

92% 99 .

81% 99 . . − .

50 0 .

00 99 .

47% 99 .

65% 98 . .

00 +0 .

50 0 .

00 99 .

10% 99 .

52% 99 . .

00 +0 . − .

50 99 .

40% 99 .

60% 99 . . − . − .

50 99 .

61% 99 .

80% 99 . .

00 +0 .

50 +0 .

50 99 .

49% 99 .

44% 99 . .

00 0 .

00 0 .

00 99 .

69% 99 .

77% 99 . .

00 0 .

00 0 .

00 99 .

84% 99 .

79% 99 . .

00 0 .

00 0 .

00 99 .

89% 99 .

73% 99 . . − .

50 0 .

00 99 .

75% 99 .

66% 99 .

05% 0061 5 .

00 +0 .

50 0 .

00 98 .

99% 98 .

81% 98 . .

00 0 .

00 0 .

00 99 .

80% 99 .

44% 98 . . − .

50 0 .

00 99 .

38% 99 .

37% 98 . .

00 +0 .

50 0 .

00 98 .

96% 98 .

38% 97 . .

00 0 .

00 0 .

00 99 .

58% 99 .

58% 99 . .

00 +0 .

50 0 .

00 99 .

85% 99 .

78% 99 . .

00 +0 .

50 0 .

00 99 .

91% 99 .

72% 99 . .

00 +0 .

50 0 .

00 99 .

90% 99 .

75% 99 . .

00 0 .

00 0 .

00 99 .

58% 99 .

57% 99 . .

00 0 .

00 0 .

00 99 .

65% 99 .

64% 99 . . − .

50 0 .

00 98 .

95% 98 .

40% 98 . .

00 0 .

00 0 .

00 99 .

60% 99 .

54% 99 . .

00 0 .

00 0 .

00 99 .

61% 99 .

60% 99 . .

50 0 .

00 0 .

00 99 .

49% 99 .

49% 99 . .

50 0 .

00 0 .

00 99 .

51% 99 .

48% 99 . . − .

50 0 .

00 99 .

55% 99 .

43% 99 . . − .

50 0 .

00 99 .

88% 99 .

86% 99 . .

00 0 .

00 0 .

00 99 .

56% 99 .

54% 99 . .

00 0 .

00 0 .

00 99 .

91% 99 .

80% 99 . . − .

50 0 .

00 99 .

45% 99 .

45% 99 . . − .

50 0 .

00 99 .

58% 99 .

71% 99 . .

00 +0 .

50 0 .

00 99 .

05% 98 .

55% 98 . . − .

50 0 .

00 99 .

59% 99 .

34% 98 . .

00 0 .

00 0 .

00 99 .

85% 99 .

77% 99 . .

00 0 .

00 0 .

00 99 .

80% 99 .

77% 99 . . − .

50 0 .

00 99 .

56% 99 .

57% 98 . . − . − .

44 99 .

13% 99 .

05% 98 . . − . − .

20 99 .

80% 99 .

78% 99 . .

00 +0 .

20 +0 .

20 99 .

82% 99 .

89% 99 . . − . − .

60 99 .

87% 99 .

83% 99 . .

00 +0 .

60 +0 .

60 99 .

48% 99 .

44% 99 . .

00 +0 .

85 +0 .

85 99 .

68% 99 .

46% 99 . . − . − .

80 99 .

72% 99 .

66% 99 . .

00 +0 .

80 +0 .

80 99 .

69% 99 .

54% 99 . .

00 +0 .

95 +0 .

95 99 .

61% 99 .

25% 98 . .

00 +0 .

97 +0 .

97 99 .

38% 98 .

77% 97 . . − . − .

90 99 .

86% 99 .

80% 99 . .

00 +0 .

90 +0 .

90 99 .

61% 99 .

36% 99 . .

00 +0 .

60 0 .

00 99 .

83% 99 .

73% 99 . .

00 0 .

00 0 .

00 99 .

70% 99 .

66% 99 . .

00 0 .

00 0 .

00 99 .

66% 99 .

84% 99 . .

00 0 .

00 0 .

00 99 .

61% 99 .

78% 99 . .

00 0 .

00 0 .

00 99 .

48% 99 .

62% 98 . .

00 +0 .

44 +0 .

44 99 .

19% 99 .

35% 99 . . − . − .

44 98 .

90% 98 .

80% 98 . .

00 +0 .

98 +0 .

98 99 .

33% 99 .

08% 98 . .

00 +0 .

50 0 .

00 99 .

58% 99 .

52% 99 . .

00 +0 .

75 +0 .

75 99 .

76% 99 .

66% 99 . .

00 +0 .

96 +0 .

96 99 .

54% 99 .

21% 98 . .

00 +0 .

99 +0 .

99 99 .

09% 98 .

85% 98 . .

00 +0 .

99 +0 .

99 98 .

80% 98 .

66% 97 . .

00 0 .

00 0 .

00 99 .

61% 99 .

59% 99 . .

00 0 .

00 0 .

00 99 .

85% 99 .

61% 99 . .

00 0 .

00 0 .

00 99 .

82% 99 .

83% 99 . .

00 0 .

00 0 .

00 99 .

71% 99 .

76% 99 . .

00 0 .

00 0 .

00 99 .

59% 99 .

64% 99 . .

99 0 .

00 0 .

00 99 .

86% 99 .

29% 98 . .

27 0 .

00 0 .

00 99 .

83% 99 .

35% 98 . .

04 0 .

00 0 .

00 99 .

76% 99 .

70% 99 . .

19 0 .

00 0 .

00 99 .

82% 99 .

45% 98 . .

17 0 .

00 0 .

00 99 .

84% 98 .

83% 97 . .

50 0 .

00 0 .

00 99 .

78% 99 .

70% 99 . .

51 0 .

00 0 .

00 99 .

49% 99 .

55% 99 . .

58 0 .

00 0 .

00 99 .

87% 99 .

52% 98 . .

50 0 .

00 0 .

00 99 .

70% 99 .

74% 99 . .

52 0 .

00 0 .

00 99 .

62% 99 .

62% 99 . .

76 0 .

00 0 .

00 99 .

88% 99 .

48% 98 . .

66 0 .

00 0 .

00 99 .

84% 99 .

22% 98 . .

52 0 .

00 0 .

00 99 .

83% 99 .

66% 99 . .

20 0 .

00 0 .

00 99 .

77% 99 .

64% 99 . .

73 0 .

00 0 .

00 99 .

82% 99 .

26% 98 . .

27 0 .

00 0 .

00 99 .

56% 99 .

69% 99 . .

32 0 .

00 0 .

00 99 .

44% 99 .

55% 99 . .

00 +0 .

60 0 .

00 99 .

04% 98 .

09% 97 . .

00 +0 .

40 0 .

00 99 .

24% 98 .

86% 98 . .

00 +0 .

40 0 .

00 99 .

23% 98 .

66% 98 . . − .

40 0 .

00 99 .

21% 99 .

55% 99 . . − .

40 0 .

00 99 .

54% 99 .

56% 98 . . − .

60 0 .

00 99 .

78% 98 .

84% 98 . . − .

90 0 .

00 99 .

26% 99 .

13% 98 . . − . − .

50 99 .

25% 99 .

17% 98 . . − .

90 0 .

00 99 .

12% 98 .

92% 98 . . − .

90 +0 .

90 99 .

86% 99 .

75% 99 . . − . − .

80 99 .

81% 99 .

70% 99 . . − .

80 +0 .

80 99 .

89% 99 .

80% 99 . . − . − .

25 99 .

12% 99 .

01% 98 . . − . − .

60 99 .

86% 99 .

84% 99 . . − .

60 0 .

00 99 .

64% 99 .

56% 99 . . − .

60 +0 .

60 99 .

80% 99 .

81% 99 . . − .

50 +0 .

50 99 .

76% 99 .

62% 99 . . − .

50 +0 .

90 99 .

72% 99 .

66% 99 . . − . − .

80 99 .

87% 99 .

82% 99 . . − .

40 +0 .

80 99 .

81% 99 .

77% 99 . . − .

30 0 .

00 99 .

67% 99 .

55% 99 . .

00 +0 .

30 0 .

00 99 .

86% 99 .

74% 99 . .

00 +0 . − .

80 99 .

83% 99 .

72% 99 . .

00 +0 .

40 +0 .

80 99 .

33% 99 .

28% 99 . .

00 +0 . − .

90 99 .

79% 99 .

69% 99 . .

00 +0 .

60 0 .

00 99 .

79% 99 .

74% 99 . .

00 +0 .

60 +0 .

60 99 .

49% 99 .

43% 99 . .

00 +0 .

65 +0 .

25 99 .

44% 99 .

39% 99 . .

00 +0 .

80 +0 .

80 99 .

72% 99 .

57% 99 . .

00 +0 .

90 0 .

00 98 .

77% 98 .

88% 98 . .

00 +0 .

90 +0 .

50 99 .

55% 99 .

49% 99 . . − .

87 +0 .

85 99 .

15% 99 .

14% 99 . . − . − .

85 99 .

23% 99 .

14% 98 . . − . − .

60 99 .

94% 99 .

88% 99 . . − .

60 0 .

00 99 .

80% 99 .

74% 99 .

57% 0237 2 . − .

60 +0 .

60 99 .

49% 99 .

44% 99 . . − . − .

50 99 .

83% 99 .

69% 99 . . − .

37 +0 .

85 99 .

77% 99 .

57% 99 . . − . − .

30 99 .

73% 99 .

61% 99 . . − .

30 0 .

00 99 .

81% 99 .

47% 99 . . − .

30 +0 .

30 99 .

79% 99 .

64% 99 . . − . − .

85 99 .

12% 98 .

97% 98 . .

00 0 . − .

60 98 .

69% 98 .

96% 98 . .

00 0 . − .

30 99 .

30% 99 .

28% 99 . .

00 0 .

00 +0 .

30 99 .

75% 99 .

53% 99 . .

00 0 .

00 +0 .

60 99 .

85% 99 .

74% 99 . .

00 +0 .

13 +0 .

85 99 .

82% 99 .

65% 99 . .

00 +0 . − .

30 99 .

58% 99 .

53% 99 . .

00 +0 .

30 0 .

00 99 .

82% 99 .

74% 99 . .

00 +0 .

30 +0 .

30 99 .

91% 99 .

74% 99 . .

00 +0 . − .

85 98 .

53% 98 .

80% 98 . .

00 +0 .

50 +0 .

50 99 .

90% 99 .

83% 99 . .

00 +0 . − .

60 99 .

93% 99 .

70% 99 . .

00 +0 .

60 0 .

00 99 .

86% 99 .

77% 99 . .

00 +0 .

60 +0 .

60 99 .

89% 99 .

83% 99 . .

00 +0 .

85 +0 .

85 99 .

76% 99 .

39% 99 . .

00 +0 . − .

85 99 .

33% 99 .

13% 98 . .

50 0 .

00 0 .

00 99 .

57% 99 .

54% 99 . . − . − .

85 98 .

40% 98 .

37% 98 . . − .

73 +0 .

85 99 .

79% 99 .

62% 99 . . − .

60 0 .

00 99 .

96% 99 .

79% 99 . . − .

60 +0 .

60 99 .

88% 99 .

75% 99 . . − . − .

60 99 .

72% 99 .

74% 99 . . − . − .

40 99 .

86% 99 .

79% 99 . . − .

60 +0 .

40 99 .

93% 99 .

81% 99 . . − . − .

50 99 .

79% 99 .

66% 99 . . − . − .

60 99 .

58% 99 .

63% 99 . . − .

40 +0 .

60 99 .

82% 99 .

62% 99 . . − . − .

30 99 .

80% 99 .

75% 99 . . − .

30 0 .

00 99 .

92% 99 .

81% 99 . . − .

30 +0 .

30 99 .

93% 99 .

85% 99 . . − . − .

85 99 .

09% 99 .

15% 99 . . − .

23 +0 .

85 99 .

82% 99 .

71% 99 . .

00 0 . − .

60 98 .

79% 99 .

04% 98 . .

00 0 . − .

30 99 .

49% 99 .

41% 99 . .

00 0 .

00 +0 .

30 99 .

84% 99 .

77% 99 . .

00 0 .

00 +0 .

60 99 .

87% 99 .

81% 99 . .

00 +0 . − .

85 98 .

53% 98 .

61% 98 . .

00 +0 .

27 +0 .

85 99 .

54% 99 .

39% 99 . .

00 +0 . − .

30 99 .

74% 99 .

71% 99 . .

00 +0 .

30 0 .

00 99 .

79% 99 .

58% 99 . .

00 +0 .

30 +0 .

30 99 .

68% 99 .

50% 99 . .

00 +0 . − .

60 99 .

60% 99 .

41% 98 . .

00 +0 .

40 +0 .

60 99 .

57% 99 .

41% 99 . .

00 +0 .

50 +0 .

50 99 .

63% 99 .

46% 99 . .

00 +0 . − .

60 99 .

79% 99 .

59% 99 . .

00 +0 . − .

40 99 .

75% 99 .

58% 99 . .

00 +0 .

60 0 .

00 99 .

71% 99 .

53% 99 . .

00 +0 .

60 +0 .

40 99 .

72% 99 .

62% 99 . .

00 +0 .

60 +0 .

60 99 .

81% 99 .

69% 99 . .

00 +0 . − .

85 99 .

85% 99 .

68% 99 . .

00 +0 .

85 +0 .

85 99 .

46% 98 .

95% 98 . .

50 0 .

00 0 .

00 99 .

77% 99 .

68% 99 . .

50 0 .

00 0 .

00 99 .

88% 99 .

80% 99 . .

50 0 .

00 0 .

00 99 .

91% 99 .

61% 99 . .

50 0 .

00 0 .

00 99 .

90% 99 .

57% 99 . .

00 0 .

00 0 .

00 99 .

90% 99 .

58% 98 . .

50 0 .

00 0 .

00 99 .

91% 99 .

38% 98 . .

50 0 .

00 0 .

00 99 .

79% 99 .

43% 98 . .

00 0 .

00 0 .

00 99 .

79% 99 .

43% 98 . .

50 0 .

00 0 .

00 99 .

77% 99 .

19% 98 . .

00 0 .

00 0 .

00 99 .

76% 99 .

04% 97 . .

00 +0 . − .

50 99 .

89% 99 .

72% 99 . .

22 +0 . − .

44 99 .

60% 99 .

40% 98 . . − .

80 +0 .

80 99 .

32% 98 .

06% 95 . . − . − .

46 99 .

27% 98 .

75% 96 . . − . − .

75 99 .

47% 99 .

35% 98 . . − . − .

80 99 .

08% 98 .

87% 96 . . − .

80 +0 .

67 99 .

35% 98 .

98% 97 . . − . − .

73 99 .

18% 99 .

27% 98 . . − . − .

70 99 .

13% 98 .

69% 96 . . − . − .

78 99 .

26% 99 .

23% 98 . .

00 +0 . − .

75 99 .

53% 99 .

41% 98 . .

00 +0 . − .

78 99 .

40% 99 .

41% 98 . . − .

74 +0 .

21 99 .

00% 98 .

86% 97 . .

37 +0 .

80 +0 .

80 99 .

13% 98 .

34% 98 . . − .

79 +0 .

07 98 .

90% 99 .

00% 98 . .

28 +0 . − .

80 99 .

09% 99 .

47% 98 . .

04 +0 .

80 +0 .

15 99 .

03% 98 .

06% 97 . .

87 +0 .

13 +0 .

80 99 .

65% 99 .

07% 98 . .

48 +0 . − .

32 99 .

40% 99 .

02% 98 . .

64 +0 .

77 +0 .

31 99 .

07% 98 .

60% 98 . .

00 +0 . − .

48 98 .

98% 98 .

22% 97 . . − . − .

18 99 .

21% 98 .

87% 97 . .

68 +0 . − .

74 99 .

18% 98 .

89% 98 . . − . − .

76 99 .

16% 99 .

53% 99 . . − .

50 +0 .

80 99 .

48% 99 .

72% 99 . . − .

80 +0 .

78 99 .

80% 99 .

62% 99 . .

16 +0 .

74 +0 .

80 99 .

88% 99 .

02% 98 . . − .

48 +0 .

52 99 .

86% 99 .

69% 99 . . − . − .

34 99 .

24% 99 .

57% 98 . . − . − .

80 98 .

99% 99 .

26% 98 . .

06 +0 . − .

80 99 .

17% 99 .

52% 99 . .

64 +0 . − .

43 99 .

64% 99 .

44% 99 . .

35 +0 . − .

78 99 .

77% 99 .

55% 99 . . − . − .

73 99 .

26% 99 .

63% 99 . . − .

40 0 .

00 99 .

31% 98 .

94% 97 . .

00 +0 .

74 +0 .

70 99 .

83% 99 .

68% 99 . .

25 +0 .

54 +0 .

80 99 .

62% 99 .

48% 99 . . − .

06 +0 .

80 99 .

51% 99 .

81% 99 . .

26 +0 .

76 +0 .

80 99 .

45% 99 .

20% 98 . .

00 +0 .

12 +0 .

11 99 .

40% 99 .

03% 98 .

16% 1461 2 . − . − .

80 99 .

15% 99 .

43% 99 . . − .

80 +0 .

51 99 .

89% 99 .

69% 99 . .

98 +0 .

61 +0 .

24 98 .

68% 98 .

47% 98 . . − . − .

32 99 .

31% 99 .

04% 98 . .

90 +0 . − .

80 99 .

93% 99 .

57% 99 . . − .

56 +0 .

80 99 .

64% 99 .

42% 99 . .

27 +0 .

51 +0 .

80 99 .

40% 99 .

81% 99 . .

85 +0 .

80 +0 .

67 99 .

74% 99 .

38% 99 . . − . − .

79 99 .

93% 99 .

72% 99 . . − . − .

80 99 .

95% 99 .

69% 99 . . − . − .

12 99 .

91% 99 .

75% 99 . .

45 +0 .

70 +0 .

79 99 .

75% 99 .

13% 99 . .

28 +0 . − .

80 99 .

60% 98 .

48% 97 . . − . − .

80 99 .

78% 99 .

64% 99 . . − .

80 +0 .

80 99 .

81% 99 .

76% 99 . .

00 +0 .

80 +0 .

80 99 .

70% 99 .

56% 99 . .

97 +0 .

80 +0 .

13 99 .

31% 99 .

35% 99 . . − . − .

80 99 .

97% 99 .

87% 99 . . − . − .

31 99 .

74% 99 .

64% 99 . .

00 +0 .

73 +0 .

79 99 .

75% 99 .

63% 99 . . − .

58 +0 .

80 99 .

77% 99 .

62% 99 . .

17 +0 . − .

19 99 .

30% 99 .

38% 99 . . − .

56 +0 .

30 99 .

84% 99 .

73% 99 . .

09 +0 . − .

40 99 .

32% 99 .

67% 99 . .

72 +0 . − .

03 99 .

24% 99 .

26% 98 . . − .

80 +0 .

51 99 .

49% 99 .

40% 99 . . − .

33 +0 .

75 99 .

86% 99 .

65% 99 . .

46 +0 . − .

17 99 .

73% 99 .

48% 99 . .

25 +0 .

41 +0 .

76 99 .

73% 99 .

24% 99 . .

66 +0 . − .

70 99 .

13% 98 .

92% 98 . . − . − .

79 99 .

67% 99 .

27% 98 . .

28 +0 .

01 +0 .

80 99 .

80% 99 .

64% 99 . . − . − .

39 99 .

72% 99 .

79% 99 . .

00 +0 .

78 +0 .

53 99 .

59% 99 .

22% 99 . .

16 +0 .

80 +0 .

03 99 .

08% 98 .

73% 98 . .

00 +0 .

68 +0 .

67 99 .

67% 99 .

59% 99 . .

03 +0 . − .

78 99 .

73% 99 .

55% 99 . . − .

75 +0 .

34 99 .

73% 99 .

65% 99 . . − . − .

20 99 .

13% 98 .

99% 98 . , 074001 (2015), URL https://doi.org/10.1088/0264-9381/32/7/074001 .[2] B. P. Abbott et al., Classical and Quan-tum Gravity , 074001 (2015), URL https://doi.org/10.1088/0264-9381/32/7/074001 .[3] B. P. Abbott et al. (LIGO Scientiﬁc Col-laboration and Virgo Collaboration), Phys.Rev. Lett. , 061102 (2016), URL https://link.aps.org/doi/10.1103/PhysRevLett.116.061102 .[4] B. P. Abbott et al. (LIGO ScientiﬁcCollaboration and Virgo Collaboration),Phys. Rev. X , 031040 (2019), URL https://link.aps.org/doi/10.1103/PhysRevX.9.031040 .[5] A. Buonanno and T. Damour, Phys. Rev. D , 084006(1999), gr-qc/9811091. [6] A. Buonanno, Y. Pan, J. G. Baker, J. Centrella,B. J. Kelly, S. T. McWilliams, and J. R. vanMeter, Phys. Rev. D , 104049 (2007), URL https://link.aps.org/doi/10.1103/PhysRevD.76.104049 .[7] A. Taracchini, Y. Pan, A. Buonanno, E. Barausse,M. Boyle, T. Chu, G. Lovelace, H. P. Pfeiﬀer, andM. A. Scheel, Phys. Rev. D , 024011 (2012), URL https://link.aps.org/doi/10.1103/PhysRevD.86.024011 .[8] A. Taracchini, A. Buonanno, Y. Pan, T. Hin-derer, M. Boyle, D. A. Hemberger, L. E. Kid-der, G. Lovelace, A. H. Mrou´e, H. P. Pfeiﬀer,et al., Phys. Rev. D , 061502 (2014), URL https://link.aps.org/doi/10.1103/PhysRevD.89.061502 .[9] S. Babak, A. Taracchini, and A. Buonanno,Phys. Rev. D , 024010 (2017), URL https://link.aps.org/doi/10.1103/PhysRevD.95.024010 . [10] A. Boh´e, L. Shao, A. Taracchini, A. Buonanno, S. Babak,I. W. Harry, I. Hinder, S. Ossokine, M. P¨urrer, V. Ray-mond, et al., Phys. Rev. D , 044028 (2017), URL https://link.aps.org/doi/10.1103/PhysRevD.95.044028 .[11] L. London, S. Khan, E. Fauchon-Jones, C. Garc´ıa,M. Hannam, S. Husa, X. Jim´enez-Forteza,C. Kalaghatgi, F. Ohme, and F. Pannarale,Phys. Rev. Lett. , 161102 (2018), URL https://link.aps.org/doi/10.1103/PhysRevLett.120.161102 .[12] R. Cotesta, A. Buonanno, A. Boh´e,A. Taracchini, I. Hinder, and S. Ossokine,Phys. Rev. D , 084028 (2018), URL https://link.aps.org/doi/10.1103/PhysRevD.98.084028 .[13] A. Nagar, G. Riemenschneider, G. Pratten, P. Rettegno,and F. Messina, Phys. Rev. D , 024077 (2020), URL https://link.aps.org/doi/10.1103/PhysRevD.102.024077 .[14] B. P. Abbott et al., The Astrophys-ical Journal , 149 (2019), URL https://doi.org/10.3847/1538-4357/ab3c2d .[15] I. M. Romero-Shaw, P. D. Lasky, and E. Thrane,Monthly Notices of the Royal Astronomi-cal Society , 5210 (2019), ISSN 0035-8711, https://academic.oup.com/mnras/article-pdf/490/4/5210/30725872/stz2996.pdf, URL https://doi.org/10.1093/mnras/stz2996 .[16] S. Wu, Z. Cao, and Z.-H. Zhu, Monthly Notices ofthe Royal Astronomical Society , 466 (2020), ISSN0035-8711, https://academic.oup.com/mnras/article-pdf/495/1/466/33387647/staa1176.pdf, URL https://doi.org/10.1093/mnras/staa1176 .[17] I. M. Romero-Shaw, N. Farrow, S. Stevenson, E. Thrane,and X.-J. Zhu, Monthly Notices of the Royal Astro-nomical Society: Letters , L64 (2020), ISSN1745-3925, https://academic.oup.com/mnrasl/article-pdf/496/1/L64/33316514/slaa084.pdf, URL https://doi.org/10.1093/mnrasl/slaa084 .[18] B. P. Abbott et al. (LIGO Scientiﬁc Col-laboration and Virgo Collaboration), Phys.Rev. Lett. , 101102 (2020), URL https://link.aps.org/doi/10.1103/PhysRevLett.125.101102 .[19] I. Romero-Shaw, P. D. Lasky, E. Thrane, and J. C.Bustillo, The Astrophysical Journal , L5 (2020), URL https://doi.org/10.3847/2041-8213/abbe26 .[20] P. C. Peters, Phys. Rev. , B1224 (1964), URL https://link.aps.org/doi/10.1103/PhysRev.136.B1224 .[21] C. K. Mishra, K. G. Arun, and B. R.Iyer, Phys. Rev. D , 084040 (2015), URL https://link.aps.org/doi/10.1103/PhysRevD.91.084040 .[22] Y. Boetzel, C. K. Mishra, G. Faye, A. Gopakumar, andB. R. Iyer, Phys. Rev. D , 044018 (2019), URL https://link.aps.org/doi/10.1103/PhysRevD.100.044018 .[23] M. Ebersold, Y. Boetzel, G. Faye, C. K.Mishra, B. R. Iyer, and P. Jetzer,Phys. Rev. D , 084043 (2019), URL https://link.aps.org/doi/10.1103/PhysRevD.100.084043 .[24] N. Yunes, K. G. Arun, E. Berti, and C. M.Will, Phys. Rev. D , 084001 (2009), URL https://link.aps.org/doi/10.1103/PhysRevD.80.084001 .[25] E. A. Huerta, P. Kumar, S. T. McWilliams,R. O’Shaughnessy, and N. Yunes,Phys. Rev. D , 084016 (2014), URL https://link.aps.org/doi/10.1103/PhysRevD.90.084016 .[26] I. Hinder, F. Herrmann, P. Laguna, and D. Shoe-maker, Phys. Rev. D , 024033 (2010), URL https://link.aps.org/doi/10.1103/PhysRevD.82.024033 .[27] E. A. Huerta, C. J. Moore, P. Kumar, D. George, A. J. K.Chua, R. Haas, E. Wessel, D. Johnson, D. Glennon,A. Rebei, et al., Phys. Rev. D , 024031 (2018), URL https://link.aps.org/doi/10.1103/PhysRevD.97.024031 .[28] T. Hinderer and S. Babak, Phys.Rev. D , 104048 (2017), URL https://link.aps.org/doi/10.1103/PhysRevD.96.104048 .[29] Z. Cao and W.-B. Han, Phys.Rev. D , 044028 (2017), URL https://link.aps.org/doi/10.1103/PhysRevD.96.044028 .[30] X. Liu, Z. Cao, and L. Shao, Phys.Rev. D , 044049 (2020), URL https://link.aps.org/doi/10.1103/PhysRevD.101.044049 .[31] D. Chiaramello and A. Nagar, Phys.Rev. D , 101501 (2020), URL https://link.aps.org/doi/10.1103/PhysRevD.101.101501 .[32] A. Nagar, P. Rettegno, R. Gamba, and S. Bernuzzi(2020), 2009.12857.[33] A. Nagar, A. Bonino, and P. Rettegno (2021),2101.08624.[34] E. Barausse”, Physical Review D , 084024 (2010).[35] E. Barausse”, Physical Review D , 104027 (2011).[36] A. Buonanno, Y. Chen, and T. Damour,Phys. Rev. D , 104005 (2006), URL https://link.aps.org/doi/10.1103/PhysRevD.74.104005 .[37] A. Buonanno and T. Damour, Phys.Rev. D , 064015 (2000), URL https://link.aps.org/doi/10.1103/PhysRevD.62.064015 .[38] T. Damour, B. R. Iyer, and A. Nagar,Phys. Rev. D , 064004 (2009), URL https://link.aps.org/doi/10.1103/PhysRevD.79.064004 .[39] Y. Pan, A. Buonanno, R. Fujita, E. Racine, andH. Tagoshi, Phys. Rev. D , 064003 (2011), URL https://link.aps.org/doi/10.1103/PhysRevD.83.064003 .[40] L. Blanchet, G. Faye, B. R. Iyer, andS. Sinha, Classical and Quantum Gravity , 165003 (2008), ISSN 1361-6382, URL http://dx.doi.org/10.1088/0264-9381/25/16/165003 .[41] K. G. Arun, A. Buonanno, G. Faye, andE. Ochsner, Phys. Rev. D , 104023 (2009), URL https://link.aps.org/doi/10.1103/PhysRevD.79.104023 .[42] A. Buonanno, G. Faye, and T. Hinderer,Phys. Rev. D , 044009 (2013), URL https://link.aps.org/doi/10.1103/PhysRevD.87.044009 .[43] Caltech-Cornell-CITA, binary black hole simulation re-sults , .[44] C. M. Will and A. G. Wiseman,Phys. Rev. D , 4813 (1996), URL https://link.aps.org/doi/10.1103/PhysRevD.54.4813 .[45] D. Bini and T. Damour, Phys.Rev. D , 124012 (2012), URL https://link.aps.org/doi/10.1103/PhysRevD.86.124012 .[46] D. Shoemaker (LIGO Scientiﬁc Collab-oration), URL https://dcc.ligo.org/cgi-bin/DocDB/ShowDocument?docid=2974 (2010).[47] A. H. Mrou´e, M. A. Scheel, B. Szil´agyi, H. P.Pfeiﬀer, M. Boyle, D. A. Hemberger, L. E. Kid-der, G. Lovelace, S. Ossokine, N. W. Taylor,et al., Phys. Rev. Lett. , 241104 (2013), URL https://link.aps.org/doi/10.1103/PhysRevLett.111.241104 .[48] M. Boyle, D. Hemberger, D. A. B. Iozzo, G. Lovelace,S. Ossokine, H. P. Pfeiﬀer, M. A. Scheel, L. C. Stein,C. J. Woodford, A. B. Zimmerman, et al., Classi- cal and Quantum Gravity , 195006 (2019), URL https://doi.org/10.1088%2F1361-6382%2Fab34e2 .[49] B. S. Sathyaprakash and S. V. Dhurand-har, Phys. Rev. D , 3819 (1991), URL https://link.aps.org/doi/10.1103/PhysRevD.44.3819 .[50] L. S. Finn and D. F. Chernoﬀ,Phys. Rev. D , 2198 (1993), URL https://link.aps.org/doi/10.1103/PhysRevD.47.2198 .[51] C. Capano, Y. Pan, and A. Buonanno,Phys. Rev. D , 102003 (2014), URL https://link.aps.org/doi/10.1103/PhysRevD.89.102003 .[52] B. Moore, T. Robson, N. Loutrel, and N. Yunes, Clas-sical and Quantum Gravity , 235006 (2018), URL https://doi.org/10.1088/1361-6382/aaea00 .[53] B. Moore and N. Yunes, Classical andQuantum Gravity , 185003 (2019), URL https://doi.org/10.1088/1361-6382/ab3778 .[54] W.-B. Han, Z. Cao, and Y.-M. Hu, Classicaland Quantum Gravity , 225010 (2017), URL https://doi.org/10.1088/1361-6382/aa891b .[55] J. C. Bustillo, N. Sanchis-Gual, A. Torres-Forn ¨¦ , J. A.Font, A. Vajpeyi, R. Smith, C. Herdeiro, E. Radu, andS. H. W. Leong, Gw190521 as a merger of proca stars:a potential new vector boson of . × − ev (2021),2009.05376. [56] J. C. Bustillo, N. Sanchis-Gual, A. Torres-Forn ¨¦ ,and J. A. Font, Confusing head-on and precessingintermediate-mass binary black hole mergers (2020),2009.01066.[57] T. Islam, V. Varma, J. Lodman, S. E. Field, G. Khanna,M. A. Scheel, H. P. Pfeiﬀer, D. Gerosa, and L. E. Kidder,

Eccentric binary black hole surrogate models for the grav-itational waveform and remnant properties: comparablemass, nonspinning case (2021), 2101.11798.[58] A. Klein, Y. Boetzel, A. Gopakumar, P. Jetzer, andL. de Vittori, Phys. Rev. D , 104043 (2018), URL https://link.aps.org/doi/10.1103/PhysRevD.98.104043 .[59] S. Tanay, A. Klein, E. Berti, and A. Nishizawa,Phys. Rev. D , 064006 (2019), URL https://link.aps.org/doi/10.1103/PhysRevD.100.064006 .[60] S. Tiwari and A. Gopakumar, Phys.Rev. D , 084042 (2020), URL https://link.aps.org/doi/10.1103/PhysRevD.102.084042 .[61] I. Hinder, L. E. Kidder, and H. P. Pfeif-fer, Phys. Rev. D , 044015 (2018), URL https://link.aps.org/doi/10.1103/PhysRevD.98.044015 .[62] S. Habib and E. A. Huerta, Phys.Rev. D , 044016 (2019), URL https://link.aps.org/doi/10.1103/PhysRevD.100.044016https://link.aps.org/doi/10.1103/PhysRevD.100.044016