Enhancing the SEOBNRv5 effective-one-body waveform model with second-order gravitational self-force fluxes

We leverage recent breakthrough calculations using second-order gravitational self-force (2GSF) theory to improve both the gravitational-mode amplitudes and radiation-reaction force in effective-one-body~(EOB) waveform models. We achieve this by introducing new calibration parameters in the SEOBNRv5HM mode amplitudes, and matching them to the newly available 2GSF energy-flux multipolar data for quasicircular nonspinning binary black holes. We find that this significantly improves the SEOBNRv5HM energy flux, when compared to numerical-relativity (NR) simulations of binary black holes with mass ratios between 1:1 and 1:20. Moreover, we find that, once the conservative part of the SEOBNRv5 dynamics is calibrated, the SEOBNRv5HM waveform model with 2GSF information reproduces the binding energy of NR simulations more accurately, providing a powerful check of the consistency and naturalness of the EOB approach. While we only include nonspinning 2GSF information, the more accurate binding energy and energy flux carry over to the SEOBNRv5 waveform models for spinning binary black holes. Thus, our results improve the latest generation of SEOBNR waveform models (i.e., SEOBNRv5), which has been recently completed for use in the upcoming fourth observing (O4) run of the LIGO-Virgo-KAGRA Collaboration.


I. INTRODUCTION
During their first, second and third observing runs [1-5], the LIGO [6] and Virgo [7] gravitational-wave (GW) observatories have detected GWs from about ninety coalescences of compact binaries, composed of black holes (BHs) and/or neutron stars.Moreover, independent confirmations of these detections, as well as claims of new ones, were obtained in Refs.[8][9][10][11][12][13].All together, these results have firmly established the field of GW astronomy.The vast majority of the observed GW signals involve binaries with comparable masses [3,10,12,14] although a few signals show evidence for binaries with greater mass asymmetry [15,16], quantified by the ratio of the component masses (q = m 1 /m 2 ≥ 1 or = m 2 /m 1 ≤ 1, m 1 and m 2 being the primary and secondary masses in the binary, respectively).
As the number of GW detections is expected to increase in upcoming observation runs [14,17,18], so likely will the number of asymmetric coalescences.It is therefore important that waveform models used for detecting, identifying, and analyzing the GW signals faithfully represent binaries in the small-mass-ratio (SMR) regime.
Effective-one-body (EOB) theory [19][20][21][22][23] provides waveform models that can be used for the analysis of GW signals by combining results from various first-principle methods for solving the two-body problem in General Relativity, such as post-Newtonian (PN) theory [24][25][26][27] and numerical relativity (NR) [28][29][30].Moreover, EOB waveform models are constructed in such a way that they reduce to test-body motion around a black hole in the limit of vanishing .Thus, to improve fidelity of EOB waveform models, it is in principle straightforward and natural to incorporate results from SMR perturbation theory or gravitational self-force (GSF) theory [31] in the EOB formalism [32].
The inclusion of GSF results in EOB waveforms has, so far, been limited to the inclusion of higher PN test-body coefficients in the energy flux and gravitational-mode amplitudes (e.g., see Refs.[33,34,36,37]), fits of EOBmode amplitudes to the inspiral Teukolsky-multipolar modes [45,46] and the calibration of EOB Hamiltonians [33,47] to match the first-order GSF correction to the nonspinning shift of the innermost stable circular orbit (ISCO) [48].Various studies [49][50][51][52][53] considered improv-ing the EOB Hamiltonian by incorporating GSF corrections to the binding energy through use of the first law of binary mechanics [54,55].In the standard gauge [19,21] used in EOB waveform models, this leads to a gauge singularity at the lightring (or photon orbit) radius [51,52].This singularity can be cured by reformulating the EOB Hamiltonian in a different gauge [52,56].So far, this has not been implemented in any fully featured EOB waveform model employed for LVK analyses.Recently, a version of the TEOBResumS model has been produced [57] that incorporates most of the previously calculated GSF corrections to the EOB Hamiltonian without, however, calibrating it to NR simulations, and focusing only on the inspiral phase to avoid issues with the lightring divergence.
Recent calculations in GSF theory have provided the second-order GSF (henceforth, 2GSF) correction to the energy flux [58] as well as corresponding post-adiabatic waveforms [59].References [60,61] carried out a detailed comparison of these results against NR waveforms using a version of the TEOBResumS waveforms; among other things, 2GSF waveforms enabled a precise assessment of the accuracy of the TEOBResumS family in the SMR regime.However, these references did not seek to improve EOB waveform models through direct use of 2GSF information.In this work we will capitalize on the 2GSF breakthrough by directly incorporating 2GSF energy flux corrections in the latest generation of SEOBNR models.We will see that, quite interestingly, including these corrections does not only improve the waveform models at small mass-ratios, but also for comparable masses.
In this paper, we employ units such that G = c = 1.The component masses of a binary are denoted m 1 and m 2 with m 1 ≥ m 2 .The total mass is M = m 1 + m 2 , the reduced mass is µ = m 1 m 2 /M , while ν denotes the symmetric mass-ratio µ/M .

