Constructing Gravitational Waves from Generic Spin-Precessing Compact Binary Inspirals

TL;DR

This paper develops the first closed-form frequency-domain gravitational waveforms for generic spin-precessing compact binary inspirals, enabling more accurate and computationally efficient detection and analysis of gravitational waves from such systems.

Contribution

It introduces a novel method to construct analytic frequency-domain waveforms for precessing binaries, combining post-Newtonian precession equations with shifted uniform asymptotics.

Findings

Waveforms are accurate for 99.2% of configurations at SNR 10.

Waveforms are suitable for unbiased parameter estimation.

Method reduces computational cost for gravitational wave data analysis.

Abstract

The coalescence of compact objects is one of the most promising sources of gravitational waves for ground-based interferometric detectors, such as advanced LIGO and Virgo. Generically, com- pact 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.

Figures

Tables

Table 1. Table 1: 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 50 50 Hz and in a frame were the orbital angular momentum is aligned with the z 𝑧 z axis.

		NSNS		BHNS		BHBH		HSNSBH
$m_1$		1.6$M_{\odot }$		10$M_{\odot }$		10$M_{\odot }$		10$M_{\odot }$
$m_2$		1.4$M_{\odot }$		1.4$M_{\odot }$		5$M_{\odot }$		1.4$M_{\odot }$
$\cos {\theta }_L$		1		1		1		1
${\varphi }_L$		0		0		0		0
$\cos {\theta }_1$		0.5		0.5		0.5		0.5
${\varphi }_1$		1.2		1.2		1.2		1.2
${\chi }_1$		0.08		0.7		0.7		0.7
$\cos {\theta }_2$		0.7		0.7		0.7		0.7
${\varphi }_2$		2.5		2.5		2.5		2.5
${\chi }_2$		0.1		0.1		0.6		0.6

Table 2. Table 2: Percentages of subthreshold systems encountered in our analysis for each type of system.

Waveform		Threshold		NSNS		BHNS		BHBH		HSNSBH
RWF		$0.965$		$0.06\%$		$0.33\%$		$1.85\%$		$1.14\%$
RWF		$0.994$		$0.3\%$		$1.6\%$		$10.4\%$		$9.2\%$
FWF		$0.965$		$0.06\%$		$0.33\%$		$1.85\%$		$1.26\%$
FWF		$0.994$		$0.3\%$		$1.6\%$		$10.7\%$		$9.9\%$

Equations

$$\dot{\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\left({\mathbf{S}}_2\times \hat{\mathbf{L}}\right)$$

$$+\mathcal{O}(v^7),$$

$${\dot{\mathbf{S}}}_1$$

$$+\frac{v^6}{2}{\mathbf{S}}_2\times {\mathbf{S}}_1+\mathcal{O}(v^7),$$

$${\dot{\mathbf{S}}}_2$$

$$+\frac{v^6}{2}{\mathbf{S}}_1\times {\mathbf{S}}_2+\mathcal{O}(v^7).$$

$$\dot{v}$$

$$T_{\text{mbox}pr}$$

$$T_{\text{mbox}rr}$$

$$\xi \equiv (1+q){\mathbf{S}}_1\cdot \hat{\mathbf{L}}+(1+q^{-1}){\mathbf{S}}_2\cdot \hat{\mathbf{L}}.$$

$$\cos {\theta }_L=\hat{\mathbf{J}}\cdot \hat{\mathbf{L}}=\frac{J^2+L^2-S^2}{2JL},$$

$$\mathbf{L}(S;\mathbf{\lambda })=L[\sin {\theta }_L,0,\cos {\theta }_L].$$

$$\mathbf{S}(S;\mathbf{\lambda })=\mathbf{J}-\mathbf{L}=[-L\sin {\theta }_L,0,J-L\cos {\theta }_L].$$

$${{\mathbf{S}}_1}^{\prime }=S_1[\sin {\theta }^{\prime }\cos {\varphi }^{\prime },\sin {\theta }^{\prime }\sin {\varphi }^{\prime },\cos {\theta }^{\prime }].$$

$$\cos {\theta }^{\prime }$$

$$\cos {\varphi }^{\prime }$$

$$-4q{S}^{2}L\xi \}/[(1-q^2)A_1{A}_2{A}_3{A}_4],$$

$$A_1$$

$$A_2$$

$$A_3$$

$$A_4$$

