# Effective-One-Body Numerical-Relativity waveform model for Eccentric spin-
precessing binary black hole coalescence

Xiaolin Liu  Institute for Frontiers in Astronomy and Astrophysics, Beijing
Normal University, Beijing 102206, China  Department of Astronomy, Beijing
Normal University, Beijing 100875, China  Zhoujian Cao 111corresponding author
<EMAIL_ADDRESS>Institute for Frontiers in Astronomy and Astrophysics, Beijing
Normal University, Beijing 102206, China  Department of Astronomy, Beijing
Normal University, Beijing 100875, China  School of Fundamental Physics and
Mathematical Sciences, Hangzhou Institute for Advanced Study, UCAS, Hangzhou
310024, China  Zong-Hong Zhu  Institute for Frontiers in Astronomy and
Astrophysics, Beijing Normal University, Beijing 102206, China  Department of
Astronomy, Beijing Normal University, Beijing 100875, China

###### Abstract

Waveform models are important to gravitational wave data analysis. People
recently pay much attention to the waveform model construction for eccentric
binary black hole coalescence. Several Effective-One-Body Numerical-Relativity
waveform models of eccentric binary black hole coalescence have been
constructed. But none of them can treat orbit eccentricity and spin-precessing
simultaneously. The current paper focuses on this problem. The authors
previously have constructed waveform model for spin-aligned eccentric binary
black hole coalescence SEOBNRE. Here we extend such waveform model to describe
eccentric spin-precessing binary black hole coalescence. We calculate the 2PN
orbital radiation-reaction forces and the instantaneous part of the decomposed
waveform for a general spinning precessing binary black hole system in
effective-one-body (EOB) coordinates. We implement these results based on our
previous SEOBNRE waveform model. We have also compared our model waveforms to
both SXS and RIT numerical relativity waveforms. We find good consistency
between our model and numerical relativity. Based on our new waveform model,
we analyze the impact of the non-perpendicular spin contributions on waveform
accuracy. We find that the non-perpendicular spin contributions primarily
affect the phase of the gravitational waveforms. For the current gravitational
wave detectors, this contribution is not significant. The future detectors may
be affected by such non-perpendicular spin contributions. More importantly our
SEOBNRE waveform model, as the first theoretical waveform model to describe
eccentric spin-precessing binary black hole coalescence, can help people to
analyze orbit eccentricity and spin precession simultaneously for
gravitational wave detection data.

##  I Introduction

After several observation runs, the LIGO and Virgo Aasi et al. (2015);
Acernese et al. (2014) ground-based gravitational wave detectors have detected
a large number of binary black hole (BBH) Abbott et al. (2016); The LIGO
Scientific Collaboration et al. (2021), binary neutron stars (BNS) Abbott et
al. (2017) and binary black hole-neutron star (BH-NS) Abbott et al. (2021)
events. Along with the increasing detector sensitivity, we expect that in
future observations, LIGO-Virgo-Kagra Akutsu et al. (2019), as well as the
under-construction ET Punturo et al. (2010), CE Reitze et al. (2019), and
space-based detectors such as LISA Amaro-Seoane et al. (2017), Taiji Hu and Wu
(2017), Tianqin Luo et al. (2016), will be able to detect more dynamic
features of binary black hole systems. This will provide us valuable physics
information about extremely strong gravitational field and provide powerful
support for the study of gravitational theory, galaxy formation theory, and
cosmology, as well as the verification of general relativity.

Merger systems of compact binary object systems are the most common
gravitational wave sources in the universe. Their dynamics and related
gravitational waves have been studied for many years Einstein et al. (1938);
Lorentz and Droste (1937); Ohta et al. (1973); Damour and Sch?fer (1985);
Blanchet (2006). The research on this yopic is still ongoing Jaranowski and
Sch?fer (2012); Damour et al. (2016); Bernard et al. (2018). The detection and
data analysis of gravitational waves depend on the construction of
gravitational waveform templates Finn (1992); Jaranowski and Krolak (1994);
Jaranowski et al. (1998). Many years ago, the designed probing targets of LIGO
are gravitational wave events from compact binary star systems in the merger
phase. Near the merger stage, people expect that the eccentricity of the
binary system has been reduced to almost zero. Therefore, many commonly used
templates now are based on circular orbit assumption for two-body systems. In
recent years, people have begun to pay attention to gravitational waves from
two-body systems along eccentric orbit, and have found some traces of
eccentricity residuals in LIGO-Virgo data Abbott et al. (2020). This situation
calls for more accurate and comprehensive gravitational wave templates for
general two-body systems.

So far, many eccentric models without spin precession and spin-precession
circular models have been established. Among them are templates based on the
effective one-body model, such as the TEOBResumS Chiaramello and Nagar (2020);
Nagar et al. (2021) model, the SEOBNRE Cao and Han (2017); Liu et al. (2020);
Yun et al. (2021); Liu et al. (2022, 2023) model and SEOBNRv4E Ramos-Buades et
al. (2022); Khalil et al. (2021) model. These models assume that the spin
direction of the binary black hole is parallel to the orbital angular
momentum, in which case the direction of the orbital angular momentum remains
unchanged. If the spin direction of one black hole is not parallel to the
direction of the orbital angular momentum, the orbital direction of the binary
black hole system will change due to the interaction between the spin of the
black hole and the orbital angular momentum, which is called spin precession
Apostolatos et al. (1994). Currently, many models have considered the spin-
precession binary black hole system on circular orbits, such as the IMRPhenomP
Pratten et al. (2021) model and the SEOBNRv4P Pan et al. (2014); Ossokine et
al. (2020) model. Recently, the next generation spin precession framework
SEOBNRv5 with higher accuracy and computational efficiency has been published
Ramos-Buades et al. (2023); Khalil et al. (2023).

Besides the afore mentioned complete inspiral-merger-ringdown waveform models,
there are other phenomenological and post-Newtonian (PN) approximated waveform
models for eccentric BBH. In Huerta et al. (2018); Habib and Huerta (2019);
Chen et al. (2021) Eccentric, Nonspinning, Inspiral-Gaussian-process Merger
Approximant $(\texttt{ENIGMA})$ model was proposed, which combines analytical
and NR results using machine learning algorithms. The authors of Hinder et al.
(2018) combine PN inspiral waveform for eccentric BBH and circular NR waveform
for merger and ringdown. The surrogate eccentric waveform model based on
numerical relativity simulations was proposed in Islam et al. (2021) for
spinless black holes. The authors of Setyawati and Ohme (2021); Wang et al.
(2023) constructed an eccentric waveform model for spinless black holes
through fitting the difference between the circular waveform and eccentric
waveform. Pioneered by Peters Peters (1964), post-Newtonian eccentric waveform
models Mishra et al. (2015); Tanay et al. (2016); Moore et al. (2016); Boetzel
et al. (2019a); Ebersold et al. (2019); Tiwari and Gopakumar (2020); Loutrel
and Yunes (2017a, b); Moore et al. (2018); Loutrel (2020, 2021) include post-
circular (PC) waveform model Yunes et al. (2009), enhanced post-circular (EPC)
waveform model Huerta et al. (2014), x-model Hinder et al. (2010) and others.

In order to aid the gravitational wave data analysis, we need complete and
accurate waveform models for generic binary systems. Especially waveform
models of eccentric spin-precessing binary black hole coalescence is needed
Romero-Shaw et al. (2020, 2023); Gamba et al. (2023); Gayathri et al. (2022).
In the current paper we extend SEOBNRE to construct such a waveform model. Our
waveform model is based on effective one body (EOB) numerical relativity
framework. The framework has been introduced in Cao and Han (2017). As the two
major blocks of the framework, we describe the extended EOB dynamics in the
next section and the extended waveform formula in Sec. III. After that we
check the performance of our waveform model by comparing the model waveform to
numerical relativity waveform in Sec. IV. Finally we conclude the paper with
some discussions in the last section.

Throughout the paper we use geometric unit system with the speed of light and
gravitational constant taking one ($c=G=1$). Some equations show $c$
explicitly for the clarity of order expansion.

##  II Effective-one-body dynamics for eccentric spin-pressing BBH

As usual we use Hamiltonian to describe the EOB dynamics. The EOB Hamiltonian
of a BBH system with mass $m_{1,2}$ reads Buonanno and Damour (1999)

$\displaystyle
H_{\text{EOB}}=M\sqrt{1+2\nu\left(\hat{H}_{\text{eff}}-1\right)},$ (1)
$\displaystyle\hat{H}_{\text{eff}}\equiv\frac{H_{\text{eff}}}{M\nu},$ (2)

where the total mass $M=m_{1}+m_{2}$ and symmetric mass ratio
$\nu=m_{1}m_{2}/M^{2}$. The effective Hamiltonian $H_{\text{eff}}$ depends on
BH’s spin $\boldsymbol{S}_{1,2}:=m_{1,2}^{2}\boldsymbol{\chi}_{1,2}$, where
$\boldsymbol{\chi}$ is dimensionless spin vector satisfying
$-1<|\boldsymbol{\chi}|<1$. We borrow the function form of $H_{\text{eff}}$
from Barausse and Buonanno (2010, 2011). Then the dynamics of eccentric spin-
pressing BBHs can be described by the following canonical equations

$\displaystyle\derivative{\boldsymbol{r}}{t}=\partialderivative{H_{\text{EOB}}}{\boldsymbol{p}},$
(3)
$\displaystyle\derivative{\boldsymbol{p}}{t}=-\partialderivative{H_{\text{EOB}}}{\boldsymbol{r}}-\boldsymbol{F},$
(4)
$\displaystyle\derivative{\boldsymbol{S}_{1,2}}{t}=\partialderivative{H_{\text{EOB}}}{\boldsymbol{S}_{1,2}}\times\boldsymbol{S}_{1,2}+\boldsymbol{F}^{s}_{1,2},$
(5)

where $\boldsymbol{r}$ is relative position vector, and
$\boldsymbol{F},\boldsymbol{F}^{s}_{1,2}$ are radiation reaction forces for
momentum $\boldsymbol{p}$ and BH spins respectively.

The spin radiation reaction forces $\boldsymbol{F}^{s}_{1,2}$ can be written
as Maia et al. (2017a, b)

$\displaystyle\boldsymbol{F}_{1}^{s}=\frac{4p_{r}}{15r^{4}}(\boldsymbol{L}\times\boldsymbol{S}_{1})\left[\nu\left(-22\frac{1}{r}+36p^{2}-60p_{r}^{2}\right)+\frac{1}{2}(1-2\nu-\delta)\left(16\frac{1}{r}-48p^{2}+75p_{r}^{2}\right)\right]$
$\displaystyle+\frac{2\nu}{15r^{6}}\left[6r^{2}p_{r}(\boldsymbol{S}_{1}\times\boldsymbol{S}_{2})\left(11\frac{1}{r}-18p^{2}+30p_{r}^{2}\right)-3r(\boldsymbol{S}_{1}\times\boldsymbol{p})\left(\boldsymbol{r}\cdot\boldsymbol{S}_{2}-\frac{1-\delta-2\nu}{2\nu}\boldsymbol{r}\cdot\boldsymbol{S}_{1}\right)\left(8\frac{1}{r}-9p^{2}+45p_{r}^{2}\right)\right.$
$\displaystyle\left.-(\boldsymbol{r}\times\boldsymbol{S}_{1})\left[15p_{r}\left(\boldsymbol{r}\cdot\boldsymbol{S}_{2}-\frac{1-\delta-2\nu}{2\nu}\boldsymbol{r}\cdot\boldsymbol{S}_{1}\right)\left(2\frac{1}{r}-3p^{2}+7p_{r}^{2}\right)+2r\left(\boldsymbol{p}\cdot\boldsymbol{S}_{2}-\frac{1-\delta-2\nu}{2\nu}\boldsymbol{p}\cdot\boldsymbol{S}_{1}\right)\left(\frac{8}{r}-9p^{2}+45p_{r}^{2}\right)\right]\right]$
$\displaystyle-\frac{2(1-\delta)}{5r^{5}}\left[4\left(\frac{1}{r}-3p^{2}+15p_{r}^{2}\right)\left[(\boldsymbol{r}\times\boldsymbol{S}_{1})\boldsymbol{r}\cdot\boldsymbol{S}_{1}+(\boldsymbol{p}\times\boldsymbol{S}_{1})\boldsymbol{r}\cdot\boldsymbol{S}_{1}\right]+15\frac{1}{r}p_{r}(\boldsymbol{r}\times\boldsymbol{S}_{1})(\boldsymbol{r}\cdot\boldsymbol{S}_{1})(3p^{2}-7p_{r}^{2})\right.$
$\displaystyle\left.-30rp_{r}(\boldsymbol{p}\times\boldsymbol{S}_{1})(\boldsymbol{p}\cdot\boldsymbol{S}_{1})\right],$
(6)

where $\boldsymbol{L}=\boldsymbol{r}\times\boldsymbol{p}$ is the orbital
angular momentum, $p_{r}=\boldsymbol{r}\cdot\boldsymbol{p}/r$, and
$\delta=(m_{1}-m_{2})/M,\boldsymbol{F}_{2}^{s}=\boldsymbol{F}_{1}^{s}(1\leftrightarrow
2)$. In Maia et al. (2017a, b), Harmonic coordinate is used. Since the
Harmonic coordinates equals to the EOB coordinates to leading PN order, we
apply the above formula of the spin radiation reaction forces
$\boldsymbol{F}^{s}_{1,2}$ to EOB coordinate directly. In the quasi-circular
limit $p_{r}=0$, the above equations indicate that
$\boldsymbol{F}^{s}_{\text{qc},1,2}\approx 0$ consequently Khalil et al.
(2023).

For general eccentric spin-precessing cases, the radiation reaction force
satisfies a balance equation with an extra Schott terms
$\dot{E}_{\text{Schott}}$ and $\dot{\boldsymbol{J}}_{\text{Schott}}$ Bini and
Damour (2012); Khalil et al. (2021)

$\displaystyle\dot{E}_{\text{sys}}+\dot{E}_{\text{Schott}}=$
$\displaystyle\mathcal{F},$
$\displaystyle\dot{\boldsymbol{J}}_{\text{sys}}+\dot{\boldsymbol{J}}_{\text{Schott}}=$
$\displaystyle\boldsymbol{\mathcal{J}},$ (7)