II. BASICS OF THE EFFECTIVE-ONE-BODY APPROACH
In the EOB formalism the dynamics of a compact binary is mapped onto that of an effective test mass (or test spin) in a deformed BH background, with the deformation parameter being the symmetric mass-ratio.The EOB approach builds semi-analytical inspiral-mergerringdown waveforms by combining analytical predictions for the inspiral and ringdown phases (from BH perturbation theory) with physically-motivated ansatzes for the plunge-merger stage.The EOB waveforms are then made highly accurate via a calibration to NR waveforms.The EOB formalism relies on three key ingredients: the EOB conservative dynamics (i.e., a two-body Hamiltonian), the EOB radiation-reaction forces (i.e., the energy flux) and the EOB gravitational modes.In this paper we shall limit to the inspiral-portion of the coalescence of nonspinning BHs.Here, we describe each of the main EOB ingredients, as necessary (e.g., see for details Ref. [40]).

A. EOB Hamiltonian
In the binary's center-of-mass frame, the motion is described by the orbital phase φ, the relative position r, the radial momentum p r and the angular momentum p φ .In the EOB formalism, the Hamiltonian H EOB , describing the conservative binary dynamics, is related to the effective Hamiltonian H eff , describing the dynamics of a test body in a deformed BH background, via the energy map [19] For nonspinning binaries, in the ν → 0 limit, H eff reduces to the Hamiltonian of a (nonspinning) test mass in a Schwarzschild background.The nonspinning EOB Hamiltonian was first derived in Refs.[19,20] with 2PN information.It was then extended to 3PN order in Ref. [21] and to 4PN order in Ref. [62].As of today, the 5PN [63][64][65] and 6PN terms [66,67] are partially known.In the non-spinning limit H eff reads: where p r * is the canonical momentum conjugate to the tortoise coordinate r * , The 5PN Taylor-expanded potential A is given by where u ≡ M/r and γ E 0.5772 is the Euler gamma constant.In Eq. ( 4), except for the log part, we replace the (partially known) coefficient of u 6 by the parameter a 6 .The latter is treated as a calibration parameter in the construction of the SEOBNRv5 waveform models [40] (and also in previous calibrations of EOB families).Furthermore, the resummed form of A, the potentials D and Q, and the spinning EOB Hamiltonian used for the calibration of the SEOBNRv5 waveform model can be found in Refs.[39,40].

B. EOB gravitational modes
The observer-frame's gravitational polarizations read where we denote with ι the binary's inclination angle (computed with respect to the direction perpendicular to the orbital plane), ϕ 0 the azimuthal direction to the observer, and Y −2 m (ι, ϕ 0 )'s the -2 spin-weighted spherical harmonics.For nonspinning binaries h m = (−1) h * −m , therefore one can restrict the discussion to ( , m) modes with m > 0.
The EOB modes for the entire coalescence (i.e., inspiral, plunge, merger and ringdown (RD)), can be written as: where t m match is defined as where t 22  peak is the peak of the (2, 2)-mode amplitude.The choice of a different attachment point for the (5,5) mode is discussed in Refs.[34,40].In the SEOBNRv5HM model, one imposes that [40] t 22 peak = t ISCO + ∆t 22  ISCO , where t ISCO is the time at which r = r ISCO , with r ISCO the ISCO radius of a Kerr BH [68] with the final mass and spin of the remnant object [69,70].The parameter ∆t 22 ISCO in Eq. ( 9) is the second EOB calibration parameter (the first being a 6 in the potential A of Eq. ( 4)), which is determined by minimizing the disagreement between EOB and NR waveforms throughout the inspiral-plunge stage.Henceforth, we focus only on the h insp−plunge m (t), and do not provide details on how the h merger−RD m (t) is constructed.
The inspiral-plunge EOB modes are written as where h F m (t)'s resum the PN-expanded GW modes for circular orbits in factorized form [71][72][73][74], while N m (t)'s are the nonquasicircular (NQC) corrections [75][76][77], aimed at incorporating relevant radial effects during the plunge (see below).The factorized inspiral modes are given by [73] The first factor, h N m , encodes the leading (Newtonian) order waveform, and its explicit expression is where d L is the luminosity distance, Y m is a scalar spherical harmonic, m is the parity of the mode, and the functions n m and c k (ν) are given by and Finally, v φ in Eq. ( 12) is given by with The second factor in Eq. ( 11), Ŝ m , is an effective source term.Depending on the parity of the mode, it is given by with The third factor in Eq. ( 11), T m , is an analytic resummation of the leading-order tail terms in the PN expansion [78], where Ω = ΩH EOB is the orbital frequency normalized by the total energy.The final two factors in Eq. ( 11) encode additional physical information included in the waveform from PN theory.The function f lm is expanded as The functions δ m , ρ m , and f S lm are fixed by requiring that the waveform at fixed orbital frequency Ω matches known analytical PN and test-body expressions.The full expressions used in SEOBNRv5HM are given in Appendix B of [40].

