{'fraction_non_alphanumeric': 0.039513082959860654, 'fraction_numerical': 0.018924790105094232, 'mean_word_length': 4.408340389500423, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Effective-one-body waveforms calibrated to numerical relativity simulations: coalescence of non-precessing, spinning, equal-mass black holes

We present the first attempt at calibrating the effective-one-body (EOB) model to accurate numerical-relativity simulations of spinning, non-precessing black-hole binaries. Aligning the EOB and numerical waveforms at low frequency over a time interval of 1000M, we first estimate the phase and amplitude errors in the numerical waveforms and then minimize the difference between numerical and EOB waveforms by calibrating a handful of EOB-adjustable parameters. In the equal-mass, spin aligned case, we find that phase and fractional amplitude differences between the numerical and EOB (2,2) mode can be reduced to 0.01 radians and 1%, respectively, over the entire inspiral waveforms. In the equal-mass, spin anti-aligned case, these differences can be reduced to 0.13 radians and 1% during inspiral and plunge, and to 0.4 radians and 10% during merger and ringdown. The waveform agreement is within numerical errors in the spin aligned case while slightly over numerical errors in the spin anti-aligned case. Using Enhanced LIGO and Advanced LIGO noise curves, we find that the overlap between the EOB and the numerical (2,2) mode, maximized over the initial phase and time of arrival, is larger than 0.999 for binaries with total mass 30-200Ms. In addition to the leading (2,2) mode, we compare four subleading modes. We find good amplitude and frequency agreements between the EOB and numerical modes for both spin configurations considered, except for the (3,2) mode in the spin anti-aligned case. We believe that the larger difference in the (3,2) mode is due to the lack of knowledge of post-Newtonian spin effects in the higher modes.


I. INTRODUCTION
Coalescing black-hole binaries are among the most promising sources for the current and future laserinterferometer gravitational-wave detectors such as LIGO/Virgo [1][2][3] and LISA [4].
In general relativity, black holes are defined only by their masses and spins; thus generically, a black-hole binary depends on eight parameters (m 1 , S 1 , m 2 , S 2 ). Hence, when black holes carry spins, it is expected that tens of thousands of waveform templates may be needed in order to extract the gravitational-wave signal from the noise using matched-filtering techniques. Considering the high computational cost of running numerical-relativity simulations of spinning binary black holes (on the order of 75, 000 CPU hours for 15 orbits, moderate spins and mild mass ratios) and the large binary parameter space, it will be impractical for numerical relativity alone to provide data analysts with a template bank. The work at the interface between analytical and numerical relativity [5][6][7][8][9][10][11][12][13][14][15] has demonstrated the possibility of modeling analytically the dynamics and the gravitational-wave emission of coalescing nonspinning black holes, thus providing data analysts with preliminary analytical template families to be used for the searches. The next important step is to extend those studies to spinning, precessing black holes. The present paper represents the first attempt in this direction, although limited to non-precessing waveforms, within the effective-one-body (EOB) formalism [16][17][18] of spinning black holes [19,20]. Recently, Ref. [21] constructed a template family of spinning, non-precessing black-hole binaries using a phenomenological approach, where the numerical-relativity waveforms are fitted to templates which resemble the post-Newtonian (PN) expansion, but in which the coefficients predicted by PN theory are replaced by many arbitrary coefficients calibrated to numerical simulations.
The first EOB Hamiltonian that included spin effects was worked out in Ref. [19]. In Ref. [22], the authors used the non-spinning EOB Hamiltonian augmented with PN spin terms to carry out the first exploratory study of the dynamics and gravitational radiation of spinning blackhole binaries during inspiral, merger and ringdown. Subsequently, Ref. [20] extended the EOB Hamiltonian of Ref. [19] to include next-to-leading order spin-orbit couplings. In those descriptions, the effective particle is endowed not only with a mass µ, but also with a spin σ. As a consequence, the effective particle interacts with the effective Kerr background (having spin S Kerr ) both via a geodesic-type interaction and via an additional spindependent interaction proportional to its spin σ. The EOB Hamiltonian developed in Refs. [19,20] (with nonspinning PN couplings through 3PN order) clarified several features of spinning two-body dynamics. However, as we shall discuss below, it is not straightforward to extend this Hamiltonian to include higher-order nonspinning PN couplings, such as the 4PN or 5PN adjustable parameters recently calibrated to numerical-relativity simulations [14,15]. Moreover, the EOB Hamiltonian of Ref. [20], based on an ad hoc test-particle limit, does not reduce to the Hamiltonian of a spinning test particle in Kerr spacetime. More recently, Ref. [23] derived the canonical Hamiltonian of a spinning test particle in a generic curved spacetime at linear order in the particle spin. The construction of an improved EOB Hamiltonian based on the results of Ref. [23] is currently under development. Despite the limitations mentioned above, the EOB Hamiltonian of Ref. [20] is an excellent starting point for exploring the calibration of numerical-relativity waveforms of spinning black holes within the EOB formalism. Thus, we have used it in this first exploratory study, augmenting it with a few adjustable parameters that we shall calibrate to two numerical-relativity simulations. For the EOB non-conservative dynamics, we use the gravitational-wave energy flux which includes spin effects and which has been computed using the factorized multipolar waveforms of Refs. [24,25].
The two numerical-relativity simulations we shall use describe the evolution of equal-mass, equal-spin, nonprecessing black-hole binaries. They are produced by the pseudospectral code SpEC of the Caltech-Cornell-CITA collaboration. In these two configurations, the spins are either aligned ("up-up", or UU) or anti-aligned ("downdown", or DD) with the orbital angular momentum, and have dimensionless magnitude χ 1 = χ 2 = 0.43655 for the UU configuration, and 0.43757 for the DD configuration. The UU simulation lasts for about 28 gravitationalwave cycles or until t = 3250M , and stops about three gravitational-wave cycles before merger. The DD simulation lasts for about 22 gravitational-wave cycles or until t = 2500M , and contains the full inspiral, merger, and ringdown waveform. Detailed information on the numerical simulation of the DD configuration can be found in Ref. [26].
For zero spin, the EOB-model considered here agrees with the waveform of the equal-mass non-spinning binary black hole [27] to a similar degree as the model constructed in our earlier work [15]. It differs from the model presented in [15] by its modeling of the energy flux of gravitational-wave radiation.
This paper is organized as follows. In Sec. II, we describe the spin EOB model adopted in this paper, including its adjustable parameters. In Sec. III, we calibrate the spin EOB model to the numerical simulations, and discuss the impact of our results on data analysis. Finally, Sec. IV summarizes our main conclusions.

II. EOB MODEL FOR SPINNING BLACK-HOLE BINARIES
In this section, we describe the spin EOB model adopted in our study and its adjustable parameters. Henceforth we use natural units G = c = 1. We use m i , X i , P i , and S i to denote the mass, the position vector, the linear momentum vector, and the spin vector of the i-th body. We work in the center-of-mass frame defined by P 1 + P 2 = 0. The two body system is described by the relative position R ≡ X = X 1 − X 2 and the relative linear momentum P ≡ P 1 = −P 2 . For convenience, we define reduced variables where M ≡ m 1 + m 2 and µ ≡ m 1 m 2 /(m 1 + m 2 ).

A. EOB conservative dynamics
Following Refs. [19,20], we assume that the effective particle in the EOB description is endowed not only with a mass µ, but also with a spin σ. As a consequence, the effective particle interacts with the effective Kerr background (having spin S Kerr and mass M ) both via a geodesic-type interaction and via an additional spindependent interaction proportional to its spin σ. We define the Kerr-like parameter a as a ≡ S Kerr /M , where S Kerr denotes the modulus of the deformed-Kerr spin vector S Kerr . Following Ref. [19], we write the effective Kerr contravariant metric components in a fixed Cartesian-like coordinate system. This is done by introducing and and writing the contravariant metric components as The Taylor approximants to the coefficients A(R) and D(R) can be written as The functions A k (r) and D k (r) all depend on the symmetric mass ratio ν ≡ µ/M through the ν-dependent coefficients a i (ν) and d i (ν). These coefficients are currently known through 3PN order (i.e. up to k = 4) and can be read off from Eqs. (47) and (48) in Ref. [13]. It is worth noticing that although through 3PN order the Padé approximant to the function ∆ t (R) of Eq. (5a) does not pose any problem [20], when including 4PN and 5PN order coefficients, the Padé approximant develops poles for several spin values a. In particular, poles are present at large separation when a > 0.7M and the 4PN and 5PN order coefficient a 5 and a 6 are included. 2 Those poles could be regularized by adding in A k (r) higherorder spin terms a 2ã 3 (ν)/r 5 , a 2ã 4 (ν)/r 6 and choosing for the coefficientsã 3 (1/4) andã 4 (1/4) negative large values (∼ −100). Since in this first exploratory study we investigate only numerical simulations of moderate spins, we do not include any regularization of the poles, and consider only the 4PN order coefficient a 5 . In the nonspinning case [9,10,14,15], the coefficient a 5 plays an important role in improving the agreement between the EOB and numerical waveforms. Here, we choose for a 5 the value obtained by taking the nonspinning limit of the spin EOB model and calibrating it to the equalmass black-hole waveform of [27], following [15]. In this way, we obtain a 5 (1/4) = 1.775; thus, a 5 is no longer an adjustable parameter in the spin EOB model. In Eq. (5a), we choose m = 1 and n = 4 so that ∆ t (R)/R 2 in the limit of a → 0 reduces to the nonspinning A(R) used in Refs. [9,10,14,15], and we choose the same 3PN D(R) function used in those references. Therefore, in the spin EOB model, we have with and where χ ≡ a/M and a 4 (ν) = 94/3 − 41/32π 2 ν. Making use of Eqs. (3) and (4), we can derive The EOB effective Hamiltonian reads [19,20] H eff (R, P , S 1 , S 2 ) = H eff Kerr (R, P , S Kerr ) + H eff part (R, P , σ) , and where Q 4 (P i ) is a quartic-momentum term at 3PN order independent of spins [18] and L ≡ R × P is the orbital angular momentum. In this paper, as a first attempt, we use the same spin coupling for the spin σ suggested in Ref. [20], even though it does not reduce to the spinning test-particle limit [23] at PN orders higher than 2.5PN.
In order for H eff to match the PN-expanded spin-orbit Hamiltonian through 2.5PN order, we need to require that the sum of the spin-orbit couplings of H eff Kerr and H eff part gives and where the two effective gyro-gravitomagnetic ratios g eff S and g eff S * read [20] g eff Here a(ν) and b(ν) are two gauge parameters related to the freedom of applying a canonical transformation involving spin variables. If we knew the exact Hamiltonian, the choice of these parameters should not affect the physics of the EOB model. However, since we start with an approximate Hamiltonian that reproduces the spin-orbit couplings only through 2.5PN order, we expect the EOB model to depend on the choice of a(ν) and b(ν).
Considering the structure of the gyro-gravitomagnetic ratios, such dependence should start at 3.5PN order as a spin-orbit coupling term. Because of this dependence, a(ν) and b(ν) can be used as adjustable parameters. Moreover, in order for H eff to match the PN-expanded spin-spin Hamiltonian through 2PN order, the simplest choice is to require that the Kerr spin [20] S Kerr = S + S * .
As a consequence, Eq. (14) implies To include higher-order spin-spin contributions in the EOB effective Hamiltonian, we introduce a 3PN spinspin term whose coefficient a 3PN SS is currently unknown and can be used as an adjustable parameter H eff (R, P , S 1 , S 2 ) = H eff Kerr (R, P , S Kerr ) + H eff part (R, P , σ) Finally, the EOB Hamiltonian is In summary, in this first exploratory study, we choose to employ only two adjustable parameters 3 : b(ν) which introduces a spin-orbit term at 3.5PN order, and a 3PN SS which introduces a 3PN spin-spin term. As we shall see, these two adjustable parameters are sufficient to reduce the phase and amplitude differences between EOB and numerical waveforms of the UU and DD configurations to (almost) the numerical error. The remaining flexibility of the spin EOB model can be exploited in the future when numerical relativity simulations of other spin configurations will become available. Thus, for the rest of the paper, we set a(ν) in Eq. (16a) to zero. Within the Hamiltonian approach, radiation-reaction effects can be incorporated into the dynamics in the following way [17,22]: Here, F i denotes the non-conservative force, which is added to the evolution equation of the (relative) momentum to take into account radiation-reaction effects. Following Ref. [22], we use 4 where Ω is the orbital frequency and L is the orbital angular momentum. The gravitational-wave energy flux dE/dt is obtained by summing over the gravitationalwave modes (l, m) as which reduces to the following expression for circular equatorial orbits in the adiabatic approximation: whereΩ is the reduced orbital frequencyΩ ≡ M Ω. We shall define the EOB waveforms h ℓm in Sec. II C. The equations of motion for the spins are simply obtained through the equations In the nonspinning case, it is useful [15,28] to replace the radial momentum P R with P R * , the conjugate momentum of the EOB tortoise radial coordinate R * : dR * /dR = √ D/A. This replacement improves the numerical stability of the EOB equations of motion because P R diverges when approaching the zero of A(r) (the EOB event horizon) but P R * does not. Therefore, in the spinning EOB Hamiltonian, we similarly choose to use the conjugate momentum to the tortoise radial coordinate of the ν-deformed Kerr geometry: In the limit a → 0, Eq. (28) reduces to the nonspinning EOB tortoise coordinate defined above. In the limit ν → 0, Eq. (28) reduces to the tortoise coordinate of the Kerr geometry: dR * /dR = (R 2 + a 2 )/∆. Since the EOB Hamiltonian and Hamilton equations are written in Cartesian coordinates, some algebra is needed to rewrite them to include this transform of the radial coordinate.
In Appendix A, we write down explicitly the transformed EOB Hamiltonian and Hamilton equations in Cartesian coordinates. In particular, Eqs. (21) and (22) should be replaced by Eqs. (A4) and (A5).
Initial conditions for the Hamilton equations are constructed following the prescription of Ref. [22], which provided post-circular initial data for quasi-spherical orbits when neglecting spin-spin and next-to-leading order spin-orbit effects. Note that exact circular orbits cease to exist in the conservative dynamics when spin-spin and next-to-leading order spin-orbit effects are present, except for special configurations in which the spins are aligned or antialigned with the orbital angular momentum. We start each evolution at a large initial separation of 50M . The EOB trajectory is sufficiently circularized when reaching a separation of ∼ 16M , where numerical waveforms start. In this way, we remove the residual eccentricity in the EOB trajectory due to imperfect initial conditions, while physical eccentricity due to spin effects is preserved.
As a final remark, the spin variable in the EOB model is the constant spin variable, i.e., its magnitude does not change during precession [29]. We identify it with the spin variable in the numerical simulation, which also remains constant during the evolution [26].

B. Characteristics of EOB orbits for spinning, non-precessing black holes
Here we investigate certain properties of the spin EOB Hamiltonian that are crucial when building the complete EOB model. Specifically, we check the existence and behavior of the innermost stable circular orbit (ISCO), the light ring (photon orbit) and the maximum of the EOB orbital frequency. This study is restricted to circular equatorial orbits in the spin aligned or anti-aligned cases. For convenience, we consider the EOB Hamiltonian written in spherical coordinates; we fix θ = π/2 and set the conjugate momentum P θ = 0.
The ISCO position is obtained by solving ∂H(R, P R * = 0, P φ )/∂R = 0, ∂ 2 H(R, P R * = 0, P φ )/∂R 2 = 0 where P R * and P φ are conjugate momenta of the tortoise radial coordinate R * and the orbital phase φ, respectively. In the nonspinning limit, we find the following ν-correction of the ISCO frequency relative to the Schwarzschild valuê where O(ν 4 ) terms contribute less than 1% to the correction. In the test-particle limit, the coefficient of the linear ν-correction term, 0.9837, should be compared to the recently available self-force result [30] (transformed to the gauge condition and mass convention used in the EOB formalism by Ref. [31]) of 1.2513. The relative difference of 21% is due to the fact that our nonspinning EOB Hamiltonian, although calibrated to equal-mass numerical simulations, does not capture all the ν-dependence correctly at 4PN order. 5 The improved spin EOB Hamiltonian [32] will incorporate consistently the self-force result (e.g., Ref. [31]) and can be better constrained by new numerical simulations. In the spin aligned or anti-aligned case, we find that the ISCO exists for all spin magnitudes. However, in the spin aligned case, when a > 0.8M , the ISCO radius (frequency) starts to increase (decrease) with increasing a. This is contrary to the monotonic dependence of the ISCO radius (frequency) on the spin magnitude in the test-particle limit. This unusual behaviour will be overcome by the improved spin EOB Hamiltonian of Ref. [32]. Nevertheless, since this problem occurs only at extreme spin magnitudes and here we have numerical waveforms of moderate spins (|a| < 0.5M ), we choose to use this spin EOB Hamiltonian in the current calibration. 5 We notice that if we used the 4PN and 5PN coefficients, a 5 and a 6 , suggested in Ref. [31], we would obtain poles in the function ∆t(R) for |a| > 0.75M . Moreover, if we adopted the values of a 5 and a 6 , suggested in Ref. [31] for the spin configurations analyzed in this paper, for which there are no poles in ∆t(R), we would obtain phase disagreements on the same order of the ones we have found.
The light ring is the unstable circular orbit of a massless particle (such as a photon) and can be computed from the deformed EOB metric or from H eff Kerr (R, P , S Kerr ). As in the nonspinning case, we do find a light ring with our spin EOB Hamiltonian. However, in contrast to the nonspinning case, for several values of the spin parameters (including the DD configuration) our spin EOB Hamiltonian does not yield a maximum orbital frequency. It is worth mentioning that if we were using only the "Kerr" part of the spin EOB Hamiltonian, i.e. we ignore H eff part (R, P , σ), then we do find a maximum of the orbital frequency and its value is quite close to the light ring position. A more detailed study has revealed that the absence of the maximum of the orbital frequency for the full spin EOB Hamiltonian is due to the spin-orbit coupling term H eff part (R, P , σ) defined in Eq. (13), which as discussed above does not reduce to the test-particle limit prediction at PN orders higher than 2.5PN. In the improved spin EOB Hamiltonian [32], H eff part (R, P , σ) will be consistent with the test-particle limit prediction at all PN orders linear in the particle spin. Analyses using the improved spin EOB Hamiltonian [32] have shown a reasonable agreement between the position of the EOB light ring and the maximum of the EOB orbital frequency.
Quite interestingly, when the numerical and EOB waveforms are aligned at low frequency, as discussed in detail in Sec. III A, we find that the EOB light ring is reached at time 0.3M before the peak of the numerical h 22 amplitude. Therefore, a nice property of the nonspinning EOB model [14,15] holds also in the spinning case, i.e. the EOB light ring position is a good approximation of the peak position of the numerical h 22 amplitude. The latter property will be a key ingredient in the EOB waveform model, as described later in Sec. II D.

C. EOB waveform: Inspiral & Plunge
Having described the inspiral dynamics, we now turn to the gravitational waveforms h ℓm . The latter can also be employed to compute consistently the inspiral dynamics through the radiation-reaction force [14]. In the nonspinning case, Refs. [14,15] have shown that the resummed, factorized PN waveforms proposed in Ref. [25] are in excellent agreement with the numerical waveforms. In Ref. [24] we have generalized the resummed factorized waveforms to include spin effects.
The resummed waveforms are written as the product of five factors, where ǫ denotes the parity of the multipolar waveform. In the circular-orbit case, ǫ is the parity of ℓ + m: ǫ = π(ℓ + m). These factors are discussed extensively in Ref. [25].
Here we simply write down the expressions used in our spin EOB model, valid for spins aligned or anti-aligned with the orbital angular momentum. Thus, we restrict ourselves to the equatorial plane (θ = π/2 and p θ = 0). The leading term h where R is the distance from the source. The n (ǫ) ℓm and c ℓ+ǫ (ν) are functions given in Eqs. (5)-(7) of Ref. [25]. The Y ℓm (θ, φ) are the scalar spherical harmonics. The tangential velocity v φ is the non-Keplerian velocity of a spherical orbit defined by v φ = r ΩΩ where and P φ,cir is the solution of the spherical orbit condition ∂H eff (R, P R = 0, P φ )/∂R = 0. As in the nonspinning case, the functionsŜ eff is an effective source term that in the test-particle, circularmotion limit contains a pole at the EOB light ring. It is given in terms of the EOB dynamics aŝ (33) is not the only possible choice; for example, one may instead chooseŜ The effect of this choice on the spin EOB model investigated in this paper is marginal, since in the equal-mass, equalspin, non-precessing binary configurations, odd parity modes contribute only a tiny fraction of the total energy flux (see Sec. III D for details). Although we choose to use the source term defined in Eq. (33), there is no evidence indicating that this choice is better or worse than others for those binary configurations in which odd parity modes are more important.
The function T ℓm in the right hand side of Eq. (30) resums leading logarithms of tail effects, and e iδ ℓm is a phase correction due to subleading logarithms. Through 2PN order, there are no tail contributions due to spin effects and T ℓm and e iδ ℓm do not differ from the nonspinning case. Their explicit expressions are given in Eqs. (19)-(29) of Ref. [25]. Finally, the functions f ℓm in the right hand side of Eq. (30) collect the remaining PN terms. We computed [24] the spin terms in f ℓm by Taylor expanding the h ℓm in Eq. (30) and comparing it to the Taylor-expanded h ℓm calculated in PN theory, including the test-particle spin effects through 4PN order. In the test-particle limit, we choose S Kerr as the spin variable of the spacetime. Expressions of f ℓm can be read from Ref. [24] 6 . Following Refs. [15,25], we resum all the nonspinning terms in f ℓm in the functional form f NS ℓm = (ρ ℓm ) ℓ that holds at known PN orders, where f NS ℓm collects the nonspinning terms in f ℓm , and ρ ℓm can be read from Appendix C of Ref. [25]. The motivation for this ρresummation is to reduce the magnitude of the 1PN coefficients in f ℓm that grow linearly with ℓ (see Sec. IID of Ref. [25]). Since at leading order we did not find such dependence on ℓ in the spinning terms [24], we do not apply the ρ-resummation to them.
Furthermore, since we are trying to reproduce effects in the numerical simulations that go beyond the quasicircular motion assumption, motivated by the PN expansion for generic orbits, we include non-quasicircular (NQC) effects in h ℓm in the form A similar expression was used in Ref. [15] except that there we usedṙ instead of p r * . For a test-particle plunging in the Kerr geometry,ṙ goes to zero at the horizon. We observe a similar behavior in the EOB ν-deformed Kerr geometry. Therefore, in contrast to the nonspinning case, the evolution ofṙ is not monotonic during the inspiral-plunge:ṙ increases during the inspiral, reaches a peak, and then starts decreasing during the plunge. By replacingṙ with p r * , we keep the NQC correction terms in Eq. (34) monotonic in time; thus, they can successfully model the monotonically increasing amplitude differences between the quasi-circular EOB and numerical waveforms. As in Ref. [15], we fix two of the four adjustable parameters a h22 i by requiring that the peaks of the numerical and EOB h 22 waveforms coincide in both time and amplitude, where the peak time of the numerical h 22 waveform is accurately predicted by the EOB light ring, as discussed above. The other two a h22 i parameters are determined by minimizing the overall amplitude difference with respect to the numerical waveform as explained in detail below. The NQC corrections in Eq. (34) also depend on spins. However, there is not enough numerical information in this work (we have only the DD configuration) to discriminate between the spinning and nonspinning contribution.

D. EOB waveform: Merger & Ringdown
The merger-ringdown waveform in the spin EOB model is built in the same way as in the nonspinning EOB model. Details on building merger-ringdown waveforms can be found in Sec. IIC of Ref. [15]. Here we briefly summarize the key points.
In the spin EOB model, the ring-down waveform is a linear combination of eight quasinormal modes (QNMs) of the final Kerr black hole. Mass and spin of the final black hole are computed from numerical data. In particular, for the numerical simulation of the DD configuration, we use M BH /M = 0.961109 ± 0.000003 and a BH /M BH = 0.54781 ± 0.00001 computed in Ref. [26]. Frequencies of the QNMs are computed by interpolating data from Ref. [33]. The eight complex coefficients of the linear combination are fixed by the hybrid comb matching described in Sec. IIC of Ref. [15]. The matching time t ℓm match is fixed to be the EOB light-ring position. The matching interval ∆t ℓm match is an adjustable parameter that is fixed by reducing the difference against numerical merger-ringdown waveforms.

III. CALIBRATING THE EOB WAVEFORMS TO NUMERICAL RELATIVITY SIMULATIONS
We now calibrate the spin EOB model against the numerical UU and DD spin configurations. We extract both the Newman-Penrose (NP) scalars Ψ ℓm 4 and the strain waveforms h ℓm from the simulations. The strain waveforms are extracted with the Regge-Wheeler-Zerilli (RWZ) formalism [34][35][36][37] (see Appendix of Ref. [15] for details of the numerical implementation used to obtain h ℓm ). We use the RWZ h ℓm to calibrate the EOB model, and use the NP Ψ ℓm 4 to check the consistency between the two wave-extraction schemes and to estimate the numerical error associated with the wave extraction schemes.
We will use the ℓ = 2, m = 2 component of the numerical waveform for tuning the EOB model. Thus, we calibrate in total the following six adjustable EOB parameters: b(ν), a 3PN SS , a h22 1 , a h22 2 , a h22 3 and ∆t 22 match .

A. Uncertainties in numerical waveforms
In this section, we compare numerical waveforms computed at different numerical resolutions and/or using different extrapolation procedures, or with different waveextraction schemes. Estimates of numerical errors in the waveforms will set our standards when calibrating the EOB model.
First, we adopt the same waveform-alignment procedure used in Ref. [15], that is we align waveforms at low frequency by minimizing the quantity over a time shift ∆t and a phase shift ∆φ, where φ 1 (t) and φ 2 (t) are the phases of the two waveforms. The range of integration (t 1 , t 2 ) is chosen to be as early as possible to maximize the length of the waveform but late enough to avoid the contamination from junk radiation present in the numerical initial data. The range of integration should also be large enough to average over numerical noise. We fix t 1 = 500M and t 2 = 1500M in Eq. (35). Using this alignment procedure, we estimate the errors on the numerical Ψ 22 4 . Figs. 1 and 2 summarize the phase errors for numerical Ψ 22 4 . The numerical waveform labeled "(N6, n=3)" [or "Ψ 4 (N6, n=3)"] is the reference numerical waveform used throughout this paper. Each waveform is extracted on a set of spheres at fixed distances from the source, and then extrapolated to future null infinity; the labels n refer to different orders of this extrapolation and are used to quantify the uncertainty in the phase due to extrapolation. The waveform labeled by N5 (as opposed to N6) is from a simulation with a lower numerical resolution and is used to quantify the uncertainty due to numerical truncation errors. The waveform labeled by "ḧ(N6, n=3)" is generated by twice differentiating the RWZ-extracted "h(N6, n=3)" waveform, and is used to quantify the uncertainty due to the systematic difference between extracting the NP scalar and extracting the strain waveform via the RWZ formalism.
The noise before t = 500M is due to spurious radiation from initial conditions. The features around t ≈ 2100M in Fig. 2 are due to a change of gauge in the numerical simulation. Extrapolation with n = 2 leads to systematic errors in the extrapolated waveform which in turn results in a systematic error in ∆t. Therefore, the green dashed lines in Figs. 1 and 2 represent a possibly overly conservative error estimate. There is a tiny frequency difference between the NP and RWZ extracted waveforms, which is magnified into a substantial time shift when the waveforms are aligned at low frequency. As a consequence, the dot-dashed brown line in Fig. 2 shows a larger phase difference which builds up during the late inspiral. It provides us with the most conservative error estimate for the DD configuration. This is better illustrated in Figs. 3 and 4, where we compare the RWZ h 22 and NP Ψ 22 4 waveforms without any time or phase shift. In blue solid lines, we show the phase and relative amplitude differences ∆φ NP and ∆A NP /A between the RWZ h 22 waveform differentiated twice with respect to time and Ψ 22 4 . In red dashed lines, we show the phase and relative amplitude differences ∆φ RWZ and ∆A RWZ /A between Ψ 22 4 integrated twice in time and the RWZ h 22 . In Fig. 4, ∆φ NP shows a slope between t = 500M and 1500M . When we apply the alignment procedure, this slope is removed through a time shift, which is transformed into a larger phase difference during late inspiral where the wave frequency is large.

B. Calibrating the EOB adjustable parameters
Here we adopt the procedure suggested in Ref. [15], augmented with the iterative scheme suggested in Ref. [14] when calibrating the adjustable parameters.
We divide the adjustable parameters into three groups and calibrate them in two steps. The first group, EOBdynamics parameters, consists of b(ν) and a 3PN SS in the EOB Hamiltonian (there is no adjustable parameter in the model of the EOB energy flux). These parameters determine the inspiral and plunge dynamics of the EOB model and affect the merger-ringdown waveform only indirectly through the waveform's phase and frequency around the matching point. The second group, EOB-NQC parameters, consists of a h ℓm i , which enter both the EOB dynamics (through the energy flux) and the EOB waveform (through the NQC correction). The third group, EOB-waveform parameters, consists of ∆t ℓm match , which affect the EOB merger-ringdown waveform but not the EOB inspiral-plunge waveform. All the EOB adjustable parameters are calibrated to the numerical RWZ h 22 . In the first step of calibration, we simultaneously reduce the difference in waveforms against the numerical UU and DD configurations by calibrating the EOB-dynamics and the EOB-NQC adjustable parameters. In the second step, using the adjustable parameters fixed in the first step, we calibrate the EOB-waveform adjustable parameters. We adopt the iterative scheme suggested in Ref. [14] to fix the EOB-dynamics and the EOB-NQC parameters in the first step of calibration.
In each iteration, we first minimize the phase difference by calibrating the EOB-dynamics adjustable parameters. To be specific, we choose to minimize the quantity max . By comparing the EOB model to both the UU and DD configurations it is possible to calibrate the parameters b(ν) and a 3PN SS separately. This is because b(ν) alters a 3.5PN spin-orbit coupling term that depends on the spin orientation, so the phases of the UU and DD waveforms change in opposite directions when varying b(ν), but a 3PN SS alters a 3PN spin-spin coupling term, so the phases of the UU and DD waveforms change in the same direction when varying a 3PN SS . The EOB-NQC adjustable parameters are calibrated only to the numerical h 22 waveform of the DD configuration, because we did not run the UU case through merger and ringdown. We first fix a h22 1 and a h22 2 by requiring that a local extremum of the EOB h 22 amplitude coincides with the peak of the numerical h 22 in time and amplitude (the peak time is predicted by the EOB light ring). We expect that in the future, the peak amplitude of the numerical h 22 will be predicted by numerical relativity with high accuracy as an interpolation function on the physical parameters. Therefore, a h22 1 and a h22 2 can be determined without a least-squares fit to the NR waveform, reducing by two the number of parameters to be determined by a least-squares fit. The other two NQC parameters, a h22 3 and a h22 4 , are calibrated to the numerical waveform to further reduce the disagreement in amplitude. The NQC parameters will enter the flux through the NQC waveform NQC h 22 in the next iteration. They are set to zeros initially to start the iteration and they usually converge within five iterations.
In the third step, we calibrate the EOB-waveform adjustable parameter ∆t 22 match by reducing the difference in the DD configuration merger-ringdown h 22 waveform.

C. Comparing the gravitational-wave modes h22
Before calibrating the EOB adjustable parameters, we investigate the phase difference for the EOB uncalibrated waveforms. For the uncalibrated model, we set a 5 = 1.775 and all six of our adjustable parameters to zero. We find that during the inspiral, the phase agreement between the numerical and spin EOB uncalibrated waveforms is already substantially better than the agreement between numerical and Taylor-expanded PN waveforms. For the latter, we consider the 3.5PN spin Taylor model (T4) of Ref. [38] with amplitude corrections through the highest PN order currently known [24,39]. In fact, using the uncalibrated spin EOB model and aligning the waveforms with the procedure defined by Eq. (35), we find that the phase differences against the numerical UU and DD waveforms, at the end of the simulation and at spin Taylor T4 model, the corresponding phase differences are 2.0 and −10.0 rads. Therefore, the spin EOB model, even uncalibrated, improves the phase agreement with numerical waveforms of Taylor-expanded PN models by resumming the PN dynamics.
When calibrating the EOB model, we find that the difference in phase and amplitude between the numerical and EOB waveforms is minimized when fixing the EOBdynamics parameters b(ν) = −1.65 and a 3PN SS = 1.5. In Fig. 5, we compare numerical and EOB h 22 waveforms for the UU configuration. The phase difference and relative amplitude difference are strictly within 0.01 rads and 1%, respectively. The systematic error in the EOB waveform in the UU configuration is therefore smaller than the numerical errors.
In Fig. 6, we compare numerical and EOB h 22 waveforms for the DD configuration. Using b(ν) = −1.65 and a 3PN SS = 1.5 again, we find that the best phase and amplitude agreement is obtained when the matching occurs at an interval of ∆t 22 match = 3.5M ended at t 22 match = 2402.0M , which is the EOB light-ring position and is 0.3M before the merger, i.e. the peak of the nu- = −4587.53. The relative amplitude difference is strictly within 1% until 2000M . After 2000M , although oscillations due to numerical gauge effects in the RWZ h 22 waveform are at the level of 2% until the merger, the average difference is still less than 1%. After the merger, the amplitude difference grows to about −5% and starts oscillating with increasing magnitude. The latter phenomenon is due to gauge effects in the RWZ h 22 waveform as discussed in Sec. III A and the Appendix of Ref. [15]. The phase difference is within 0.01 rads until about 1800M and grows to −0.28 rads until merger and settles to about 0.1 rads before the exponentially decaying amplitude results in increased errors in the extracted gravitational-wave phase.
In the spin DD configuration, the error in the EOB waveform (especially in the phase) is not within the numerical errors. The phase difference of ∼ 0.13 rads at late inspiral around t = 2300M can be reduced to within the numerical errors of ∼ 0.01 rads by calibrating the EOBdynamics adjustable parameters, i.e. b(ν) and a 3PN SS . However, this leads to an increase of the phase difference around the merger. Since we choose to minimize the span of the phase difference over the range of the full inspiral, merger and ringdown waveform, the phase difference at late inspiral is larger than what it could have been. The largest phase difference around merger can not be removed by calibrating the chosen adjustable parameters. Nevertheless, we can substantially reduce the phase difference if we allow one of the EOB-dynamics parameters b(ν) and a 3PN SS to be different in the UU and DD cases, or if we add one more spin-independent adjustable parameter. For instance, there can be a NQC correction factor on the right-hand side of Eq. (34) that contributes solely to the phase of the waveform in the form of [28] h insp−plunge We can reduce the phase difference at merger to < 0.05 rads by calibrating this extra EOB-NQC adjustable parameter b h ℓm 1 . However, until we study a larger number of waveforms, we will not over-tune parameters, since the main purpose of this preliminary and exploratory study on the spin EOB model is to show that by making a very simple and minimal choice of adjustable parameters, we can achieve a quite fair agreement with the numerical simulations.
We shall emphasize that, despite the small phase difference that exceeds the numerical errors in the DD configuration, the faithfulness of the EOB waveforms with the numerical waveforms is very good. Using the noise curves of Enhanced LIGO and Advanced LIGO 8 , for both the UU and DD configurations, we find that the faithfulness is always better than 0.999 for black-hole binaries with a total mass of 30-200M ⊙ . Note that the numerical waveforms start roughly at 40Hz for binaries with total mass 30M ⊙ and at 10Hz for binaries with total mass 100M ⊙ . Since the Advanced LIGO noise curve has a low frequency cutoff at 10Hz, the numerical waveforms are not long enough to cover the Advanced LIGO sensitivity band for binaries with total mass smaller than 100M ⊙ . When computing overlaps for these lower mass binaries using Advanced LIGO noise curve, we start the integration at the starting frequency of the numerical waveforms instead of 10 Hz. For the available numerical waveforms, the overlaps with EOB waveforms are well above the requirement on the accuracy of binary blackhole waveforms for detection and measurement purposes in gravitational-wave observations [40].
In Fig. 7, we compare the amplitude and frequency of numerical and EOB h 22 waveforms together with the orbital frequency of the EOB model, for the DD configuration. Unlike the nonspinning case [15], the orbital frequency Ω continues to grow during the plunge. However, the EOB light-ring is very close to the peak of the numerical h 22 , as discussed in Sec. II A. Note that during the ringdown, the frequency computed from the numerical h 22 shows increasingly large oscillations. We also plot the frequency computed from the numerical Ψ 22 4 model. This frequency shows much smaller, and bounded, oscillations deep into the ringdown regime. Advanced LIGO, we use the broadband configuration power spectral density given at http://www.ligo.caltech.edu/advLIGO/scripts/ref_des.shtml. D. Comparing the gravitational-wave modes hℓm Here we generate inspiral higher-order modes, h ℓm , using the same dynamics-adjustable parameters calibrated to the numerical h 22 mode in the previous section. The EOB-NQC parameters and the EOB-waveform parameters for these modes are not calibrated, since higherorder numerical waveforms show large numerical errors before reaching their peaks. For this reason, we constrain the comparison between numerical and EOB higher-order modes to the inspiral stage. The higher-order modes are aligned at low frequencies using the same time and phase shifts (modulo a factor of m/2 in the phase shifts) applied to the EOB h 22 mode.
In Fig. 8, we compare the EOB (dashed lines) and numerical (solid lines) amplitudes of the first five modes that dominate the signal power. In the DD case we show results only until t = 2000M because at later times the numerical data are affected by large oscillations, likely due to gauge effects. Except for the h 32 mode in the DD configuration, the agreement is very good for all the subdominant modes, as well as for the h 22 mode, in both the UU and DD waveforms. We believe that the difference seen for the DD h 32 mode is due to the lack of knowledge of PN spin couplings in the amplitude of the higher modes. In fact, only the leading-order PN spin term is known in the amplitude of the h 32 mode, and no PN spin terms are known in the amplitudes of the other subdominant modes shown in Fig. 8. Being resummed in the form of Eq. (30), the leading spin term in h 32 leaves a residual term in f 32 at the leading order. We test two choices of the odd-parity source term in Eq. (33),L eff andĤ eff , and find this residual term always dominating over all the nonspinning terms and causing f 32 to decrease and cross zero at high frequency, thus showing the odd behavior in the DD (3,2) mode of Fig. 8. We also try to apply the ρ-resummation discussed in Sec. II C on the spin terms of the f 32 . Although when applying the ρ-resummation, the leading order residual term in ρ 32 is reduced by 1/ℓ = 1/3 with respect to the residual term in f 32 , it still dominates over other terms and causes (ρ 32 ) 3 to cross zero at high frequency.
In Fig. 8, we also show the Taylor-expanded PN amplitudes (dotted lines). Their expressions can be read from Refs. [39] and [24], and they depend on dynamical variables only through the orbital velocity. We calculate these amplitudes using the non-Keplerian orbital velocity defined in Eq. (32) for the leading term and the Keplerian orbital velocity for all the next-to-leading terms. We calculate the non-Keplerian and Keplerian velocities using the EOB dynamics. That is to say, these amplitudes and the resummed amplitudes are calculated using exactly the same dynamical evolutions. In particular, the energy flux in the dynamics is always modeled by resummed waveforms, even when we calculate the Taylorexpanded PN amplitudes. These Taylor-expanded PN amplitudes are not to be confused with the amplitudes of the adiabatic PN approximants, such as the TaylorT1 and TaylorT4 approximants [41], because the underlying dynamics of the latter is completely different. In Fig. 8, although the Taylor-expanded PN amplitudes work reasonably well for the h 22 mode during inspiral, and probably by chance also for the h 32 mode in the UU configuration, their performance is not as good as that of the resummed amplitudes in general. Especially, for the h 44 and h 42 modes, the Taylor amplitudes are not monotonic. This unpleasant behavior is caused by their 1PN order non-spinning terms. Furthermore, the insert of Fig. 8 shows that the performance of Taylor-expanded PN amplitudes becomes worse for the h 22 mode during the late inspiral and plunge in the DD configuration.
Given the current information from PN theory and numerical simulations, we consider the agreement in Fig. 8 reasonable and do not dwell further on the choice of the waveform modeling options. The differences have little impact on the EOB model since the largest difference in h 32 affects the energy flux by less than 10 −4 , which is overwhelmed by other uncertainties in the EOB dynamics.
For the five dominant modes, the relative differences between the numerical and EOB h ℓm frequencies are within 0.5%, except for the (3, 2) mode in the DD configuration where the difference is within 1%. Since the h ℓm frequency depends on both the orbital frequency and its amplitude, the larger amplitude difference in the (3, 2) mode affects its gravitational-wave frequency. Except for the (3,2) mode in the DD configuration, all frequency agreement is within the numerical errors.

IV. CONCLUSIONS
In this paper, we carried out the first calibration of the spin EOB model to accurate numerical-relativity simulations of spinning, non-precessing black-hole binaries. We focused on two equal-mass black-hole binaries having spins both aligned, or both anti-aligned with the or-bital angular momentum, and dimensionless magnitude ∼ 0.44 [26].
For the EOB conservative dynamics, we adopted the spin EOB Hamiltonian suggested in Refs. [19,20], augmented with the 4PN order non-spinning parameter a 5 and two adjustable parameters. For the EOB nonconservative dynamics, we employed the gravitationalwave energy flux which includes spin effects and which has been computed using the factorized multipolar waveforms of Ref. [24].
As in previous cases [13,15], we aligned the EOB and numerical waveforms at low frequency over a time interval of 1000M , and minimized the difference between numerical and EOB waveforms by calibrating a handful of EOB-adjustable parameters. In particular, in this first exploration, we calibrated two EOB-dynamics adjustable parameters [b(ν) in Eq. (16b) which introduces a spin-orbit term at 3.5PN order, and a 3PN SS in Eq. (19) which introduces a 3PN spin-spin term], and three EOB-NQC adjustable parameters [see Eq. (34)] which enter the gravitational-wave energy flux and the EOB gravitational-wave (2,2) mode. Finally, we also calibrated the EOB-waveform adjustable parameter ∆t 22 match . Quite interestingly, similar to the case of non-spinning waveforms, we found that for spinning waveforms, once the EOB-dynamics adjustable parameters are calibrated at low frequency, the EOB light-ring coincides with the peak of the numerical-relativity waveform. Thus, for both spinning and non-spinning binary black holes, the EOB light-ring marks the most natural point at which to match the EOB inspiral-plunge waveform to the EOB merger-ringdown waveform.
In the equal-mass, spin aligned case, we found that phase and fractional amplitude differences between the numerical and EOB (2, 2) mode can be reduced to 0.01 rads and 1%, respectively, over the entire inspiral waveforms. In the equal-mass, spin anti-aligned case, these differences can be reduced to 0.13 rads and 1% during inspiral, and to 0.4 rads and 10% during merger and ringdown. The waveform agreement is within numerical errors in the spin aligned case while slightly over numerical errors in the spin anti-aligned case. Despite this difference, we found that using Enhanced LIGO and Advanced LIGO noise curves, the overlap maximized with respect to reference time and phase between the EOB and the numerical (2,2) mode, is larger than 0.999 for binaries with total mass 30-200M ⊙ . This is well above the accuracy requirement of binary black-hole waveforms for detection and measurement purposes in gravitationalwave observations [40].
In addition to comparing the numerical and EOB waveforms for the leading (2,2) mode, we also compared them for the next four subleading modes. Except for the h 32 mode in the DD configuration, the amplitude and frequency agreements are very good for all the subdominant modes, as well as for the h 22 mode, in both the UU and DD waveforms. We believe that the difference seen for the DD h 32 mode is due to the lack of knowledge of PN spin couplings in the amplitude of the subleading modes.
The spin EOB Hamiltonian [19,20] adopted in this paper was an excellent starting point to explore the calibration of the EOB model against spinning numerical simulations; however, as discussed above, and in particular in Sec. II B, the spin EOB Hamiltonian we used exhibits some unusual behaviour. Especially when extended at 4PN and 5PN, in some regions of the parameter space the Hamiltonian does not have an ISCO or the ISCO radius grows as the spin magnitude increases. This is opposite to the result in the test-particle limit case. Moreover, although the spin EOB Hamiltonian has a light ring, in some regions of the parameter space (including the anti-aligned case discussed in this paper) the orbital frequency does not reach a maximum. Those features turned out to be crucial when calibrating the EOB model to nonspinning numerical waveforms, and we believe will be crucial also when modeling spinning numerical waveforms. We found that the lack of those features in the current EOB Hamiltonian is due to the ad hoc spin coupling term H eff part (R, P , σ), defined in Eq. (13). This spin coupling term does not reproduce the results of a spinning test-particle at PN orders higher than 2.5PN. Analyses using an improved spin EOB Hamiltonian [32] obtained by building on Ref. [23] have shown that those features can be recovered.
variables R, P * , S 1 and S 2 . The transform from P to P * is a coordinate transform, not a canonical transform. In this section, we derive explicitly the transform to tortoise coordinate for the Hamiltonian and Hamilton equations of motion implemented in Cartesian coordinates.
The transform between P and P * is determined by the invariance in their tangential components and the rescaling in their radial components, that is R × P = R × P * , ξ a (R) R · P = R · P * . (A1) Choosing three independent equations out of the four above, we can write the transform in components as or explicitly as