$${\mathbf{S}}_1(S;\mathbf{\lambda })=\mathbb{R}(\hat{\mathbf{y}},{\theta }_S){{\mathbf{S}}_1}^{\prime },$$

$$\cos {\theta }_S=\hat{\mathbf{S}}\cdot \hat{\mathbf{J}}=\frac{J^2+S^2-L^2}{2JS}.$$

$${\mathbf{S}}_2(S;\mathbf{\lambda })=\mathbf{J}-\mathbf{L}-{\mathbf{S}}_1.$$

$${\left(\frac{d{S}^2}{dt}\right)}^2=-A^2\left(S^6+B{S}^4+C{S}^2+D\right),$$

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

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

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

$$m=\frac{S_{+}^2-S_{-}^2}{S_{+}^2-S_3^2}.$$

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

$$\text{sign}(\sin {\varphi }^{\prime })=\text{sign}({\mathbf{S}}_1\cdot {\mathbf{y}}^{\prime }),$$

$$\text{sign}(\sin {\varphi }^{\prime })$$

$$=\text{sign}\left[\mathrm{sn}(\psi ,\mathrm{m})\mathrm{cn}(\psi ,\mathrm{m})\right],$$

$$\frac{d{\varphi }_z}{dt}$$

$$-\frac{3(1+q)}{2q{A}_1^2{A}_2^2}\left(1-\xi v\right)[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 )]\}.$$

Full text

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 Spruit and Phinney (1998). 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 Spruit and Phinney (1998) and tilting the orbital plane, since the kick’s direction is typically correlated with the spin of the exploding star Wang et al. (2006, 2007); the various angular momenta become misaligned Kalogera (2000); Gerosa et al. (2014). 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 Kalogera (2000); Gerosa et al. (2014).

Few mechanism 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 Kalogera (2000). This possibility is remote, with models and data suggesting that the kick is actually aligned with the spin angular momentum Wang et al. (2006, 2007); Kaplan et al. (2008). A mechanism for realignment is through torques exerted by a circumbinary accretion disk on BH binaries Bogdanovic et al. (2007). 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 Rodriguez et al. (2016). Finally, we note that the binary might undergo spin-orbit resonances during its evolution Schnittman (2004). These resonances, however, do not align the spins with the orbital angular moment, but rather maintain certain special (“resonant”) precessing configurations Schnittman (2004); Gerosa et al. (2014); Kesden et al. (2014); Gerosa et al. (2015a, b).

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 Chatziioannou et al. (2015a); Chatziioannou et al. (2014).

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 approximation111The SPA is the leading-order term in the asymptotic expansion of a Fourier integral through the method of steepest descent Droz et al. (1999); Bender and Orszag (1999). (SPA) and rendering it nonapplicable Klein et al. (2014).

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 Blanchet (2014). 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 order222A term is of $A$PN order is proportional to $(u/c)^{2A}$ relative to its controlling factor, where $u$ is some characteristic velocity and $c$ the speed of light. 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 Apostolatos et al. (1994); Apostolatos (1995, 1996a, 1996b); Hannam et al. (2014). The ensuing motion is simple precession and the resulting waveform is ideal for BHNS systems Lundgren and O’Shaughnessy (2013). The second representation is based on the effective-one-body formulation of the general relativistic two-body problem Buonanno and Damour (1999a); Buonanno and Damour (1999b); Damour (2001). The resulting waveform is ever improving through fits of its nonprecessing part to numerical relativity simulations Taracchini et al. (2014); Pan et al. (2014); PhysRevD.93.044007 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 Boyle et al. (2011); Schmidt et al. (2011) to compute a simpler waveform PhysRevD.93.044007 and map it back to the source frame Schmidt et al. (2012, 2014). 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 Hannam et al. (2014)333It is worth emphasizing, thought, that the first GW detection Abbott et al. (2016a); Abbott et al. (2016b) did not suffer from such systematics due to its orientation and minimal precession Abbott et al. (2016a); Abbott et al. (2016c)..

The final representation of GWs from inspiraling spin-precessing systems was through a multiple scale analysis (MSA) Bender and Orszag (1999), 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 Klein et al. (2013) and slowly spinning Chatziioannou et al. (2013a) systems. The latter are accurate representations of NSNS inspirals, both for detection and parameter estimation Chatziioannou et al. (2014); Chatziioannou et al. (2015a, b). 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. Kesden et al. (2014); Gerosa et al. (2015b, a). 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 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. Klein et al. (2014), 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. Chatziioannou et al. (2016), 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 Klein et al. (2014), 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:

• •

Vectors are written in boldface, with components $\bm{A}=[A_{x},A_{y},A_{z}]$ and magnitude $A$. Unit vectors are denoted with a hat, e.g. $\bm{\hat{A}}$.

• •

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}$.

• •

The Newtonian orbital angular momentum of the system is $\bm{L}$, the spin angular momentum of each body is $\bm{S}_{A}$, and the total angular momentum is $\bm{J}=\bm{L}+\bm{S}_{1}+\bm{S}_{2}$. The dimensionless spin parameter of each object is $\chi_{A}\equiv S_{A}/m_{A}^{2}$ with $A\in\{1,2\}$.

• •

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}$.

• •

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 1 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 $\bm{L}$, while $\theta_{A}$ and $\phi_{A}$ are the polar angles of $\bm{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 orbit444Orbit-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., the precession equations governing the conservative evolution of the orbital and spin angular momenta are Thorne and Hartle (1984); Barker and O’Connell (1979); Buonanno et al. (2006)555The precession equations used here are only strictly valid for BHs, as for NSs they acquire additional terms describing the quadrupole moment of the bodies Poisson (1998). However, the extra terms are degenerate with the spins Krishnendu et al. (2017) and difficult to measure.

[TABLE]

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 Poisson (2004); Chatziioannou et al. (2013b). The magnitude $L$ is related to the evolution of the orbital frequency $\omega$, and the PN expansion parameter $v$, leading to

[TABLE]

The coefficients $\{g_{n},g^{\ell}_{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 $\bm{L},\bm{S_{1}}$ and $\bm{S_{2}}$. We use only the leading PN order expressions in each interaction666Spin-orbit corrections can be found in Faye et al. (2006); Bohe et al. (2013a), spin-spin in Bohé et al. (2015), and spin-cubed in Marsat (2015).. 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 Bohé et al. (2015); Racine et al. (2009). In principle, we could have included the spin-orbit corrections. However, as explained later, our solution makes use of a certain quantity Racine (2008) 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.5$PN order in all spin interactions Racine et al. (2009), $4$PN in linear-in-spin terms Racine et al. (2009); Favata (2009); Blanchet et al. (2011); Bohe et al. (2013a, b); Marsat et al. (2014) and $22$PN order in the point particle limit, neglecting spins and BH absorption effects Tanaka et al. (1996); Shibata et al. (1995); Mino et al. (1997); Fujita (2012). In our analysis we keep terms in Eq. (4) to $3.5$PN 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 Bohé et al. (2015) we can include it in our model.

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

[TABLE]

while the later evolve on the radiation reaction time scale

[TABLE]

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

Recently, Kesden et al. Kesden et al. (2014) 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 Gerosa et al. (2015a).

Here we take a different approach. Rather that precession-averaging Eqs. (1)-(3) and numerically accounting for Eq. (4), we make explicit use of the fact that $T_{{\mbox{\tiny pr}}}/T_{{\mbox{\tiny 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_{{\mbox{\tiny pr}}}/T_{{\mbox{\tiny 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 Kesden et al. (2014).

A precessing binary has a total of 9 degrees of freedom arising from the $3$ components of $3$ Newtonian vectors $(\bm{L},\bm{S_{1}},\bm{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=|\bf{S_{1}}+\bf{S_{2}}|$.

The conserved quantities are $\bm{\lambda}\equiv(S_{1},S_{2},L,J,\bm{\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 Racine (2008)

[TABLE]

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 $\bm{\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)$.

III.1 Precession in a noninertial frame

The identification of $\bm{\hat{J}}$ as a conserved quantity suggests a coordinate frame where $\bm{\hat{z}}=\bm{\hat{J}}$ (see Fig. 1). We further pick the $x$ and $y$ axes to be precessing around $\bm{\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

[TABLE]

from the $\bf{\hat{z}}$ axis.

This allows us to express $\bm{L}$ as

[TABLE]

The total spin angular momentum then is

[TABLE]

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

[TABLE]

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

[TABLE]

where

[TABLE]

In the original unprimed system

[TABLE]

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

[TABLE]

Once we have $\bm{S_{1}}$ in the original unprimed system, then

[TABLE]

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^{\prime}$ 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):

[TABLE]

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 systems777In 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..

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

[TABLE]

The solution to this equation is

[TABLE]

where $\rm{sn}$ is a Jacobi elliptic function (see Sec. 16 of Abramowitz and Stegun (1972) for a detailed introduction to the Jacobi elliptic functions, and Brizard (2007) 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

[TABLE]

and

[TABLE]

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^{2}_{3}=S^{2}_{+}$ is when $S^{2}_{+}=S^{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

[TABLE]

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 $\bm{\hat{z}}$ is the sign of $\sin{\phi^{\prime}}$. Equation (11) implies that

[TABLE]

which after some algebra can be shown to be equivalent to

[TABLE]

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

III.2 Precession in an inertial frame

All angular momenta so far have been expressed in a noninertial frame that precesses around $\bm{\hat{J}}$. An Euler rotation of $\bm{L},\bm{S_{1}}$, and $\bm{S_{2}}$ around $\bm{\hat{z}}$ by some angle $\phi_{z}$ and substitution into the precession equations yields Kesden et al. (2014)

[TABLE]

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

[TABLE]

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

[TABLE]

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 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 $\bm{\hat{J}}=\bm{\hat{z}}$. In summary, at some initial time:

• •

The orbital angular momentum $\bm{L}$ is given by Eq. (9), which depends on the angle $\theta_{L}$ given in Eq. (8). The latter depends on $S$, which varies on the precession timescale as described in Eq. (23);

• •

The spin angular momentum of the heavier body $\bm{S_{1}}$ is given in Eq. (18), which depends on the angle $\theta_{S}$ given in Eq. (19) as well as on $S_{1}^{\prime}$ given in Eq. (11) in terms of the angles $(\theta^{\prime},\phi^{\prime})$ of Eqs. (12) and (13). All of these depend on $S$, which again is described by Eq. (23);

• •

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 $\bm{\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 Bender and Orszag (1999).

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_{{\mbox{\tiny pr}}}$ to the radiation reaction time scale $T_{{\mbox{\tiny 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.

IV.1 Choice of an inertial frame

The precession solution of Sec. III was built around the assumption that $\bm{\hat{J}}$ is conserved and aligned with $\bm{\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 Apostolatos et al. (1994) that in the context of simple precession $(\bm{S_{2}}=0)$ the variation of $\bm{\hat{J}}$ in a precession cycle averages out. Here we show that this is approximately true for generic precession as well Buonanno et al. (2003).

Equations (1)-(3) imply

[TABLE]

and after some algebra we can rewrite this as

[TABLE]

Averaging over $\phi_{z}$ we find

[TABLE]

This averaging induces an error in $\dot{\bm{\hat{J}}}$ that is ${\cal{O}}(T_{{\mbox{\tiny pr}}}/T_{{\mbox{\tiny rr}}})$, or \dot{\bm{\hat{J}}}-\langle\dot{\bm{\hat{J}}}\rangle_{\text{\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 ${\cal{O}}(v^{3})$. They are therefore of the same order as the averaging error, and can be neglected.

Working to this order we have

[TABLE]

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

[TABLE]

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{\bm{\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{\bm{\hat{J}}}\rangle_{\text{\phi_{z}}}=0; radiation reaction changes the magnitude of $\bm{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 1) and plot the components of $\bm{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 that $J_{z}$ and oscillate around [math], 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 $\bm{\hat{J}}$ and align it with $\bm{\hat{z}}$. This 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.

IV.2 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 Kesden et al. (2014)

[TABLE]

This can be integrated exactly to yield

[TABLE]

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, $\left\langle S^{2}\right\rangle_{\text{{\mbox{\tiny pr}}}}$ is constant when ignoring high-order PN effects, and the integral of Eq. (40) can be calculated to give

[TABLE]

The quantity $\left\langle S^{2}\right\rangle_{\text{{\mbox{\tiny pr}}}}$ can be computed from Eq. (23):

[TABLE]

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

[TABLE]

which together with Eq. (25) yields

[TABLE]

and

[TABLE]

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

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

[TABLE]

To verify that this approximate $\bm{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_{x}$ and $J_{y}$ are identically zero, so we omit them. The bottom panel shows the fractional error in the magnitude of the total angular momentum when approximated by Eq. (41). The maximum discrepancy in the magnitude $J$ is of $\mathcal{O}\left(10^{-2}\right)$ in those particular examples indicating both that Eq. (41) is accurate and that setting $J_{x}$ and $J_{y}$ equal to zero is justified.

IV.3 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^{2}(t)=S^{2}(t_{{\mbox{\tiny pr}}},t_{{\mbox{\tiny rr}}})$, where $t_{{\mbox{\tiny pr}}}$ denotes variation on the precession time scale, while $t_{{\mbox{\tiny rr}}}=\epsilon t_{{\mbox{\tiny pr}}}$ denotes variations on the radiation reaction time scale, with $\epsilon$ a bookkeeping parameter.

Expanding $S^{2}$ as

[TABLE]

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

[TABLE]

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_{{\mbox{\tiny rr}}}$:

[TABLE]

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

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

[TABLE]

where we keep terms of $\mathcal{O}\left(\epsilon\right)$ by taking the derivative with respect to $t$ rather than $t_{{\mbox{\tiny 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

[TABLE]

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^{2}_{+}$ and $S^{2}_{-}$. For the NSNS and BHBH systems, the analytic phase $\psi$ 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 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 future888For example, if and when LIGO’s sensitivity increases, so will its

requirement for more accurate waveforms. we can obtain it by carrying out the expansion of Eq. (51) to higher order.

IV.4 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

[TABLE]

and expand the precession angle as

[TABLE]

The reason $\phi_{z}$ includes a term of $\mathcal{O}\left(\epsilon^{-1}\right)$ 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}\left(\epsilon^{-1}\right)$

[TABLE]

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

[TABLE]

and averaging over $t_{{\mbox{\tiny pr}}}$ we find

[TABLE]

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

[TABLE]

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

[TABLE]

where

[TABLE]

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.

IV.4.1 Leading order MSA

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

[TABLE]

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 $\left\langle\dot{\phi}_{z}\right\rangle_{\text{{\mbox{\tiny pr}}}}$. We write

[TABLE]

where on the third line we average over precession. Unfortunately, no closed form expressions exist for $\left\langle\,\text{sn}^{2}(\psi,m)\right\rangle_{\text{{\mbox{\tiny pr}}}}$ and $\left\langle\,\text{sn}^{4}(\psi,m)\right\rangle_{\text{{\mbox{\tiny pr}}}}$ for arbitrary m. We can, however, calculate these averages as an expansion in $m\ll 1$ since, as already discussed, $m\sim{\cal{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

[TABLE]

where we have defined

[TABLE]

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 $\left\langle\Omega_{z}\right\rangle_{\text{{\mbox{\tiny 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

[TABLE]

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

[TABLE]

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

IV.4.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 $\left\langle\Omega_{z}\right\rangle_{\text{{\mbox{\tiny pr}}}}$ does not depend on the precession time scale $t_{{\mbox{\tiny pr}}}$, and the result is $\left\langle\Omega_{z}\right\rangle_{\text{{\mbox{\tiny pr}}}}t_{{\mbox{\tiny pr}}}$. In that expression, we choose for convenience to substitute $t_{{\mbox{\tiny pr}}}=\psi/\dot{\psi}$.

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

[TABLE]

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.

IV.4.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 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 Blanchet (2014); Apostolatos et al. (1994); Arun et al. (2009); Boyle et al. (2011); Lundgren and O’Shaughnessy (2013):

[TABLE]

where

[TABLE]

are the antenna pattern functions, $h_{+,\times}$ are the GW polarization states, $(\theta_{N}^{\prime},\phi_{N}^{\prime})$ are the polar angles of $\bm{\hat{N}}$ in a frame tied to the arms of the detector with $\bm{\hat{z}}^{\prime}$ the normal to the detector plane, and $\psi_{p}$ is given by

[TABLE]

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

The polarization states can be decomposed into a spin-weighted spherical harmonic basis Arun et al. (2009); Boyle et al. (2011); Lundgren and O’Shaughnessy (2013)

[TABLE]

where

[TABLE]

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

[TABLE]

where ${}_{s}Y_{lm}$ are the spin-weighted spherical harmonics, the amplitudes $h^{lm}$ are in Blanchet (2014), $D^{l}_{m,m^{\prime}}$ are the Wigner D-matrices, the angles $\theta_{L}$ and $\phi_{z}$ are the spherical angles of $\bm{\hat{L}}$ in the frame where $\bm{\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 Klein et al. (2014)

[TABLE]

where $t_{m}$ and $T_{m}$ are defined by

[TABLE]

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

[TABLE]

for $p\in\{0,\ldots,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 1 using only the leading $(\ell=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.

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\pi]$, 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

[TABLE]

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 factors999A fitting factor is the faithfulness maximized over all model parameters, quantifying a model’s suitability for detection., selecting a value for the faithfulness that is “good enough” is not straightforward. The nominal fitting factor threshold of $0.965$ corresponds to a $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 Abbott et al. (2016d), 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 Klein et al. (2014) 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

[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 Shoemaker (2010).

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 $1M_{\odot}$ to $20M_{\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 $4\text{th}$ 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 2 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. (150) and (151)] 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.

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

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

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

Finally, in Sec IV.4 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:

• •

The black solid line corresponds to the fully analytic waveform also studied in Fig. 6, with approximately $10\%$ of the systems below $0.994$.

• •

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 $\sim 4\%$.

• •

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 $\sim 2\%$ of subthreshold systems.

• •

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 $\sim 1\%$.

• •

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 $\sim 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 $\sim 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 Apostolatos et al. (1994). 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 $\bm{L}\cdot\bm{N}$ on our results. It is well known that binaries observed approximately edge-on ($\bm{L}\cdot\bm{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 Abbott et al. (2016c). We should note, however, that edge-on systems are less likely to be detected by aLIGO due to selection effects Littenberg et al. (2015).

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. Kesden et al. (2014) 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) Klein et al. (2014). 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 Smith et al. (2016). 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 Hannam et al. (2014); Schmidt et al. (2012, 2014). 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 Pan et al. (2014) due to the analytic treatment of the inspiral dynamics.

Acknowledgements.

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 PHY-1306702. N. C. and N. Y. acknowledge 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 $\dot{v}$

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

[TABLE]

with all other terms vanishing. The coefficients $\{a_{i},b_{i}\}$ are given in Appendix A of Chatziioannou et al. (2013a). 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 that the order to which we know the $\dot{v}$ expansion. Explicitly, we use

[TABLE]

where $v_{0}$ corresponds to the value of $v$ at the initial time, $\xi,S^{2}_{{\mbox{\tiny av}}},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

[TABLE]

The coefficients of Eq. (30) are

[TABLE]

The coefficients of Eq. (31) are

[TABLE]

where

[TABLE]

Appendix C Coefficients of $\psi$

The coefficients in Eq. (51) are

[TABLE]

where

[TABLE]

Appendix D Coefficients of $\phi_{z}$

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

[TABLE]

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

[TABLE]

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

[TABLE]

where we have defined

[TABLE]

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{{\mbox{\tiny 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 $\left\langle{\phi}_{z}\right\rangle_{\text{{\mbox{\tiny 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 $\bm{S}$ is approximately (anti)aligned with $\bm{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 $\bm{S}$ and $\bm{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 situation101010We have verified that this is the case both analytically and numerically..

We argue that even though Eq. (62) is well behaved in both the small mass ratio and the (anti)alignment between $\bm{S}$ and $\bm{L}$ limit, the same need not be true for its PN expansion. Consider the following function111111The similarity between our toy function and J given in Eq. (41) is not accidental.

[TABLE]

and its expansion around $x=0$

[TABLE]

Clearly, $h(x)$ is finite as $h_{2}\rightarrow 0$. However, $h^{{}_{\text{exp}}}(x)$ is not, and the $h_{2}\rightarrow 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 $\bm{S}$ and $\bm{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 $\bm{S}$ and $\bm{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 $\left\langle\phi_{z}\right\rangle_{\text{{\mbox{\tiny 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 $\left\langle{\phi}_{z}\right\rangle_{\text{{\mbox{\tiny 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:

[TABLE]

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 $\bm{S}$ and $\bm{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 $\bm{L}$ can be calculated by solving

[TABLE]

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

[TABLE]

where

[TABLE]

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

[TABLE]

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

[TABLE]

and

[TABLE]

where the $\Omega_{z,i}$’s, $f_{d}$ and $h_{d}$ are given in Appendix D and

[TABLE]

The first correction to MSA is given by

[TABLE]

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

[TABLE]

Appendix G Faithfulness requirement

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

[TABLE]

In the high SNR regime, a typical waveform sample from the posterior distribution function Veitch et al. (2015) is given by

[TABLE]

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

[TABLE]

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

[TABLE]

where

[TABLE]

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

[TABLE]

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}\simeq\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. (178) 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 an 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

[TABLE]

Thus

[TABLE]

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

[TABLE]

where $d$ is the data. The peak of the posterior, evaluated at the best-fit parameters is

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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 Lindblom et al. (2008) that were written in terms of requirements on the GW amplitude and phase to requirement on the faithfulness.