C. EOB equations of motion and radiation-reaction force
The EOB equations of motion read where the radiation-reaction (RR) forces F r and F φ are given by [23] where F EOB is the energy flux, which is given as a sum over ( , m) modes [73], Each of the ( , m)-mode contributions to the energy flux is obtained from the inspiral-plunge waveform modes, assuming quasi-circular orbits As discussed in Sec.II B, in the SEOBNRv5 model, as in other EOB variants, the inspiral-plunge quasi-circular waveform is enhanced by the NQC corrections N m (see Eq. ( 10)), which for the SEOBNR models take the form with the constants a h m i and b h m i chosen such that the EOB modes agree with NR modes at the point where the inspiral-plunge modes are matched to the merger-RD modes (see Eq. ( 7) above and Ref. [40]).
In the initial SEOBNR models, the NQC corrections were included in the energy flux (and RR forces) through Eq. ( 25) via an iterative procedure or fits (e.g., see Refs.[47,79]).However, starting from the SEOBNRv4 model [35], the inspiral-plunge modes in Eq. ( 25) only contain the factorized modes.The NQC corrections are included only in the gravitational polarization modes.The initial TEOBResumS models also included the NQC corrections in the energy flux, but then in subsequent versions the energy flux did not take them in.Recently, the TEOBResumsS model of Ref. [80] incorporates fits to the NQC corrections in the energy flux (and RR force) through Eq. (25).Finally, in the SEOBNRv5 waveform model used in this paper, the NQC corrections do not enter the RR forces.As we shall discuss below, in the nonspinning case, the calibration to 2GSF results improves the EOB energy fluxes considerably, thus the NQC corrections play a subdominant role.Furthermore, the inclusion of higher-order PN spin terms in the gravitational EOB modes, which were absent in the previous SEOBNRv4 model, reduces the disagreement between EOB and NR fluxes and pushes it mostly to the very late inspiral, where the effective test-body motion is almost unaffected by dissipative effects (see for details Ref. [40]).

III. BASICS OF THE GRAVITATIONAL SELF-FORCE APPROACH
The development of the GSF approach has traditionally been driven by the need to model GW emission from extreme-mass-ratio inspirals.This approach expands the metric of the binary around the metric of the primary in powers of = m 2 /m 1 (see Eq. ( 27) below).It is well known [81] that in order to get a waveform phase error that scales as O( ), the expansion of the metric perturbation must be carried through O(2 ) (i.e., at 2GSF).That is to say, the waveform phase error scales as O( 0 ) if second-order (in the mass-ratio) corrections are not included in the metric.
Practical 2GSF calculations have recently been carried out using a two-timescale framework [82].Within this approach the (multipolar) flux for a quasicircular, nonspinning binary was recently computed [58].This has since been combined with a calculation of the binding energy [83] to compute the associated inspiral dynamics and waveforms [59].Important additional details are available in Ref. [60].Here we will briefly review the calculation of the 2GSF flux.
Restricting to quasicircular orbits, we expand the metric of the binary as where g αβ is the Schwarzschild metric of the primary and φ is the orbital azimuthal angle of the secondary.The metric amplitudes, h n,m αβ , depend on the binary's slowly evolving orbital frequency Ω ≡ dφ/dt 2 .During the inspiral (i.e., sufficiently far from the ISCO), the orbital frequency, and thus the metric amplitudes, evolve on the slow RR timescale t RR ∼ 1/( Ω), whereas φ evolves on the fast orbital timescale t orb ∼ 1/Ω.The two-timescale framework treats t RR and t orb as independent.By substituting Eq. ( 27) into the Einstein field equations, we can split them into (i) a set of Fourier-domain partial differential equations for the amplitudes at fixed Ω, and (ii) evolution equations that determine Ω and φ as functions of time.Decomposing the metric perturbation onto a basis of tensor spherical harmonics and working in the Lorenz gauge [84,85], we arrive at the field equations which can be found explicitly in Ref. [82].Constructing the secondorder source, applying appropriate boundary conditions, and integrating the field equations required the development of a raft of new techniques and codes [86][87][88][89][90][91].
In our two-timescale scheme we assume the secondary follows a quasicircular orbit in which Ω = O( ) = 0.In order to satisfy the Einstein field equations through second order in the mass ratio, we consistently account for the nonzero Ω everywhere that it appears (which is in numerous places in the second-order field equation and second-order flux).However, the assumption Ω ∼ breaks down near the ISCO, and the expansions based on that assumption cause Ω to unphysically diverge at the ISCO.Due to the presence of Ω terms, the 2GSF flux computed using the above two-timescale expansion also diverges at the ISCO.This non-physical divergence can be removed by transitioning to a new expansion in the vicinity of the ISCO [92].The location where this transition occurs provides an estimate for where the inspiral two-timescale expansion breaks down.Reference [60] estimated that, for small ν, this breakdown occurs around where v ISCO Ω = 1/ √ 6 ≈ 0.408.The GSF energy flux is calculated from the ( , m) modes of h 1,m αβ + 2 h 2,m αβ at null infinity [93].With the polarizations expanded in -2 spin-weighted spherical harmonics (see Eq. ( 6)), the flux is given by Defining y = (m 1 Ω) 2/3 we can write the flux as an expansion in at fixed y as The symmetry of the physical binary system under the interchange of the labels m 1 ↔ m 2 suggests that, for comparable-mass binaries, it is natural to re-expand in the symmetric mass-ratio ν at fixed total mass M .This is known to improve agreement of perturbative results with NR simulations of comparable-mass binaries [94,95] and also it is natural for comparing with PN and EOB models.Defining x = (M Ω) 2/3 = v 2 Ω and using Figure 1.Comparison of the (5, 5)-mode flux extracted from NR simulations and GSF calculations at r/GM ≡ 1/v 2 Ω = 9.5 as function of the symmetric mass-ratio, ν.After subtracting the 1GSF (resp.2GSF) flux from the NR flux, we see the residual scales as ν 3 (ν 4 ), as expected.We further observe the significant improvement in the agreement with NR when using the re-expanded normalized GSF flux defined in Eq. (36).Similar results for the other mode fluxes considered in this work are given in Appendix B.
we can re-expand the flux as where Hereafter we shall use For interfacing with the EOB model, it is useful to define a re-expanded (Newtonian-)normalized flux with We hereafter define FGSF Using the re-expanded normalized flux also results in a significant improvement in the agreement between the GSF and NR fluxes for all modes other than the (2,2) mode.
The first calculation of the 2GSF flux was presented in Ref. [58], where remarkable agreement was found between the total GSF flux (first plus second order) and the flux computed from NR simulations of comparablemass binaries.Figure 4 of Ref. [58] showed that for the (2,2) and (3,3) mode the difference between the NR and GSF flux scaled as O(ν 4 ), as expected, over mass ratios ranging from 10:1 to 1:1.In Fig. 1 we show that this scaling also holds for the (5,5) mode across a wider range of mass ratios from 20:1 to 3:1.In Appendix B we show similar plots for the other modes considered in this work.This excellent agreement with NR gives us further confidence that the calculated 2GSF flux is capturing all contributions to the flux through O(ν 3 ).For this work we also computed the 2GSF flux at many more orbital frequencies than were previously presented in Ref. [58].
Figure 1 (and the related plots in Appendix B) also shows the improvement gained by using the re-expanded normalized GSF flux.We see that factoring out the leading Newtonian behaviour brings the GSF flux much closer to the NR flux, especially in the case of comparable masses.This factorization is also a key part of the resummation of the EOB flux as described in Sec.II.This further motivates us to incorporate the new 2GSF flux information in the EOB flux.