where $E_{\text{sys}},\boldsymbol{J}_{\text{sys}}$ are the energy and total
angular momentum of the binary system respectively,
$\mathcal{F}=\differential{E_{\text{gw}}}/\differential{t}$ and
$\boldsymbol{\mathcal{J}}$ respectively represent the energy flux and the
angular momentum flux of gravitational waves at infinity. The post-Newtonian
approximation of the energy flux and the angular momentum flux of
gravitational waves has been calculated in Kidder et al. (1993); Buonanno et
al. (2013). For our usage we transform the Harmonic coordinate in these
calculation results to EOB coordinate. We follow the detail transformation
given in Khalil et al. (2021).

The radiation reaction forces act on the binary system and cause the energy
and angular momentum loss Wang and Will (2007)

$\displaystyle\dot{E}_{\text{sys}}=\dot{\boldsymbol{r}}\cdot\boldsymbol{F},$
$\displaystyle\dot{\boldsymbol{J}}_{\text{sys}}=\boldsymbol{r}\times\boldsymbol{F}+\boldsymbol{F}^{s}_{1}+\boldsymbol{F}^{s}_{2}.$
(8)

Substituting Eq. (8) into Eq. (7), we can solve the radiation reaction force
and get

$\displaystyle\boldsymbol{F}=\frac{1}{\dot{\boldsymbol{r}}\cdot\boldsymbol{r}}\left[\boldsymbol{r}(\mathcal{F}-\dot{E}_{\text{Sch}})-\dot{\boldsymbol{r}}\times\left[\boldsymbol{\mathcal{J}}-\dot{\boldsymbol{J}}_{\text{Sch}}-(\boldsymbol{F}^{s}_{1}+\boldsymbol{F}^{s}_{2})\right]\right].$
(9)

Regarding to the Schott term $E_{\text{Sch}}$, we borrow the instant term
shown in the first equation of Eq. (31) of Khalil et al. (2021). And more we
borrow the SO term shown in the first equation of Eq. (28) of Khalil et al.
(2021) and replace $\chi_{1,2}$ with
$\boldsymbol{L}\cdot\boldsymbol{S}_{1,2}$. For the SS term we adopt form

$\displaystyle\qquad\qquad
E^{\text{SS}}_{\text{Sch}}=\frac{\eta^{2}p_{r}}{r^{4}c^{4}}\left[\boldsymbol{S}_{1}^{2}(a_{1}p_{r}^{2}+a_{2}p^{2}+\frac{a_{3}}{r})+(\boldsymbol{S}_{1}\cdot\boldsymbol{S}_{2})(a_{4}p_{r}^{2}+a_{5}p^{2}+\frac{a_{6}}{r})+\boldsymbol{S}_{2}^{2}(a_{7}p_{r}^{2}+a_{8}p^{2}+\frac{a_{9}}{r})\right]+$
$\displaystyle\frac{1}{r^{7}c^{4}}\left[(\boldsymbol{p}\cdot\boldsymbol{S}_{1})(\boldsymbol{r}\cdot\boldsymbol{S}_{1})\left(\frac{h_{1}}{r}+h_{2}p^{2}+h_{3}p_{r}^{2}\right)+(\boldsymbol{r}\cdot\boldsymbol{S}_{1})(\boldsymbol{p}\cdot\boldsymbol{S}_{2})\left(\frac{h_{4}}{r}+h_{5}p^{2}+h_{6}p_{r}^{2}\right)+1\leftrightarrow
2\right]+$
$\displaystyle\frac{p_{r}}{r^{7}c^{4}}\left[d_{1}(\boldsymbol{r}\cdot\boldsymbol{S}_{1})^{2}+d_{2}(\boldsymbol{r}\cdot\boldsymbol{S}_{1})(\boldsymbol{r}\cdot\boldsymbol{S}_{2})+1\leftrightarrow
2\right],$ (10)

Regarding to the Schott term $\boldsymbol{J}_{\text{Sch}}$, we borrow the
instant term for $\boldsymbol{J}_{\text{Sch}}$ shown in the second equation of
Eq. (31) of Khalil et al. (2021) and replace the magnitude of the angular
momentum $L$. For the SO term, we need to introduce components besides the one
along the angular momentum direction $\boldsymbol{L}$

$\displaystyle\boldsymbol{J}^{\text{SO}}_{\text{Sch}}=\frac{\eta^{2}p_{r}}{r^{2}c^{3}}\left[\frac{\boldsymbol{L}}{r}\left(\boldsymbol{L}\cdot\boldsymbol{S}_{1}(\beta_{1}p_{r}^{2}+\beta_{2}p^{2}+\frac{\beta_{3}}{r})\right)+\right.$
$\displaystyle\left.\boldsymbol{p}\left(\frac{\beta_{4}}{r}\boldsymbol{r}\cdot\boldsymbol{S}_{1}+\beta_{5}\boldsymbol{p}\cdot\boldsymbol{S}_{1}\right)+\boldsymbol{S}_{1}(\beta_{6}p_{r}^{2}+\beta_{7}p^{2}+\frac{\beta_{8}}{r})\right]+$
$\displaystyle\frac{\beta_{9}\eta^{2}\boldsymbol{r}}{r^{4}c^{3}}\boldsymbol{p}\cdot\boldsymbol{S}_{1}+1\leftrightarrow
2.$ (11)

For the SS term we adopt form

$\displaystyle\qquad\qquad\boldsymbol{J}=\frac{\eta^{2}p_{r}\boldsymbol{L}}{r^{4}c^{4}}\left[a_{10}\boldsymbol{S}_{1}^{2}+a_{11}\boldsymbol{S}_{1}\cdot\boldsymbol{S}_{2}+a_{12}\boldsymbol{S}_{2}^{2}+\frac{1}{r^{2}}\left[f_{1}(\boldsymbol{r}\cdot\boldsymbol{S}_{1})(\boldsymbol{r}\cdot\boldsymbol{S}_{1})+f_{2}(\boldsymbol{r}\cdot\boldsymbol{S}_{1})(\boldsymbol{r}\cdot\boldsymbol{S}_{2})+f_{3}(\boldsymbol{r}\cdot\boldsymbol{S}_{2})(\boldsymbol{r}\cdot\boldsymbol{S}_{2})\right]+\right.$
$\displaystyle\left.r\left[f_{4}(\boldsymbol{p}\cdot\boldsymbol{S}_{1})(\boldsymbol{p}\cdot\boldsymbol{S}_{1})+f_{5}(\boldsymbol{p}\cdot\boldsymbol{S}_{1})(\boldsymbol{p}\cdot\boldsymbol{S}_{2})+f_{6}(\boldsymbol{p}\cdot\boldsymbol{S}_{2})(\boldsymbol{p}\cdot\boldsymbol{S}_{2})\right]\right]+\frac{\eta^{2}p_{r}}{r^{4}c^{4}}\left[\boldsymbol{p}\times\boldsymbol{S}_{1}\left(b_{1}(\boldsymbol{r}\cdot\boldsymbol{S}_{1})+b_{2}(\boldsymbol{r}\cdot\boldsymbol{S}_{2})\right)+\frac{\boldsymbol{r}\times\boldsymbol{S}_{1}}{r}\right.$
$\displaystyle\left.\left(\left(\frac{b_{3}}{r}+b_{4}p_{r}^{2}+b_{5}p^{2}\right)(\boldsymbol{p}\cdot\boldsymbol{S}_{1})+\left(\frac{b_{6}}{r}+b_{7}p_{r}^{2}+b_{8}p^{2}\right)(\boldsymbol{p}\cdot\boldsymbol{S}_{2})\right)+1\leftrightarrow
2\right]+\frac{\eta^{2}}{r^{2}c^{4}}\left[\boldsymbol{p}\times\boldsymbol{S}_{1}\left((\boldsymbol{p}\cdot\boldsymbol{S}_{1})\left(c_{1}p^{2}+c_{2}p_{r}^{2}+\frac{c_{3}}{r}\right)+\right.\right.$
$\displaystyle\left.\left.(\boldsymbol{p}\cdot\boldsymbol{S}_{2})\left(c_{4}p^{2}+c_{5}p_{r}^{2}+\frac{c_{6}}{r}\right)\right)+\frac{\boldsymbol{r}\times\boldsymbol{S}_{1}}{r^{3}}\left((\boldsymbol{r}\cdot\boldsymbol{S}_{1})\left(d_{1}p^{2}+d_{2}p_{r}^{2}+\frac{d_{3}}{r}\right)+(\boldsymbol{r}\cdot\boldsymbol{S}_{2})\left(d_{4}p^{2}+d_{5}p_{r}^{2}+\frac{d_{6}}{r}\right)\right)+1\leftrightarrow
2\right]+$
$\displaystyle\frac{\eta^{2}\boldsymbol{L}}{r^{6}}\left[(\boldsymbol{r}\cdot\boldsymbol{S}_{1})(\boldsymbol{p}\cdot\boldsymbol{S}_{1})\left(\frac{e_{1}}{r}+e_{2}p^{2}+e_{3}p_{r}^{2}\right)+(\boldsymbol{r}\cdot\boldsymbol{S}_{1})(\boldsymbol{p}\cdot\boldsymbol{S}_{2})\left(\frac{e_{4}}{r}+e_{5}p^{2}+e_{6}p_{r}^{2}\right)+1\leftrightarrow
2\right].$ (12)

In the above equations about Schott terms we have introduced a set of to-be-
determined coefficients. In order to fix them, we tried condition

$\displaystyle\frac{\boldsymbol{F}\cdot\boldsymbol{r}L}{|\boldsymbol{r}\times\boldsymbol{F}|rp_{r}}=$
$\displaystyle 1+\order{p_{r}^{2}}$ (13)

used in Pan et al. (2014). But we find this condition is not enough to
determine the set of coefficients mentioned above.

In fact the condition (13) can be looked as a nearly circular approximation.
Keeping this fact in mind and we note that the radiation reaction force in
circular limit should take the form Pan et al. (2014)

$\displaystyle\boldsymbol{F}|_{p_{r}=0,\dot{p}_{r}=0}=\frac{\mathcal{F}_{\text{qc}}}{\Omega_{\text{qc}}|\boldsymbol{L}|}\boldsymbol{p},$
(14)

where $\Omega_{\text{qc}}$ means the orbital angular velocity. The above
equation implies

$\displaystyle\frac{\Omega_{\text{qc}}L}{\mathcal{F}_{\text{qc}}}\boldsymbol{F}=\boldsymbol{p}\left[\order{1}+\order{p_{r}}+\order{\dot{p}_{r}}\right]+$
$\displaystyle\boldsymbol{r}\left[\order{p_{r}}+\order{\dot{p}_{r}}\right]+\boldsymbol{L}\left[\order{p_{r}}+\order{\dot{p}_{r}}\right]+$
$\displaystyle\boldsymbol{S}_{1}\left[\order{p_{r}}+\order{\dot{p}_{r}}\right]+\boldsymbol{S}_{2}\left[\order{p_{r}}+\order{\dot{p}_{r}}\right].$
(15)

We plug the Schott terms into (9) and expand the vector $\boldsymbol{F}$
according to the form shown in the above equation. The condition (15) tells us
that in the $\boldsymbol{p}$ direction, the lowest order in the expansions
using $p_{r}$ and $\dot{p}_{r}$ is $\order{1}$, while in other directions,
such as $\boldsymbol{r}$ and $\boldsymbol{L}$, the lowest order is at least
$\order{p_{r}}$ or $\order{\dot{p}_{r}}$. These facts give us a set of
equations. We solve these equations and fix the to-be-determined coefficients
introduced in the Schott terms. After that combining the energy flux, angular
momentum flux, Schott terms and $\boldsymbol{F}^{s}_{1,2}$ we archive the form
of the radiation reaction force

$\displaystyle\boldsymbol{F}=\left(f^{r}_{\text{ns}}+\frac{1}{c^{3}}f^{r}_{\text{so}}+\frac{1}{c^{4}}f^{r}_{\text{ss}}\right)\boldsymbol{r}+\left(f^{p}_{\text{ns}}+\frac{1}{c^{3}}f^{p}_{\text{so}}+\frac{1}{c^{4}}f^{p}_{\text{ss}}\right)\boldsymbol{p}$
$\displaystyle+\frac{1}{c^{3}}f^{L}_{\text{so}}\boldsymbol{L}+\frac{1}{c^{4}}f^{1}_{\text{ss}}\boldsymbol{S}_{1}+\frac{1}{c^{4}}f^{2}_{\text{ss}}\boldsymbol{S}_{2},$
(16)

The coefficients $f$s are listed in the Appendix A.

Next, we re-sum this radiation reaction force into a factorized form,

$\displaystyle\boldsymbol{F}=\frac{\mathcal{F}}{\Omega
L}(\mathcal{F}^{e}_{p}\boldsymbol{p}+\mathcal{F}^{e}_{r}\boldsymbol{r}+\mathcal{F}^{e}_{L}\boldsymbol{L}+\mathcal{F}^{e}_{1}\boldsymbol{S}_{1}+\mathcal{F}^{e}_{2}\boldsymbol{S}_{2}),$
(17)
$\displaystyle\Omega\equiv|\boldsymbol{r}\times\dot{\boldsymbol{r}}|/r^{2},$
(18) $\displaystyle L\equiv|\boldsymbol{r}\times\boldsymbol{p}|.$ (19)

The factors $\mathcal{F}^{e}_{...}$ are listed in the supplementary materiel.
Eq. (17) is used in the updated SEOBNRE code to calculate the radiation
reaction force. Numerically we calculate $\mathcal{F}_{\text{qc}}$ with the
waveform got in the code according Cao and Han (2017)

$\displaystyle\mathcal{F}=\frac{1}{16\pi}\sum_{\ell}\sum_{m=-\ell}^{\ell}|\dot{h}_{\ell
m}|^{2},$ (20)

where $h_{\ell m}$ are spherical harmonic modes of the gravitational waves.

Regarding to numerical implementation, we need note one fact. Near merger,
$\dot{p}_{r}$ becomes very large which strongly affects the numerical
accuracy. To treat this problem we follow the trick used by usual EOB code Pan
et al. (2010) to replace the original canonical momentum $p_{r}$ with the
momentum of the tortoise coordinates $p_{r*}$. The detail replacement
relations are

