<figure><img src="image_1.png" /></figure>

HAL
open science

# Constructing Gravitational Waves from Generic Spin-Precessing Compact Binary Inspirals

Katerina Chatziioannou, Antoine Klein, Nicolás Yunes, Neil Cornish

► To cite this version:

Katerina Chatziioannou, Antoine Klein, Nicolás Yunes, Neil Cornish. Constructing Gravitational Waves from Generic Spin-Precessing Compact Binary Inspirals. Physical Review D, 2017, 95 (10), pp.104004. 10.1103/PhysRevD.95.104004. hal-01554567

HAL Id: hal-01554567

https://hal.science/hal-01554567v1

Submitted on 4 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Constructing gravitational waves from generic spin-precessing compact
binary inspirals

Katerina Chatziioannou,¹,² Antoine Klein,³,⁴,⁵ Nicolás Yunes,¹ and Neil Cornish¹

¹eXtreme Gravity Institute, Department of Physics, Montana State University,
Bozeman, Montana 59717, USA

²Canadian Institute for Theoretical Astrophysics, 60 St. George Street,
University of Toronto, Toronto, ON M5S 3H8, Canada

³Department of Physics and Astronomy, The University of Mississippi University, Mississippi 38677, USA

⁴CENTRA, Departamento de Física, Instituto Superior Técnico, Universidade de Lisboa,
Avenida Rovisco Pais 1, 1049 Lisboa, Portugal

⁵CNRS, UMR 7095, Institut d'Astrophysique de Paris, 98 bis Bd Arago, 75014 Paris, France
(Received 11 March 2017; published 5 May 2017)

The coalescence of compact objects is one of the most promising sources, as well as the source of the first detections, of gravitational waves for ground-based interferometric detectors, such as advanced LIGO and Virgo. Generically, compact objects in binaries are expected to be spinning with spin angular momenta misaligned with the orbital angular momentum, causing the orbital plane to precess. This precession adds rich structure to the gravitational waves, introducing such complexity that an analytic closed-form description has been unavailable until now. We here construct the first closed-form frequency-domain gravitational waveforms that are valid for generic spin-precessing quasicircular compact binary inspirals. We first construct time-domain gravitational waves by solving the post-Newtonian precession equations of motion with radiation reaction through multiple scale analysis. We then Fourier transform these time-domain waveforms with the method of shifted uniform asymptotics to obtain closed-form expressions for frequency-domain waveforms. We study the accuracy of these analytic, frequency-domain waveforms relative to waveforms obtained by numerically evolving the post-Newtonian equations of motion and find that they are suitable for unbiased parameter estimation for 99.2%(94.6%) of the binary configurations we studied at a signal-to-noise ratio of 10(25). These new frequency-domain waveforms could be used for detection and parameter estimation studies due to their accuracy and low computational cost.

DOI: 10.1103/PhysRevD.95.104004

# I. INTRODUCTION

Spin is ubiquitous in Nature; under the influence of a generic perturbation any astrophysical system will rotate even if the initial state was perfectly spherically symmetric. When massive stars give birth to neutron stars (NSs) or black holes (BHs) through supernova explosions, the newly-born remnant typically spins rapidly even if the progenitor was spinning slowly. This is possibly due to asymmetries in the supernova explosion inducing a “kick” on the remnant causing it to rotate [1]. At a fundamental level, the physics at play here is the same as that which causes a soccer ball to spin after kicked.

The spin angular momenta of the components of a compact binary system will not necessarily be aligned with the orbital angular momentum. For example, consider an isolated binary system of two stars with spins aligned with the orbital angular momentum. The most massive star will first fill its Roche lobe and transfer mass to the companion before going supernova. The compact remnant (BH or NS) receives a kick spinning it up [1] and tilting the orbital plane, since the kick's direction is typically correlated with the spin of the exploding star [2,3]; the various angular momenta become misaligned [4,5]. Eventually, the

second star also fills its Roche lobe and the binary enters a common envelope phase. In this phase, the angular momenta could partially align through tidal effects. But, again, the phase ends with the star going supernova and endowing the binary with a kick that spins the remnant up and typically tilts the orbital plane, misaligning the angular momenta yet again [4,5].

Few mechanisms exist that could prevent misalignment or realign the spins of compact binary components with the orbital angular momentum. One possibility is if the supernova kicks are in the orbital plane, i.e. perpendicular to the spin angular momentum, such that the orbital plane is not tilted [4]. This possibility is remote, with models and data suggesting that the kick is actually aligned with the spin angular momentum [2,3,6]. A mechanism for realignment is through torques exerted by a circumbinary accretion disk on BH binaries [7]. Such disks, however, are only expected in galaxy mergers of supermassive BHs, whose gravitational waves (GWs) would be outside of the sensitivity band of the ground-based detectors advanced LIGO (aLIGO) and advanced Virgo (AdV). Besides the formation mechanism described above, dynamical formation channels are expected to result in binaries with arbitrarily
distributed spin directions [8]. Finally, we note that the binary might undergo spin-orbit resonances during its evolution [9]. These resonances, however, do not align the spins with the orbital angular moment, but rather maintain certain special (“resonant”) precessing configurations [5,9–12].

Compact binaries emitting GWs in the sensitivity band of aLIGO and AdV will thus have spins with arbitrary magnitudes and directions, though the spin of NSs will typically be smaller than that of BHs. Such GWs are very different from those emitted by nonspinning binaries or binaries with spins aligned with the orbital angular momentum for one main reason: *precession*. When spins are misaligned with the orbital angular momentum, all angular momenta precess about the total angular momentum, causing the orbital plane to precess too. This induces amplitude and phase modulations in the GWs, for example through the changing inclination angle of the system relative to the line of sight. Accurately modeling these modulations can be important for detection and parameter estimation of GW sources [13,14].

GWs with such rich and complex structure are double-edged swords: on the one hand, this structure can encode new information about the source and break degeneracies in parameter estimation; on the other hand, this intrinsic complexity comes at the cost of an increased difficulty to model these waves. First, the temporal evolution of the orbital phase depends on the angular momenta, which themselves satisfy certain precession equations, increasing the overall complexity of the differential system. Second, precession introduces *mathematical catastrophes* when computing the Fourier transform of the GWs: degenerate critical points in the orbital phase, where the first and second time derivatives vanish, violating the assumptions of the standard *stationary phase approximation*.<sup>1</sup> (SPA) and rendering it nonapplicable [17].

Accurate modeling of GWs is important for detection and crucial for parameter estimation with ground-based detectors, since the expected signals might be deeply buried in the detector noise. The optimal strategy for extracting known signals from noise is by fitting *waveform models* to the interferometric data and minimizing the residual. The efficiency of this method relies on the accuracy of the waveform model, with parameter estimation placing more stringent requirements on the accuracy of the models used. Inaccuracies in the models can lead to missed signals or systematic errors in the extracted parameters.

This has motivated the construction of waveform models for coalescing compact binaries. During the inspiral, the binary can be modeled with the *post-Newtonian* (PN) formalism, an expansion in small characteristic velocities

and weak gravitational fields [18]. When the binary components are not spinning or when their spin is aligned/anti-aligned with the orbital angular momentum, the equations of motion have been derived and solved up to 3.5PN order<sup>2</sup> including radiation reaction due to GW energy loss. When the spins are misaligned with the orbital angular momentum, the orbital equations of motion and the precession equations have been derived to 2.5PN order. In this case, a closed form solution has not been obtained due to the complexity of the differential system.

To this day, four main representations of GWs from spin-precessing compact binaries exist. The first representation is based on the fact that the precession equations admit a closed-form analytic solution when only one object is spinning [19–23]. The ensuing motion is *simple precession* and the resulting waveform is ideal for BHNS systems [24]. The second representation is based on the *effective-one-body* formulation of the general relativistic two-body problem [25–27]. The resulting waveform is ever improving through fits of its nonprecessing part to numerical relativity simulations [28–30] and describes the full coalescence, albeit at the expense of prohibitive computational cost. The third representation utilizes a coordinate frame in which precessional effects are minimized [31,32] to compute a simpler waveform [30] and map it back to the source frame [33,34]. This approach applies to the full binary coalescence and the waveform was found to be sufficiently good for detection, but could introduce biases in parameter estimation [23].<sup>3</sup>

The final representation of GWs from inspiraling spin-precessing systems was through a multiple scale analysis (MSA) [16], a well-known mathematical technique to solve differential systems that have distinct characteristic scales by expanding in the ratio of these scales. For the problem at hand, the orbital time scale is much shorter than the precession time scale, which in turn is much shorter than the radiation reaction time scale. This technique has already been applied successfully to nearly aligned [38] and slowly spinning [39] systems. The latter are accurate representations of NSNS inspirals, both for detection and parameter estimation [13,14,40]. In this paper we utilize two recent breakthroughs to construct waveforms for spin-precessing systems with arbitrary spin magnitudes and orientations with MSA.

The first breakthrough in the modeling of generic spin-precessing binaries was by Kesden *et al.* [10–12]. Neglecting radiation reaction, the authors found an *exact* solution to the precession equations that govern the evolution of the orbital and the spin angular momenta of

<footnote n="1"><sup>1</sup>The SPA is the leading-order term in the asymptotic expansion of a Fourier integral through the method of steepest descent [15,16].</footnote>

<footnote n="2"><sup>2</sup>A term is of APN order is proportional to $(u/c)^{2A}$ relative to its controlling factor, where *u* is some characteristic velocity and *c* the speed of light.</footnote>

<footnote n="3"><sup>3</sup>It is worth emphasizing, thought, that the first GW detection [35,36] did not suffer from such systematics due to its orientation and minimal precession [35,37].</footnote>
the binary. By identifying certain constants of the precessional motion, they were able to express all angular momenta as functions of the total spin magnitude, which satisfies an ordinary differential equation. We here solve this differential equation analytically and obtain an exact solution to the precession equations in the absence of radiation reaction. We then use MSA to introduce radiation reaction perturbatively as an expansion in the ratio of the precession to the radiation reaction time scale. With this at hand, we obtain time-domain waveforms in terms of the parameters of the system only.

The second breakthrough in the modeling of generic spin-precessing binaries was by Klein et al. [17], and tackles the failure of the SPA. The authors introduced the *shifted uniform asymptotics* (SUA) method where the waveform is decomposed into Bessel functions, the Fourier integral is evaluated term by term in the SPA, and then resumed using the exponential shift theorem. The result is a closed-form analytic expression for the gravitational wave in the frequency domain as a series of time-domain waveforms evaluated at shifted stationary times. Unlike previous approaches, both the time- and frequency-domain waveforms that we obtain are valid for arbitrary mass ratios, arbitrary spin magnitudes and arbitrary spin orientations. This waveform was first presented in Ref. [41], while this paper provides the details of its derivation.

Closed-form expressions for the waveforms have several advantages. From a theoretical standpoint, analytic solutions shed light on the physical processes at play, the structure of the resultant signal, and the transition through different resonant states. From a practical standpoint, analytic solutions are in general faster to evaluate, avoiding costly numerical integrations and discrete Fourier transforms. Estimating the computational gain from closed-form, analytic expressions relative to numerical ones is not straightforward since it depends heavily on the implementation. However, we estimate that in the restricted waveform case (when only one harmonic is used) the closed-form, analytic frequency-domain waveforms computed here can be an order of magnitude faster than the implementation of [17], and are at worst comparable.

The remainder of the paper provides the details of the waveform construction described above. Throughout, we use geometric units where $G = c = 1$ and use the following conventions:

(i) Vectors are written in boldface, with components $\mathbf{A} = [A_x, A_y, A_z]$ and magnitude $A$. Unit vectors are denoted with a hat, e.g. $\hat{\mathbf{A}}$.

(ii) The masses of the two binary components are $m_A$, with $A \in \{1, 2\}$, the total mass $M \equiv m_1 + m_2$ is set equal to 1, the mass ratio is $q \equiv m_2/m_1 < 1$, the symmetric mass ratio is $\eta \equiv m_1 m_2$ and the mass difference is $\delta m = m_1 - m_2$.

(iii) The Newtonian orbital angular momentum of the system is $\mathbf{L}$, the spin angular momentum of each

<table><caption>TABLE I. Parameters of the systems we use for comparisons of our analytic solution to the numerical solution to the PN precession equations. All parameters are defined at 50 Hz and in a frame were the orbital angular momentum is aligned with the z axis.</caption><thead><tr><th></th><th>NSNS</th><th>BHNS</th><th>BHBH</th><th>HSNSBH</th></tr></thead><tbody><tr><td>$m_1$</td><td>1.6 $M_\odot$</td><td>10 $M_\odot$</td><td>10 $M_\odot$</td><td>10 $M_\odot$</td></tr><tr><td>$m_2$</td><td>1.4 $M_\odot$</td><td>1.4 $M_\odot$</td><td>5 $M_\odot$</td><td>1.4 $M_\odot$</td></tr><tr><td>$\cos \theta_L$</td><td>1</td><td>1</td><td>1</td><td>1</td></tr><tr><td>$\phi_L$</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>$\cos \theta_1$</td><td>0.5</td><td>0.5</td><td>0.5</td><td>0.5</td></tr><tr><td>$\phi_1$</td><td>1.2</td><td>1.2</td><td>1.2</td><td>1.2</td></tr><tr><td>$\chi_1$</td><td>0.08</td><td>0.7</td><td>0.7</td><td>0.7</td></tr><tr><td>$\cos \theta_2$</td><td>0.7</td><td>0.7</td><td>0.7</td><td>0.7</td></tr><tr><td>$\phi_2$</td><td>2.5</td><td>2.5</td><td>2.5</td><td>2.5</td></tr><tr><td>$\chi_2$</td><td>0.1</td><td>0.1</td><td>0.6</td><td>0.6</td></tr></tbody></table>

body is $S_A$, and the total angular momentum is $\mathbf{J} = \mathbf{L} + \mathbf{S}_1 + \mathbf{S}_2$. The dimensionless spin parameter of each object is $\chi_A \equiv S_A/m_A^2$ with $A \in \{1, 2\}$.

(iv) The orbital angular frequency in a frame fixed to the orbital plane is $\omega$, while the PN expansion parameter we use is $v \equiv \omega^{1/3} = \eta L^{-1}$.

(v) We test our analytic solution by comparing it to the numerical solution to the precession equations for certain systems. We select a NSNS, a BHNS, a BHBH, and a highly spinning (HS)NSBH system with parameters given in Table I in a frame were the z axis is aligned with the orbital angular momentum. The angles $\theta_L$ and $\phi_L$ are the polar angles of $\mathbf{L}$, while $\theta_A$ and $\phi_A$ are the polar angles of $S_A$.

## II. SPIN AND ANGULAR MOMENTUM EVOLUTION