IV. MATCHING EOB AND GSF MULTIPOLAR FLUXES
In order to incorporate information from the 2GSF flux into the EOB flux, we need to compare the two in a gauge-invariant manner.In both cases, the energy flux is decomposed into −2 spin-weighted spherical modes.Consequently, we can compare the ( , m)-mode fluxes individually at a fixed value of the orbital frequency M Ω.By its nature, the GSF result is given as an expansion in powers of ν.We need to do the same with the EOB energy flux modes F EOB m .Combining Eqs. ( 11) and ( 25) gives Instead of expanding the full flux, it is more convenient to expand the flux normalized by its leading Newtonian contribution, This allows us to preserve the nonpolynomial dependence of F N m on ν in the final EOB result, and ensure its correct behaviour in the equal-mass limit.One potential complication in doing this is that F N m is not written directly in terms of Ω, but instead depends indirectly on Ω through v φ .However, expanding the dependence of v φ in powers of ν and v Ω (see Eqs. ( 16) and ( 19)), we find that Since in this comparison we are only interested in next-toleading order corrections in ν, we can thus safely replace v φ by v Ω everywhere in the expansion.Next, we need to expand the individual factors in Eq. ( 39) in powers of ν at fixed values of Ω. Starting with the effective source Ŝ m given in Eq. ( 18), we write At leading order, this is simply given by the well-known test-body result which is reproduced exactly by the EOB Hamiltonian (by construction).At next-to-leading order, the effective source for quasicircular inspirals depends on the linearin-ν correction to the EOB A potential.In principle, this contribution depends on the exact details of the implementation of the A-potential in the SEOBNRv5HM model, including any calibration of a 6 to NR results.(This is one of the routes through which the calibration of the A-potential becomes degenerate with calibration of the energy flux in the EOB RR force.)However, through use of the first law of binary mechanics [54], it is possible to directly compute the linear-in-ν correction to the A potential [49,50] in terms of the Detweiler redshift invariant z (1) [96], the exact value of which can be computed numerically in the GSF context [49,51,97].This gives eff , + m is even with We use interpolated data for the redshift z (1) generated with the code [97] from a previous work [52].The next step is to expand the tail term T m , given in Eq. ( 20).This has a hidden dependence on ν through Ω = ΩH EOB .We start by expanding H EOB : which straightforwardly gives Following [73], we note that the modulus square of the tail term can be written as We thus find the expansion of this term as with and We now write the expansion of the final factor in Eq. (38) as and compare to the re-expanded normalized GSF flux in Eq. (36).We can find the exact values of ρ (0) m and ρ (1) m in terms of the GSF flux by matching the two expressions for the normalized flux at fixed Ω order-by-order in ν, yielding and Equation ( 52), of course, matches the expression previously found in [73].The expression for ρ (1),GSF m is the new expression needed to incorporate the 2GSF flux into the EOB flux.
For including 2GSF information in the EOB mode amplitudes and energy flux, we focus on the 7 dominant ( , m) modes that are inlcuded in the SEOBNRv5HM model.For these modes we first determine the contributions to ρ (1), GSF 22 from applying Eq. ( 53) to the 2GSF fluxes, and the base EOB ρ (1), EOB 22 given by (54).In addition, we show the corrected ρ (1) 22 after adding the fitted correction (55a).The bottom panel shows the absolute difference between the GSF and EOB values with and without ∆ρ expanding ρ m in powers of ν,3 ρ (1),EOB 22 (1),EOB 21 (1),EOB 33 (1),EOB 32 (1),EOB 44 (1),EOB 43 (1),EOB 55 We augment the ρ (1),EOB m by adding an additional polynomial ∆ρ Ω starting at the lowest order in v 2 Ω not already included.The ∆ρ m are determined by fitting to the numerical ρ (1),GSF m results.While these extra terms take the form of higher-order PN terms, we emphasize that the goal here is not to estimate the next-order terms in the PN series (which would in general also contain log v Ω contributions).Instead, the goal is to capture Figure 3. Same as in Fig. 2, but now for the (3,3), (3,2), (4,4), (4,3), (5,5), and (2,1) modes.The panels on the left show the even parity modes, while the panels on the right show the odd parity modes.We note a different behaviour of the even and odd parity 2GSF modes close to the divergence induced by the two-time scale expansion in the proximity of the ISCO.
as much of the behaviour of the numerical ρ (1),GSF m data as possible.
To see how this fit works in practice, let us focus on the case of the (2, 2)-mode shown in Fig. 2.There are two complicating factors in performing the fit.The first is that the GSF data has finite numerical accuracy.This causes issues in the weak-field regime, where the EOB approximation is more accurate than the GSF data, and we thus run the risk of overfitting the numerical noise.The second complication is that the GSF data diverges at the Schwarzschild ISCO at v Ω = 1/ √ 6, where the inspiral two-timescale expansion breaks down (see Sec. III).This feature is not physical and should not be reproduced by the EOB flux.The combination of these complications with the small number of fitting parameters makes it most practical to manually fit the parameters using the residuals to judge both the location of the noise floor and the divergence.In the case of the (2, 2)-mode this produces ∆ρ Repeating the process for the six remaining modes (shown in Fig. 3) yields the following fits: Thus, in the GSF-augmented EOB model, ρ m in Eq. ( 21) is replaced as both when computing the EOB gravitational polarizations and RR force (taking ∆ρ (1) m = 0 for modes for which it has not been calculated, yet).
For two modes (the (3,2) and (4,3)) these fits contain higher PN terms than included in the corresponding ρ (0) m in previous SEOBNR models.For the sake of consistency, SEOBNRv5HM augments the corresponding ρ (0) m -terms to the appropriate PN order with corresponding 1GSF flux terms, which are known up to very high PN order [98].