$\displaystyle
p_{r}=p_{r*}\left(1+\frac{2}{r^{2}c^{2}}+\frac{4-3\nu}{r^{2}c^{4}}\right),$
(21)
$\displaystyle\dot{p}_{r}=\dot{p}_{r*}-\frac{2(p_{r*}^{2}-r\dot{p}_{r*})}{r^{2}c^{2}}+\frac{1}{r^{3}c^{4}}\left[r(r\dot{p}_{r*}(4-3\nu)\right.$
$\displaystyle\left.+rp^{4}_{r*}(1+\nu)+p^{2}_{r*}(5(\nu-1)+r^{2}\dot{p}_{r*}(1+\nu))\right].$
(22)

Besides the dynamical equations we have described above, another important
issue about dynamics is the initial conditions. Adiabatically the spin-
precession orbits can be looked as evolved Keplerian orbits. So we need only
to determine the initial Keplerian orbit to quantify the initial condition.
The initial Keplerian orbit can be described by mean orbital angular velocity
$f_{0}$ and eccentricity $e_{0}$. Similar to what we have done in Cao and Han
(2017), we solve

$\displaystyle\partialderivative{H_{\text{EOB}}}{p_{\phi}}$
$\displaystyle=\frac{\pi}{f_{0}},$ (23)
$\displaystyle\partialderivative{H_{\text{EOB}}}{r}$
$\displaystyle=-\frac{e_{0}}{r^{2}},$ (24)

for the initial conditions. Here we call $f_{0}$ the initial orbital angular
velocity and $e_{0}$ the initial eccentricity. The remaining steps follow
those proposed in Buonanno et al. (2006). There are several different
definitions of eccentricity for binary system in general relativity Loutrel et
al. (2019); Knee et al. (2022). Many definitions are gauge dependent. In the
current work, we do not pay attention to the concept subtlety of eccentricity.
In stead we just take it as a parameter to describe the waveform.

##  III Factorized waveform

Following Pan et al. (2014) we calculate gravitational waveform for spin-
precessing BBHs in the precession frame. Similar to Khalil et al. (2023), we
define three mutually orthogonal normal vectors
$(\boldsymbol{n},\boldsymbol{\lambda},\boldsymbol{e})$ via

$\displaystyle\boldsymbol{n}:=\frac{\boldsymbol{r}}{r},$ (25)
$\displaystyle\boldsymbol{e}:=\frac{\boldsymbol{L}_{N}}{L_{N}},$ (26)
$\displaystyle\boldsymbol{\lambda}:=\boldsymbol{e}\times\boldsymbol{n},$ (27)

where $\boldsymbol{L}_{N}:=\boldsymbol{r}\times\boldsymbol{v}$, with
$\boldsymbol{v}:=\dot{\boldsymbol{r}}$ is velocity. These vectors constitute a
coordinate system in the precession frame.

Up to 2PN order, the gravitational wave waveform for eccentric spin-precessing
BBHs in harmonic coordinates has been calculated in Will and Wiseman (1996)
for orbital motion part and in Buonanno et al. (2013) for spin-orbit part. We
then transform the harmonic coordinate to EOB coordinate according Khalil et
al. (2021). Then we apply a rotation to the waveform to get the waveform in
comoving frame. The rotation is the one to let $\boldsymbol{r}$ and
$\boldsymbol{p}$ lie in x-y plane. After that we decompose the waveform in
comoving frame into spin-weighted -2 spherical harmonic modes. Finally we
express the results in terms of $r$, $p_{r}$, $L$ and

$\displaystyle\chi_{(n,\lambda,e)}\equiv\boldsymbol{\chi}\cdot(\boldsymbol{n},\boldsymbol{\lambda},\boldsymbol{e}).$
(28)

The detail expressions are presented in the supplementary material. Our
waveform can recover the results of Khalil et al. (2021) in the spin aligned
binary cases.

Next we re-sum the waveform into the factorized form Pan et al. (2011); Damour
et al. (2009); Ossokine et al. (2020); Khalil et al. (2021)

$\displaystyle h^{F}_{\ell m}=h^{N}_{\ell
m}\hat{S}_{\text{eff}}T^{\text{qc}}_{\ell m}e^{i\delta_{\ell m}}f_{\ell m},$
(29)

where $h^{N}_{\ell m}$ is the Newtonian order of the waveform. We adopt the
form used in Ossokine et al. (2020) for $T^{\text{qc}}_{\ell
m}e^{i\delta_{\ell m}}$ to count the tail contributions. The effective source
term $\hat{S}_{\text{eff}}$ is given by

$\displaystyle\hat{S}_{\text{eff}}=\begin{cases}\hat{H}_{\text{eff}}&\text{$\ell+m$
is even}\\\ v_{\Omega}L&\text{$\ell+m$ is odd}\end{cases}.$ (30)

In Barausse and Buonanno (2010, 2011), $\hat{H}_{\text{eff}}$ has been given
in terms of $r$, $p_{r}$, $L$, $\chi_{1,2n}$, $\chi_{1,2\lambda}$ and
$\chi_{1,2e}$. In the current work we adopt this function form but replace $L$
with $\beta\equiv\sqrt{1+r^{2}\dot{p}_{r}}$. In addition we expand
$v_{\Omega}\equiv\Omega^{1/3}$ and $L$ in terms of $r$, $p_{r}$, and $\beta$
through relations (Eq. (20) of Khalil et al. (2021))

$\displaystyle\frac{L^{2}}{r}=\beta^{2}+\frac{1}{2rc^{2}}\left[\beta^{4}(\nu+1)-3\beta^{2}(\nu-1)+2(\nu+1)+r\left(\beta^{2}(\nu+1)-\nu+3\right)p_{r}^{2}\right]+\frac{3\beta\left(\chi_{se}(\nu-2)-2\chi_{ae}\sqrt{1-4\nu}\right)}{c^{3}r^{3/2}}$
(31)
$\displaystyle+\frac{1}{8r^{2}c^{4}}\left[\beta^{6}\nu^{2}+5\beta^{6}\nu+\beta^{6}-3\beta^{4}\nu^{2}-5\beta^{4}\nu+13\beta^{4}+2\beta^{2}\nu^{2}-12\beta^{2}\nu+34\beta^{2}+48\delta\chi_{a\lambda}\chi_{s\lambda}\nu-64\beta^{2}\chi_{ae}^{2}\nu\right.$
$\displaystyle\left.+16\beta^{2}\chi_{ae}^{2}+2rp_{r}^{2}\left(3\beta^{2}+\nu\left(3\beta^{4}+\beta^{2}(\nu-5)+16\chi_{a\lambda}^{2}-16\chi_{s\lambda}^{2}(\nu-1)-\nu\right)+4\chi_{ae}^{2}(4\nu-1)-4\chi_{se}^{2}(1-2\nu)^{2}\right.\right.$
$\displaystyle\left.\left.-4\chi_{a\lambda}^{2}-4\chi_{s\lambda}^{2}-22\nu+5\right)-64\beta^{2}\delta\chi_{ae}\chi_{se}\nu+32\beta^{2}\delta\chi_{ae}\chi_{se}+48\delta\chi_{ae}\chi_{se}\nu+32\delta\chi_{ae}\chi_{se}\nu
rp_{r}^{2}-16\delta\chi_{ae}\chi_{se}rp_{r}^{2}\right.$
$\displaystyle\left.+64\beta^{2}\chi_{se}^{2}\nu^{2}-64\beta^{2}\chi_{se}^{2}\nu+16\beta^{2}\chi_{se}^{2}-48\chi_{se}^{2}\nu^{2}+48\chi_{se}^{2}\nu-48\chi_{s\lambda}^{2}\nu^{2}+48\chi_{s\lambda}^{2}\nu-12\nu-8\left(2\beta^{2}-3\right)(4\nu-1)\chi_{an}^{2}\right.$
$\displaystyle\left.-16\delta\chi_{an}\chi_{sn}\left(\beta^{2}(4\nu-2)+3\right)+8\chi_{sn}^{2}\left(2\beta^{2}(1-2\nu)^{2}-3\right)-\left(\beta^{2}-1\right)((\nu-1)\nu+1)r^{2}p_{r}^{4}+32\delta\chi_{a\lambda}\chi_{s\lambda}\nu
rp_{r}^{2}\right.$
$\displaystyle\left.-16\delta\chi_{a\lambda}\chi_{s\lambda}rp_{r}^{2}+24\right],$
(32)
$\displaystyle\sqrt{r}v_{\Omega}=\beta^{1/3}+\frac{1}{12r\beta^{5/3}c^{2}}\left[-\left(\beta^{4}(\nu+1)\right)+\beta^{2}(\nu-1)+2(\nu+1)-r\left(\beta^{2}(\nu+1)+\nu-3\right)p_{r}^{2}\right]+\frac{\chi_{se}(\nu-2)-2\delta\chi_{ae}}{6\beta^{2/3}r^{3/2}c^{3}}$
(33)
$\displaystyle+\frac{1}{288\beta^{11/3}r^{2}}\left[\beta^{8}((\nu-16)\nu+1)-2\beta^{6}(\nu(4\nu-9)+14)-7\beta^{4}(\nu+3)(3\nu+1)+2rp_{r}^{2}\left(\beta^{6}(\nu(7\nu-4)+7)\right.\right.$
$\displaystyle\left.\left.+\beta^{4}\left(-6\nu^{2}+4\nu-8\right)+\beta^{2}\left(96\delta\chi_{a\lambda}\chi_{s\lambda}\nu-48\delta\chi_{a\lambda}\chi_{s\lambda}+24\chi_{ae}^{2}(4\nu-1)+48\delta\chi_{ae}\chi_{se}(2\nu-1)-24\chi_{se}^{2}(1-2\nu)^{2}\right.\right.\right.$
$\displaystyle\left.\left.\left.+24\chi_{a\lambda}^{2}(4\nu-1)-96\chi_{s\lambda}^{2}\nu^{2}+96\chi_{s\lambda}^{2}\nu-24\chi_{s\lambda}^{2}-39\nu^{2}-84\nu-41\right)+10(\nu-3)(\nu+1)\right)\right.$
$\displaystyle\left.+4\beta^{2}\left(\nu\left(72\left(\delta\chi_{a\lambda}\chi_{s\lambda}+\delta\chi_{ae}\chi_{se}-\nu\left(\chi_{se}^{2}+\chi_{s\lambda}^{2}\right)+\chi_{se}^{2}+\chi_{s\lambda}^{2}\right)+19\nu-18\right)+36(4\nu-1)\chi_{an}^{2}-72\delta\chi_{an}\chi_{sn}\right.\right.$
$\displaystyle\left.\left.-36\chi_{sn}^{2}+17\right)-20(\nu+1)^{2}+r^{2}\left(\beta^{4}(\nu(13\nu+8)+13)+\beta^{2}\left(20\nu^{2}-34\nu-36\right)-5(\nu-3)^{2}\right)p_{r}^{4}\right].$
(34)

So based on the re-summation relation we can determine $f_{\ell m}$. And more
we divide the got $f_{\ell m}$ into two parts. One is the quasi-circular part
corresponding to $p_{r}=\dot{p}_{r}=0$. The another is the elliptical
correction part $f^{e}_{\ell m}$ corresponding to the left part other than
$p_{r}=\dot{p}_{r}=0$. For odd $m$ cases, the quasi-circular part is divided
into spinless part $\rho^{\ell}_{\ell m}$ and spin part $f^{s}_{\ell m}$ more
Khalil et al. (2021). Consequently we have

$\displaystyle\qquad f_{\ell m}=\begin{cases}\rho^{\ell}_{\ell
m}(r,\beta,\chi_{1,2n},\chi_{1,2\lambda},\chi_{1,2e})+f^{e}_{\ell
m}(r,p_{r},\beta,\chi_{1,2n},\chi_{1,2\lambda},\chi_{1,2e})&\text{$m$ is
even}\\\ \rho^{\ell}_{\ell m}(r,\beta)+f^{s}_{\ell
m}(r,\beta,\chi_{1,2n},\chi_{1,2\lambda},\chi_{1,2e})+f^{e}_{\ell
m}(r,p_{r},\beta,\chi_{1,2n},\chi_{1,2\lambda},\chi_{1,2e})&\text{$m$ is
odd}\end{cases}.$ (35)

We rearrange the auto-variables for the quasi-circular part to form
$v_{\Omega}$. So we can express $f_{\ell m}$ as

$\displaystyle\qquad f_{\ell m}=\begin{cases}\rho^{\ell}_{\ell
m}(v_{\Omega},\chi_{1,2n},\chi_{1,2\lambda},\chi_{1,2e})+f^{e}_{\ell
m}(r,p_{r},\beta,\chi_{1,2n},\chi_{1,2\lambda},\chi_{1,2e})&\text{$m$ is
even}\\\ \rho^{\ell}_{\ell m}(v_{\Omega})+f^{e}_{\ell
m}(r,p_{r},\beta,\chi_{1,2n},\chi_{1,2\lambda},\chi_{1,2e})+f^{s}_{\ell
m}(v_{\Omega},\chi_{1,2n},\chi_{1,2\lambda},\chi_{1,2e})&\text{$m$ is
odd}\end{cases}.$ (36)

---
Figure 1: Waveform comparison for $(2,2)$ mode between numerical relativity
and EOB models including newly got SEOBNRE model in the current paper and
SEOBNRv4PHM. The top row is the real part of the $(2,2)$ mode. The bottom row
is the imaginary part of the $(2,2)$ mode. SXS:BBH:1106 (mass ratio $q=1.68$,
BH’s spin $\vec{\chi}_{1}=(0.72,-0.00,0.36)$ and
$\vec{\chi}_{2}=(-0.14,-0.07,0.05)$) Caltech-Cornell-CITA  case is used in
this plot.
---
Figure 2: Waveform comparison for the real parts of higher modes $(2,1)$,
$(3,3)$, $(4,4)$ and $(5,5)$ between numerical relativity and EOB models
including newly got SEOBNRE model in the current paper and SEOBNRv4PHM.
Similar to Fig. 1, SXS:BBH:1106 (mass ratio $q=1.68$, BH’s spin
$\vec{\chi}_{1}=(0.72,-0.00,0.36)$ and $\vec{\chi}_{2}=(-0.14,-0.07,0.05)$)
Caltech-Cornell-CITA  case is used in this plot.

Comparing to the waveform used in Khalil et al. (2021), we introduce the extra
contribution of the BH’s spin component perpendicular to the orbital angular
momentum in the circular orbital part $\rho_{\ell m}$ for even-parity modes
and spin contribution $f^{s}_{\ell m}$ for odd-parity modes. We list the
detail expressions for these extra terms in Appendix B. The elliptical
corrections $f^{e}_{\ell m}$ are more complicated. We present them in the
supplementary material.

Note that the waveform (29) is respect to the precession frame. In order to
get the waveform in inertial frame we firstly calculate the state variables
according to the dynamical equations described in the last section. From the
calculated state variables we can get $r$, $\beta$, $\chi_{1,2n}$,
$\chi_{1,2\lambda}$ and $\chi_{1,2e}$. Consequently we can calculate the
waveform in the precession frame using (29). After that we apply an inverse
rotation operation on the waveforms $h_{\ell m}$ Pan et al. (2014). The
rotation operation is the one lets $\boldsymbol{r}$ and $\boldsymbol{p}$ lie
in x-y plane.

##  IV Performance of the spin-precession SEOBNRE waveform model

Along our series work about SEOBNRE waveform models Cao and Han (2017); Liu et
al. (2020, 2022, 2023), the current work is the first one describing spin-
precession. Before considering the spin-precession effect, we have checked the
performance of our new SEOBNRE waveform model for the spin-aligned cases. We
find our new model can achieve as good performance as the last version SEOBNRE
waveform model Liu et al. (2022). Compared to the last version SEOBNRE
waveform model Liu et al. (2022), the (5,5) mode is newly constructed. In
Appendix C, we show the test of the spin-precession SEOBNRE waveform model
against the spin-aligned BBHs. Regarding to spin-precession, we would like to
firstly check the performance of our waveform model for circular spin-
precessing cases. After that we check generic eccentric spin-precessing cases.
Then we quantify the performance of our waveform model with matching factor
which has been widely used to estimate the accuracy of a given waveform model.
Specifically SXS numerical relativity waveform catalog and RIT numerical
relativity waveform catalog are used to do the waveform comparison and
matching factor calculation.

###  IV.1 Circular spin-precessing cases

---
Figure 3: Corresponding to Fig. 1 and 2, the orbital plan precession behavior
(bottom panel) and the waveform correction factor $f_{22}$ (top panel:
amplitude of $f_{22}$, middle panel: phase of $f_{22}$) are plotted.
Comparison between SEOBNRv4PHM waveform model and SEOBNRE waveform model is
shown.

SEOBNRv4PHM waveform model is valid for circular spin-precession cases. In
order to check the performance of our new waveform model, we can compare our
resulted waveform to both the SEOBNRv4PHM waveform and the numerical
relativity waveform. As an example, we randomly choose a circular spin-
precession waveform from SXS catalog Caltech-Cornell-CITA , SXS:BBH:1106 which
has mass ratio $q=1.68$ and spin $\vec{\chi}_{1}=(0.72,-0.00,0.36)$
$\vec{\chi}_{2}=(-0.14,-0.07,0.05)$. For the EOB waveforms we set the
parameters of mass ratio and BH’s spin accordingly. For comparison, we align
the time of the waveforms with setting the corresponding time of maximal
$|h_{22}|$ as 0. After that we adjust the initial phase of EOB waveforms to
match the NR waveform. Then we plot the results in Fig. 1. From this
comparison we can see that both SEOBNRv4PHM waveform model and our newly
developed SEOBNRE waveform model can recover the NR results quite well. And
more we can see that the real part and imaginary part show similar behavior.
Such similar behavior happens for all cases we have checked and for all
waveform modes. So in the following we just plot real part. In Fig. 2 we
compare the higher modes for $(2,1)$, $(3,3)$, $(4,4)$ and $(5,5)$. The
consistency between NR waveforms and EOB models can be seen there.

---
Figure 4: Mismatch factor of $(2,2)$ spherical mode between NR waveform and
EOB model waveforms. The left panel is for SEOBNRv4PHM model. The right panel
is for SEOBNRE model. The Advanced LIGO designed sensitivity is used in this
figure.

Phenomenologically spin-precession makes the binary orbit plane precess. If we
initially set the orbit plane coincide with x-y plane, such precession
behavior shows variation of EOB spherical coordinate $\theta$. Using
SXS:BBH:1106 as an example, we investigate such plane precession behavior
respect to time. The result is plotted in bottom panel of Fig. 3. The
oscillation behavior of $\theta$ is clear. Along with this $\theta$ we also
plot the model waveform correction factor $f_{22}$ described in previous
section. For comparison, both factors for SEOBNRv4PHM waveform model and
SEOBNRE waveform model. Interestingly, we find that the factor of SEOBNRE
waveform model can closely catch up the oscillation behavior of $\theta$. In
contrast, the factor of SEOBNRv4PHM waveform model does not count such
oscillation at all. We understand this fact as following. SEOBNRv4PHM waveform
model neglects the contribution of the non-perpendicular spin by assuming that
such contribution is small. Differently SEOBNRE waveform model faithfully
counts the contribution of the non-perpendicular spin. It is this non-
perpendicular spin that introduces the oscillation behavior of $f_{22}$. In
fact other correction factors including $f_{\ell m}$ all admit this
oscillation behavior in the SEOBNRE waveform model.

In order to quantitatively check the consistency between the NR waveform and
EOB model waveforms, we calculate the mismatch factors. Due to initial
oscillations in NR simulations, when computing the best-matching EOB waveform,
we search for a waveform that matches well in the vicinity of the initial
frequency given by NR, serving as the best-matching waveform. We define the
mismatch factor $\mathcal{M}$ between two waveforms, $h_{1}$ and $h_{2}$ via

$\displaystyle\mathcal{M}=1-\max_{t_{c},\phi_{c}}{\frac{\innerproduct{h_{1}}{h_{2}}}{\sqrt{\innerproduct{h_{1}}{h_{1}}\innerproduct{h_{2}}{h_{2}}}}},$
(37)

where the inner product $\innerproduct{\cdot}{\cdot}$ is given by

$\displaystyle\innerproduct{h_{1}}{h_{2}}=\int_{f_{\text{min}}}^{f_{\text{max}}}\frac{\tilde{h}_{1}(f)\tilde{h}_{2}^{*}(f)}{S_{n}(f)}\differential{f}.$
(38)

We use the power spectrum density (PSD) of the high power designed Advanced
LIGO D. Shoemaker (2010) (LIGO Scientific Collaboration) as the one sided
power spectrum of detectors noise $S_{n}$. We compute the mismatch factor
$\mathcal{M}$ between the (2,2) mode of the numerical relativity waveforms and
the EOB waveforms. This mismatch factor $\mathcal{M}$ depends on the total
mass $M$ of the binary black hole system. The chosen range for the total mass
is from 20 to 200 solar masses, which roughly corresponds to the detection
range of ground-based detectors. The best-matching waveform is the one that
achieves the highest overall matching in the selected mass range.

In Fig. 4, we use $\mathcal{M}_{0}$ to denote the mismatch factor between the
NR waveform and the SEOBNRv4PHM model waveform and use $\mathcal{M}_{1}$ to
denote the mismatch factor between the NR waveform and the SEOBNRE model
waveform. The plot shows the mismatch factor as a function of the BBH total
mass. In all 63 NR waveforms are used here. We can see that SEOBNRE model
behaves as good as SEOBNRv4PHM model in most cases. In some cases, SEOBNRE
model can even show a little bit better behavior. This small improvement,
about $10^{-2}$, is due to the non-perpendicular spin contribution afore
mentioned.

###  IV.2 Eccentric spin-precessing cases

In the previous subsection, we have investigated the behavior of SEOBNRE
waveform model in the circular orbit spin-precessing cases. In this subsection
we move on to check the behavior of SEOBNRE waveform model for general cases
with both spin precession and orbit eccentricity. To the best of our
knowledge, our SEOBNRE waveform model is the first theoretical waveform model
to treat such generic BBHs.

---
Figure 5: Waveform comparison for $(2,2)$ mode between numerical relativity
and SEOBNRE model. The top row is the real part of the $(2,2)$ mode. The
bottom row is the imaginary part of the $(2,2)$ mode. RIT:eBBH:1632 (mass
ratio $q=1$, BH’s spin $\vec{\chi}_{1}=(0.7,0,0)$ and
$\vec{\chi}_{2}=(0.7,0,0)$) Healy and Lousto (2022) case is used in this plot.
---
Figure 6: Waveform comparison for the real parts of higher modes $(2,1)$,
$(3,3)$, $(4,4)$ and $(5,5)$ between numerical relativity and SEOBNRE model.
Similar to Fig. 5, RIT:eBBH:1632 (mass ratio $q=1$, BH’s spin
$\vec{\chi}_{1}=(0.7,0,0)$ and $\vec{\chi}_{2}=(0.7,0,0)$) Healy and Lousto
(2022) case is used in this plot.

Since SXS catalog has not simulation results for general BBH with both spin
precession and orbit eccentricity, we adopt RIT catalog instead Healy and
Lousto (2022). As an example, we chose RIT:eBBH:1632 to compare the waveform
between the NR and the SEOBNRE model. The mass ratio is 1, the initial spin
configuration is $\vec{\chi}_{1}=(0.7,0,0),\vec{\chi}_{2}=(0.7,0,0)$. The
initial eccentricity is $0.2775$ at $Mf_{22}=0.002463$. Numerical relativistic
simulations exhibit oscillations in the initial stage. After drop initial part
NR waveform, we simultaneously search for the initial frequency $f_{0}$ and
initial eccentricity $e_{0}$ through comparing the NR waveform to the SEOBNRE
model waveforms.

---
Figure 7: Similar to Fig. 3 and corresponding to Figs. 5 and 6, the SEOBNRE
model waveform correction factor $f_{22}$ (top panel: amplitude of $f_{22}$,
middle panel: phase of $f_{22}$) are plotted. Along with that, the orbit plane
precession behavior of $\theta$ is plotted in the third panel. The fourth
panel is the model waveform frequency of (2,2) mode $\omega_{22}$
corresponding to the co-precession frame. The fifth panel if the calculated
eccentricity according to $\omega_{22}$ shown in the fourth panel.

We extend the method proposed in Arif Shaikh et al. (2023) to estimate the
eccentricity $e_{\rm gw}(t)$ based on the gravitational waveform. Since there
are spin-precession effects in our case, we calculate the eccentricity in the
co-precession frame. The eccentricity of gravitational waves is defined by

$\displaystyle e_{\rm gw}=\cos{(\Psi/3)}-\sqrt{3}\sin{(\Psi/3)},$ (39)

where

$\displaystyle\Psi:=\arctan{\left(\frac{1-e^{2}_{\omega_{22}}}{2e_{\omega_{22}}}\right)}.$
(40)

The $e_{\omega_{22}}$ represent the generalized eccentricity that defined by
the frequency of the dominant $(2,2)$ mode $\omega_{22}$,

$\displaystyle e_{\omega_{22}}$
$\displaystyle=\frac{\sqrt{\omega^{p}_{22}}-\sqrt{\omega^{a}_{22}}}{\sqrt{\omega^{p}_{22}}+\sqrt{\omega^{a}_{22}}},$
(41)

where the superscripts $p,a$, denotes the frequency that evaluated at
pericenters and apocenters.

In Fig. 5, we present the comparison results between the best-matching SEOBNRE
model waveform and the NR waveform. The consistency between these two
waveforms is clear. This figure is similar to Fig. 1 but for general cases
with both spin precession and orbit eccentricity. Similar to Fig. 2 we plot
the higher modes in Fig. 6. Similar to Fig. 3, we check the model waveform
correction factor $f_{22}$ and correspondingly the orbit plane precession and
the evolution behavior of eccentricity in Fig. 7. Besides the behavior of
$f_{22}$, we have also plotted the orbit plane precession behavior of
$\theta$, the model waveform frequency of (2,2) mode $\omega_{22}$
corresponding to the co-precession frame and the calculated eccentricity
according to the got $\omega_{22}$. In addition to the oscillation behavior of
$f_{22}$ resulted from the orbit plan precession, we can also see another
oscillation behavior which happens when the orbit passes by the pericenters.
Interestingly there is a local maximal of $|f_{22}|$ happening exactly when
the orbit locates at the pericenters. Not like the spin-aligned cases, the
eccentricity may be effected by the orbit precession. During the evolution
process, when the orbit precession becomes strongly enough, the orbit
eccentricity shows oscillation behavior as shown in the bottom panel of Fig.
7, near merger.

---
Figure 8: Mismatch factor of $(2,2)$ spherical mode between the NR waveform
and the SEOBNRE model waveform. Similar to Fig. 5, RIT:eBBH:1632 (mass ratio
$q=1$, BH’s spin $\vec{\chi}_{1}=(0.7,0,0)$ and $\vec{\chi}_{2}=(0.7,0,0)$)
Healy and Lousto (2022) case is investigated in this plot.

In order to quantitatively estimate the accuracy of the SEOBNRE model
waveform, we again calculate the mismatch factor between the model waveform
and the NR waveform. We plot the result in Fig. 8. Using RIT:eBBH:1632 (mass
ratio $q=1$, BH’s spin $\vec{\chi}_{1}=(0.7,0,0)$ and
$\vec{\chi}_{2}=(0.7,0,0)$) Healy and Lousto (2022) as an example of generic
eccentric spin-precession BBH, we find that the matching factor achieved
within the mass range of 20 to 200 solar masses better than 96.2%.

Comparing Fig. 6 to Fig. 2 and 12, we can find that the higher modes for
eccentric spin-precession case perform worse than individual eccentric case
and spin-precession circular case. This means our model has not treated the
interaction between eccentricity and spin-precession well. But the good news
is the matching factor of dominating (2,2) mode between NR waveform and
SEOBNRE model waveform is better than 96% as we have seen in Fig. 8.

---
Figure 9: Waveforms comparison for different spin amplitudes. Corresponding to
RIT:eBBH:1632, mass ratio $q=1$ and initial eccentricity $e_{0}=0.37$ at
$Mf_{22}=0.003$ are used in this figure. Different lines correspond to
different BHs’ spin amplitude $\chi_{x}$. The BHs’ spins are set as
$\vec{\chi}_{1}=(\chi_{x},0,0)$ and $\vec{\chi}_{2}=(\chi_{x},0,0)$.
---
Figure 10: Waveforms comparison for different initial orbit eccentricity
$e_{0}$. Corresponding to RIT:eBBH:1632, mass ratio $q=1$ and BH’s spin
$\vec{\chi}_{1}=(0.7,0,0)$ and $\vec{\chi}_{2}=(0.7,0,0)$ are used in this
figure. Different lines correspond to different initial orbit eccentricity at
$Mf_{22}=0.003$.

Now we can use SEOBNRE waveform model as tool to investigate the effect of
spin precession and orbit eccentricity individually on the waveform. We use
RIT:eBBH:1632 as a reference but change BHs’ spin or change orbit eccentricity
to check the waveform behavior.

In Fig. 9 we investigate the effect of BH spin amplitude on the waveform. The
effect of spin on waveform comes in at relative high post-Newtonian order.
Consequently the waveforms corresponding to different spin amplitude are
similar to each other. In the current spin-precession cases, larger spin
results in larger precession. We can consequently see precession correction to
the waveform amplitude at about time $t=4000M$ and $t=8000M$. Interestingly we
find that the correction to the waveform phase is small.

In Fig. 10 we investigate the effect of orbit eccentricity on the waveform. As
ones expected, larger eccentricity results in faster inspiral and merger. In
addition, larger eccentricity results in stronger waveform burst behavior.
Consequently the waveform phase for larger eccentricity shows faster
increasing behavior. From Fig. 10 we can see different initial eccentricity
gives clearly different behavior.

|
---|---
Figure 11: Mismatch factors between the spin amplitude affected and orbit
eccentricity affected waveforms and the referenced SEOBNRE waveform modeling
RIT:eBBH:1632. The left panel is for the spin amplitude affected waveforms
corresponding to Fig. 9. The right panel is for the orbit eccentricity
affected waveforms corresponding to Fig. 10.

For LIGO-like ground-based detectors, the detected waveform is short. Only
merger part is available. So it is some hard to distinguish spin precession
and orbit eccentricity. Although we do not intend to investigate this issue in
detail in the current work, we would like to mimic such detection through
aligning waveforms at merger time and calculating corresponding mismatch
factor. After time alignment, if the two waveforms admit different length, we
just take the common time duration for the mismatch factor calculation. We
plot the resulted mismatch factors between the spin amplitude affected
(corresponding to Fig. 9) and orbit eccentricity (corresponding to Fig. 10)
affected waveforms and the referenced SEOBNRE waveform modeling RIT:eBBH:1632
in Fig. 11. Just as we analyzed above, the orbit eccentricity results in much
larger mismatch factor. So we expect our SEOBNRE waveform model can help
people to distinguish spin precession and orbit eccentricity in real data
analysis. Off course, much more studies on the parameters degeneracy are
needed.

##  V conclusion

In this work, we apply the method proposed in Khalil et al. (2021) to
calculate the 2PN orbital dynamics of a general spinning precessing binary
black hole system along eccentric orbit in EOB coordinates. We have also
calculated the instantaneous part of the factorized waveform. Specifically, we
have quantitatively analyzed the impact of non-perpendicular spin
contributions on the waveforms.

We extend our previous constructed SEOBNRE waveform model by incorporating the
instantaneous part of general spin contributions and radiation reaction force
correction. By comparing SEOBNRE waveform model and SEOBNRv4PHM waveform model
with NR waveforms with the PSD of ground-based detectors, we found that the
inclusion of non-perpendicular spin contributions led to a small improvement
in waveform accuracy. This improvement mainly stems from the modification of
the gravitational wave phase due to the non-perpendicular spin contributions,
which first appears at the 1PN order but is not prominent. For ground-based
detectors, since stellar-mass binary black hole systems are close to the
merger phase when they enter the detection frequency band, the phase
correction caused by non-perpendicular spin contributions is not significant,
and the aligned spin waveform exhibits high accuracy at this stage. For future
next-generation ground-based detectors and space-based detectors with lower
detection frequency bands, the influence of this contribution may become
important in waveform modeling.

As the first theoretical waveform model for generic spin-precession BBH with
orbit eccentricity, we have also investigated the accuracy of the SEOBNRE
waveform model through comparing the model waveform to NR waveform. Both SXS
Caltech-Cornell-CITA  and RIT Healy and Lousto (2022) NR waveform catalogs are
used in the current work. We have compared the new model with a NR waveform
featuring significant spin precession and a relatively large initial
eccentricity. We find that the new model approximately matches the numerical
relativistic waveform, indicating that the corrections computed in this study
can effectively describe this most general type of spin precession orbit with
eccentricity. Good consistency between our model waveform and NR waveform is
shown. In the near future our SEOBNRE waveform model can help people as a tool
to distinguish spin precession and orbit eccentricity in real data analysis.

In the following work we plan to improve our SEOBNRE waveform model in several
respects. Future efforts should involve constructing more comprehensive
models, such as the 2PN complete EOB decomposed waveform presented in this
paper and incorporating radiation-reaction forces in the dynamical system,
followed by hyper-parameter tuning. Furthermore, calculating the contributions
of the hereditary part of waveforms for general spinning precessing orbits is
another topic worth discussing. Currently, the hereditary part of waveforms
for generic equatorial orbits can be computed using quasi-Keplerian methods
and methods based on small eccentricity expansions Boetzel et al. (2019b).
However, for orbits with spin precession, constructing periodicity-based
descriptions to estimate the contributions of their hereditary part requires
further investigation.

## Acknowledgments

We would like to thank Alessandra Buonanno, Lijing Shao and Antoni Ramos
Buades for many helpful discussions. This work was supported in part by the
National Key Research and Development Program of China Grant No.
2021YFC2203001, the Strategic Priority Research Program of the Chinese Academy
of Sciences (No. XDB23000000) and in part by the NSFC (No. 11920101003 and No.
12021003). Z. Cao and Z.-H. Zhu were supported by “the Interdiscipline
Research Funds of Beijing Normal University”.

##  Appendix A The radiation reaction force

The radiation reaction force for general orbits up to 2PN order we get in this
paper reads

$\displaystyle\boldsymbol{F}=\left(f^{r}_{\text{ns}}+\frac{1}{c^{3}}f^{r}_{\text{so}}+\frac{1}{c^{4}}f^{r}_{\text{ss}}\right)\boldsymbol{r}+\left(f^{p}_{\text{ns}}+\frac{1}{c^{3}}f^{p}_{\text{so}}+\frac{1}{c^{4}}f^{p}_{\text{ss}}\right)\boldsymbol{p}+\frac{1}{c^{3}}f^{L}_{\text{so}}\boldsymbol{L}+\frac{1}{c^{4}}f^{1}_{\text{ss}}\boldsymbol{S}_{1}+\frac{1}{c^{4}}f^{2}_{\text{ss}}\boldsymbol{S}_{2},$
(42)

where the no-spin coefficients $f^{r}_{\text{ns}},f^{p}_{\text{ns}}$ are

$\displaystyle
f^{r}_{\text{ns}}=-\frac{8\nu^{2}p_{r}^{3}}{r^{4}}\left\\{1+\frac{1}{c^{2}}\left[-\frac{3}{14}(7\nu-5)p^{2}-\frac{3508\nu+1683}{420r}\right]+\frac{1}{c^{4}}\left[\left(\frac{185\nu^{2}}{504}-\frac{2171\nu}{504}+\frac{269}{504}\right)p^{4}\right.\right.$
$\displaystyle\left.\left.+\left(\frac{19\nu^{2}}{9}+\frac{155\nu}{36}-\frac{13}{6}\right)p^{2}p_{r}^{2}+\left(\frac{2617\nu^{2}}{180}-\frac{5557\nu}{72}-\frac{77837}{7560}\right)\frac{p^{2}}{r}+\left(-\frac{202\nu^{2}}{45}+\frac{3194\nu}{315}+\frac{1151}{630}\right)\frac{p_{r}^{2}}{r}\right.\right.$
$\displaystyle\left.\left.+\left(\frac{281\nu^{2}}{54}-\frac{102167\nu}{3780}+\frac{359}{1260}\right)\frac{1}{r^{2}}\right]\right\\},$
(43) $\displaystyle
f^{p}_{\text{ns}}=-\frac{8\nu^{2}}{15r^{3}}\left\\{10p^{2}-39p_{r}^{2}-\frac{22}{r}+\frac{1}{c^{2}}\left[\left(-\frac{3\nu}{7}-\frac{279}{28}\right)p^{4}+\left(\frac{9\nu}{7}+\frac{837}{28}\right)p^{2}p_{r}^{2}+\left(-\frac{121\nu}{14}-\frac{3833}{56}\right)\frac{p^{2}}{r}\right.\right.$
$\displaystyle\left.\left.+\left(\frac{45\nu}{2}-\frac{225}{14}\right)p_{r}^{4}+\left(\frac{2445\nu}{14}+\frac{9599}{28}\right)\frac{p_{r}^{2}}{r}+\frac{421\nu}{14}+\frac{6213}{56}\frac{1}{r^{2}}\right]+\frac{1}{c^{4}}\left[\left(\frac{1319\nu^{2}}{42}+\frac{109609\nu}{336}-\frac{45875}{1008}\right)p^{6}\right.\right.$
$\displaystyle\left.\left.+\left(-\frac{1319\nu^{2}}{14}-\frac{109609\nu}{112}+\frac{45875}{336}\right)p^{4}p_{r}^{2}+\left(-\frac{2491\nu^{2}}{42}-\frac{100847\nu}{336}+\frac{7595}{144}\right)\frac{p^{4}}{r}+\left(\frac{439\nu^{2}}{21}-\frac{2117\nu}{12}+\frac{16775}{56}\right)\frac{p^{2}}{r^{2}}\right.\right.$
$\displaystyle\left.\left.+\left(-\frac{925\nu^{2}}{168}+\frac{10855\nu}{168}-\frac{1345}{168}\right)p^{2}p_{r}^{4}+\left(-\frac{15121\nu^{2}}{56}-\frac{63683\nu}{42}+\frac{13037}{63}\right)\frac{p^{2}p_{r}^{2}}{r}+\left(-\frac{401\nu^{2}}{42}-\frac{15405\nu}{28}-\frac{44015}{56}\right)\frac{p_{r}^{2}}{r^{2}}\right.\right.$
$\displaystyle\left.\left.+\left(-\frac{95\nu^{2}}{3}-\frac{775\nu}{12}+\frac{65}{2}\right)p_{r}^{6}+\left(-\frac{139\nu^{2}}{12}+\frac{129\nu}{7}+\frac{1542}{7}\right)\frac{p_{r}^{4}}{r}+\left(\nu^{2}+\frac{1216\nu}{21}-\frac{47561}{378}\right)\frac{1}{r^{3}}\right]\right\\}.$
(44)

the spin-orbit coefficients
$f^{r}_{\text{so}},f^{p}_{\text{so}},f^{L}_{\text{so}}$ are

$\displaystyle
f^{r}_{\text{so}}=\frac{2p_{r}^{3}\nu}{15r^{5}L^{2}}\left[\boldsymbol{L}\cdot\boldsymbol{S}_{1}\left(-19\delta+10\nu+90p^{2}r(3\delta-2\nu-3)-105r(3\delta-2\nu-3)p_{r}^{2}+19\right)\right.$
$\displaystyle\qquad\left.+\boldsymbol{L}\cdot\boldsymbol{S}_{2}\left(19\delta+10\nu-90p^{2}r(3\delta+2\nu+3)+105r(3\delta+2\nu+3)p_{r}^{2}+19\right)\right],$
(45) $\displaystyle
f^{p}_{\text{so}}=\frac{2\nu}{15r^{4}L^{2}}\left\\{\boldsymbol{L}\cdot\boldsymbol{S}_{1}\left[2p^{4}r(\delta+11\nu-1)+p^{2}\left(2(53\delta-73\nu-53)+3r(41\delta-74\nu-41)p_{r}^{2}\right)\right.\right.$
$\displaystyle\qquad\left.\left.+p_{r}^{2}\left(-125\delta+156\nu+10r(-17\delta+23\nu+17)p_{r}^{2}+125\right]\right)\right.$
$\displaystyle\qquad\left.+\boldsymbol{L}\cdot\boldsymbol{S}_{2}\left[-2p^{4}r(\delta-11\nu+1)-p^{2}\left(2(53(\delta+1)+73\nu)+3r(41(\delta+1)+74\nu)p_{r}^{2}\right)\right.\right.$
$\displaystyle\qquad\left.\left.+p_{r}^{2}\left(125(\delta+1)+156\nu+10r(17(\delta+1)+23\nu)p_{r}^{2}\right)\right]\right\\},$
(46) $\displaystyle
f^{L}_{\text{so}}=\frac{2p_{r}^{2}}{15r^{4}L^{2}}\left\\{\frac{p_{r}\nu}{r}\boldsymbol{r}\cdot\boldsymbol{S}_{1}\left[-19\delta+10\nu+15p^{2}r(-3\delta+2\nu+3)+19\right]+\frac{p_{r}\nu}{r}\boldsymbol{r}\cdot\boldsymbol{S}_{2}\left[19(\delta+1)+10\nu+15p^{2}r(3\delta+2\nu+3)\right]\right.$
$\displaystyle\qquad\left.+\boldsymbol{p}\cdot\boldsymbol{S}_{1}\left[\nu(19\delta-10\nu-19)-15p^{2}r\left(5(\delta-3)\nu-3\delta-2\nu^{2}+3\right)+15r\left(8\delta\nu-3\delta-4\nu^{2}-18\nu+3\right)p_{r}^{2}\right]\right.$
$\displaystyle\qquad\left.+\boldsymbol{p}\cdot\boldsymbol{S}_{2}\left[-\nu(19\delta+10\nu+19)+15p^{2}r\left(5(\delta+3)\nu-3(\delta+1)+2\nu^{2}\right)-15r\left(8\delta\nu-3\delta+4\nu^{2}+18\nu-3\right)p_{r}^{2}\right]\right\\}$
(47)

the spin-spin coefficients
$f^{r}_{\text{ss}},f^{p}_{\text{ss}},f^{1}_{\text{ss}},f^{2}_{\text{ss}}$ are

$\displaystyle
f^{r}_{\text{ss}}=\frac{2p_{r}^{2}}{15r^{7}}\left[15rp_{r}\left(\boldsymbol{S}_{1}\cdot\boldsymbol{S}_{1}\left(7\delta+2\nu^{2}+14\nu-7\right)+4(\nu-7)\nu\boldsymbol{S}_{1}\cdot\boldsymbol{S}_{2}+\boldsymbol{S}_{2}\cdot\boldsymbol{S}_{2}\left(-7\delta+2\nu^{2}+14\nu-7\right)\right)\right.$
$\displaystyle\qquad\left.+\nu(244\nu\boldsymbol{p}\cdot\boldsymbol{S}_{1}\boldsymbol{r}\cdot\boldsymbol{S}_{1}+488\nu\boldsymbol{p}\cdot\boldsymbol{S}_{1}\boldsymbol{r}\cdot\boldsymbol{S}_{2}-105\boldsymbol{p}\cdot\boldsymbol{S}_{1}\boldsymbol{r}\cdot\boldsymbol{S}_{2}+105\boldsymbol{p}\cdot\boldsymbol{S}_{2}\boldsymbol{r}\cdot\boldsymbol{S}_{1}+244\nu\boldsymbol{p}\cdot\boldsymbol{S}_{2}\boldsymbol{r}\cdot\boldsymbol{S}_{2})\right],$
(48) $\displaystyle f^{p}_{\text{ss}}=-\frac{1}{30r^{8}}\left[-171\delta
p^{2}r^{3}\boldsymbol{S}_{1}\cdot\boldsymbol{S}_{1}-400\nu^{2}p^{2}r^{3}\boldsymbol{S}_{1}\cdot\boldsymbol{S}_{1}-342\nu
p^{2}r^{3}\boldsymbol{S}_{1}\cdot\boldsymbol{S}_{1}+171p^{2}r^{3}\boldsymbol{S}_{1}\cdot\boldsymbol{S}_{1}-800\nu^{2}p^{2}r^{3}\boldsymbol{S}_{1}\cdot\boldsymbol{S}_{2}\right.$
$\displaystyle\qquad\left.+756\nu
p^{2}r^{3}\boldsymbol{S}_{1}\cdot\boldsymbol{S}_{2}+171\delta
p^{2}r^{3}\boldsymbol{S}_{2}\cdot\boldsymbol{S}_{2}-400\nu^{2}p^{2}r^{3}\boldsymbol{S}_{2}\cdot\boldsymbol{S}_{2}-342\nu
p^{2}r^{3}\boldsymbol{S}_{2}\cdot\boldsymbol{S}_{2}+171p^{2}r^{3}\boldsymbol{S}_{2}\cdot\boldsymbol{S}_{2}\right.$
$\displaystyle\qquad\left.+1035\delta
r^{3}\boldsymbol{S}_{1}\cdot\boldsymbol{S}_{1}p_{r}^{2}+632\nu^{2}r^{3}\boldsymbol{S}_{1}\cdot\boldsymbol{S}_{1}p_{r}^{2}+2070\nu
r^{3}\boldsymbol{S}_{1}\cdot\boldsymbol{S}_{1}p_{r}^{2}-1035r^{3}\boldsymbol{S}_{1}\cdot\boldsymbol{S}_{1}p_{r}^{2}+1264\nu^{2}r^{3}\boldsymbol{S}_{1}\cdot\boldsymbol{S}_{2}p_{r}^{2}\right.$
$\displaystyle\qquad\left.-4500\nu
r^{3}\boldsymbol{S}_{1}\cdot\boldsymbol{S}_{2}p_{r}^{2}-1035\delta
r^{3}\boldsymbol{S}_{2}\cdot\boldsymbol{S}_{2}p_{r}^{2}+632\nu^{2}r^{3}\boldsymbol{S}_{2}\cdot\boldsymbol{S}_{2}p_{r}^{2}+2070\nu
r^{3}\boldsymbol{S}_{2}\cdot\boldsymbol{S}_{2}p_{r}^{2}-1035r^{3}\boldsymbol{S}_{2}\cdot\boldsymbol{S}_{2}p_{r}^{2}\right.$
$\displaystyle\qquad\left.+480\nu^{2}r(\boldsymbol{r}\cdot\boldsymbol{S}_{1})^{2}p_{r}^{2}+960\nu^{2}r(\boldsymbol{r}\cdot\boldsymbol{S}_{1})(\boldsymbol{r}\cdot\boldsymbol{S}_{2})p_{r}^{2}+480\nu^{2}r(\boldsymbol{r}\cdot\boldsymbol{S}_{2})^{2}p_{r}^{2}+48(\boldsymbol{p}\cdot\boldsymbol{S}_{1})^{2}r^{3}\left(9\delta+14\nu^{2}+3\nu-9\right)\right.$
$\displaystyle\qquad\left.+96\nu(14\nu-3)(\boldsymbol{p}\cdot\boldsymbol{S}_{1})(\boldsymbol{p}\cdot\boldsymbol{S}_{2})r^{3}-48(\boldsymbol{p}\cdot\boldsymbol{S}_{2})^{2}r^{3}\left(9\delta-14\nu^{2}-3\nu+9\right)+585\delta
r^{2}\boldsymbol{S}_{1}\cdot\boldsymbol{S}_{1}+704\nu^{2}r^{2}\boldsymbol{S}_{1}\cdot\boldsymbol{S}_{1}\right.$
$\displaystyle\qquad\left.+1170\nu
r^{2}\boldsymbol{S}_{1}\cdot\boldsymbol{S}_{1}-585r^{2}\boldsymbol{S}_{1}\cdot\boldsymbol{S}_{1}+1408\nu^{2}r^{2}\boldsymbol{S}_{1}\cdot\boldsymbol{S}_{2}-2364\nu
r^{2}\boldsymbol{S}_{1}\cdot\boldsymbol{S}_{2}-585\delta
r^{2}\boldsymbol{S}_{2}\cdot\boldsymbol{S}_{2}+704\nu^{2}r^{2}\boldsymbol{S}_{2}\cdot\boldsymbol{S}_{2}\right.$
$\displaystyle\qquad\left.+1170\nu
r^{2}\boldsymbol{S}_{2}\cdot\boldsymbol{S}_{2}-585r^{2}\boldsymbol{S}_{2}\cdot\boldsymbol{S}_{2}-1368\delta(\boldsymbol{r}\cdot\boldsymbol{S}_{1})^{2}-624\nu^{2}(\boldsymbol{r}\cdot\boldsymbol{S}_{1})^{2}-2736\nu(\boldsymbol{r}\cdot\boldsymbol{S}_{1})^{2}+1368(\boldsymbol{r}\cdot\boldsymbol{S}_{1})^{2}\right.$
$\displaystyle\qquad\left.-1248\nu^{2}(\boldsymbol{r}\cdot\boldsymbol{S}_{1})(\boldsymbol{r}\cdot\boldsymbol{S}_{2})+5472\nu(\boldsymbol{r}\cdot\boldsymbol{S}_{1})(\boldsymbol{r}\cdot\boldsymbol{S}_{2})+1368\delta(\boldsymbol{r}\cdot\boldsymbol{S}_{2})^{2}-624\nu^{2}(\boldsymbol{r}\cdot\boldsymbol{S}_{2})^{2}-2736\nu(\boldsymbol{r}\cdot\boldsymbol{S}_{2})^{2}+1368(\boldsymbol{r}\cdot\boldsymbol{S}_{2})^{2}\right],$
(49) $\displaystyle
f^{1}_{\text{ss}}=\frac{2p_{r}^{2}\nu}{15r^{6}}\left[244\nu\boldsymbol{r}\cdot\boldsymbol{S}_{1}p_{r}+105\boldsymbol{r}\cdot\boldsymbol{S}_{2}p_{r}-244\
\nu\boldsymbol{p}\cdot\boldsymbol{S}_{1}r-105\boldsymbol{p}\cdot\boldsymbol{S}_{2}r\right],$
(50) $\displaystyle
f^{2}_{\text{ss}}=\frac{2p_{r}^{2}\nu}{15r^{6}}\left[p_{r}(488\nu\boldsymbol{r}\cdot\boldsymbol{S}_{1}-105\boldsymbol{r}\cdot\boldsymbol{S}_{1}+244\nu\boldsymbol{r}\cdot\boldsymbol{S}_{2})+(105-488\nu)\boldsymbol{p}\cdot\boldsymbol{S}_{1}r-244\nu\boldsymbol{p}\cdot\boldsymbol{S}_{2}r\right].$
(51)

##  Appendix B The quasi-circular waveform correction

The quasi-circular corrections resulted by the spin miss-aligned part are
given by

$\displaystyle\rho_{22}=1+\frac{1}{84}v_{\Omega}^{2}(21\chi_{a\lambda}+55\nu+21i\chi_{an}+21\delta(\chi_{s\lambda}+i\chi_{sn})-86)-\frac{2}{3}v_{\Omega}^{3}(\delta\chi_{ae}-\nu\chi_{se}+\chi_{se})$
$\displaystyle+v_{\Omega}^{4}\left(\delta\chi_{a\lambda}\chi_{s\lambda}\nu+\frac{7\delta\chi_{a\lambda}\chi_{s\lambda}}{16}+\frac{11\delta\chi_{s\lambda}\nu}{336}+\frac{13\delta\chi_{s\lambda}}{28}-2\chi_{ae}^{2}\nu+\frac{\chi_{ae}^{2}}{2}+\delta\chi_{ae}\chi_{se}+\frac{\chi_{se}^{2}}{2}-\chi_{a\lambda}^{2}\nu+\frac{7\chi_{a\lambda}^{2}}{32}\right.$
$\displaystyle\left.-\frac{121\chi_{a\lambda}\nu}{336}+\frac{13\chi_{a\lambda}}{28}-\chi_{s\lambda}^{2}\nu^{2}+\frac{9\chi_{s\lambda}^{2}\nu}{8}+\frac{7\chi_{s\lambda}^{2}}{32}+\frac{19583\nu^{2}}{42336}-\frac{33025\nu}{21168}+\frac{\chi_{an}^{2}}{32}-\frac{25}{16}i\delta\chi_{s\lambda}\chi_{an}+6i\chi_{a\lambda}\nu\chi_{an}\right.$
$\displaystyle\left.-\frac{25i\chi_{a\lambda}\chi_{an}}{16}-\frac{317i\nu\chi_{an}}{112}-\delta\nu\chi_{an}\chi_{sn}+\frac{\delta\chi_{an}\chi_{sn}}{16}+\frac{37i\chi_{an}}{42}+\nu^{2}\chi_{sn}^{2}-\frac{9\nu\chi_{sn}^{2}}{8}+\frac{\chi_{sn}^{2}}{32}-\frac{25}{16}i\delta\chi_{a\lambda}\chi_{sn}\right.$
$\displaystyle\left.-\frac{167}{336}i\delta\nu\chi_{sn}+\frac{37i\delta\chi_{sn}}{42}+\frac{1}{4}i\chi_{s\lambda}\nu\chi_{sn}-\frac{25i\chi_{s\lambda}\chi_{sn}}{16}-\frac{20555}{10584}\right),$
(52) $\displaystyle
f^{s}_{21}=v_{\Omega}\left(-\frac{3\chi_{ae}}{2\delta}-\frac{3\chi_{se}}{2}\right)+v_{\Omega}^{2}\left(\frac{\chi_{s\lambda}\nu}{\delta}-6i\chi_{an}+\frac{2i\nu\chi_{sn}}{\delta}-\frac{6i\chi_{sn}}{\delta}\right)+v_{\Omega}^{3}\left(-\frac{6\chi_{ae}\chi_{a\lambda}\nu}{\delta}+\frac{3\chi_{ae}\chi_{a\lambda}}{2\delta}\right.$
$\displaystyle\left.+\frac{131\chi_{ae}\nu}{84\delta}+\frac{61\chi_{ae}}{12\delta}-3\chi_{ae}\chi_{s\lambda}\nu+\frac{3\chi_{ae}\chi_{s\lambda}}{2}-\frac{30i\chi_{ae}\nu\chi_{an}}{\delta}+\frac{15i\chi_{ae}\chi_{an}}{2\delta}+3i\chi_{ae}\nu\chi_{sn}+\frac{15i\chi_{ae}\chi_{sn}}{2}\right.$
$\displaystyle\left.+\frac{6\chi_{se}\chi_{s\lambda}\nu^{2}}{\delta}-\frac{6\chi_{se}\chi_{s\lambda}\nu}{\delta}+\frac{3\chi_{se}\chi_{s\lambda}}{2\delta}-3\chi_{se}\chi_{a\lambda}\nu+\frac{3\chi_{se}\chi_{a\lambda}}{2}+\frac{79\chi_{se}\nu}{84}+3i\chi_{se}\nu\chi_{an}+\frac{15i\chi_{se}\chi_{an}}{2}\right.$
$\displaystyle\left.-\frac{6i\chi_{se}\nu^{2}\chi_{sn}}{\delta}+\frac{6i\chi_{se}\nu\chi_{sn}}{\delta}+\frac{15i\chi_{se}\chi_{sn}}{2\delta}+\frac{61\chi_{se}}{12}\right)$
(53) $\displaystyle
f^{s}_{33}=v_{\Omega}^{2}\left(-\frac{16\chi_{s\lambda}\nu}{9\delta}-\frac{16i\nu\chi_{sn}}{9\delta}\right)+v_{\Omega}^{3}\left(\frac{19\chi_{ae}\nu}{2\delta}-\frac{2\chi_{ae}}{\delta}+\frac{5\chi_{se}\nu}{2}-2\chi_{se}\right)$
(54)
$\displaystyle\rho_{32}=1+\frac{4v_{\Omega}\chi_{se}\nu}{3-9\nu}+v_{\Omega}^{2}\left[-\frac{7\delta\chi_{s\lambda}\nu^{2}}{8(1-3\nu)^{2}}-\frac{17\delta\chi_{s\lambda}\nu}{24(1-3\nu)^{2}}+\frac{\delta\chi_{s\lambda}}{3(1-3\nu)^{2}}-\frac{16\chi_{se}^{2}\nu^{2}}{9(1-3\nu)^{2}}+\frac{17\chi_{a\lambda}\nu^{2}}{8(1-3\nu)^{2}}\right.$
$\displaystyle\left.-\frac{41\chi_{a\lambda}\nu}{24(1-3\nu)^{2}}+\frac{\chi_{a\lambda}}{3(1-3\nu)^{2}}-\frac{45i\nu^{2}\chi_{an}}{2(1-3\nu)^{2}}+\frac{27i\nu\chi_{an}}{2(1-3\nu)^{2}}-\frac{2i\chi_{an}}{(1-3\nu)^{2}}-\frac{i\delta\nu^{2}\chi_{sn}}{(1-3\nu)^{2}}+\frac{19i\delta\nu\chi_{sn}}{3(1-3\nu)^{2}}-\frac{2i\delta\chi_{sn}}{(1-3\nu)^{2}}\right]$
(55) $\displaystyle
f^{s}_{31}=v_{\Omega}^{2}\left(\frac{16\chi_{s\lambda}\nu}{\delta}-\frac{16i\nu\chi_{sn}}{\delta}\right)+v_{\Omega}^{3}\left(\frac{11\chi_{ae}\nu}{2\delta}-\frac{2\chi_{ae}}{\delta}+\frac{13\chi_{se}\nu}{2}-2\chi_{se}\right)$
(56)
$\displaystyle\rho_{44}=1+v_{\Omega}^{2}\left[\frac{81\delta\chi_{s\lambda}\nu}{256(3\nu-1)}-\frac{81\chi_{a\lambda}\nu}{256(3\nu-1)}+\frac{175\nu^{2}}{88(3\nu-1)}-\frac{587\nu}{132(3\nu-1)}+\frac{269}{220(3\nu-1)}-\frac{81i\nu\chi_{an}}{256(3\nu-1)}+\frac{81i\delta\nu\chi_{sn}}{256(3\nu-1)}\right]$
(57) $\displaystyle
f^{s}_{43}=\frac{5v_{\Omega}\eta}{2(1-2\eta)}\left(\chi_{se}-\frac{\chi_{ae}}{\delta}\right)$
(58)
$\displaystyle\rho_{42}=1+v_{\Omega}^{2}\left[-\frac{39\delta\chi_{s\lambda}\nu}{32(3\nu-1)}+\frac{39\chi_{a\lambda}\nu}{32(3\nu-1)}+\frac{19\nu^{2}}{88(3\nu-1)}-\frac{353\nu}{132(3\nu-1)}+\frac{191}{220(3\nu-1)}-\frac{21i\nu\chi_{an}}{16(3\nu-1)}+\frac{21i\delta\nu\chi_{sn}}{16(3\nu-1)}\right]$
(59) $\displaystyle
f^{s}_{41}=\frac{5v_{\Omega}\eta}{2(1-2\eta)}\left(\chi_{se}-\frac{\chi_{ae}}{\delta}\right),$
(60)

##  Appendix C test the spin-precession SEOBNRE waveform model against the
spin-aligned BBH

In Liu et al. (2022) we have developed higher modes for SEOBNRE waveform
model, including (2,1), (3,3) and (4,4). In the current paper we construct
(5,5) mode also. Although the SEOBNRE waveform model proposed here is for
generic eccentric spin-precession BBHs, it also works for spin-aligned BBHs.
As an example we show the comparison between the NR waveform SXS:BBH:1374 and
the corresponding SEOBNRE waveform in Fig. 12.

---
Figure 12: Waveforms comparison between the NR waveform SXS:BBH:1374 and the
corresponding SEOBNRE waveform. This is a spinless BBH with mass ratio $q=3$.

##  Appendix D SXS waveforms used in the current work

Table 1: The NR waveforms used in this work. Here, we list the parameters for each waveform, including mass ratio $q$, initial frequency $Mf_{\text{ini}}$, and initial spin configuration. Additionally, we provide the maximum mismatch between the two EOB models and these waveforms. The subscription ‘old’ means SEOBNRv4PHM waveform model, and ‘new’ means SEOBNRE waveform model. SXS ID | $q$ | $\boldsymbol{\chi}_{1}$ | $\boldsymbol{\chi}_{2}$ | $10^{3}Mf_{\text{ini}}$ | $\max\mathcal{M}_{\text{old}}$ | $\max\mathcal{M}_{\text{new}}$
---|---|---|---|---|---|---
0401 | 1.50 | (0.38,0.42,-0.57) | (-0.03,0.03,-0.80) | 0.45 | 0.16% | 0.19%
0622 | 1.20 | (0.42,0.42,0.60) | (0.54,0.54,0.37) | 0.53 | 2.17% | 2.15%
0623 | 1.11 | (0.89,0.00,0.11) | (0.89,0.00,-0.12) | 0.31 | 2.48% | 2.31%
0628 | 2.06 | (-0.38,0.75,-0.10) | (0.08,-0.40,-0.30) | 0.61 | 0.60% | 0.71%
0632 | 1.08 | (0.72,0.24,-0.09) | (0.63,-0.49,-0.39) | 0.51 | 1.21% | 1.42%
0633 | 1.34 | (0.38,-0.16,-0.66) | (-0.79,-0.06,0.26) | 0.50 | 3.80% | 1.44%
0651 | 1.33 | (-0.00,0.00,-0.00) | (0.58,0.55,0.00) | 0.48 | 0.27% | 0.25%
0715 | 2.00 | (-0.79,-0.10,-0.04) | (0.30,-0.74,0.08) | 0.50 | 0.92% | 0.85%
0721 | 2.00 | (0.31,0.74,-0.03) | (-0.79,0.11,0.08) | 0.49 | 1.23% | 1.00%
0753 | 1.00 | (0.44,-0.35,0.57) | (-0.00,-0.00,-0.00) | 0.49 | 0.20% | 0.18%
0787 | 2.00 | (-0.48,-0.25,-0.59) | (0.09,-0.79,0.06) | 0.48 | 2.49% | 1.26%
0888 | 2.00 | (-0.77,-0.21,-0.01) | (-0.50,-0.21,0.59) | 0.49 | 2.51% | 1.78%
0891 | 1.00 | (0.01,-0.04,-0.80) | (-0.56,-0.07,0.56) | 0.46 | 0.62% | 0.45%
0904 | 2.00 | (-0.00,-0.03,0.80) | (0.56,-0.08,0.57) | 0.52 | 0.27% | 0.24%
0906 | 2.00 | (0.57,-0.56,0.01) | (0.42,-0.33,0.59) | 0.50 | 1.22% | 1.26%
0940 | 2.00 | (0.00,-0.00,-0.00) | (-0.54,-0.16,-0.57) | 0.49 | 0.70% | 0.66%
0951 | 2.00 | (0.38,0.70,-0.01) | (0.07,0.59,-0.54) | 0.49 | 6.16% | 2.11%
0961 | 2.00 | (0.00,0.00,-0.00) | (-0.28,0.75,0.00) | 0.50 | 0.45% | 0.41%
1106 | 1.68 | (0.72,-0.00,0.36) | (-0.14,-0.07,0.05) | 0.29 | 6.14% | 0.92%
1184 | 1.00 | (0.06,-0.04,0.85) | (0.40,0.70,0.01) | 0.49 | 1.72% | 1.48%
1185 | 1.00 | (-0.73,0.43,-0.00) | (-0.03,-0.05,-0.80) | 0.45 | 0.35% | 0.33%
1188 | 1.00 | (0.62,0.58,0.05) | (0.19,-0.78,-0.05) | 0.48 | 0.10% | 0.13%
1189 | 1.00 | (0.61,0.59,-0.05) | (-0.77,0.23,0.06) | 0.47 | 0.31% | 0.23%
1190 | 1.00 | (0.74,0.42,-0.00) | (-0.03,0.05,-0.80) | 0.45 | 0.21% | 0.31%
1191 | 1.00 | (0.20,-0.82,0.06) | (-0.77,0.23,-0.06) | 0.48 | 0.35% | 0.41%
1192 | 1.00 | (-0.01,-0.85,0.00) | (0.06,-0.00,-0.80) | 0.46 | 0.28% | 0.36%
1193 | 2.00 | (0.54,-0.66,0.03) | (0.35,-0.78,0.00) | 0.50 | 2.54% | 2.86%
1194 | 2.00 | (0.51,-0.68,0.05) | (-0.84,0.12,-0.08) | 0.49 | 1.30% | 1.14%
1195 | 1.00 | (-0.81,0.24,-0.06) | (0.20,-0.82,0.06) | 0.48 | 0.59% | 0.52%
1196 | 1.00 | (0.23,-0.82,0.01) | (0.23,-0.82,0.01) | 0.49 | 2.08% | 1.84%
1197 | 2.00 | (-0.84,-0.11,-0.04) | (0.32,-0.78,0.09) | 0.50 | 1.25% | 1.13%
1198 | 1.00 | (0.43,-0.73,0.01) | (-0.07,-0.03,0.85) | 0.49 | 1.68% | 1.49%
1199 | 2.00 | (0.39,-0.75,0.01) | (0.09,0.03,-0.84) | 0.48 | 1.87% | 1.67%
1200 | 2.00 | (-0.84,-0.14,0.01) | (-0.85,0.09,-0.01) | 0.48 | 2.02% | 2.68%
1201 | 1.00 | (-0.03,-0.05,-0.85) | (-0.73,0.44,-0.00) | 0.45 | 0.25% | 0.35%
1202 | 2.00 | (-0.01,-0.04,-0.85) | (-0.73,0.44,0.00) | 0.46 | 0.06% | 0.08%
1203 | 2.00 | (-0.85,0.04,-0.00) | (-0.02,-0.09,-0.85) | 0.47 | 2.51% | 2.31%
1204 | 2.00 | (0.33,0.78,0.04) | (0.32,-0.78,-0.08) | 0.50 | 1.33% | 0.88%
1205 | 1.00 | (-0.81,0.24,0.06) | (0.61,0.59,-0.06) | 0.47 | 0.35% | 0.24%
1206 | 1.00 | (0.20,-0.82,-0.05) | (0.62,0.58,0.06) | 0.48 | 0.32% | 0.20%
1208 | 1.00 | (-0.03,0.05,-0.85) | (0.75,0.41,-0.00) | 0.45 | 0.51% | 0.38%
1209 | 2.00 | (0.06,-0.01,0.85) | (0.18,0.83,0.01) | 0.53 | 1.20% | 1.37%
1211 | 2.00 | (-0.03,0.03,-0.85) | (0.74,0.41,-0.00) | 0.46 | 0.26% | 0.13%
1212 | 2.00 | (0.17,0.83,0.01) | (0.10,-0.04,0.84) | 0.50 | 2.49% | 2.24%
1213 | 2.00 | (0.46,0.72,0.00) | (-0.07,0.06,-0.84) | 0.48 | 1.57% | 1.49%
1214 | 1.00 | (0.60,0.60,0.00) | (0.60,0.60,0.00) | 0.48 | 2.89% | 1.97%
1215 | 2.00 | (0.30,0.79,0.02) | (0.50,0.69,-0.01) | 0.51 | 1.52% | 1.59%
1216 | 2.00 | (0.40,-0.75,0.01) | (0.09,0.03,-0.79) | 0.48 | 1.97% | 1.78%
1217 | 1.00 | (0.23,-0.82,0.01) | (0.21,-0.77,0.01) | 0.49 | 1.53% | 1.79%
1218 | 1.00 | (0.20,-0.82,-0.05) | (0.58,0.55,0.05) | 0.48 | 0.28% | 0.19%
1219 | 1.00 | (0.42,-0.74,0.01) | (-0.06,-0.02,0.80) | 0.49 | 1.58% | 1.34%
1909 | 4.00 | (0.00,-0.00,-0.00) | (-0.32,-0.73,0.01) | 0.56 | 0.07% | 0.07%
1918 | 4.00 | (-0.80,-0.07,0.02) | (-0.05,-0.10,-0.79) | 0.54 | 2.25% | 1.40%
1930 | 4.00 | (0.34,0.72,0.00) | (-0.09,0.09,-0.79) | 0.54 | 1.86% | 1.32%
1968 | 4.00 | (-0.34,0.43,0.58) | (0.11,0.03,0.79) | 0.58 | 1.34% | 1.08%
1969 | 4.00 | (0.15,-0.53,0.58) | (0.11,-0.00,-0.79) | 0.57 | 1.01% | 1.11%
1971 | 4.00 | (-0.19,-0.53,0.57) | (-0.08,0.05,0.79) | 0.58 | 0.63% | 0.71%
1981 | 4.00 | (0.08,-0.59,-0.53) | (-0.79,-0.06,-0.08) | 0.52 | 3.93% | 3.81%
1999 | 4.00 | (0.20,-0.78,0.00) | (0.69,-0.38,0.14) | 0.55 | 2.11% | 2.42%
2009 | 4.00 | (0.29,0.74,0.01) | (0.09,-0.08,0.79) | 0.55 | 2.07% | 1.67%
2012 | 4.00 | (-0.45,0.66,0.04) | (-0.14,-0.01,-0.79) | 0.54 | 2.70% | 2.15%
2060 | 4.00 | (0.00,-0.00,-0.00) | (-0.23,-0.52,-0.56) | 0.55 | 0.30% | 0.31%
2080 | 4.00 | (0.78,-0.19,0.05) | (-0.79,0.06,-0.08) | 0.54 | 1.39% |

## References

  * Aasi et al. (2015) J. Aasi, B. Abbott, R. Abbott, T. Abbott, M. Abernathy, K. Ackley, C. Adams, T. Adams, P. Addesso, R. Adhikari, et al., Classical and quantum gravity 32, 074001 (2015).
  * Acernese et al. (2014) F. Acernese et al., Classical and Quantum Gravity 32, 024001 (2014).
  * Abbott et al. (2016) B. P. Abbott et al. (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. Lett. 116, 061102 (2016), URL https://link.aps.org/doi/10.1103/PhysRevLett.116.061102.
  * The LIGO Scientific Collaboration et al. (2021) The LIGO Scientific Collaboration, The Virgo Collaboration, The KAGRA Collaboration, R. Abbott, et al. (2021).
  * Abbott et al. (2017) B. P. Abbott et al. (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. Lett. 119, 161101 (2017), URL https://link.aps.org/doi/10.1103/PhysRevLett.119.161101.
  * Abbott et al. (2021) R. Abbott et al., The Astrophysical Journal Letters 915, L5 (2021).
  * Akutsu et al. (2019) T. Akutsu et al., Nature Astronomy 3, 35 (2019).
  * Punturo et al. (2010) M. Punturo et al., Classical and Quantum Gravity 27, 194002 (2010).
  * Reitze et al. (2019) D. Reitze, R. X. Adhikari, S. Ballmer, B. Barish, L. Barsotti, G. Billingsley, D. A. Brown, Y. Chen, D. Coyne, R. Eisenstein, et al. (2019).
  * Amaro-Seoane et al. (2017) P. Amaro-Seoane et al., arXiv e-prints arXiv:1702.00786 (2017), eprint 1702.00786.
  * Hu and Wu (2017) W.-R. Hu and Y.-L. Wu, National Science Review 4, 685 (2017).
  * Luo et al. (2016) J. Luo, L.-S. Chen, H.-Z. Duan, Y.-G. Gong, S. Hu, J. Ji, Q. Liu, J. Mei, V. Milyukov, M. Sazhin, et al., Classical and Quantum Gravity 33, 035010 (2016).
  * Einstein et al. (1938) A. Einstein, L. Infeld, and B. Hoffmann, The Annals of Mathematics 39, 65 (1938).
  * Lorentz and Droste (1937) H. A. Lorentz and J. Droste, pp. 330–355 (1937).
  * Ohta et al. (1973) T. Ohta, H. Okamura, T. Kimura, and K. Hiida, Progress of Theoretical Physics 50, 492 (1973).
  * Damour and Sch?fer (1985) T. Damour and G. Sch?fer, General Relativity and Gravitation 17, 879 (1985).
  * Blanchet (2006) L. Blanchet, Living Reviews in Relativity 9 (2006).
  * Jaranowski and Sch?fer (2012) P. Jaranowski and G. Sch?fer, Physical Review D 86, 061503 (2012).
  * Damour et al. (2016) T. Damour, P. Jaranowski, and G. Sch?fer, Physical Review D 93, 084014 (2016).
  * Bernard et al. (2018) L. Bernard, L. Blanchet, G. Faye, and T. Marchand, Physical Review D 97, 044037 (2018).
  * Finn (1992) L. S. Finn, Phys. Rev. D 46, 5236 (1992), URL https://link.aps.org/doi/10.1103/PhysRevD.46.5236.
  * Jaranowski and Krolak (1994) P. Jaranowski and A. Krolak, Phys. Rev. D 49, 1723 (1994), URL https://link.aps.org/doi/10.1103/PhysRevD.49.1723.
  * Jaranowski et al. (1998) P. Jaranowski, A. Królak, and B. F. Schutz, Phys. Rev. D 58, 063001 (1998), URL https://link.aps.org/doi/10.1103/PhysRevD.58.063001.
  * Abbott et al. (2020) R. Abbott et al. (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. Lett. 125, 101102 (2020), URL https://link.aps.org/doi/10.1103/PhysRevLett.125.101102.
  * Chiaramello and Nagar (2020) D. Chiaramello and A. Nagar, Phys. Rev. D 101, 101501 (2020), URL https://link.aps.org/doi/10.1103/PhysRevD.101.101501.
  * Nagar et al. (2021) A. Nagar, A. Bonino, and P. Rettegno, Phys. Rev. D 103, 104021 (2021), URL https://link.aps.org/doi/10.1103/PhysRevD.103.104021.
  * Cao and Han (2017) Z. Cao and W.-B. Han, Phys. Rev. D 96, 044028 (2017), URL https://link.aps.org/doi/10.1103/PhysRevD.96.044028.
  * Liu et al. (2020) X. Liu, Z. Cao, and L. Shao, Phys. Rev. D 101, 044049 (2020), URL https://link.aps.org/doi/10.1103/PhysRevD.101.044049.
  * Yun et al. (2021) Q. Yun, W.-B. Han, X. Zhong, and C. A. Benavides-Gallego, Phys. Rev. D 103, 124053 (2021), URL https://link.aps.org/doi/10.1103/PhysRevD.103.124053.
  * Liu et al. (2022) X. Liu, Z. Cao, and Z.-H. Zhu, Classical and Quantum Gravity 39, 035009 (2022).
  * Liu et al. (2023) X. Liu, Z. Cao, and L. Shao, International Journal of Modern Physics D 32, 2350015 (2023), eprint 2306.15277.
  * Ramos-Buades et al. (2022) A. Ramos-Buades, A. Buonanno, M. Khalil, and S. Ossokine, Phys. Rev. D 105, 044035 (2022), URL https://link.aps.org/doi/10.1103/PhysRevD.105.044035.
  * Khalil et al. (2021) M. Khalil, A. Buonanno, J. Steinhoff, and J. Vines, Phys. Rev. D 104, 024046 (2021), URL https://link.aps.org/doi/10.1103/PhysRevD.104.024046.
  * Apostolatos et al. (1994) T. A. Apostolatos, C. Cutler, G. J. Sussman, and K. S. Thorne, Phys. Rev. D 49, 6274 (1994), URL https://link.aps.org/doi/10.1103/PhysRevD.49.6274.
  * Pratten et al. (2021) G. Pratten, C. García-Quirós, M. Colleoni, A. Ramos-Buades, H. Estellés, M. Mateu-Lucena, R. Jaume, M. Haney, D. Keitel, J. E. Thompson, et al., Phys. Rev. D 103, 104056 (2021), URL https://link.aps.org/doi/10.1103/PhysRevD.103.104056.
  * Pan et al. (2014) Y. Pan, A. Buonanno, A. Taracchini, L. E. Kidder, A. H. Mroué, H. P. Pfeiffer, M. A. Scheel, and B. Szilágyi, Phys. Rev. D 89, 084006 (2014), URL https://link.aps.org/doi/10.1103/PhysRevD.89.084006.
  * Ossokine et al. (2020) S. Ossokine, A. Buonanno, S. Marsat, R. Cotesta, S. Babak, T. Dietrich, R. Haas, I. Hinder, H. P. Pfeiffer, M. Pürrer, et al., Phys. Rev. D 102, 044055 (2020), URL https://link.aps.org/doi/10.1103/PhysRevD.102.044055.
  * Ramos-Buades et al. (2023) A. Ramos-Buades, A. Buonanno, H. Estellés, M. Khalil, D. P. Mihaylov, S. Ossokine, L. Pompili, and M. Shiferaw (2023).
  * Khalil et al. (2023) M. Khalil, A. Buonanno, H. Estellés, D. P. Mihaylov, S. Ossokine, L. Pompili, and A. Ramos-Buades (2023).
  * Huerta et al. (2018) E. A. Huerta, C. J. Moore, P. Kumar, D. George, A. J. K. Chua, R. Haas, E. Wessel, D. Johnson, D. Glennon, A. Rebei, et al., Phys. Rev. D 97, 024031 (2018), URL https://link.aps.org/doi/10.1103/PhysRevD.97.024031.
  * Habib and Huerta (2019) S. Habib and E. A. Huerta, Phys. Rev. D 100, 044016 (2019), URL https://link.aps.org/doi/10.1103/PhysRevD.100.044016.
  * Chen et al. (2021) Z. Chen, E. A. Huerta, J. Adamo, R. Haas, E. O’Shea, P. Kumar, and C. Moore, Phys. Rev. D 103, 084018 (2021), URL https://link.aps.org/doi/10.1103/PhysRevD.103.084018.
  * Hinder et al. (2018) I. Hinder, L. E. Kidder, and H. P. Pfeiffer, Phys. Rev. D 98, 044015 (2018), URL https://link.aps.org/doi/10.1103/PhysRevD.98.044015.
  * Islam et al. (2021) T. Islam, V. Varma, J. Lodman, S. E. Field, G. Khanna, M. A. Scheel, H. P. Pfeiffer, D. Gerosa, and L. E. Kidder, Phys. Rev. D 103, 064022 (2021), URL https://link.aps.org/doi/10.1103/PhysRevD.103.064022.
  * Setyawati and Ohme (2021) Y. Setyawati and F. Ohme, Phys. Rev. D 103, 124011 (2021), URL https://link.aps.org/doi/10.1103/PhysRevD.103.124011.
  * Wang et al. (2023) H. Wang, Y.-C. Zou, and Y. Liu, Phys. Rev. D 107, 124061 (2023), URL https://link.aps.org/doi/10.1103/PhysRevD.107.124061.
  * Peters (1964) P. C. Peters, Phys. Rev. 136, B1224 (1964), URL https://link.aps.org/doi/10.1103/PhysRev.136.B1224.
  * Mishra et al. (2015) C. K. Mishra, K. G. Arun, and B. R. Iyer, Phys. Rev. D 91, 084040 (2015), URL https://link.aps.org/doi/10.1103/PhysRevD.91.084040.
  * Tanay et al. (2016) S. Tanay, M. Haney, and A. Gopakumar, Phys. Rev. D 93, 064031 (2016), URL https://link.aps.org/doi/10.1103/PhysRevD.93.064031.
  * Moore et al. (2016) B. Moore, M. Favata, K. G. Arun, and C. K. Mishra, Phys. Rev. D 93, 124061 (2016), URL https://link.aps.org/doi/10.1103/PhysRevD.93.124061.
  * Boetzel et al. (2019a) Y. Boetzel, C. K. Mishra, G. Faye, A. Gopakumar, and B. R. Iyer, Phys. Rev. D 100, 044018 (2019a), URL https://link.aps.org/doi/10.1103/PhysRevD.100.044018.
  * Ebersold et al. (2019) M. Ebersold, Y. Boetzel, G. Faye, C. K. Mishra, B. R. Iyer, and P. Jetzer, Phys. Rev. D 100, 084043 (2019), URL https://link.aps.org/doi/10.1103/PhysRevD.100.084043.
  * Tiwari and Gopakumar (2020) S. Tiwari and A. Gopakumar, Phys. Rev. D 102, 084042 (2020), URL https://link.aps.org/doi/10.1103/PhysRevD.102.084042.
  * Loutrel and Yunes (2017a) N. Loutrel and N. Yunes, Classical and Quantum Gravity 34, 044003 (2017a), eprint 1607.05409.
  * Loutrel and Yunes (2017b) N. Loutrel and N. Yunes, Classical and Quantum Gravity 34, 135011 (2017b), eprint 1702.01818.
  * Moore et al. (2018) B. Moore, T. Robson, N. Loutrel, and N. Yunes, Classical and Quantum Gravity 35, 235006 (2018), eprint 1807.07163.
  * Loutrel (2020) N. Loutrel, Classical and Quantum Gravity 37, 075008 (2020), eprint 1909.02143.
  * Loutrel (2021) N. Loutrel, Classical and Quantum Gravity 38, 015005 (2021), eprint 2003.13673.
  * Yunes et al. (2009) N. Yunes, K. G. Arun, E. Berti, and C. M. Will, Phys. Rev. D 80, 084001 (2009), URL https://link.aps.org/doi/10.1103/PhysRevD.80.084001.
  * Huerta et al. (2014) E. A. Huerta, P. Kumar, S. T. McWilliams, R. O’Shaughnessy, and N. Yunes, Phys. Rev. D 90, 084016 (2014), URL https://link.aps.org/doi/10.1103/PhysRevD.90.084016.
  * Hinder et al. (2010) I. Hinder, F. Herrmann, P. Laguna, and D. Shoemaker, Phys. Rev. D 82, 024033 (2010), URL https://link.aps.org/doi/10.1103/PhysRevD.82.024033.
  * Romero-Shaw et al. (2020) I. Romero-Shaw, P. D. Lasky, E. Thrane, and J. Calderón Bustillo, Astrophys. J. Letters 903, L5 (2020), eprint 2009.04771.
  * Romero-Shaw et al. (2023) I. M. Romero-Shaw, D. Gerosa, and N. Loutrel, Monthly Notices of the Royal Astronomical Society 519, 5352 (2023), eprint 2211.07528.
  * Gamba et al. (2023) R. Gamba, M. Breschi, G. Carullo, S. Albanesi, P. Rettegno, S. Bernuzzi, and A. Nagar, Nature Astronomy 7, 11 (2023), eprint 2106.05575.
  * Gayathri et al. (2022) V. Gayathri, J. Healy, J. Lange, B. O’Brien, M. Szczepańczyk, I. Bartos, M. Campanelli, S. Klimenko, C. O. Lousto, and R. O’Shaughnessy, Nature Astronomy 6, 344 (2022).
  * Buonanno and Damour (1999) A. Buonanno and T. Damour, Phys. Rev. D 59, 084006 (1999), URL https://link.aps.org/doi/10.1103/PhysRevD.59.084006.
  * Barausse and Buonanno (2010) E. Barausse and A. Buonanno, Phys. Rev. D 81, 084024 (2010), URL https://link.aps.org/doi/10.1103/PhysRevD.81.084024.
  * Barausse and Buonanno (2011) E. Barausse and A. Buonanno, Phys. Rev. D 84, 104027 (2011), URL https://link.aps.org/doi/10.1103/PhysRevD.84.104027.
  * Maia et al. (2017a) N. T. Maia, C. R. Galley, A. K. Leibovich, and R. A. Porto, Phys. Rev. D 96, 084064 (2017a), URL https://link.aps.org/doi/10.1103/PhysRevD.96.084064.
  * Maia et al. (2017b) N. T. Maia, C. R. Galley, A. K. Leibovich, and R. A. Porto, Phys. Rev. D 96, 084065 (2017b), URL https://link.aps.org/doi/10.1103/PhysRevD.96.084065.
  * Bini and Damour (2012) D. Bini and T. Damour, Phys. Rev. D 86, 124012 (2012), URL https://link.aps.org/doi/10.1103/PhysRevD.86.124012.
  * Kidder et al. (1993) L. E. Kidder, C. M. Will, and A. G. Wiseman, Phys. Rev. D 47, R4183 (1993), URL https://link.aps.org/doi/10.1103/PhysRevD.47.R4183.
  * Buonanno et al. (2013) A. Buonanno, G. Faye, and T. Hinderer, Phys. Rev. D 87, 044009 (2013), URL https://link.aps.org/doi/10.1103/PhysRevD.87.044009.
  * Wang and Will (2007) H. Wang and C. M. Will, Phys. Rev. D 75, 064017 (2007), URL https://link.aps.org/doi/10.1103/PhysRevD.75.064017.
  * Pan et al. (2010) Y. Pan, A. Buonanno, L. T. Buchman, T. Chu, L. E. Kidder, H. P. Pfeiffer, and M. A. Scheel, Phys. Rev. D 81, 084041 (2010), URL https://link.aps.org/doi/10.1103/PhysRevD.81.084041.
  * Buonanno et al. (2006) A. Buonanno, Y. Chen, and T. Damour, Phys. Rev. D 74, 104005 (2006), URL https://link.aps.org/doi/10.1103/PhysRevD.74.104005.
  * Loutrel et al. (2019) N. Loutrel, S. Liebersbach, N. Yunes, and N. Cornish, Classical and Quantum Gravity 36, 025004 (2019), eprint 1810.03521.
  * Knee et al. (2022) A. M. Knee, I. M. Romero-Shaw, P. D. Lasky, J. McIver, and E. Thrane, Astrophys. J. 936, 172 (2022), eprint 2207.14346.
  * Will and Wiseman (1996) C. M. Will and A. G. Wiseman, Phys. Rev. D 54, 4813 (1996), URL https://link.aps.org/doi/10.1103/PhysRevD.54.4813.
  * Pan et al. (2011) Y. Pan, A. Buonanno, R. Fujita, E. Racine, and H. Tagoshi, Phys. Rev. D 83, 064003 (2011), URL https://link.aps.org/doi/10.1103/PhysRevD.83.064003.
  * Damour et al. (2009) T. Damour, B. R. Iyer, and A. Nagar, Phys. Rev. D 79, 064004 (2009), URL https://link.aps.org/doi/10.1103/PhysRevD.79.064004.
  * (82) Caltech-Cornell-CITA, _binary black hole simulation results_ , http://www.black-holes.org/waveforms.
  * D. Shoemaker (2010) (LIGO Scientific Collaboration) D. Shoemaker (LIGO Scientific Collaboration), _Advanced LIGO designed high power PSD_ (2010), URL https://dcc.ligo.org/cgi-bin/DocDB/ShowDocument?docid=2974.
  * Healy and Lousto (2022) J. Healy and C. O. Lousto, Phys. Rev. D 105, 124010 (2022), URL https://link.aps.org/doi/10.1103/PhysRevD.105.124010.
  * Arif Shaikh et al. (2023) M. Arif Shaikh, V. Varma, H. P. Pfeiffer, A. Ramos-Buades, and M. van de Meent, arXiv e-prints arXiv:2302.11257 (2023), eprint 2302.11257.
  * Boetzel et al. (2019b) Y. Boetzel, C. K. Mishra, G. Faye, A. Gopakumar, and B. R. Iyer, Phys. Rev. D 100, 044018 (2019b), URL https://link.aps.org/doi/10.1103/PhysRevD.100.044018.