A quasicircular binary system consisting of generic spinning compact objects is subject to spin-orbit and spin-spin interactions that force all angular momenta to precess. Averaging over one orbit,<fnref n="4" /> the precession equations governing the conservative evolution of the orbital and spin angular momenta are [42–44]<fnref n="5" />

$$ \dot{\mathbf{L}} = \left\{ \left(2 + \frac{3}{2}q\right) - \frac{3}{2}\frac{v}{\eta} \left[ (\mathbf{S}_2 + q\mathbf{S}_1) \cdot \hat{\mathbf{L}} \right] \right\} v^6 (\mathbf{S}_1 \times \hat{\mathbf{L}}) \\ + \left\{ \left(2 + \frac{3}{2q}\right) - \frac{3}{2}\frac{v}{\eta} \left[ \left(\mathbf{S}_1 + \frac{1}{q}\mathbf{S}_2\right) \cdot \hat{\mathbf{L}} \right] \right\} v^6 (\mathbf{S}_2 \times \hat{\mathbf{L}}) \\ + \mathcal{O}(v^7), \qquad (1) $$

<footnote n="4">⁴Orbit-averaging should be well justified provided there is a clean separation between the orbital and the precessional time scales, as is the case in the early inspiral.</footnote>

<footnote n="5">⁵The precession equations used here are only strictly valid for BHs, as for NSs they acquire additional terms describing the quadrupole moment of the bodies [45]. However, the extra terms are degenerate with the spins [46] and difficult to measure.</footnote>
$$ \dot{\mathbf{S}}_1 = \left\{ \eta \left( 2 + \frac{3}{2}q \right) - \frac{3v}{2} \left[ (q\mathbf{S}_1 + \mathbf{S}_2) \cdot \hat{\mathbf{L}} \right] \right\} v^5 (\hat{\mathbf{L}} \times \mathbf{S}_1) + \frac{v^6}{2} \mathbf{S}_2 \times \mathbf{S}_1 + \mathcal{O}(v^7), \quad (2) $$

$$ \dot{\mathbf{S}}_2 = \left\{ \eta \left( 2 + \frac{3}{2q} \right) - \frac{3v}{2} \left[ \left( \frac{1}{q} \mathbf{S}_2 + \mathbf{S}_1 \right) \cdot \hat{\mathbf{L}} \right] \right\} v^5 (\hat{\mathbf{L}} \times \mathbf{S}_2) + \frac{v^6}{2} \mathbf{S}_1 \times \mathbf{S}_2 + \mathcal{O}(v^7). \quad (3) $$

Radiation reaction drives the evolution of the magnitude of the orbital angular momentum, leaving the magnitude of the spin angular momenta unaltered to our current knowledge of the PN expansion and ignoring all energy and angular momentum flux through BH horizons [47,48]. The magnitude $L$ is related to the evolution of the orbital frequency $\omega$, and the PN expansion parameter $v$, leading to