V. IMPACT OF GSF INFORMATION ON THE SEOBNRV5HM MODEL ACCURACY
In this section we study the impact of including the 2GSF information in the energy flux and mode amplitudes on the overall faithfulness of the SEOBNRv5HM model developed in Ref. [40].We start with comparing the energy flux of the SEOBNRv5HM model to NR simulations from the SXS collaboration [99,100].Details of the simulations used can be found in Appendix C. In Fig. 4 we compare the energy flux of an NR simulation with mass-ratio q = m 1 /m 2 = 4 to the SEOBNRv5HM flux (24) with and without the 2GSF corrections.We see that even at the modest mass-ratio, the 2GSF corrections improve the agreement with the NR flux by a factor of a few across the frequencies spanned.This improvement is much more substantial than that obtained by including the NQCs in the energy flux [40].Moreover, we see that adding the NQCs to the flux on top of the 2GSF corrections leads to no substantial improvement except in the last fraction of a GW cycle before merger.
To understand how the improvement of the SEOBNRv5HM flux due to the 2GSF corrections scales with mass-ratio, in Fig. 5 we plot the same quantities as in Fig. 4, but now for different NR simulations at varying mass-ratio and fixed value of v Ω = 0.37.We again see that the 2GSF calibration improves the agreement with the NR flux by a factor of a few across the range of sampled mass-ratios, and provides a much more substantial improvement than merely including NQC corrections in the energy flux.At low ν, adding the NQC corrections on top of the 2GSF corrections provides an additional small improvement, while near equal masses the NQC corrections actually make the agreement with NR slightly worse.Naively, one might expect the relative error of the SEOBNRv5HM flux with the 2GSF calibration to scale with ν 2 .However, instead we see a relative error which is almost constant.This suggests that insufficient accuracy in the ρ (0),EOB m (i.e. the test-body flux) is the dominant source of error.However, note that while v Ω = 0.37 is smaller than v break Ω (28) for all mass-ratios, it is still close enough to the ISCO for corrections to the flux from the transition to plunge to be relevant.Such contributions would lead to an almost flat relative error scaling as ν 2/5 (see e.g.Fig. 15 in Appendix B).
It thus appears that the calibration of the SEOBNRv5HM flux against 2GSF results is successful in bringing the EOB flux closer to the NR flux.Ultimately, the true measure of the model is the waveforms that it produces.In Fig. 6 we compare the (2,2)-mode of the SEOBNRv5HM model (including 2GSF calibration) to waveforms from the pure GSF-based waveform from Ref. [59] (in its 1PAT1 form) 4 , and NR waveforms at three different mass ratios.The 2GSF waveforms are shown until the binary velocity reaches v break Ω .For the first part, the waveforms are visually indistinguishable, and only in the last ∼ 900M before merger we start to see differences, especially with the (inspiral only!) 2GSF waveforms with mass ratios 4 and 1.Indeed, when we look at the dephasing in the lower panel, the 2GSF waveforms for mass ratios 9.99 and 4 (1) stay below 0.1 (0.3) radians up to ∼ 900M before merger.Remarkably, the dephasing of the 2GSF waveform for mass ratio 9.99 is still below 0.1 radians when the velocity reaches v break Ω , highlighting the importance of including large mass-ratio corrections, while for mass ratios 4 and 1, the dephasing reaches 1 radian and ∼ 6 radians, respectively, at v break Ω .The dephasing of the SEOBNRv5HM waveforms is shown throughout the coalescence (i.e., during the inspiral, merger and ringdown stages) and it is much smaller than that of the GSF waveforms.This is expected since the SEOBNRv5HM waveforms have been calibrated to NR simulations [40].
To provide a more quantitative assessment of the impact of including the 2GSF information in the SEOBNRv5HM model, we calculate the mismatch (or unfaithfulness) between (2,2)-modes of the SEOBNRv5HM waveforms and of a set of NR waveforms with varying mass ratios using (e.g., see Ref. [40]) where we maximize over the relative shift in phase (δφ) and time (δt) between the two waveforms, while (•|•) denotes the noise weighted inner product [101,102] given by where S n (f ) is the one-sided power spectral density (PSD) of the detector noise, which we assume to be the design zero-detuned high-power noise of Advanced LIGO [103].
For each NR simulation, we calculate the mismatches for a range of total masses between 10M and 290M .In Fig. 7, we show a histogram of the mismatches for the three cases: the SEOBNRv5HM model without 2GSF corrections, the SEOBNRv5HM model with 2GSF corrections in the RR force (i.e., in the energy flux) and polarization modes, and the latter with the addition of NQC corrections in the RR force.Each waveform model is calibrated to NR by tuning the two calibration parameters introduced in Sec.II: a 6 , which appears in the Hamiltonian, and ∆t 22  ISCO , which determines the time at which the inspiral-plunge waveform is matched to the merger-ringdown one (for details see Ref. [40]).We stress that the NR calibration is done by demanding that the SEOBNRv5 (2,2) mode has mismacthes with NR below 10 −3 throughout the inspiral, merger and ringdown stages.
After calibration to the NR simulations, the histograms of the three models in Fig. 7 are very similar.To gain insight into this, it is instructive to compare the mismatches between the models and NR to the NR error.Generally, there are several contributions to the NR-error budget, including truncation error and error in extrapolating the waveforms to infinity. 5For mismatch studies, the former is more dominant [99], so we compute the mismatch between the highest and second highest NR resolution as a conservative measure of the NR error.From Fig. 7 one can see that the mismatches between the models and NR are close to the NR error.This suggests that part of the reason for minor differences between the models is that they are hitting the limits due to NR error.However, a potentially more important factor is that there is a large degree of degeneracy in the mismatch between changes in the RR force and changes in the Hamiltonian controlling the conservative dynamics (see Eqs. ( 22) and ( 53)).Consequently, the calibration of the Hamiltonian (through the calibration of the waveform modes) can largely compensate for imperfections in the dissipative sector of the EOB approach.In Fig. 8, we see the values of the two main calibration parameters, a 6 and ∆t 22  ISCO of the SEOBNRv5HM models with and without 2GSF corrections in the RR force and polariza-tion modes.The presence of the 2GSF corrections has clear impact on the calibration coefficients.This implies that the two calibrated models have somewhat different dynamics; however, those differences do not lead to appreciable differences in the corresponding waveforms, as can be seen in Fig. 7.
Since the calibration parameter a 6 controls part of the EOB A-potential, the two SEOBNRv5HM models with and without 2GSF information have different Hamiltonians, and therefore differ in their binding energy.The latter is given by In Fig. 9 we compare the SEOBNRv5HM binding energy to the one extracted from NR simulations from Ref. [104].The SEOBNRv5HM with 2GSF corrections reproduces the NR binding energy much more faithfully, staying within the NR error estimates until roughly one GW cycle before merger.This improvement persists even when modeling aligned-spin binaries despite only adding 2GSF corrections to the nonspinning part of the waveform and RR force.In Fig. 10, we compare the SEOBNRv5HM binding energy for models with and without the 2GSF information to that of a set of spin-aligned NR waveforms [104] at a fixed value v Ω = 0.447.Without the 2GSF calibration, the binding energy can be off by as much as 2.5% especially at high values of the effective spin χ eff = (χ 1 m 1 + χ 2 m 2 )/M .However, with the 2GSF calibration, we find deviations from the NR binding energy to be at the sub-percent level, with an average difference of just 0.24%.bind inferred from NR simulations, and E EOB bind from the SEOBNRv5HM models with and without 2GSF corrections.The shaded area shows the estimated error on the NR binding energy in the case q = 1, which is taken as representative for the general NR error.The ticks on the top x-axis show the number of GW cycles before merger for the q = 1 NR simulation.