$$ \dot{v} = \frac{v^9}{3} \frac{1}{\sum_{n=0}^{7} [g_n + 3g_n' \ln(v)] v^n}. \quad (4) $$

The coefficients $\{g_n, g_n'\}$ are functions of the symmetric mass ratio and inner products of the angular momenta, given in Appendix A.

Equations (1)–(3) describe the conservative dynamics, while Eq. (4) describes the dissipative dynamics. The former models the spin-spin and spin-orbit interactions, that change only the direction of $L, S_1$ and $S_2$. We use only the leading PN order expressions in each interaction.<fnref n="6" /> We do not use higher PN order corrections because the spin-spin and spin-cubed terms have not been fully calculated for generic precessing orbits yet [51,53]. In principle, we could have included the spin-orbit corrections. However, as explained later, our solution makes use of a certain quantity [54] that is conserved by the leading-order in spin-orbit and spin-spin interactions precession equations. If we use partial precession equations (including spin-orbit but not spin-spin corrections) it is not clear if we can modify this quantity so that it remains conserved. Once the spin-spin and spin-cubed terms have been fully calculated we can revisit this problem.

The dissipative dynamics govern the GW frequency evolution by changing the magnitude of the Newtonian orbital angular momentum $L = \eta/v$. This equation is known to 2.5PN order in all spin interactions [53], 4PN in linear-in-spin terms [50,53,55–58] and 22PN order in the point particle limit, neglecting spins and BH absorption effects [59–62]. In our analysis we keep terms in Eq. (4) to 3.5PN order since this is the highest complete PN order, ignoring spin-spin terms. In this case, we can easily include partial PN terms in radiation reaction to make the evolution

more accurate. When the 3PN spin-spin term has been fully calculated for precessing orbits [51] we can include it in our model.

Conservative and dissipative equations evolve on distinct time scales. The former evolve on the precession time scale

$$ T_{\text{pr}} \equiv \frac{|\mathbf{S}_1|}{|\dot{\mathbf{S}}_1|} \sim v^{-5}, \quad (5) $$

while the later evolve on the radiation reaction time scale

$$ T_{\text{rr}} \equiv \frac{v}{\dot{v}} \sim v^{-8}. \quad (6) $$

The ratio $T_{\text{pr}}/T_{\text{rr}} \sim v^3$ is a small quantity in the inspiral and thus a natural expansion parameter.

Recently, Kesden et al. [10] found an exact solution to the precession equations [Eqs. (1)–(3)] ignoring radiation reaction [Eq. (4)]. This solution can be used to “precession-average” the full precession equations with radiation reaction (analogously to orbit-averaging). The final precession-averaged equations depend only on quantities that vary on the radiation reaction time scale, and can be numerically integrated with a larger step size [11].

Here we take a different approach. Rather than precession-averaging Eqs. (1)–(3) and numerically accounting for Eq. (4), we make explicit use of the fact that $T_{\text{pr}}/T_{\text{rr}} \sim v^3$ to solve the precession equations analytically. We use a perturbation theory technique known as multiple scale analysis (MSA) and treat radiation reaction as a slowly-evolving perturbation on top of precession. This approach allows us to find a solution to the full set of Eqs. (1)–(4) as an expansion in $T_{\text{pr}}/T_{\text{rr}}$.

### III. ANALYTIC SOLUTION TO THE PRECESSION EQUATIONS WITHOUT RADIATION REACTION

Ignoring radiation reaction, the precession equations can be solved analytically by making use of certain conserved quantities of the system. Below we review and complete the solution first presented in [10].

A precessing binary has a total of 9 degrees of freedom arising from the 3 components of 3 Newtonian vectors $(L, S_1, S_2)$. The precession equations lead to 7 conserved quantities, reducing the degrees of freedom to 2. Of the remaining degrees of freedom, one is associated with the choice of a coordinate system, while the other corresponds to a dynamical quantity that changes with time. This dynamical quantity is chosen to be the magnitude of the total spin angular momentum $S = |\mathbf{S}_1 + \mathbf{S}_2|$.

The conserved quantities are $\lambda \equiv (S_1, S_2, L, J, \hat{J}, \xi)$: the magnitudes of the spin angular momenta, the magnitude of the orbital angular momentum, the magnitude and direction of the total angular momentum, and the mass weighted effective spin [54]

<footnote n="6">⁶Spin-orbit corrections can be found in [49,50], spin-spin in [51], and spin-cubed in [52].</footnote>
$$ \xi \equiv (1 + q)S_1 \cdot \hat{L} + (1 + q^{-1})S_2 \cdot \hat{L}. \quad (7) $$

In the effective-one-body formalism, $\xi$ corresponds to the projection of the spin angular momentum of the body at the center of mass onto the orbital angular momentum. Once the system is allowed to evolve under radiation reaction, $S_1$, $S_2$, and $\xi$ are still conserved, while $L$, $J$ and $\hat{J}$ evolve on the radiation reaction time scale.

In the remainder of this section we use these 7 conserved quantities to geometrically solve for the 9 components of the angular momenta as a function of $S$ in a specific coordinate system. We then complete the solution for the angular momenta as a function of time by solving a differential equation to determine $S(t)$.

### A. Precession in a noninertial frame

The identification of $\hat{J}$ as a conserved quantity suggests a coordinate frame where $\hat{z} = \hat{J}$ (see Fig. 1). We further pick the $x$ and $y$ axes to be precessing around $\hat{z}$ (a noninertial frame), following the precession of the orbital angular momentum which is chosen to be in the $x-z$ plane, at an angle

$$ \cos \theta_L = \hat{J} \cdot \hat{L} = \frac{J^2 + L^2 - S^2}{2JL}, \quad (8) $$

from the $\hat{z}$ axis. This allows us to express $L$ as

$$ L(S; \lambda) = L[\sin \theta_L, 0, \cos \theta_L]. \quad (9) $$

The total spin angular momentum then is

$$ S(S; \lambda) = J - L = [-L \sin \theta_L, 0, J - L \cos \theta_L]. \quad (10) $$

In another frame with $\hat{z}' = \hat{S}, \hat{y}' = \hat{J}$, and $\hat{x}' = \hat{y}' \times \hat{z}'$, we define angles $(\theta', \phi')$ (see Fig. 1) such that

$$ S_1' = S_1[\sin \theta' \cos \phi', \sin \theta' \sin \phi', \cos \theta'.] \quad (11) $$

<figure><img src="image_3.png" /><figcaption>FIG. 1. Initial configuration of the angular momenta in a noninertial frame precessing around $\hat{z}$.</figcaption></figure>

Using the definition of $\xi$ given in Eq. (7) we get

$$ \cos \theta' = \hat{S}_1 \cdot \hat{S} = \frac{S_1^2 + S_2^2 - S_3^2}{2SS_1}, \quad (12) $$

$$ \cos \phi' = \frac{\{ (J^2 - L^2 - S^2)[S^2(1+q)^2 - (S_1^2 - S_2^2)(1-q^2)] - 4qS^2L\xi \}}{[(1-q^2)A_1A_2A_3A_4]}, \quad (13) $$

where

$$ A_1 = \sqrt{J^2 - (L - S)^2}, \quad (14) $$

$$ A_2 = \sqrt{(L + S)^2 - J^2}, \quad (15) $$

$$ A_3 = \sqrt{S^2 - (S_1 - S_2)^2}, \quad (16) $$

$$ A_4 = \sqrt{(S_1 + S_2)^2 - S^2}. \quad (17) $$

In the original unprimed system

$$ S_1(S; \lambda) = \mathbb{R}(\hat{y}, \theta_S)S_1', \quad (18) $$

where $\mathbb{R}(\hat{y}, \theta_S)$ is a rotation around $\hat{y}$ by an angle $\theta_S$ and

$$ \cos \theta_S = \hat{S} \cdot \hat{J} = \frac{J^2 + S^2 - L^2}{2JS}. \quad (19) $$

Once we have $S_1$ in the original unprimed system, then

$$ S_2(S; \lambda) = J - L - S_1. \quad (20) $$

Equations (9), (18), and (20) determine the angular momenta in a noninertial frame as a function of $S$ up to the sign of $\sin \phi'$ in Eq. (18), which we will tackle shortly.

At this point, the various orbital angular momenta have been written in a noninertial frame in terms of $S$ using purely geometrical arguments. The evolution equation of $S$ can be derived from Eqs. (1)–(3):

$$ \left(\frac{dS^2}{dt}\right)^2 = -A^2(S^6 + BS^4 + CS^2 + D), \quad (21) $$

where the coefficients $A, B, C, D$ depend only on quantities that change on the radiation reaction time scale. Their explicit form is given in Appendix B. The roots of the polynomial on the right-hand side of Eq. (21) have a simple interpretation. When $S^2$ is equal to one of the roots, its derivative is zero. Therefore, two of the roots are the maximum $S_+^2$ and the minimum $S_-^2$ of $S^2$. The third root $S_3^2$
does not correspond to any physically interesting scenario; in fact, it is negative for most systems.<fnref n="7" />

Making explicit use of the roots of the polynomial, we can rewrite Eq. (21) as

$$ \left(\frac{dS^2}{dt}\right)^2 = -A^2(S^2 - S_+^2)(S^2 - S_-^2)(S^2 - S_3^2). \quad (22) $$

The solution to this equation is

$$ S^2 = S_+^2 + (S_-^2 - S_+^2)\mathrm{sn}^2(\psi, m) \quad (23) $$

where sn is a Jacobi elliptic function (see Sec. 16 of [63] for a detailed introduction to the Jacobi elliptic functions, and [64] for a physics-oriented approach), $\psi$ is its phase, and $m \in [0, 1]$. When $m = 0$, sn reduces to a sine, while for $m = 1$ it gives a hyperbolic tangent. The period of $S^2$ is $2K(m)$, where $K(m)$ is the complete elliptic integral of the first kind. The phase and the parameter $m$ are given by

$$ \frac{d\psi}{dt} = \frac{A}{2} \sqrt{S_+^2 - S_3^2} \quad (24) $$

and

$$ m = \frac{S_+^2 - S_-^2}{S_+^2 - S_3^2}. \quad (25) $$

Clearly, this solution requires that $S_+^2 \neq S_3^2$, which is almost always the case because $S_+^2$ and $S_3^2$ are defined to be the largest and smallest roots respectively. The only possible case when $S_3^2 = S_+^2$ is when $S_+^2 = S_3^2$, but then $S^2$ is constant in the first place and there is no precession. The phase $\psi$ can be obtained by noticing that $\dot{\psi}$ is constant if we ignore radiation reaction, so that

$$ \psi = \frac{A}{2} \sqrt{S_+^2 - S_3^2} t. \quad (26) $$

The final ingredient we need in order to have a complete expression for all angular momenta as function of time in a non-inertial frame precessing around $\hat{z}$ is the sign of $\sin \phi'$. Equation (11) implies that

$$ \mathrm{sign}(\sin \phi') = \mathrm{sign}(\mathbf{S}_1 \cdot \mathbf{y}'), \quad (27) $$

which after some algebra can be shown to be equivalent to

$$ \begin{aligned} \mathrm{sign}(\sin \phi') &= \mathrm{sign}[(\hat{\mathbf{L}} \times \mathbf{S}_1) \cdot \mathbf{S}_2] = \mathrm{sign}\left(-\frac{dS^2}{dt}\right) \\ &= \mathrm{sign}[\mathrm{sn}(\psi, m)\mathrm{cn}(\psi, m)], \end{aligned} \quad (28) $$

where $\mathrm{cn}(\psi, m)$ is another Jacobi elliptic function and in the last equality we have used Eq. (23).

## B. Precession in an inertial frame

All angular momenta so far have been expressed in a noninertial frame that precesses around $\hat{\mathbf{J}}$. An Euler rotation of $\mathbf{L}$, $\mathbf{S}_1$, and $\mathbf{S}_2$ around $\hat{z}$ by some angle $\phi_z$ and substitution into the precession equations yields [10]

$$ \frac{d\phi_z}{dt} = \Omega_z = \frac{J}{2} v^6 \left\{ 1 + \frac{3}{2\eta}(1 - \xi v) - \frac{3(1+q)}{2qA_1^2A_2^2}(1 - \xi v)[4(1-q)L^2(S_1^2 - S_2^2) - (1+q)(J^2 - L^2 - S^2)(J^2 - L^2 - S^2 - 4\eta L \xi)] \right\}. \quad (29) $$

The precession angle $\phi_z$ changes on the precession time scale through $S$ and on the radiation-reaction time scale through $J$ and $L$. We recast it in the form

$$ \frac{\dot{\phi}_z}{J} = a + \frac{c_0 + c_2 \mathrm{sn}^2(\psi, m) + c_4 \mathrm{sn}^4(\psi, m)}{d_0 + d_2 \mathrm{sn}^2(\psi, m) + d_4 \mathrm{sn}^4(\psi, m)}, \quad (30) $$

where $a$, the $d_i$'s and the $c_i$'s are quantities that evolve on the radiation reaction time scale only. Their explicit form is given in Appendix B. Now $\dot{\phi}_z$ can be integrated exactly in the absence of radiation reaction to give

$$ \frac{\phi_z}{J} = A_{\phi} \frac{\psi}{\dot{\psi}} + iB_{\phi} \frac{F[\mathrm{isinh}^{-1}(\mathrm{sc}(\psi, m)), 1-m]}{\dot{\psi}} + iC_{\phi} \frac{\Pi[n_c, \mathrm{isinh}^{-1}(\mathrm{sc}(\psi, m)), 1-m]}{\dot{\psi}} + iD_{\phi} \frac{\Pi[n_d, \mathrm{isinh}^{-1}(\mathrm{sc}(\psi, m)), 1-m]}{\dot{\psi}}, \quad (31) $$

where $\dot{\psi}$ is given by Eq. (24), $F$ is the elliptic integral of the first kind, $\Pi$ is the elliptic integral of the third kind, and $\mathrm{sc}$ is a Jacobi elliptic function. The quantities $A_{\phi}$, $B_{\phi}$, $C_{\phi}$, $D_{\phi}$, $n_c$, $n_d$ are functions of $\{a, c_i, d_i\}$, and they are constant in the absence of radiation reaction. They are given in Appendix B.

This concludes the solution to the precession equations in the absence of radiation reaction in a frame where $\hat{\mathbf{J}} = \hat{\mathbf{z}}$. In summary, at some initial time:

(i) The orbital angular momentum $\mathbf{L}$ is given by Eq. (9), which depends on the angle $\theta_L$ given in Eq. (8). The

<footnote n="7">⁷In the most generic case, a third order polynomial with real coefficients can have complex roots. However, we argue that this is an unphysical scenario. Unless two of the roots are real, $S^2$ will increase or decrease with no bound. If two roots of a third order polynomial with real coefficients are real, then the third root must be real too.</footnote>
latter depends on $S$, which varies on the precession timescale as described in Eq. (23);

(ii) The spin angular momentum of the heavier body $S_1$ is given in Eq. (18), which depends on the angle $\theta_S$ given in Eq. (19) as well as on $S'_1$ given in Eq. (11) in terms of the angles $(\theta', \phi')$ of Eqs. (12) and (13). All of these depend on $S$, which again is described by Eq. (23);

(iii) The spin angular momentum of the lighter body $S_2$ is given by Eq. (20), which depends on $L$ and $S_1$ described above.

The full precessional motion of these angular momenta in an inertial frame is obtained by rotating them around $\hat{z}$ by $\phi_z$, given in Eq. (31).

## IV. ADDITION OF RADIATION REACTION

The exact solution to the precession equations obtained in the previous section is valid only in the absence of radiation reaction. The problem of including radiation reaction admits a perturbative solution owing to its two distinct time scales: radiation reaction unfolds on a much longer time scale than precession. This natural separation of timescales allows us to treat radiation reaction as a slow perturbation of the more rapid precession, a technique formally known as multiple scale analysis [16].

In MSA, every quantity is expanded in the ratio of the two distinct timescales. In our case, we expand in the ratio of the precessional time scale $T_{pr}$ to the radiation reaction time scale $T_{rr}$; radiation reaction is a 1.5PN effect on top of precession. This is not the first application of MSA to the precession problem. In fact, the precession equations we started with are orbit-averaged, which would be the first term in an MSA expansion about the ratio of the fast orbital time scale to the precession time scale.

### A. Choice of an inertial frame

The precession solution of Sec. III was built around the assumption that $\hat{J}$ is conserved and aligned with $\hat{z}$. Our first task when adding radiation reaction is to check whether this remains true. If it does, then the functional form of Eqs. (9), (18), and (20) holds, since they were derived solely on geometrical arguments.

Radiation reaction does not strictly conserve the direction of the total angular momentum. However, it has been argued [19] that in the context of simple precession ($S_2 = 0$) the variation of $\hat{J}$ in a precession cycle averages out. Here we show that this is approximately true for generic precession as well [65].

Equations (1)–(3) imply

$$ \dot{J} = \dot{L} \hat{L}, \qquad (32) $$

and after some algebra we can rewrite this as

$$ \dot{\hat{J}} = \frac{\dot{L}}{JL} \hat{L} - \frac{\dot{J}}{J^2} \hat{J}. \qquad (33) $$

Averaging over $\phi_z$ we find

$$ \langle \dot{J}_x \rangle_{\phi_z} = \left\langle \frac{\dot{L}}{J} \sin \theta_L \cos \phi_z \right\rangle_{\phi_z}, \qquad (34) $$

$$ \langle \dot{J}_y \rangle_{\phi_z} = \left\langle \frac{\dot{L}}{J} \sin \theta_L \sin \phi_z \right\rangle_{\phi_z}, \qquad (35) $$

$$ \langle \dot{J}_z \rangle_{\phi_z} = 0. \qquad (36) $$

This averaging induces an error in $\dot{\hat{J}}$ that is $O(T_{pr}/T_{rr})$, or $\dot{\hat{J}} - \langle \dot{\hat{J}} \rangle_{\phi_z} \sim v^3$. At this order, we can treat $\dot{L}$ as a constant, since the spin couplings in Eq. (4) first enter at $O(v^3)$. They are therefore of the same order as the averaging error, and can be neglected.

Working to this order we have

$$ \langle \dot{J}_{x,y} \rangle_{\phi_z} \sim \langle \sin \theta_L \cos \phi_z \rangle_{\phi_z} \\ \sim \left\langle \sqrt{1 - \left( \frac{J^2 + L^2 - S^2}{2JL} \right)^2} \cos \phi_z \right\rangle_{\phi_z}. \qquad (37) $$

Since $L \sim O(v^{-1})$, $J \sim O(v^{-1})$, and $S \sim O(v^0)$, a PN expansion yields schematically

$$ \langle \dot{J}_{x,y} \rangle_{\phi_z} \sim \langle \cos \phi_z \rangle_{\phi_z} + v^2 \langle S^2 \cos \phi_z \rangle_{\phi_z} + O(v^4). \qquad (38) $$

The first term vanishes, while the second is of higher PN order and we neglect it. This situation is different from simple precession. In the latter the averaging out of $\dot{\hat{J}}$ is exact, while here it requires a PN expansion. We therefore expect this result to become less and less accurate as the binary approaches merger.

The above calculation implies that $\langle \dot{\hat{J}} \rangle_{\phi_z} = 0$; radiation reaction changes the magnitude of $\hat{J}$ while leaving its direction approximately constant. The components $J_x$ and $J_y$ are expected to oscillate with an amplitude much smaller than $J_z$ without exhibiting any secular growth. Figure 2 tests the validity of this statement. We select 4 systems with typical parameters as expected for NSNS, BHNS, BHNS, and HSNSBH binaries (see Table I) and plot the components of $\hat{J}$ obtained by numerically solving Eqs. (1)–(4) as a function of the GW frequency $f$. In all cases $J_x$ and $J_y$ are at least 2 orders of magnitude smaller than $J_z$ and oscillate around 0, with no signs of secular growth.

Based on this result we can build a solution to the precession equations including radiation reaction in the inertial frame introduced in Sec. III. That is, we neglect any variation in the direction of $\hat{J}$ and align it with $\hat{z}$. This
<figure><img src="image_4.png" /><figcaption>FIG. 2. (Top Panel) Comparison between the numerical PN and the analytic components of the total angular momentum as a function of the GW frequency for the NSNS (Top Left), the BHNS (Top Right), the BHBH (Bottom Left), and the HSNSBH (Bottom Right) system of Table I. (Bottom Panel) Fractional error between the magnitude of the total angular momentum obtained numerically and analytically.</figcaption></figure>

choice of frame automatically means that the functional form of Eqs. (9), (18), and (20) for the orbital and spin angular momenta respectively is still valid, since they were derived on purely geometric arguments. On the contrary, any quantity that was derived based on Eqs. (1)–(3) needs to be revisited and recalculated by taking Eq. (4) into account. This involves the remaining 5 conserved quantities of precession ($S_1, S_2, L, J, \xi$), Eq. (22) for the magnitude of the total spin angular momentum, and Eq. (30) for the precession angle.

## B. Constants of the precessional motion

In principle, the constants of the precessional motion need not remain constant when radiation reaction is invoked. The magnitudes of the two spin angular momenta $S_1$ and $S_2$, and the mass weighted effective spin $\xi$ remain constant under radiation reaction to the PN order we work here and ignoring horizon absorption. The magnitude of the

orbital angular momentum $L$ is updated by definition through $L = \eta/v$. The magnitude of the total angular momentum $J$ depends on $L$ and also changes under radiation reaction. The evolution equation for $J$ averaged over one period of $S(t)$ is [10]

$$ \left\langle \frac{dJ}{dL} \right\rangle_{\text{pr}} = \frac{J^2 + L^2 - \langle S^2 \rangle_{\text{pr}}}{2JL} \qquad (39) $$

This can be integrated exactly to yield

$$ J^2 = L^2 + \frac{2c_1}{\eta} L - L \int \frac{\langle S^2 \rangle_{\text{pr}}}{L^2} dL, \qquad (40) $$

where here and in what follows $J$ is approximated by its precession average, and $c_1$ is an integration constant. As we will show below, $\langle S^2 \rangle_{\text{pr}}$ is constant when ignoring
high-order PN effects, and the integral of Eq. (40) can be
calculated to give

$$
J^2 = L^2 + \frac{2c_1}{v} + \langle S^2 \rangle_{\text{pr}} + \mathcal{O}(v). \quad (41)
$$

The quantity 〈S²〉pr can be computed from Eq. (23):

$$
S_{\text{av}}^2 \equiv \langle S^2 \rangle_{\text{pr}} = \frac{1}{m} \left[ (m-1)S_+^2 + S_-^2 + \frac{E(m)}{K(m)} (S_+^2 - S_-^2) \right], \quad (42)
$$

where $K(m)$ and $E(m)$ are the complete elliptic integrals of
the first and second kind respectively. PN expanding $S_+^2$
and $S_-^2$ around their initial value we find

$$
S_{\pm}^{2} = S_{\pm,0}^{2} + \mathcal{O}(v), \quad S_{3}^{2} = \mathcal{O}(v^{-2}), \quad (43)
$$

which together with Eq. (25) yields

$$
m = \mathcal{O}(v^2),
\quad (44)
$$

and

$$
S_{\text{av}}^2 = \frac{1}{2} (S_{+,0}^2 + S_{-,0}^2) + \mathcal{O}(v^2). \quad (45)
$$

In the above expressions $S_{\pm,0}^2$ are the roots computed from
the initial conditions.

Combining the result for $J$ obtained here and Sec. IVA
where we justified keeping $\hat{J}$ aligned with $\hat{z}$, our analytic
approximation for the total angular momentum is

$$
\mathbf{J} = [0, 0, J]. \tag{46}
$$

To verify that this approximate **J** stays close to the
numerical PN solution we plot it in Fig. 2 as a function
of the GW frequency for our three study systems. The
analytic *J*<sub>x</sub> and *J*<sub>y</sub> are identically zero, so we omit them.
The bottom panel shows the fractional error in the
magnitude of the total angular momentum when approxi-
mated by Eq. (41). The maximum discrepancy in the
magnitude *J* is of *O*(10<sup>-2</sup>) in those particular examples
indicating both that Eq. (41) is accurate and that setting *J*<sub>x</sub>
and *J*<sub>y</sub> equal to zero is justified.

C. Magnitude of the total spin angular momentum

Once radiation reaction is included, Eq. (22) for the
magnitude of the total spin angular momentum needs to be
solved with MSA. We first explicitly separate the time
scales by writing S²(t) = S²(tpr, tt), where tpr denotes
variation on the precession time scale, while tt = etpr

denotes variations on the radiation reaction time scale, with
$\epsilon$ a bookkeeping parameter.

Expanding $S^2$ as

$$
S^2(t_{\text{pr}}, t_{\text{rr}}) = \sum_{n \ge 0} \epsilon^n S_n^2(t_{\text{pr}}, t_{\text{rr}}) \quad (47)
$$

and substituting this expression into Eq. (22), at leading
order in $\epsilon$, we recover Eq. (22) for $S_0^2(t_{\text{pr}}, t_{\text{rr}})$ with the time
derivative taken on the precession time scale $t_{\text{pr}}$:

$$
\left(\frac{\partial S_0^2}{\partial t_{\text{pr}}}\right)^2 = -A^2(t_{\text{rr}})[S_0^2(t_{\text{pr}}, t_{\text{rr}}) - S_+^2(t_{\text{rr}})] \\
\times [S_0^2(t_{\text{pr}}, t_{\text{rr}}) - S_-^2(t_{\text{rr}})][S_0^2(t_{\text{pr}}, t_{\text{rr}}) - S_3^2(t_{\text{rr}})]. \tag{48}
$$

The solution to this differential equation is similar to
Eq. (23), except that quantities that were previously
constant are now promoted to functions of $t_{\text{rr}}$:

$$
S_0^2 = S_+^2(t_{\text{rr}}) + [S_-^2(t_{\text{rr}}) - S_+^2(t_{\text{rr}})] \text{sn}[\psi(t_{\text{pr}}, t_{\text{rr}}), m(t_{\text{rr}})], \quad (49)
$$

where $S_+^2(t_{\text{rr}})$, $S_-^2(t_{\text{rr}})$, and $m(t_{\text{rr}})$ now depend on time
through $L(t_{\text{rr}})$ and $J(t_{\text{rr}})$.

The angle $\psi(t_{\text{pr}}, t_{\text{rr}})$ satisfies

$$
\frac{d\psi}{dt} = \frac{A(t_{\text{rr}})}{2} \sqrt{S_{+}^{2}(t_{\text{rr}}) - S_{3}^{2}(t_{\text{rr}})}. \quad (50)
$$

where we keep terms of $\mathcal{O}(\epsilon)$ by taking the derivative with
respect to $t$ rather than $t_{\text{pr}}$. We can integrate this equation
using a PN integration, i.e. expanding it in powers of $v$ and
integrating term by term. The result is

$$
\psi = \psi_0 - \frac{3g_0}{4} \delta m v^{-3} (1 + \psi_1 v + \psi_2 v^2), \quad (51)
$$

where $\psi_0$ is an integration constant, and the constants $\psi_1$, $\psi_2$ are given in Appendix C. We find that expanding Eq. (51) to relative 1PN order suffices.

We test this solution for *S* in Fig. 3 by plotting the
numerical PN, analytic, and hybrid magnitude of the total
spin angular momentum *S* as a function of the GW
frequency for the 4 systems we study. The hybrid *S* is
obtained through Eq. (49) but with a numerical solution to
Eq. (50). For all systems, the amplitude of *S* shows
excellent agreement with the numerical PN results, which
is controlled by the roots *S*<sub>+</sub><sup>2</sup> and *S*<sub>-</sub><sup>2</sup>. For the NSNS and
BHBH systems, the analytic phase *ψ* also shows very good
agreement with the numerical PN result, although the
dephasing for the BHNS and HSBHNS systems is about
2 cycles. However, both systems are dominated by the spin
of the BH, making the motion close to that of simple
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>FIG. 3. Comparison between the numerical PN (black solid), the analytic (red dashed), and the hybrid (blue dot-dashed) magnitude of the total spin angular momentum as a function of the GW frequency for the NSNS (Top Left), the BHNS (Top Right), the BHBH (Bottom Left), and the HSNSBH (Bottom Right) system of Table I.</figcaption></figure>

precession; the variation in S is very small as demonstrated by the scale of the y axis of the right panels of Fig. 3 and this dephasing should not affect the emitted waveform considerably.

On the other hand, the phase of the hybrid S is always in excellent agreement with the numerical solution, indicating that if we do indeed need an improved solution in the future<fnref n="8" /> we can obtain it by carrying out the expansion of Eq. (51) to higher order.

### D. Precession angle

The final quantity that needs to be recalculated to account for radiation reaction is the precession angle. Its derivative, given in Eq. (30), depends both on the precession and the radiation reaction time scale, so it requires a MSA treatment.

We write

$$ \frac{d\phi_z}{dt} = \Omega_z[S(t), L(t), J(t)] = \Omega_z[S(t_{pr}, t_{rr}), L(t_{rr}), J(t_{rr})], \quad (52) $$

and expand the precession angle as

$$ \phi_z(t_{pr}, t_{rr}) = \epsilon^{-1} \phi_{z,-1}(t_{pr}, t_{rr}) + \phi_{z,0}(t_{pr}, t_{rr}) + \mathcal{O}(\epsilon). \quad (53) $$

The reason $\phi_z$ includes a term of $\mathcal{O}(\epsilon^{-1})$ is because the binary precesses even in the absence of radiation reaction.

Solving Eq. (52) order by order in $\epsilon$, we find to $\mathcal{O}(\epsilon^{-1})$

$$ \frac{1}{\epsilon} \frac{\partial \phi_{z,-1}}{\partial t_{pr}} = 0, \quad (54) $$

which means $\phi_{z,-1} = \phi_{z,-1}(t_{rr})$. To next order, we find

<footnote n="8">⁸For example, if and when LIGO’s sensitivity increases, so will its requirement for more accurate waveforms.</footnote>
$$ \frac{\partial \phi_{z,-1}}{\partial t_{rr}} + \frac{\partial \phi_{z,0}}{\partial t_{pr}} = \Omega_z(t_{pr}, t_{rr}), \qquad (55) $$

and averaging over $t_{pr}$ we find

$$ \frac{d\phi_{z,-1}}{dt_{rr}} = \langle \Omega_z \rangle_{\text{pr}}(t_{rr}), \qquad (56) $$

where we set $\langle \partial \phi_{z,0} / \partial t_{pr} \rangle_{\text{pr}} = 0$ to cancel secular terms. Equation (56) can be solved with a PN integration. Going back to Eq. (55) we get

$$ \frac{\partial \phi_{z,0}}{\partial t_{pr}} = \Omega_z(t_{pr}, t_{rr}) - \langle \Omega_z \rangle_{\text{pr}}(t_{rr}). \qquad (57) $$

Integrating the first term on the right-hand side of Eq. (57) we recover Eq. (31) for $\phi_z$ in the absence of radiation reaction. Integrating the second term is straightforward.

The full solution for $\phi_z$ is then

$$ \phi_z = \phi_{z,-1} + \phi_{z,0} + \mathcal{O}(\epsilon), \qquad (58) $$

where

$$ \phi_{z,-1} = \int \langle \Omega_z \rangle_{\text{pr}}(t_{rr}) dt_{rr}, \qquad (59) $$

$$ \phi_{z,0} = \int \Omega_z(t_{pr}, t_{rr}) dt_{pr} - \int \langle \Omega_z \rangle_{\text{pr}}(t_{rr}) dt_{pr}. \qquad (60) $$

The meaning of each term in the MSA expansion is clear. The first term $\phi_{z,-1}$ is averaged over the fast (relative to radiation reaction) precession time scale, and then integrated over radiation reaction. The next term $\phi_{z,0}$ is a first order correction to this precession averaging.

## 1. Leading order MSA

The leading order MSA term is defined in Eq. (56) which to first order in $\epsilon$ is equivalent to

$$ \left\langle \frac{d\phi_z}{dt} \right\rangle_{\text{pr}} = \langle \Omega_z \rangle_{\text{pr}}. \qquad (61) $$

The average of $\Omega_z$ can be obtained by taking the difference between Eq. (31) evaluated at $\psi = 0$ and at $\psi = 2K(m)$, where recall that $K(m)$ is the complete elliptic integral of the first kind. However, for reasons explained in Appendix E, we prefer to use Eq. (30) and find an alternative way of calculating $\langle \dot{\phi}_z \rangle_{\text{pr}}$. We write

$$ \begin{aligned} \frac{\dot{\phi}_z}{J} - a &\equiv \dot{\phi}_z^{\text{red}} = \frac{c_0 + c_2 \text{sn}^2(\psi, m) + c_4 \text{sn}^4(\psi, m)}{d_0 + d_2 \text{sn}^2(\psi, m) + d_4 \text{sn}^4(\psi, m)} \Rightarrow \\ & \quad [d_0 + d_2 \text{sn}^2(\psi, m) + d_4 \text{sn}^4(\psi, m)] \dot{\phi}_z^{\text{red}} \\ &= c_0 + c_2 \text{sn}^2(\psi, m) + c_4 \text{sn}^4(\psi, m) \Rightarrow \\ & \quad d_0 \langle \dot{\phi}_z^{\text{red}} \rangle_{\text{pr}} + d_2 \langle \text{sn}^2(\psi, m) \dot{\phi}_z^{\text{red}} \rangle_{\text{pr}} + d_4 \langle \text{sn}^4(\psi, m) \dot{\phi}_z^{\text{red}} \rangle_{\text{pr}} \\ &= c_0 + c_2 \langle \text{sn}^2(\psi, m) \rangle_{\text{pr}} + c_4 \langle \text{sn}^2(\psi, m) \rangle_{\text{pr}}, \end{aligned} $$

where on the third line we average over precession. Unfortunately, no closed form expressions exist for $\langle \text{sn}^2(\psi, m) \rangle_{\text{pr}}$ and $\langle \text{sn}^4(\psi, m) \rangle_{\text{pr}}$ for arbitrary $m$. We can, however, calculate these averages as an expansion in $m \ll 1$ since, as already discussed, $m \sim \mathcal{O}(v^2)$. We could in principle retain high order in $m$ terms in this expansion, but in practice we find that working to leading order in $m$ suffices. Expanding the above expression to leading order in $m \ll 1$, we find

$$ \begin{aligned} & d_0 \langle \dot{\phi}_z^{\text{red}} \rangle_{\text{pr}} + d_2 \langle \dot{\phi}_z^{\text{red}} \sin^2 \psi \rangle_{\text{pr}} + d_4 \langle \dot{\phi}_z^{\text{red}} \sin^4 \psi \rangle_{\text{pr}} \\ &= c_0 + \frac{1}{2}c_2 + \frac{3}{8}c_4 \Rightarrow \\ & d_0 \langle \dot{\phi}_z^{\text{red}} \rangle_{\text{pr}} + d_2 D_2 \langle \dot{\phi}_z^{\text{red}} \rangle_{\text{pr}} + d_4 D_4 \langle \dot{\phi}_z^{\text{red}} \rangle_{\text{pr}} \\ &= c_0 + \frac{1}{2}c_2 + \frac{3}{8}c_4 \Rightarrow \\ & \langle \dot{\phi}_z \rangle_{\text{pr}} = J \left( a + \frac{c_0 + \frac{1}{2}c_2 + \frac{3}{8}c_4}{d_0 + d_2 D_2 + d_4 D_4} \right), \end{aligned} \qquad (62) $$

where we have defined

$$ D_2 \equiv \frac{\langle \dot{\phi}_z^{\text{red}} \sin^2 \psi \rangle_{\text{pr}}}{\langle \dot{\phi}_z^{\text{red}} \rangle_{\text{pr}}} = \frac{\frac{c_0+c_2 \sin^2 \psi+c_4 \sin^4 \psi}{d_0+d_2 \sin^2 \psi+d_4 \sin^4 \psi} \sin^2 \psi}{\langle c_0+c_2 \sin^2 \psi+c_4 \sin^4 \psi \rangle_{\text{pr}}}, \qquad (63) $$

$$ D_4 \equiv \frac{\langle \dot{\phi}_z^{\text{red}} \sin^4 \psi \rangle_{\text{pr}}}{\langle \dot{\phi}_z^{\text{red}} \rangle_{\text{pr}}} = \frac{\frac{c_0+c_2 \sin^2 \psi+c_4 \sin^4 \psi}{d_0+d_2 \sin^2 \psi+d_4 \sin^4 \psi} \sin^4 \psi}{\langle c_0+c_2 \sin^2 \psi+c_4 \sin^4 \psi \rangle_{\text{pr}}}. \qquad (64) $$

The quantities $D_2$ and $D_4$ are functions of $v$ and can be calculated exactly. For reasons explained in Appendix E we do not wish to use these full expressions, but rather we keep the quantities $D_2$ and $D_4$ constant and set them equal to their leading PN order expressions.

We can now integrate the right-hand side of Eq. (62) by first PN expanding it. However, we find it more convenient to factor $J$ out of $\langle \Omega_z \rangle_{\text{pr}}$ and PN expand the remaining terms. We do so to avoid artificial divergences in the small mass ratio limit arising from expanding around essentially $\eta/v$; see Appendix E. We, then, have to perform an integral of the form
$$ \phi_{z,-1} = \int \frac{J}{\xi^3} \sum_{n=0}^{5} \langle \Omega_z \rangle^{(n)} v^n d\xi, \qquad (65) $$

where the coefficients $\langle \Omega_z \rangle^{(n)}$ are given in Appendix D. This integral can be directly calculated to give

$$ \phi_{z,-1} = \sum_{n=0}^{5} \langle \Omega_z \rangle^{(n)} \phi_z^{(n)} + \phi_{z,-1}^{0}, \qquad (66) $$

where $\phi_z^{(n)}$ are functions given in Appendix D and $\phi_{z,-1}^0$ is an integration constant.

## 2. Correction to MSA

The first-order correction to MSA is given in Eq. (60). The solution to the first integral is Eq. (31) where we set $m = 0$. The second integral is trivial since $\langle \Omega_z \rangle_{pr}$ does not depend on the precession time scale $t_{pr}$, and the result is

$\langle \Omega_z \rangle_{pr} t_{pr}$. In that expression, we choose for convenience to substitute $t_{pr} = \psi/\dot{\psi}$.

Collecting all the elements together, the correction to the precession phase is given by

$$ \phi_{z,0} = \frac{C_{\phi} \sqrt{n_c}}{\dot{\psi} (n_c - 1)} \arctan \left[ \frac{(1 - \sqrt{n_c}) \tan \psi}{1 + \sqrt{n_c} \tan^2 \psi} \right] + \frac{D_{\phi} \sqrt{n_d}}{\dot{\psi} (n_d - 1)} \arctan \left[ \frac{(1 - \sqrt{n_d}) \tan \psi}{1 + \sqrt{n_d} \tan^2 \psi} \right], \qquad (67) $$

where $\dot{\psi}$ is given in Eq. (24), $\psi$ is given in Eq. (51) and $C_{\phi}$, $D_{\phi}$, $n_c$ and $n_d$ are functions of $v$ given in Appendix B.

## 3. Comparisons

In Fig. 4 we plot the numerical PN and analytic solutions for $\phi_z$ with and without the MSA corrections. The small oscillations of the numerical PN phase are reproduced by the analytic phase with MSA corrections. These oscillations are more pronounced for the NSNS and BHBH

<figure><img src="image_7.png" /><figcaption>FIG. 4. (Top Panel) Comparison between the numerical PN and the analytic precession phase as a function of the GW frequency for the NSNS (Top Left), the BHNS (Top Right), the BHBH (Bottom Left), and the HSNSBH (Bottom Right) system of Table I. (Bottom Panel) Error in $\phi_z$ with and without the MSA corrections.</figcaption></figure>
systems where both spins contribute significantly to the dynamics. The bottom panel shows the error in the precession phase with and without MSA corrections.

V. BUILDING THE WAVEFORM

Using the solution for the angular momenta described above, we calculate an analytic time-domain waveform for generic precessing binaries. The gravitational wave signal emitted by a precessing binary system as observed in an interferometric detector is [18,19,24,31,66]:

$$h(t) = F_{+}h_{+} + F_{\times}h_{\times}, \qquad (68)$$

where

$$F_{+} = \frac{1}{2}(1 + \cos^2\theta'_{N}) \cos 2\phi'_{N} \cos 2\psi_{p} - \cos\theta'_{N} \sin 2\phi'_{N} \sin 2\psi_{p}, \qquad (69)$$

$$F_{\times} = \frac{1}{2}(1 + \cos^2\theta'_{N}) \cos 2\phi'_{N} \sin 2\psi_{p} + \cos\theta'_{N} \sin 2\phi'_{N} \cos 2\psi_{p}, \qquad (70)$$

are the antenna pattern functions, $h_{+,×}$ are the GW polarization states, $(\theta'_{N}, \phi'_{N})$ are the polar angles of $\hat{N}$ in a frame tied to the arms of the detector with $\hat{z}'$ the normal to the detector plane, and $\psi_{p}$ is given by

$$\psi_p = \arctan \left[ \frac{(P_N \hat{J}) \cdot \hat{z}'}{(\hat{N} \times \hat{J}) \cdot \hat{z}'} \right], \qquad (71)$$

where $P_N$ acts as a projection along $\hat{N}$.

The polarization states can be decomposed into a spin-weighted spherical harmonic basis [24,31,66]

$$h_{+} - i h_{\times} = \sum_{l \ge 2} \sum_{m = -l}^{l} H^{lm}(\theta_s, \phi_s) e^{-im\Phi}, \qquad (72)$$

where

$$\Phi = \phi_{\text{orb}} - 3v^3(2 - \eta v^2) \ln v, \qquad (73)$$

and $(\theta_s, \phi_s)$ are the spherical angles of $\hat{N}$ in a frame where $\hat{J}$ is aligned with the z-axis, $\phi_{\text{orb}}$ is the orbital phase, and

$$H^{lm} = h^{lm} \sum_{m'=-l}^{l} D_{m',m}^{l}(\phi_z, \theta_L, \zeta)_{-2} Y_{lm'}(\theta_s, \phi_s), \qquad (74)$$

where $s Y_{lm}$ are the spin-weighted spherical harmonics, the amplitudes $h^{lm}$ are in [18], $D_{m,m'}^{l}$ are the Wigner D-matrices, the angles $\theta_L$ and $\phi_z$ are the spherical angles of $\hat{L}$ in the frame where $\hat{J}$ is aligned with the z-axis, and $\zeta$

satisfies $\dot{\zeta} = \dot{\phi}_z \cos \theta_L$. In order to solve for $\zeta$ we can employ the same techniques as for $\phi_z$, namely MSA. An explicit expression for $\zeta$ is given in Appendix F.

The above prescribe a waveform $h(t)$ in the time domain. To compute its Fourier transform, we use the shifted uniform asymptotics method of [17]

$$\tilde{h}(f) = \sqrt{2\pi} \sum_{m \ge 1} T_m e^{i(2\pi f t_m - m\Phi - \pi/4)} \\ \times \sum_{l \ge 2} \sum_{k=-k_{\text{max}}}^{k_{\text{max}}} \frac{a_{k,k_{\text{max}}}}{2 - \delta_{k,0}} \mathcal{H}_{lm}(t_m + kT_m), \qquad (75)$$

where $t_m$ and $T_m$ are defined by

$$2\pi f = m\dot{\Phi}(t_m), \qquad (76)$$

$$T_m = \frac{1}{\sqrt{m\dot{\Phi}(t_m)}}, \qquad (77)$$

$$\mathcal{H}_{lm} = \frac{1}{2}(F_{+} + iF_{\times}) \\ \times \sum_{m'=-l}^{l} h^{lm} D_{m',m}^{l}(\phi_z, \theta_L, \zeta)_{-2} Y_{lm'}(\theta_s, \phi_s) \\ + \frac{1}{2}(F_{+} - iF_{\times}) \\ \times \sum_{m'=-l}^{l} h^{l,-m} D_{m',-m}^{l}(\phi_z, \theta_L, \zeta)_{-2} Y_{lm'}(\theta_s, \phi_s), \qquad (78)$$

and the constants $a_{k,k_{\text{max}}}$ satisfy the linear system

$$\frac{(-i)^p}{2^p p!} = \sum_{k=0}^{k_{\text{max}}} a_{k,k_{\text{max}}} \frac{k^{2p}}{(2p)!}, \qquad (79)$$

for $p \in \{0, \dots, k_{\text{max}}\}$. In this expression, Eq. (76) expresses the stationary time $t_m$ as a function of the frequency $f$. For a LIGO-type detector, $\mathcal{H}_{lm}$ depends on time through $\phi_z$, $\theta_L$, and $\zeta$.

Figure 5 compares the frequency domain GWs for the 4 systems of Table I using only the leading ($l = 2, m = 2$) harmonics of Eq. (75). The two waveforms are computed with the numerical solution to the PN precession equations and with the analytic solution described in Secs. III and IV. Both waveforms are Fourier-transformed with SUA, allowing us to assess the effect of our new analytic solution to the GW amplitude and phase. The agreement between the wave amplitudes is excellent over a wide range of frequencies, while the dephasing between the two waveforms never exceeds 0.3 radians, even for our BHBH system. This figure serves as a first indication of the accuracy of our model to accurately capture generic precessing features in GWs.
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>FIG. 5. (Top Panel) Amplitude of the GW including only the dominant ($\ell = 2, m = 2$) harmonic as a function of the GW frequency for the numerical PN and analytic SUA waveforms for the NSNS (Top Left), the BHNS (Top Right), the BHBH (Bottom Left), and the HSNSBH (Bottom Right) system of Table I. The reference amplitude $|\tilde{h}_{\text{ref}}|$ is the numerical PN SUA amplitude at 50 Hz. (Bottom Panel) GW dephasing between the numerical PN and analytic SUA waveforms.</figcaption></figure>

## VI. WAVEFORM COMPARISON

In order to have a more complete picture of our waveform’s ability to model generic systems, we carry out a Monte Carlo study randomizing over the 15 parameters describing a quasicircular compact binary waveform. For the randomization, we draw the components’ masses form a flat distribution in log space between [1, 2.5] $M_{\odot}$ for NSs and [2.5, 20] $M_{\odot}$ for BHs, while the components’ spin magnitudes are uniformly distributed in [0, 0.1] for NSs and [0, 1] for BHs. We selected seemingly low BH masses in order to focus on systems for which the inspiral part is the most important. Indeed, those are the ones for which the accurate modeling of the precession effects are the most challenging, due to the increased number of precession cycles that low masses entail. All directions (spin, sky location, orbital angular momentum) are drawn uniformly on a unit sphere. The phase of coalescence is assumed to be uniform in [0, 2π], while the time of coalescence and the distance are fixed at $10^5$ seconds and 100 Mpc respectively.

The large number of systems simulated can only be analyzed through some appropriate and efficient statistic; we use the faithfulness (or match) defined as

$$F \equiv \max_{t_c, \phi_c} \frac{(h_1|h_2)}{\sqrt{(h_1|h_1)(h_2|h_2)}}. \qquad (80)$$

The faithfulness is calculated between two waveforms $h_1$ and $h_2$ with the same physical parameters, but maximized over any unphysical parameters: the time $t_c$ and phase $\phi_c$ of coalescence. As such, it is a good estimator of a model’s suitability for parameter estimation. The faithfulness always falls between -1 and 1, with the latter indicating perfect agreement between the waveforms.

Unlike fitting factors,<fnref n="9" /> selecting a value for the faithfulness that is “good enough” is not straightforward. The nominal fitting factor threshold of 0.965 corresponds to a

<footnote n="9">⁹A fitting factor is the faithfulness maximized over all model parameters, quantifying a model’s suitability for detection.</footnote>
10% drop in detection rates. On the other hand, a faithfulness threshold should be translatable to a requirement about parameter estimation accuracy: the systematic mismodeling error should be smaller than the statistical measurement error. The latter depends on the signal-to-noise ratio (SNR) of the signal, while the former does not, meaning that any faithfulness threshold should take the strength of the signal into account. In Appendix G we calculate the faithfulness threshold as a function of the SNR and find that for an SNR of 10(25)[50], a faithfulness of 0.96(0.9936)[0.9984] suffices for accurate parameter estimation. Led by the SNR of the first detected GW [67], we set our faithfulness threshold to 0.994.

In our study $h_1$ is a waveform calculated by numerically solving the precession equations, while $h_2$ uses our new analytic solution. Both waveforms are Fourier-transformed with the SUA method, justified by [17] where it was shown that SUA induces a negligible loss of faithfulness compared to a discrete Fourier transform. The use of SUA in both waveforms allows us to isolate the effect of our new solution: any mismatch is solely caused by the solution to the precession equations described in this paper.

The inner product in Eq. (80) is defined in the usual way

$$ (h_1|h_2) \equiv 4\Re \int_{f_{\min}}^{f_{\max}} \frac{\tilde{h}_1(f)\tilde{h}_2^*(f)}{S_n(f)} df, \quad (81) $$

<figure><img src="image_12.png" /><figcaption>FIG. 6. Distribution of $1-F$ for NSNS (Top Left), BHNS (Top Right), BHBH (Bottom Left) and the 4th generic set containing all masses and spins (Bottom Right) for waveforms with full harmonic content (solid black line) and waveforms restricted to the leading ($\ell = 2, m = 2$) mode (dashed red line). Top panels show cumulative distribution functions, and bottom panels give the corresponding probability distribution function. The leftmost vertical line denotes a faithfulness $F = 0.994$, corresponding to 0.3% of systems for NSNS (both waveforms), 1.6% for BHNS (both waveforms), 10.4% (RWF) and 10.7% (FWF) for BHBH, and 9.2% (RWF) and 9.9% (FWF) for the 4th set. The other vertical lines correspond to the medians of the distributions, which are $1-F = 2.7 \times 10^{-6}$ (RWF), and $1-F = 6 \times 10^{-6}$ (FWF) for NSNS, $1-F = 7.6 \times 10^{-5}$ (RWF), and $1-F = 8.3 \times 10^{-5}$ (FWF) for NSBH, $1-F = 7.1 \times 10^{-4}$ (RWF), and $1-F = 7.4 \times 10^{-4}$ (FWF) for BHBH, and $1-F = 9.3 \times 10^{-4}$ (RWF), and $1-F = 8.6 \times 10^{-4}$ (FWF) for the 4th set.</figcaption></figure>
<table><caption>TABLE II. Percentages of subthreshold systems encountered in our analysis for each type of system.</caption><thead><tr><th>Waveform</th><th>Threshold</th><th>NSNS</th><th>BHNS</th><th>BHBH</th><th>HSNSBH</th></tr></thead><tbody><tr><td>RWF</td><td>0.965</td><td>0.06%</td><td>0.33%</td><td>1.85%</td><td>1.14%</td></tr><tr><td>RWF</td><td>0.994</td><td>0.3%</td><td>1.6%</td><td>10.4%</td><td>9.2%</td></tr><tr><td>FWF</td><td>0.965</td><td>0.06%</td><td>0.33%</td><td>1.85%</td><td>1.26%</td></tr><tr><td>FWF</td><td>0.994</td><td>0.3%</td><td>1.6%</td><td>10.7%</td><td>9.9%</td></tr></tbody></table>

where $f_{min} = 10$ Hz is aLIGO's lower frequency cutoff, $f_{max}$ is the frequency that corresponds to an orbital separation of $6M$, and $S_n(f)$ is aLIGO's design zero-detuning, high power noise spectral density [68].

Figure 6 shows the distributions of $1-F$ for 4 sets, each containing 10,000 systems. The first 3 sets contain systems with masses and spins corresponding to NSNS, BHNS, and BHBH systems, respectively. The fourth set contains an additional, less astrophysically motivated but useful to test our model in the most challenging setting, type of system where both masses were drawn from a log-flat distribution ranging from $1 M_{\odot}$ to $20 M_{\odot}$, and both spin magnitudes uniformly distributed in $[0, 1]$. We study 2 different types of waveforms: full waveforms (FWF) contain all the known harmonics in Eq. (75), while restricted waveforms (RWF) contain only the dominant ($\ell = 2, m = 2$) harmonic.

The agreement between our analytical waveform and the numerical PN one is excellent for a wide range of parameters. In the NSNS case we find that only 0.06% (0.3%) of the systems have a faithfulness below the 0.965 (0.994) for both waveforms, while for BHNS systems, this number is 0.33%(1.6%) for both waveforms. The percentage of systems below the nominal faithfulness threshold is increased to 1.85%(10.4%) (RWF) and 1.85%(10.7%) (FWF) in the case of BHBHs and 1.14%(9.2%) (RWF) and 1.26%(9.9%) (FWF) for the 4th generic set. This increase is not unexpected, since precessional features are more pronounced, and hence more difficult to model, when the spins are large and the masses different. In the next section we study the various sources of error in our analytical waveform and quantify their effect. Table II summarizes these results.

## VII. SOURCE OF ERROR

The subthreshold systems of Fig. 6 can be split into two rough categories: systems for which the faithfulness is very low $F \lesssim 0.8$, and systems for which the faithfulness is high, but not high enough $0.8 \lesssim F \lesssim 0.994$. Systems falling into the first category can mainly be explained by the effect described in Appendix E. For them the orbital angular momentum becomes approximately (anti)aligned with the total spin angular momentum at some point in the evolution of the systems. In this case the PN expansion of Eq. (62) becomes ill-defined. Our specific choice for the values of $D_2$ and $D_4$ in Eqs. (63) and (64) to some extent ameliorates

this problem, yet it does not fully solve it. We have explored many choices for $D_2$ and $D_4$, some even leading up to 8% of systems with faithfulnesses below 0.96 in the BHBH case. The particular values for $D_2$ and $D_4$ we employ in our model [Eqs. (E3) and (E4)] yield the best results among all the expressions we tested.

Systems falling in the second category can be modeled accurately only for low SNR signals. The unfaithfulness of these systems can be attributed to the various approximations we have used in our model construction. In order to quantify the effect of each approximation, we retrace the steps we followed in Sec. IV to add radiation reaction effects to the exact precession solution of Sec III.

(1) Our first task when adding radiation reaction effects is to specify a coordinate system. In Sec. IVA we assume that $\hat{J}$ is constant and identify it with the $\hat{z}$ axis of our system.

(2) In Sec. IV B we use MSA to solve for the magnitude of the total angular momentum $J$.

(3) In Sec. IV C we use MSA and a PN approximation to solve for the total spin magnitude $S$.

(4) Finally, in Sec IV D we use the $J$ and $S$ obtained above to solve for the precession angle $\phi_z$.

Overall, the addition of radiation reaction effects requires the identification of a coordinate system and the solution to 3 coupled differential equations. Below we perform each of these steps numerically and each time compute matches for BHBH systems in order to quantify the improvement. The unfaithfulness distributions are given in Fig. 7, while the inset focuses on the region of interest $F \in [0.9, 0.999]$ with the vertical line denoting $F = 0.994$. The different curves in this figure represent the following:

<figure><img src="image_13.png" /><figcaption>FIG. 7. Cumulative distribution of $1-F$ for waveforms including different degrees of numerical and analytical calculations. The difference between these distributions is a quantification of the error from each analytic approximation we have made. The inset shows the faithfulness range of interest $F \in [0.9, 0.999]$, while the vertical line denotes the faithfulness threshold $F = 0.994$. See the text for more details and discussion.</figcaption></figure>
(i) The black solid line corresponds to the fully analytic waveform also studied in Fig. 6, with approximately 10% of the systems below 0.994.

(ii) The maroon dashed line was created with a waveform that evaluates Eq. (52) for $\phi_z$ with the analytical $J$ and $S$, and then solves this equation numerically, bringing the percentage of subthreshold systems down to ~4%.

(iii) The green dot-dashed line again uses a numerical solution to Eq. (52) but where now $S$ and $J$ are obtained and substituted by numerically solving their corresponding differential equations. We stress that in this version of the waveform, $S$ and $J$ are solved for numerically only when used in Eq. (52). This corresponds to the most accurate solution for $\phi_z$ possible and results in ~2% of subthreshold systems.

(iv) The blue dot-double-dashed line again uses the most accurate $\phi_z$ from the previous waveform, and now the numerical $S$ is also used for the entire waveform. This brings the percentage of subthreshold systems down to ~1%.

(v) Finally, the red dotted line is produced with a waveform that solves for all $\phi_z$, $S$, and $J$ numerically. The faithfulness distribution for this waveform is almost indistinguishable from the previous waveform with ~1% of subthreshold systems.

The above results suggest a clear-cut way to improve our model if a more faithful waveform model is required in the future. The largest improvement would be obtained if we found a more accurate solution to Eq. (52) for the precession angle $\phi_z$, either by improving the PN solution to Eq. (61), or by taking Eq. (67) to higher order in MSA. We have studied those two error sources separately, and have found the former to dominate over the latter. The second step would be to improve the solution for the total spin magnitude by employing elements of MSA to solve Eq. (22), or a more accurate PN prescription for Eq. (50) for the phase of $S$. Finally, improving the solution for the magnitude of the total angular momentum $J$ will not affect the waveform considerably.

Even with these improvements, we still have ~1% of our systems with $F < 0.994$. The only approximation we have made for these systems is that the total angular momentum has a fixed direction, and to verify that this approximation indeed breaks, we studied 100 of these systems explicitly. We evolved $\vec{J}$ and found that 2 of these systems undergo the effect of transitional precession [19]. The remaining 98 systems do not undergo transitional precession but (i) have parameters that are statistically consistent with what is required for transitional precession (small mass ratios, large and misaligned spin for the larger body, everything else random), and (ii) the components $J_x$ and $J_y$ always remain smaller than $J_z$, but by about an order of magnitude only. Compare this with Fig. 2 where $J_x$ and $J_y$ remain at least 3 orders of magnitude below $J_z$ for all frequencies.

As a concluding remark, we should mention the effect of the inclination angle $L \cdot N$ on our results. It is well known that binaries observed approximately edge-on ($L \cdot N \sim 0$) exhibit the largest precessional effects. We indeed find that all other things being equal, the more edge-on the binary is, the lower its faithfulness. This is because when precessional effects become more pronounced, a better and better modeling of them is required in order to achieve a certain goodness of fit. In other words, inclination does not cause problems on its own, but rather it amplifies preexisting ones, an effect also observed in [37]. We should note, however, that edge-on systems are less likely to be detected by aLIGO due to selection effects [69].

VIII. DISCUSSION

We have constructed the first closed-form fully analytic GW template in the frequency domain that can accurately model quasicircular systems of generic masses, spin magnitudes, and spin orientations in the inspiral phase. We expand the exact solution to the precession equations in the absence of radiation reaction derived by Kesden et al. [10] to include radiation reaction using elements from multiple scale analysis. This allows us to derive the first closed-form time-domain GW model valid for generic inspirals. We then use the method of shifted uniform asymptotics to transform this waveform from the time domain to the frequency domain.

The resulting waveform is ideal for extracting parameters from generically precessing quasicircular inspirals as demonstrated by a Monte Carlo study of 40,000 system; only 0.8%(5.4%) of them had a faithfulness with a numerical PN waveform that solves the precession equations numerically below 0.965(0.994). The remaining inaccuracies of our model can be mapped back to specific assumptions we made while solving the spin-precession equations including radiation reaction. Analytical understanding of all these assumptions and the elements that enter our waveform construction enable us to improve the accuracy of our model if deemed necessary when more sensitive GW detector networks become available. This is, perhaps, the most attractive feature of having analytic control over complicated processes like spin-precession.

Finally, analytic methods have the potential to be much faster than numerical ones, while still encompassing all precessional effects. We estimate that our analytic SUA waveform can be up to 15 times faster to evaluate than the numerical SUA waveform in certain regions of the parameter space. Interestingly, the region of the parameter space where the analytic SUA waveform presents the maximum improvement over the numerical SUA waveforms (the BHBH case) is distinct from the region where the numerical SUA waveform is much faster than fully numerical PN time-domain models (the NSNS case) [17]. This suggests that a hybrid model where the numerical or the analytical SUA is called depending on the system's mass can achieve both
high accuracy and numerical efficiency. Further improvement could be obtained through reduced order modeling and reduced order quadrature integration in data analysis implementations [70]. We leave such studies to future work.

As a final remark, we note that our results lay the framework for the construction of full inspiral-merger-ringdown (IMR) waveforms following similar procedures as IMRPhenomP [23,33,34]. This is a promising avenue for future research since such an IMR waveform has the potential to be more accurate than the IMRPhenomP due to more accurate description of precessional dynamics, as well as faster than SEOBNRv3 [29] due to the analytic treatment of the inspiral dynamics.

ACKNOWLEDGMENTS

We would like to thank Emanuele Berti, Mike Kesden, Sylvain Marsat, and Frank Ohme for helpful discussions. K. C. acknowledges support from the Onassis Foundation. N. Y. acknowledges support NSF CAREER Grant No. PHY-1250636. N. C. acknowledges support from the NSF Award No. PHY-1306702. N. C. and N. Y. acknowledges support from NASA Grant No. NNX16AB98G. A. K. is supported by NSF CAREER Grant No. PHY-1055103, by FCT contract IF/00797/2014/CP1214/CT0012 under the IF2014 Programme, and by the H2020-MSCA-RISE-2015 Grant No. StronGrHEP-690904. This work was supported by the Centre National d'Études Spatiales.

APPENDIX A: COEFFICIENTS OF v

The coefficients of the evolution of the PN parameter v as defined in Eq. (4) are

$$g_0 = \frac{1}{a_0}, \qquad (A1)$$

$$g_2 = -\frac{a_2}{a_0}, \qquad (A2)$$

$$g_3 = -\frac{a_3}{a_0}, \qquad (A3)$$

$$g_4 = -\frac{a_4 - a_2^2}{a_0}, \qquad (A4)$$

$$g_5 = -\frac{a_5 - 2a_3a_2}{a_0}, \qquad (A5)$$

$$g_6 = -\frac{a_6 - 2a_4a_2 - a_3^2 + a_2^3}{a_0}, \qquad (A6)$$

$$g_6' = -\frac{3b_6}{a_0}, \qquad (A7)$$

$$g_7 = -\frac{a_7 - 2a_5a_2 - 2a_4a_3 + 3a_3a_2^2}{a_0}, \qquad (A8)$$

with all other terms vanishing. The coefficients {$a_i, b_i$} are given in Appendix A of [39]. The spin couplings in the above expressions are evaluated with all angular momenta averaged over one precession cycle using the solution of Secs. III and IV. The error induced by this is of 4PN order, higher than the order to which we know the v expansion. Explicitly, we use

$$S_1 \cdot \hat{L} \rightarrow \langle S_1 \cdot \hat{L} \rangle_{\text{pr}} = \frac{c_1(1+q) - q\eta\xi}{\eta(1-q^2)}, \qquad (A9)$$

$$S_2 \cdot \hat{L} \rightarrow \langle S_2 \cdot \hat{L} \rangle_{\text{pr}} = -q \frac{c_1(1+q) - \eta\xi}{\eta(1-q^2)}, \qquad (A10)$$

$$S_1 \cdot S_2 \rightarrow \langle S_1 \cdot S_2 \rangle_{\text{pr}} = \frac{S_{\text{av}}^2}{2} - \frac{S_1^2 + S_2^2}{2}, \qquad (A11)$$

$$(S_1 \cdot \hat{L})^2 \rightarrow \langle (S_1 \cdot \hat{L})^2 \rangle_{\text{pr}} = \langle S_1 \cdot \hat{L} \rangle_{\text{pr}}^2 + \frac{(S_+^2 - S_-^2)^2 v_0^2}{32\eta^2(1-q)^2}, \qquad (A12)$$

$$(S_2 \cdot \hat{L})^2 \rightarrow \langle (S_2 \cdot \hat{L})^2 \rangle_{\text{pr}} = \langle S_2 \cdot \hat{L} \rangle_{\text{pr}}^2 + \frac{q^2(S_+^2 - S_-^2)^2 v_0^2}{32\eta^2(1-q)^2}, \qquad (A13)$$

$$(S_1 \cdot \hat{L})(S_2 \cdot \hat{L}) \rightarrow \langle (S_1 \cdot \hat{L})(S_2 \cdot \hat{L}) \rangle_{\text{pr}} = \langle S_1 \cdot \hat{L} \rangle_{\text{pr}} \langle S_2 \cdot \hat{L} \rangle_{\text{pr}} - \frac{q(S_+^2 - S_-^2)^2 v_0^2}{32\eta^2(1-q)^2}, \qquad (A14)$$

where $v_0$ corresponds to the value of v at the initial time, $\xi$, $S_{\text{av}}^2$, $c_1$ are defined in Eqs. (7), (45), and (41) respectively, and $S_+^2, S_-^2$ are the roots of the right-hand side of Eq. (21).

APPENDIX B: COEFFICIENTS OF THE PRECESSION SOLUTION

Below, we provide explicit expressions for the coefficients appearing in the precession solution of Sec. III. The coefficients of Eq. (21) are

$$A = -\frac{3}{2\sqrt{\eta}} v^6 (1 - \xi v), \qquad (B1)$$

$$B = (L^2 + S_1^2)q + 2L\xi - 2J^2 - S_1^2 - S_2^2 + \frac{L^2 + S_2^2}{q}, \qquad (B2)$$

$$C = (J^2 - L^2)^2 - 2L\xi(J^2 - L^2) - 2\frac{1-q}{q}(S_1^2 - qS_2^2)L^2 + 4\eta L^2 \xi^2 - 2\delta m(S_1^2 - S_2^2)\xi L + 2\frac{1-q}{q}(qS_1^2 - S_2^2)J^2, \qquad (B3)$$
$$D = \frac{1-q}{q} (S_2^2 - qS_1^2)(J^2 - L^2)^2 + \frac{\delta m^2}{\eta} (S_1^2 - S_2^2)^2 L^2 + 2\delta mL\xi(S_1^2 - S_2^2)(J^2 - L^2). \quad (B4)$$

The coefficients of Eq. (30) are

$$a = \frac{1}{2} v^6 \left\{ 1 + \frac{3}{2\eta} (1 - \xi v) \right\}, \quad (B5)$$

$$c_0 = \frac{3}{4}(1 - \xi v)v^2 \left\{ \eta^3 + 4\eta^3\xi v - 2\eta[J^2 - S_+^2 + 2(S_1^2 - S_2^2)\delta m]v^2 - 4\eta\xi(J^2 - S_+^2)v^3 + \frac{(J^2 - S_+^2)^2}{\eta}v^4 \right\}, \quad (B6)$$

$$c_2 = -\frac{3\eta}{2}(S_+^2 - S_-^2)\left(1 + 2\xi v - \frac{J^2 - S_+^2}{\eta^2}v^2\right)(1 - \xi v)v^4, \quad (B7)$$

$$c_4 = \frac{3}{4\eta}(S_+^2 - S_-^2)^2(1 - \xi v)v^6, \quad (B8)$$

$$d_0 = -[J^2 - (L + S_+)^2][J^2 - (L - S_+)^2], \quad (B9)$$

$$d_2 = -2(S_+^2 - S_-^2)(J^2 + L^2 - S_+^2), \quad (B10)$$

$$d_4 = -(S_+^2 - S_-^2)^2. \quad (B11)$$

The coefficients of Eq. (31) are

$$A_{\phi} = A + \frac{c_4}{d_4}, \quad (B12)$$

$$B_{\phi} = \left( \frac{c_4}{d_4} - \frac{c_0 + c_2 + c_4}{d_0 + d_2 + d_4} \right), \quad (B13)$$

$$C_{\phi} = C_1 + C_2, \quad (B14)$$

$$D_{\phi} = C_1 - C_2, \quad (B15)$$

$$n_c = 2 \frac{d_0 + d_2 + d_4}{2d_0 + d_2 + s_d}, \quad (B16)$$

$$n_d = \frac{2d_0 + d_2 + s_d}{2d_0}, \quad (B17)$$

where

$$C_1 = -\frac{1}{2} \left( \frac{c_0}{d_0} - \frac{c_0 + c_2 + c_4}{d_0 + d_2 + d_4} \right), \quad (B18)$$

$$C_2 = \frac{c_0(-2d_0d_4 + d_2^2 + d_2d_4) - c_2d_0(d_2 + 2d_4)}{2d_0(d_0 + d_2 + d_4)s_d} + \frac{c_4d_0(2d_0 + d_2)}{2d_0(d_0 + d_2 + d_4)s_d}, \quad (B19)$$

$$s_d = \sqrt{d_2^2 - 4d_0d_4}. \quad (B20)$$

## APPENDIX C: COEFFICIENTS OF $\psi$

The coefficients in Eq. (51) are

$$\psi_1 = 3 \frac{2\xi\eta^2 - c_1}{\eta\delta m^2}, \quad (C1)$$

$$\psi_2 = \frac{3g_2}{g_0} + \frac{3}{2\eta^3} \left\{ 2\Delta - 2\frac{\eta^2}{\delta m^2} S_{\text{av}}^2 - 10\frac{\eta}{\delta m^4} c_1^2 + 2\frac{\eta^2}{\delta m^2} \frac{7+6q+7q^2}{(1-q)^2} c_1\xi - \frac{\eta^3}{\delta m^2} \frac{3+4q+3q^2}{(1-q)^2} \xi^2 + \frac{\eta}{(1-q)^2} [q(2+q)S_1^2 + (1+2q)S_2^2] \right\}, \quad (C2)$$

where

$$\Delta = \left\{ \left\{ \frac{c_1^2 \eta}{q \delta m^4} - \frac{2 c_1 \eta^3 (1+q)}{q \delta m^4} \xi - \frac{\eta^2}{\delta m^4} [\delta m^2 S_1^2 - \eta^2 \xi^2] \right\} \times \left\{ \frac{c_1^2 \eta^2}{\delta m^4} - \frac{2 c_1 \eta^3 (1+q)}{\delta m^4} \xi - \frac{\eta^2}{\delta m^4} [\delta m^2 S_2^2 - \eta^2 \xi^2] \right\} \right\}^{1/2}. \quad (C3)$$

## APPENDIX D: COEFFICIENTS OF $\phi_z$

For the PN expansions of the coefficients of Eq. (65) we define

$$R_m = S_+^2 - S_-^2, \quad (D1)$$

$$c_p = (S_+^2 \eta^2 - c_1^2), \quad (D2)$$

$$c_m = (S_-^2 \eta^2 - c_1^2), \quad (D3)$$

$$a_1 = \frac{1}{2} + \frac{3}{4}\eta, \quad (D4)$$

$$a_2 = -\frac{3}{4\eta}\xi, \quad (D5)$$

$$a_d = \frac{-3(S_1^2 - S_2^2) \frac{\eta}{\delta m} + 3 \frac{c_1}{\eta} (c_1 - 2\xi\eta^2)}{4\sqrt{c_p c_m}}, \quad (D6)$$
$$c_d = -\frac{3}{128} \frac{R_m}{\eta \sqrt{c_p c_m}}, \qquad (D7)$$

$$h_d = \frac{c_1}{\eta^2} \left( 1 - \frac{c_p + c_m}{2\sqrt{c_p c_m}} \right), \qquad (D8)$$

$$\Omega_{z,0} = a_1 + a_d, \qquad (D9)$$

$$\Omega_{z,1} = a_2 - a_d \xi - a_d h_d, \qquad (D10)$$

$$\Omega_{z,2} = a_d h_d \xi + c_d - a_d f_d + a_d h_d^2, \qquad (D11)$$

$$\Omega_{z,3} = (a_d f_d - c_d - a_d h_d^2)(\xi + h_d) + a_d f_d h_d, \qquad (D12)$$

$$\Omega_{z,4} = (c_d + a_d h_d^2 - 2a_d f_d)(h_d \xi + h_d^2 - f_d) - a_d f_d^2, \qquad (D13)$$

$$\Omega_{z,5} = (c_d - a_d f_d + a_d h_d^2) f_d (\xi + 2h_d) - (c_d + a_d h_d^2 - 2a_d f_d) h_d^2 (\xi + h_d) - a_d f_d^2 h_d, \qquad (D14)$$

With these definitions, the coefficients of Eq. (65) are

$$\langle \Omega_z \rangle^{(0)} = 3g_0 \Omega_{z,0}, \qquad (D15)$$

$$\langle \Omega_z \rangle^{(1)} = 3g_0 \Omega_{z,1}, \qquad (D16)$$

$$\langle \Omega_z \rangle^{(2)} = 3(g_0 \Omega_{z,2} + g_2 \Omega_{z,0}), \qquad (D17)$$

$$\langle \Omega_z \rangle^{(3)} = 3(g_0 \Omega_{z,3} + g_2 \Omega_{z,1} + g_3 \Omega_{z,0}), \qquad (D18)$$

$$\langle \Omega_z \rangle^{(4)} = 3(g_0 \Omega_{z,4} + g_2 \Omega_{z,2} + g_3 \Omega_{z,1} + g_4 \Omega_{z,0}), \qquad (D19)$$

$$\langle \Omega_z \rangle^{(5)} = 3(g_0 \Omega_{z,5} + g_2 \Omega_{z,3} + g_3 \Omega_{z,2} + g_4 \Omega_{z,1} + g_5 \Omega_{z,0}). \qquad (D20)$$

The functions $\phi_z^{(n)}$ in Eq. (66) are

$$\phi_z^{(0)} = \frac{J}{\eta^4} \left( \frac{c_1^2}{2} - \frac{c_1 \eta^2}{6v} - \frac{S_{av}^2 \eta^2}{3} - \frac{\eta^4}{3v^2} \right) - \frac{c_1}{2\eta} \left( \frac{c_1^2}{\eta^4} - \frac{S_{av}^2}{\eta^2} \right) l_1, \qquad (D21)$$

$$\phi_z^{(1)} = -\frac{J}{2\eta^2} (c_1 + \eta L) + \frac{1}{2\eta^3} (c_1 - \eta^2 S_{av}^2) l_1, \qquad (D22)$$

$$\phi_z^{(2)} = -J + \sqrt{S_{av}^2} l_2 - \frac{c_1}{\eta} l_1, \qquad (D23)$$

$$\phi_z^{(3)} = Jv - \eta l_1 + \frac{c_1}{\sqrt{S_{av}^2}} l_2, \qquad (D24)$$

$$\phi_z^{(4)} = \frac{J}{2S_{av}^2} v(c_1 + vS_{av}^2) - \frac{1}{2(S_{av}^2)^{3/2}} (c_1^2 - \eta^2 S_{av}^2) l_2, \qquad (D25)$$

$$\phi_z^{(5)} = -Jv \left( \frac{c_1^2}{2(S_{av}^2)^2} - \frac{c_1 v}{6S_{av}^2} - \frac{v^2}{3} - \frac{\eta^2}{3S_{av}^2} \right) + \frac{c_1}{2(S_{av}^2)^{5/2}} (c_1^2 - \eta^2 S_{av}^2) l_2, \qquad (D26)$$

where we have defined

$$l_1 = \ln(c_1 + J\eta + L\eta), \qquad (D27)$$

$$l_2 = \ln(c_1 + J\sqrt{S_{av}^2 v + S_{av}^2 v}), \qquad (D28)$$

In the above expressions we keep the roots $S_+^2$ and $S_-^2$ constant and equal to their initial value. The complexity of the roots' PN expansion makes its use prohibitive. Note that we do not expand the roots at all, but rather use their initial value as a form of partial resummation to increase the accuracy of our results. We find that this approximation does not affect our final result for the GW significantly.

## APPENDIX E: JUSTIFICATION OF THE $\phi_z$ CALCULATION

The precession-averaged $\langle \dot{\phi}_z \rangle_{\text{pr}}$ given in Eqs. (62)–(64) is exact. In principle, we could calculate $D_2$ and $D_4$ as functions of $v$, substitute them in Eq. (62), and carry out a PN expansion and integration to obtain $\langle \phi_z \rangle_{\text{pr}} = \phi_{z,-1}$.

Though this approach should work, in practice we run into 2 considerable problems. First, the resulting $\phi_{z,-1}$ is ill-behaved in the small mass ratio limit, despite never having assumed comparable masses. Second, $\phi_{z,-1}$ diverges when, at any point in the evolution of a precessing system, the total spin angular momentum is (anti)aligned with the orbital angular momentum. We stress that this does not mean that $S$ is approximately (anti)aligned with $L$ all the time; a brief moment of (anti)alignment suffices.

Both issues are not caused by real physical divergences in $\Omega_z$. First, at no point did we assume comparable masses. The second issue is more subtle. It might be true that the denominator of $\Omega_z$ vanishes if $S$ and $L$ are (anti)aligned. However, the binary (and $\phi_z$) is well behaved at the moment of (anti)alignment since the numerator of $\Omega_z$ vanishes too, leading to a 0/0 type situation<fnref n="10" />

We argue that even though Eq. (62) is well behaved in both the small mass ratio and the (anti)alignment between $S$ and $L$ limit, the same need not be true for its PN expansion. Consider the following function<fnref n="11" />

<footnote n="10">¹⁰We have verified that this is the case both analytically and numerically.</footnote>

<footnote n="11">¹¹The similarity between our toy function and $J$ given in Eq. (41) is not accidental.</footnote>
$$h(x; h_2, h_1, h_0) = \sqrt{\frac{h_2^2}{x^2} + \frac{h_1}{x} + h_0}, \quad (E1)$$

and its expansion around $x = 0$

$$h^{\text{exp}}(x; h_2, h_1, h_0) = \frac{h_2}{x} + \frac{h_1}{2h_2} - \frac{h_1^2 - 4h_0h_4}{8h_2^3} x + \mathcal{O}(x^2). \quad (E2)$$

Clearly, $h(x)$ is finite as $h_2 \to 0$. However, $h^{\text{exp}}(x)$ is not, and the $h_2 \to 0$ limit is worse and worse as we keep more terms in the $x$ expansion.

This is exactly the situation we encounter with Eq. (62) both in the small mass ratio limit, and in the approximate S and L (anti)alignment limit. Fixing the small mass ratio limit is straightforward: we identify the problem as originating from expanding the J multiplying the entire right-hand side of Eq. (62), and factor it out. This is the reason behind the form Eq. (65) has.

The second problem is more complicated. We can still identify the terms that, when expanded, cause the limit when S and L are (anti)aligned to be problematic. However, if we do not expand them, we can no longer perform the integral of Eq. (65). Using this fully expanded $\langle\phi_z\rangle_{\text{pr}}$ causes 5% of the systems studied here to have faithfulnesses below threshold (see Sec. VI).

In light of this, we tried a number of alternative, approximate methods for calculating $\langle\phi_z\rangle_{\text{pr}}$. We discovered that if we keep the terms $D_2$ and $D_4$ in Eqs. (63) and (64), our results are greatly improved by about an order of magnitude: only 0.8% of the systems are below the faithfulness threshold (see Sec. VI). We examined a number of different definitions for $D_2$ and $D_4$, from using their initial value as given directly from Eqs. (63) and (64) to retaining different orders in a PN expansion, but evaluated at the initial time. We found out that these methods give comparable results, so we choose, for simplicity, to set $D_2$ and $D_4$ equal to their leading PN order:

$$D_2 \rightarrow \frac{c_p - \sqrt{c_p c_m}}{R_m \eta^2}, \quad (E3)$$

$$D_4 \rightarrow \frac{c_p(c_p - \sqrt{c_p c_m})}{R_m^2 \eta^4} - \frac{\sqrt{c_p c_m}}{2R_m \eta^2}, \quad (E4)$$

where $c_p, c_m, R_m$ are defined in Appendix D.

We expect this problem to be solved if we consistently PN expand both ($D_2, D_4$) and the roots $S_+^2, S_-^2$ (see Appendix D). The complexity of the roots' expansion poses some serious problems in this calculation and we here opt for the approach described above and the partial resummation of the roots explained in Appendix D. This approach yields satisfactory results for the waveform precision required for aLIGO (see Fig. 6), but can be

improved if need be through expansions appropriate for these systems, like a small misalignment between S and L expansion.

## APPENDIX F: COEFFICIENTS OF $\zeta$

The angle $\zeta$ that enters in the transformation to the waveform to the frame corotating with the precession of L can be calculated by solving

$$\dot{\zeta} = \dot{\phi}_z \cos \theta_L = \Omega_z \cos \phi_L \equiv \Omega_{\zeta}. \quad (F1)$$

The solution to this equation can be obtained through MSA and it is very similar to the solution to Eq. (30)

$$\zeta = \zeta_{-1} + \zeta_0, \quad (F2)$$

where

$$\zeta_{-1} = \int \langle \Omega_{\zeta} \rangle_{\text{pr}}(t_{\text{rr}}) dt_{\text{rr}}, \quad (F3)$$

$$\zeta_0 = \int [\Omega_{\zeta}(t_{\text{pr}}, t_{\text{rr}}) - \langle \Omega_{\zeta} \rangle_{\text{pr}}(t_{\text{rr}})] dt_{\text{pr}}. \quad (F4)$$

Following the same steps as in Sec. IV D 1 and for reasons explained in Appendix E we find

$$\zeta_{-1} = \eta v^{-3} \sum_{i=0}^{5} \langle \Omega_{\zeta} \rangle^{(n)} v^n + \zeta_{-1}^{0}, \quad (F5)$$

where $\zeta_{-1}^0$ is a constant of integration and we have defined

$$\langle \Omega_{\zeta} \rangle^{(0)} = -g_0 \Omega_{\zeta,0}, \quad (F6)$$

$$\langle \Omega_{\zeta} \rangle^{(1)} = -\frac{3}{2} g_0 \Omega_{\zeta,1}, \quad (F7)$$

$$\langle \Omega_{\zeta} \rangle^{(2)} = -3(g_0 \Omega_{\zeta,2} + g_2 \Omega_{\zeta,0}), \quad (F8)$$

$$\langle \Omega_{\zeta} \rangle^{(3)} = 3(g_0 \Omega_{\zeta,3} + g_2 \Omega_{\zeta,1} + g_3 \Omega_{\zeta,0}), \quad (F9)$$

$$\langle \Omega_{\zeta} \rangle^{(4)} = 3(g_0 \Omega_{\zeta,4} + g_2 \Omega_{\zeta,2} + g_3 \Omega_{\zeta,1} + g_4 \Omega_{\zeta,0}), \quad (F10)$$

$$\langle \Omega_{\zeta} \rangle^{(5)} = \frac{3}{2}(g_0 \Omega_{\zeta,5} + g_2 \Omega_{\zeta,3} + g_3 \Omega_{\zeta,2} + g_4 \Omega_{\zeta,1} + g_5 \Omega_{\zeta,0}), \quad (F11)$$

and

$$\Omega_{\zeta,0} = \Omega_{\zeta,0}, \quad (F12)$$

$$\Omega_{\zeta,1} = \Omega_{\zeta,1} + \frac{c_1}{\eta^2} \Omega_{\zeta,0}, \quad (F13)$$
$$ \Omega_{\zeta,2} = \Omega_{\zeta,2} + \frac{c_1}{\eta^2} \Omega_{\zeta,1}, \qquad (\text{F14}) $$

$$ \Omega_{\zeta,3} = \Omega_{\zeta,3} + \frac{c_1}{\eta^2} \Omega_{\zeta,2} + g_d, \qquad (\text{F15}) $$

$$ \Omega_{\zeta,4} = \Omega_{\zeta,4} + \frac{c_1}{\eta^2} \Omega_{\zeta,3} - g_d \xi - g_d h_d, \qquad (\text{F16}) $$

$$ \Omega_{\zeta,5} = \Omega_{\zeta,5} + \frac{c_1}{\eta^2} \Omega_{\zeta,4} + g_d h_d \xi + g_d (h_d^2 - f_d), \qquad (\text{F17}) $$

where the $\Omega_{\zeta,i}$'s, $f_d$ and $h_d$ are given in Appendix D and

$$ g_d = \frac{3 R_m^2 c_1 - \eta^2 \xi}{64 \eta^3 \sqrt{c_p c_m}}. \qquad (\text{F18}) $$

The first correction to MSA is given by

$$ \zeta_0 = \frac{A_{\theta_L}}{\dot{\psi}} (C_{\phi} + D_{\phi}) + 2d_0 \frac{B_{\theta_L}}{\dot{\psi}} \left( \frac{C_{\phi}}{s_d - d_2} - \frac{D_{\phi}}{s_d + d_2} \right), \qquad (\text{F19}) $$

where $C_{\phi}$, $D_{\phi}$, $d_0$ and $d_2$ are given in Appendix D, $\dot{\psi}$ is Eq. (51), and

$$ A_{\theta_L} = \frac{J^2 + L^2 - S_+^2}{2JL}, \qquad (\text{F20}) $$

$$ B_{\theta_L} = \frac{S_+^2 - S_-^2}{2JL}. \qquad (\text{F21}) $$

## APPENDIX G: FAITHFULNESS REQUIREMENT

The agreement between two waveforms $h, \bar{h}$ with parameters $\vec{\lambda}$ is measured in terms of the faithfulness $F$

$$ F(h, \bar{h}) = \frac{(h|\bar{h})}{\sqrt{(h|h)(\bar{h}|\bar{h})}}. \qquad (\text{G1}) $$

In the high SNR regime, a typical waveform sample from the posterior distribution function [71] is given by

$$ \bar{h} = h + h_{,i} \Delta \lambda^i + \frac{1}{2} h_{,ij} \Delta \lambda^i \Delta \lambda^j + \dots \qquad (\text{G2}) $$

where the $\Delta\vec{\lambda} = \vec{\lambda} - \vec{\lambda}_0$ are described by the multivariate normal distribution

$$ p(\Delta\vec{\lambda}) = \sqrt{\det(\Gamma/2\pi)} e^{-\Gamma_{ij}\Delta\lambda^i \Delta\lambda^j/2}, \qquad (\text{G3}) $$

with $\Gamma_{ij} = (h_{,i}|h_{,j})$ and $\vec{\lambda}_0$ are the true parameters. Treating the $\Delta\lambda^i$ as small and expanding we get

$$ F = 1 - \frac{1}{2} g_{ij} \Delta \lambda^i \Delta \lambda^j + \dots \qquad (\text{G4}) $$

where

$$ g_{ij} = \frac{(h_i|h_j)}{(h|h)} - \frac{(h|h_{,i})(h|h_{,j})}{(h|h)^2}. \qquad (\text{G5}) $$

Using $E[\Delta\lambda^i \Delta\lambda^j] = C^{ij} \approx \Gamma_{ij}^{-1}$, we find

$$ E[F] \approx 1 - \frac{(D-1)}{2\text{SNR}^2}, \qquad (\text{G6}) $$

for the expectation value of the faithfulness, where $D$ is the dimension of $\vec{\lambda}$. The factor of $D$ comes from $C^{ij}\Gamma_{ij} \approx \delta_i^i = D$ and the factor of $-1$ from the $(h|h_{,i})(h|h_{,j})/(h|h)^2$ removing the dependence on the overall amplitude of the waveform, thus reducing the dimensions count by one.

The expected value of the faithfulness in Eq. (G6) describes the impact of statistical errors. In deciding how accurate a waveform model needs to be, we should at a minimum demand that the systematic errors from mismodeling are smaller than the statistical errors. If we wish to model spin-precessing binaries with $D=8$ intrinsic parameters for systems with SNRs up to 50 then the modeling unfaithfulness should be below 8/5000 = 0.0016 (there is no $-1$ for just intrinsic parameters, the amplitude is extrinsic). For a SNR of 25 we obtain the faithfulness requirement of 0.994 that we used in Sec. VI.

To calculate the variance, it is easier to work with the unfaithfulness, $1-F$. The expectation of the square is given by

$$ \begin{aligned} E[(1-F)^2] &= \frac{1}{4} g_{ij} g_{kl} E[\Delta\lambda^i \Delta\lambda^j \Delta\lambda^k \Delta\lambda^l] \\ &= \frac{1}{4} g_{ij} g_{kl} (C^{ij} C^{kl} + C^{ik} C^{jl} + C^{il} C^{jk}) \\ &\approx \frac{3(D-1)^2}{4\text{SNR}^4}. \end{aligned} \qquad (\text{G7}) $$

Thus

$$ \text{var}[1-F] = \frac{2(D-1)^2}{4\text{SNR}^4}. \qquad (\text{G8}) $$

This shows that the average faithfulness is slightly less than $1-\sigma$ from a perfect faithfulness ($\sigma/\sqrt{2}$ to be precise). This agrees with what we see when computing the distribution of the match from MCMC waveform samples. The distribution is not Gaussian, and has a larger tail toward small values of the match.

An alternative derivation of the faithfulness requirement makes direct use of the posterior distribution function in the case of uniform priors

$$ p(\vec{\lambda}) \sim e^{-\frac{(d-h|d-h)}{2}}, \qquad (\text{G9}) $$
where d is the data. The peak of the posterior, evaluated at the best-fit parameters is

$$p(\vec{\lambda}_{bf}) \sim e^{-\frac{(d-h_{bf})(d-h_{bf})}{2}} \sim e^{-\frac{(d)(h_{bf})(h_{bf})-2(d-h_{bf})}{2}} \\ \sim e^{-\frac{\text{SNR}^2+\text{SNR}^2-2\text{SNR}^2\text{FF}}{2}} \sim e^{-\text{SNR}^2(1-\text{FF})}, \quad (G10)$$

where FF is the fitting factor, or the faithfulness maximized over all model parameters. The posterior on the true parameters is

$$p(\vec{\lambda}_0) \sim e^{-\frac{(d-h_0)(d-h_0)}{2}} \sim e^{-\text{SNR}^2(1-F)}, \quad (G11)$$

From Eq. (G3) we can calculate the value of the multi-dimensional posterior $1 - \sigma$ away from the best-fit parameters

$$p(\vec{\lambda}_{1-\sigma}) \sim e^{-\frac{\Gamma_{ij} \Delta^{ij} \Delta_{ij}}{2}} \sim e^{-\frac{\Gamma_{ij} c^{ij}}{2}} \sim e^{-\frac{D}{2}}, \quad (G12)$$

Assuming that the model can fit the data perfectly for some parameters (an assumption that will lead to a conservative faithfulness threshold) we set FF = 1 and requiring that the true parameters are less than $1 - \sigma$ away from the best-fit ones we find

$$1-F < \frac{D}{2\text{SNR}^2}, \quad (G13)$$

where D is the number of parameters whose measurability is affected by the model inaccuracy. For spin-precessing models with 8 intrinsic parameters, $D = 8$.

This derivation translates the results of [72] that were written in terms of requirements on the GW amplitude and phase to requirement on the faithfulness.

[1] H. Spruit and E. S. Phinney, Nature (London) 393, 139 (1998).

[2] C. Wang, D. Lai, and J. Han, Astrophys. J. 639, 1007 (2006).

[3] C. Wang, D. Lai, and J. Han, Astrophys. J. 656, 399 (2007).

[4] V. Kalogera, Astrophys. J. 541, 319 (2000).

[5] D. Gerosa, R. O'Shaughnessy, M. Kesden, E. Berti, and U. Sperhake, Phys. Rev. D 89, 124025 (2014).

[6] D. L. Kaplan, S. Chatterjee, B. M. Gaensler, and J. Anderson, Astrophys. J. 677, 1201 (2008).

[7] T. Bogdanovic, C. S. Reynolds, and M. C. Miller, Astrophys. J. Lett. 661, L147 (2007).

[8] C. L. Rodriguez, M. Zevin, C. Pankow, V. Kalogera, and F. A. Rasio, Astrophys. J. 832, L2 (2016).

[9] J. D. Schnittman, Phys. Rev. D 70, 124020 (2004).

[10] M. Kesden, D. Gerosa, R. O'Shaughnessy, E. Berti, and U. Sperhake, Phys. Rev. Lett. 114, 081103 (2015).

[11] D. Gerosa, M. Kesden, U. Sperhake, E. Berti, and R. O'Shaughnessy, Phys. Rev. D 92, 064016 (2015).

[12] D. Gerosa, M. Kesden, R. O'Shaughnessy, A. Klein, E. Berti, U. Sperhake, and D. Trifiro, Phys. Rev. Lett. 115, 141102 (2015).

[13] K. Chatziioannou, N. Cornish, A. Klein, and N. Yunes, Astrophys. J. Lett. 798, L17 (2015).

[14] K. Chatziioannou, N. Cornish, A. Klein, and N. Yunes, Phys. Rev. D 89, 104023 (2014).

[15] S. Droz, D. J. Knapp, E. Poisson, and B. J. Owen, Phys. Rev. D 59, 124016 (1999).

[16] C. M. Bender and S. A. Orszag, Advanced Mathematical Methods for Scientists and Engineers 1, Asymptotic Methods and Perturbation Theory (Springer, New York, 1999).

[17] A. Klein, N. Cornish, and N. Yunes, Phys. Rev. D 90, 124029 (2014).

[18] L. Blanchet, Living Rev. Relativ. 17, 2 (2014).

[19] T. A. Apostolatos, C. Cutler, G. J. Sussman, and K. S. Thorne, Phys. Rev. D 49, 6274 (1994).

[20] T. A. Apostolatos, Phys. Rev. D 52, 605 (1995).

[21] T. A. Apostolatos, Phys. Rev. D 54, 2421 (1996).

[22] T. A. Apostolatos, Phys. Rev. D 54, 2438 (1996).

[23] M. Hannam, P. Schmidt, A. Bohé, L. Haegel, S. Husa, F. Ohme, G. Pratten, and M. Pürrer, Phys. Rev. Lett. 113, 151101 (2014).

[24] A. Lundgren and R. O'Shaughnessy, Phys. Rev. D 89, 044021 (2014).

[25] A. Buonanno and T. Damour, Phys. Rev. D 59, 084006 (1999).

[26] A. Buonanno and T. Damour, Phys. Rev. D 59, 084006 (1999).

[27] T. Damour, Phys. Rev. D 64, 124013 (2001).

[28] A. Taracchini et al., Phys. Rev. D 89, 061502 (2014).

[29] 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).

[30] S. Khan, S. Husa, M. Hannam, F. Ohme, M. Pürrer, X. J. Forteza, and A. Bohé, Phys. Rev. D 93, 044007 (2016).

[31] M. Boyle, R. Owen, and H. P. Pfeiffer, Phys. Rev. D 84, 124011 (2011).

[32] P. Schmidt, M. Hannam, S. Husa, and P. Ajith, Phys. Rev. D 84, 024046 (2011).

[33] P. Schmidt, M. Hannam, and S. Husa, Phys. Rev. D 86, 104063 (2012).

[34] P. Schmidt, F. Ohme, and M. Hannam, Phys. Rev. D 91, 024043 (2015).

[35] B. P. Abbott et al. (Virgo, LIGO Scientific Collaboration), Phys. Rev. Lett. 116, 241102 (2016).

[36] B. Abbott et al. (Virgo, LIGO Scientific Collaboration), Phys. Rev. X 6, 041014 (2016).

[37] B. P. Abbott et al. (Virgo, LIGO Scientific Collaboration), arXiv:1611.07531.
[38] A. Klein, N. Cornish, and N. Yunes, Phys. Rev. D **88**, 124015 (2013).

[39] K. Chatzioannou, A. Klein, N. Yunes, and N. Cornish, Phys. Rev. D **88**, 063011 (2013).

[40] K. Chatzioannou, K. Yagi, A. Klein, N. Cornish, and N. Yunes, Phys. Rev. D **92**, 104008 (2015).

[41] K. Chatzioannou, A. Klein, N. Cornish, and N. Yunes, Phys. Rev. Lett. **118**, 051101 (2017).

[42] K. S. Thorne and J. B. Hartle, Phys. Rev. D **31**, 1815 (1985).

[43] B. Barker and R. O'Connell, Gen. Relativ. Gravit. **11**, 149 (1979).

[44] A. Buonanno, Y. Chen, and T. Damour, Phys. Rev. D **74**, 104005 (2006).

[45] E. Poisson, Phys. Rev. D **57**, 5287 (1998).

[46] N. V. Krishnendu, K. G. Arun, and C. K. Mishra, arXiv: 1701.06318.

[47] E. Poisson, Phys. Rev. D **70**, 084044 (2004).

[48] K. Chatzioannou, E. Poisson, and N. Yunes, Phys. Rev. D **87**, 044022 (2013).

[49] G. Faye, L. Blanchet, and A. Buonanno, Phys. Rev. D **74**, 104033 (2006).

[50] A. Bohe, S. Marsat, G. Faye, and L. Blanchet, Classical Quantum Gravity **30**, 075017 (2013).

[51] A. Bohé, G. Faye, S. Marsat, and E. K. Porter, Classical Quantum Gravity **32**, 195010 (2015).

[52] S. Marsat, Classical Quantum Gravity **32**, 085008 (2015).

[53] E. Racine, A. Buonanno, and L. Kidder, Phys. Rev. D **80**, 044010 (2009).

[54] E. Racine, Phys. Rev. D **78**, 044021 (2008).

[55] M. Favata, Phys. Rev. D **80**, 024002 (2009).

[56] L. Blanchet, A. Buonanno, and G. Faye, Phys. Rev. D **84**, 064041 (2011).

[57] A. Bohe, S. Marsat, and L. Blanchet, Classical Quantum Gravity **30**, 135009 (2013).

[58] S. Marsat, A. Bohé, L. Blanchet, and A. Buonanno, Classical Quantum Gravity **31**, 025023 (2014).

[59] T. Tanaka, H. Tagoshi, and M. Sasaki, Prog. Theor. Phys. **96**, 1087 (1996).

[60] M. Shibata, M. Sasaki, H. Tagoshi, and T. Tanaka, Phys. Rev. D **51**, 1646 (1995).

[61] Y. Mino, M. Sasaki, M. Shibata, H. Tagoshi, and T. Tanaka, Prog. Theor. Phys. Suppl. **128**, 1 (1997).

[62] R. Fujita, Prog. Theor. Phys. **128**, 971 (2012).

[63] M. Abramowitz and I. A. Stegun, *Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables* (Dover, New York, 1972).

[64] A. J. Brizard, arXiv:0711.4064.

[65] A. Buonanno, Y.-b. Chen, and M. Vallisneri, Phys. Rev. D **67**, 104025 (2003).

[66] K. G. Arun, A. Buonanno, G. Faye, and E. Ochsner, Phys. Rev. D **79**, 104023 (2009).

[67] B. P. Abbott et al. (Virgo, LIGO Scientific Collaboration), Phys. Rev. Lett. **116**, 061102 (2016).

[68] D. Shoemaker, Tech. Report No. LIGO-T0900288-v3, 2010.

[69] T. B. Littenberg, B. Farr, S. Coughlin, V. Kalogera, and D. E. Holz, Astrophys. J. **807**, L24 (2015).

[70] R. Smith, S. E. Field, K. Blackburn, C.-J. Haster, M. Pürrer, V. Raymond, and P. Schmidt, Phys. Rev. D **94**, 044031 (2016).

[71] J. Veitch et al., Phys. Rev. D **91**, 042003 (2015).

[72] L. Lindblom, B. J. Owen, and D. A. Brown, Phys. Rev. D **78**, 124020 (2008).