VI. DISCUSSION
In this paper, we have enhanced the accuracy of the (factorized) gravitational modes used in the SEORBNRv5 models of Refs.[39][40][41] by calibrating them to nonspinning, quasi-circular 2GSF multipolar data of Ref. [58].This calibration affects also the EOB radiation reaction (RR) force driving the dynamical evolution of the binary black holes in the SEORBNRv5 model.
By direct comparison of the energy flux in the SEORBNRv5HM model to that extracted from NR simulations, we have confirmed in Figs. 4 and 5 that the 2GSF calibration of the flux leads to a significant improvement in the faithfulness of the SEORBNRv5HM flux.In particular, the improvement seems to make the inclusion of NQC corrections in the RR force subdominant, and limited to the very late inspiral (plunge), where the effective testbody motion is almost unaffected by dissipative effects.
Furthermore, when looking at the mismatches between the SEORBNRv5HM and NR waveforms in Fig. 7, the inclusion of the 2GSF calibration seems to only have a marginal impact on the waveform mismatches after calibration to NR.If anything, this is a testament to the effectiveness of the SEORBNRv5 Hamiltonian's calibration to NR, which is obtained by demanding that the mismatches of the SEORBNRv5 inspiral-merger-ringdown (2,2) modes are below 10 −3 .Since the waveforms are computed using the EOB equations of motion, which depend on the conservative and dissipative dynamics, the mismatches have significant degeneracy between the calibra- tion terms in the EOB Hamiltonian and in the RR force (notably the 2GSF terms in the energy flux).This is one reason why the flux calibration terms that we have added to the EOB flux could not have been added through the NR calibration performed in Ref. [40] 6 .This degeneracy also means that calibrating SEOBNRv5HM to NR with and without the 2GSF calibration of the flux leads to a different EOB Hamiltonian.The Hamiltonian itself however is supposed to correspond to a gauge invariant observable of the binary, the binding energy.Comparing the EOB binding energy to results extracted from NR simulations in Fig. 9, we find that the Hamiltonian calibrated with the 2GSF information included reproduces the NR bind- 6 One could explore in the future the possibility of calibrating directly the Hamiltonian from the binding energy extracted from NR simulations instead of doing it indirectly using the waveforms [104].However, this procedure would require the computation of the binding energy for the entire set of 442 aligned-spin SXS NR waveforms used to calibrate the SEOBNRv5 model.
ing energy much more faithfully than the Hamiltonian calibrated without.This is a significant consistency test of the SEOBNRv5HM model, and one that extends to binary BHs with spins, as well (see Fig. 10).So, while adding the 2GSF information to the SEOBNRv5HM model does not necessarily improve the faithfulness of the corresponding waveforms in the regime where they are calibrated to NR, it does improve the overall consistency and naturalness of the model.This gives us greater confidence that the SEOBNRv5HM model will remain (somewhat) faithful to NR when extrapolated beyond the calibration region, in particular for higher mass ratios.
In this work we focused on adding 2GSF corrections to the nonspinning sector of the SEOBNRv5HM waveforms.However, numerical results are also available for corrections to the 2GSF flux linear in either the primary or secondary spin [58,105,106].In principle, the matching procedure employed here can also be used to calibrate the SEOBNRv5HM modes to these data.We will leave the implementation of this to future work.
A limiting factor in this work has been that the 2GSF multipolar flux data we used do not include corrections due to the transition from inspiral to plunge, causing it to diverge at the ISCO.This has limited our ability to calibrate the RR force and gravitational modes in the strong-field regime.Inclusion of these terms could lead to further improvements of our results, and will be addressed once new 2GSF data becomes available.The accurate mismatches and small binding-energy disagreements with NR found in Sec.V have been obtained using the SEOBNRv5 nonspinning Hamiltonian of Ref. [39,40].Here, we want to understand whether those results are somehow tied to the particular structure (or resummation) of the Hamiltonian and the PN content.Thus, we repeat some of the analyses using the non-spinning Hamiltonian from the previous SEOBNR family, SEOBNRv4 [33,47].
Figure 11 is similar to Fig. 7, but now we compute the mismatches between the SEOBNRv4.5HMand NR (2,2) waveforms, where the SEOBNRv4.5HMmodel is constructed calibrating the SEOBNRv4 Hamiltonian, and using the SEOBNRv5HM RR force and gravitational modes with and without the 2GSF information.We find that although the mismatches are a bit higher than for the SEOBNRv5 ones, there is actually a noticeable improvement when including the 2GSF information.
In Fig. 12 we revisit Fig. 9 with the SEOBNRv4.5HMmodel.We again see that calibrating the model using the 2GSF information in the RR force and gravitational modes leads to a more accurate recovery of the binding energy, albeit less striking than in the case of the SEOBNRv5HM model.

Appendix B: Extended comparison between GSF and NR multipolar fluxes
In this appendix we show further comparisons between GSF and NR multipolar energy fluxes.The comparisons for the (3,2), (3,3), (4,3) and (4,4) modes are shown in Fig. 13.The residual after subtracting the GSF flux from the NR flux clearly shows the expected scaling.The comparisons for the (2,1) and (2,2) modes are presented in Fig. 14.The scaling of the residuals for these modes is less clear for the reasons given below and in the caption of the figure.The SXS datasets used to make Figs. 1, 13, and 14 are given in Table I.
Figure 14 shows that for the (2,2) mode, the flux's scaling with ν at v Ω = 0.324 is likely affected by the transition to plunge.Such a transition occurs over a frequency interval ∼ ν 2/5 /M on a time scale ∼ M/ν 1/5 , during which the small parameter ν 1/5 enters into the SMR expansion [92].The behavior of the energy flux in this case can be obtained by combining Eq. (10) of Ref. [108] with Eqs. ( 22) and (23) of Ref. [60].In those equations, we define R = (r − 6M )/ν 2/5 and ∆Ω = (Ω − Ω ISCO )/ν 2/5 , where r is the orbital separation, such that R ∼ M and ∆Ω ∼ 1/M in the transition regime.The cited equations then give us Here E is the specific binding energy, meaning it is related to the flux by F = −νdE/dt, which implies F ∼ ν 2 +ν 12/5 +. ... (In all of these schematic equations, the reader should understand that powers of ν come with ∆Ω-dependent coefficients, which we omit to avoid introducing additional notation.)Figure 15 repeats the comparison in Fig. 14 at a frequency closer to the ISCO (v Ω = 0.370), and there we see clear evidence of the ν 12/5 term appearing in the flux.After the first-order flux is subtracted from the NR flux, the residual scales as ν 3 .For the (2,1) mode after the second-order flux is also subtracted, the residual is within the magnitude of the oscillations in the NR data and so the scaling of the residual is less clear.For the (2,2) mode the residual does not clearly follow the dash-dotted (blue) ν 4 curve as it is likely contaminated by effects related to the transition to plunge.The effect of this transition is clear for orbital radii close to the ISCO, as one can see in Fig. 15.
Figure 15.The same as the left panel of Fig. 14, but now for r/GM ≡ 1/v 2 Ω = 7.3.The dashed (yellow) curve is a ν 12/5 reference, which is the expected scaling for the flux near the transition to plunge.

mFigure 2 .
Figure 2. The top panel shows the numerical values of ρ

Figure 4 .Figure 5 .
Figure 4.The normalized energy flux F compared between a quasicircular nonspinning NR simulation at mass-ratio q = 4, and the SEOBNRv5HM flux with and without 2GSF calibration.In addition we show what happens when the NQC corrections are included in the SEOBNRv5HM flux.The top panel shows the full flux as a function of vΩ.The circles indicate the merger (peak of |h insp−plunge 22 |), and the other markers indicate 1 (diamond), 2 (square), 4 (triangle), and 10 (hexagon) GW cycles before merger.The bottom panels show the relative difference between NR and the SEOBNRv5HM fluxes.The vertical dashed line indicates the fixed frequency used for Fig. 5.

Figure 6 .
Figure 6.Top panel: Comparison of the (2,2)-mode waveforms from NR, the GSF 1PAT1 model and the SEOBNRv5HM model at three different mass-ratios q.The waveforms at each mass-ratio are aligned at the start of the NR waveforms using the procedure described in Ref.[40].The last −900M before the peak of the (2,2)-mode are shown magnified.Bottom panel: Dephasing of the 1PAT1 and SEOBNRv5HM models relative to the NR waveforms.The GSF waveforms are truncated at v break Ω

Figure 7 . 50 Figure 8 .
Figure 7.A histogram of the mismatches of NR versus SEOBNRv5HM models with and without 2GSF and NQC corrections.As an indicator of the NR error, the mismatch of the highest resolution NR waveforms against the next higher resolution is shown in gray.

Figure 9 .
Figure 9. Relative difference between the binding energy E NRbind inferred from NR simulations, and E EOB bind from the SEOBNRv5HM models with and without 2GSF corrections.The shaded area shows the estimated error on the NR binding energy in the case q = 1, which is taken as representative for the general NR error.The ticks on the top x-axis show the number of GW cycles before merger for the q = 1 NR simulation.

Figure 10 .
Figure 10.Comparison of the SEOBNRv5HM binding energy to NR data for spin-aligned waveforms at a fixed value vΩ = 0.447.The top (bottom) panel shows the binding energy from the SEOBNRv5HM model with (without) 2GSF information.

Figure 11 .Figure 12 .
Figure 11.Similar to Fig.7, but now we show the mismatches for the model built by calibrating the SEOBNRv4 Hamiltonian, and using the SEOBNRv5HM RR force and gravitational modes with and without the 2GSF information.We label this model SEOBNRv4.5HM.
Appendix A: Independence of results on the specific EOB Hamiltonian

Figure 14 .
Figure 14.Comparison of the fluxes extracted from NR simulations and GSF calculations for the l = 2 modes at r/GM ≡ 1/v 2 Ω = 9.5 as a function of the symmetric mass-ratio, ν.After the first-order flux is subtracted from the NR flux, the residual scales as ν 3 .For the (2,1) mode after the second-order flux is also subtracted, the residual is within the magnitude of the oscillations in the NR data and so the scaling of the residual is less clear.For the (2,2) mode the residual does not clearly follow the dash-dotted (blue) ν 4 curve as it is likely contaminated by effects related to the transition to plunge.The effect of this transition is clear for orbital radii close to the ISCO, as one can see in Fig.15.

Table I .
Details of the SXS simulations used in Figures throughout the paper.