Analysis of spin precession in binary black hole systems including quadrupole-monopole interaction

# Analysis of spin precession in binary black hole systems including quadrupole-monopole interaction

## Abstract

We analyze in detail the spin precession equations in binary black hole systems, when the tidal torque on a Kerr black hole due to quadrupole-monopole coupling is taken into account. We show that completing the precession equations with this term reveals the existence of a conserved quantity at 2PN order when averaging over orbital motion. This quantity allows one to solve the (orbit-averaged) precession equations exactly in the case of equal masses and arbitrary spins, neglecting radiation reaction. For unequal masses, an exact solution does not exist in closed form, but we are still able to derive accurate approximate analytic solutions. We also show how to incorporate radiation reaction effects into our analytic solutions adiabatically, and compare the results to solutions obtained numerically. For various configurations of the binary, the relative difference in the accumulated orbital phase computed using our analytic solutions versus a full numerical solution vary from to over orbital cycles accumulated while sweeping over the orbital frequency range . This typically corresponds to a discrepancy of order radians. While this may not be accurate enough for implementation in LIGO template banks, we still believe that our new solutions are potentially quite useful for comparing numerical relativity simulations of spinning binary black hole systems with post-Newtonian theory. They can also be used to gain more understanding of precession effects, with potential application to the gravitational recoil problem, and to provide semi-analytical templates for spinning, precessing binaries.

## I Introduction and Summary

Coalescing black hole binary systems are a class of sources of gravitational radiation that can potentially be detected by observatories such as LIGO LIGO (), VIRGO VIRGO (), GEO GEO () and TAMA TAMA (). Data analysis algorithms searching for such sources are based on the method of matched filtering, which essentially consists of correlating the detector output with a bank of (theoretical) template signals. For detection to be at all possible, these templates must be accurate to 1 radian over hundreds to thousands of cycles. In the regime where the black holes are separated by a large enough distance, the post-Newtonian (PN) expansion is expected to be accurate enough to be useful in generating gravitational waveform templates, provided one carries the expansion to high enough (3.5 PN) order (e.g. see BlanchetLiving () and references therein).

While the distribution of spins for black holes in binaries is still a matter of ongoing research, it is generally thought very plausible that these spins may be large. It is therefore important to include in template banks waveforms that account for black hole spins, as they leave significant imprints on the signal, e.g. strong modulation of the waveform amplitude and phase Apostolatos (); BCV2 (); Pan (). Ignoring effects of black hole spins could very well lead to missing such binaries when analyzing interferometer data.

In order to predict the contribution of the spins to the gravitational waveform to the accuracy required by the data analysis process, the time evolution of the spins must be computed at least up to 2.5 post-Newtonian order (e.g. see Kidder (); BBF () and references therein). However in this paper we restrict attention, for simplicity, to the precession equations accurate at 2PN order, i.e. including leading order spin-orbit and spin-spin terms. To that order the spin evolution equations, which will also be loosely referred to as ”precession equations” throughout, take the functional form

 dS1,2dt = f0(t)J×S1,2+f1(t)S2,1×S1,2 (1) +termsofhigherPNorder,

where are the spin 3-vectors, is the system’s total angular momentum, and are known, implicit functions of coordinate time. In the literature the functions are generally computed assuming a model of spinning point-particles, i.e. the precession equations are derived from solving the PN expanded Einstein equations assuming a distributional stress-energy tensor that contains mass monopole and spin dipole degrees of freedom. However this procedure ignores the precession due to the coupling of the mass quadrupole moment of each spinning hole to the tidal field of the companion1. Since the mass quadrupole moment of a Kerr black hole scales in order of magnitude as , where is here the black hole mass and the magnitude of its angular momentum, the precession of the hole’s spin due its mass quadrupole contributes at the same PN order as the spin-spin term in Eq.(1), but is nevertheless often omitted in the literature, with the following notable exceptions.

Let us first point out the works of Barker and O’Connell barker () and Damour EOB (), who consider quadrupole-monopole coupling in the precession problem. Barker and O’Connell derive generic precession equations from a Hamiltonian point of view, while Damour focuses to the problem of spinning binary black holes, which allows a very compact Hamiltonian formulation at 2PN order. Next we note the contributions of Poisson Poisson () and Pan et al. Pan (), who have looked at some contributions of this term in the context of gravitational wave data analysis. In Ref.Poisson () Poisson argues that since the quadrupole-monopole precession contributes at the same PN order as spin-spin precession, the main features of the work of Apostolatos et al., who study extensively spin precession in binary black holes at leading order in post-Newtonian theory, should remain qualitatively unchanged by the addition of the quadrupole-monopole precession term, as it is a sub-leading term. In Ref.Pan (), Pan et al. devote a subsection to a comparison of gravitational wave templates with and without the quadrupole-monopole precession, their conclusion being that the contribution of the quadrupole-monopole precession can be neglected. Lastly it is worth mentioning the works of Gergely, Keresztes and Mikóczi g1 (); g2 (), who consider the quadrupole-monopole contributions to both conservative and dissipative pieces of orbital dynamics of compact binaries, and in particular provide a complete solution to the binary’s radial motion.

However in light of the analyses of PoissonPoisson () and Pan et al.Pan (), it is understandable that a complete study of the precession equations including the quadrupole-monopole contribution for black hole binaries has not yet been completed. The main purpose of this paper is to perform this analysis carefully, being motivated by the following two considerations. The first element is consistency: in black hole binaries the quadrupole-monopole term contributes at the same PN order than spin-spin coupling. Thus if spin-spin coupling is taken into account in constructing template waveforms, then so should quadrupole-monopole coupling in principle. The second motivating element is essentially academic interest in an unsolved problem, coupled with the fact that sometimes order of magnitude estimates miss important facts and more careful investigations reveal hidden features in a problem. In the analysis we present below, such a feature will be unravelled, namely the existence of a new conserved quantity in the orbit-averaged precession equations, when the quadrupole-monopole precession is taken into account.

Our analysis is organized as follows. We first review in section II the derivation of the spin precession term due to the mass quadrupole of a Kerr black hole. In section III we show that including the mass quadrupole term in the orbit-averaged precession equations leads to the discovery of a new conserved quantity2. We then show in section IV that this conserved quantity allows one to solve the orbit-averaged precession equations exactly in the case of equal masses and arbitrary spins, provided one neglects radiation reaction. In the case of unequal masses, we derive an accurate perturbative analytic solution. These solutions allow one to identify all the relevant frequencies characterizing the precession and nutation motions of each spin relative to the total angular momentum axis. Finally in section V, we derive approximate analytic solutions to the precession equations taking into account radiation damping, and compare these solutions to numerically integrated solutions.

## Ii Quadrupole-monopole interaction

In Newtonian gravity, it is well-known that a body possessing a mass quadrupole moment (or any higher order mass multipole moment) placed into a prescribed gravitational tidal field (e.g created by a companion object) experiences a torque, which leads to a non-trivial time evolution of its spin. This effect is, for example, what causes the precession of the Earth’s spin axis (precession of the equinoxes). A simple computation of the torque in Newtonian gravity, using a coordinate system that is mass-centered on the body, yields

 τi = ∫d3xρ(x,t)ϵijkxj[−∂kΦext(x,t)] (2) = ∫d3xρ(x,t)ϵijkxj ×∂k∞∑l=01l!G(t)x,

where the brackets denote the symmetric trace-free projection, is the mass density of the body, is the Newtonian potential created by an external source, and where3

 G(t)=−∂i1∂i2...∂ilΦext(x,t)∣∣x=0. (3)

The right-hand side of (3) is clearly symmetric, and Laplace’s equation ensures that is also trace-free on all indices. Using the Newtonian definition of mass multipole moments

 M(t)=∫d3xρ(x,t)x, (4)

we obtain

 τi = ϵijk∞∑l=01l!Gki1...il∫d3xρ(x,t)xjx (5) = ϵijk∞∑l=01l!Gki1...il∫d3xρ(x,t)x = ϵijk∞∑l=01l!Mji1...ilGki1...il,

where the second line follows from the fact that and differ by trace terms, which do not contribute since . The form (5) of the torque follows the syntax of Damour, Soffel and Xu DSX (). When specializing to a binary system and restricting attention solely to the quadrupole-monopole interaction, we find

 dSi1dt=3ϵijkQ1M2nr3, (6)

where is the mass quadrupole of body 1, is the mass of body 2, is a unit vector pointing from body 1 to body 2 and is orbital separation. While Eq.(5) has been revisited here in the context of Newtonian gravity, its validity in the case of a binary black hole system in the regime where the post-Newtonian expansion is valid has been proven rigorously by means of a surface-integral method by Racine in Racine (). For that situation, the mass multipole moments of each object are defined through their imprints in the far-field metric. Thus the mass quadrupole moment entering Eq.(6) in the case of a spinning black hole is well-known and is given by (e.g. see Poisson Poisson ())

 Q1=−1M1(Si1Sj1−13δijSk1Sk1), (7)

where , being the dimensionless magnitude of the black hole spin, and the unit vector along the spin of body 1. Substituting (7) into (6) we obtain

 dS1dt=3r3M2M1(^n⋅S1)^n×S1 (8)

for the quadrupole-monopole precession term. When adding this term to the precession equations typically found in the literature (see e.g. Kidder (); ThorneHartle ()), the full precession equations are

 dS1dt = 1r3{[2+32q]LN−S2+3μM1[^n⋅S0]^n}×S1, (9a) dS2dt = 1r3{[2+3q2]LN−S1+3μM2[^n⋅S0]^n}×S2,

where , , and where

 S0=(1+M2M1)S1+(1+M1M2)S2. (10)

The vector was originally introduced by Damour in the 2PN Hamiltonian for spinning binary black holes EOB (). To close the system (neglecting radiation-reaction for the time being), one uses conservation of total angular momentum4 to derive the precession equation for the (Newtonian) orbital angular , which is

 dLNdt=1r3[(12S+32S0)×LN+3μM(^n⋅S0)^n×S0], (11)

where . One may also recognize the often-encountered combination

 2Seff=12S+32S0. (12)

## Iii A new conserved quantity

From Eqs.(9), one can estimate the timescale for precession of the spins to be , being the orbital period. Thus it is customary to average the precession equations (9)-(11) over an orbital period if one is interested in the secular evolution of the spins, i.e. neglecting small fluctuations occurring over the orbital timescale. This average over orbital motion is performed using (see e.g. Schnittman JS ())

 ⟨1r3⟩ = 1a3(1−e2)3/2, (13a) ⟨ninjr3⟩ = 12a3(1−e2)3/2(δij−^LiN^LjN), (13b)

where is the semi-major axis and is the eccentricity. Defining , we then obtain the following system of evolution equations

 dS1dt = 12d3{[4+3M2M1−3μM1λ]LN+S2}×S1, (14a) dS2dt = 12d3{[4+3M1M2−3μM2λ]LN+S1}×S2, (14b) dLNdt = (14c)

where

 λ=LN⋅S0|LN|2. (15)

A key property of system (14) is that it possesses a previously unknown conserved quantity, whose identification was only possible by completing the precession equations with the quadrupole-monopole interaction. This conserved quantity turns to be , defined in Eq.(15). This can be shown by a direct computation. First, clearly is conserved from (14c), giving

 dλdt ∝ {LN⋅dS0dt+S0⋅dLNdt} ∝ {M2LN⋅(S2×S1)+M1LN⋅(S1×S2) +μS0⋅(S×LN)} = {M2LN⋅(S2×S1)+M1LN⋅(S1×S2) +(M2S1+M1S2)⋅[(S1+S2)×LN)]} = {M2[LN⋅(S2×S1)+S1⋅(S2×LN)] +M1[LN⋅(S1×S2)+S2⋅(S1×LN)]}=0.

It is important to note here that the conservation of is a property of the orbit-averaged precession equations only. Repeating the same computation using (9) and (11) instead of (14) yields . This implies, for example, that the conservation of (when averaging over orbital motion) cannot be deduced directly from the 2PN Hamiltonian of spinning binary black holes EOB (), unless the orbital degrees of freedom are averaged out.

In addition, proof (LABEL:proof) is only valid over the precession timescale. Over the radiation reaction timescale, the orbital angular momentum evolution equation must be supplemented by the following dissipative term [restricting henceforth attention to orbits with negligible eccentricy, i.e. ]

 [˙LN]rr = −325μM2d4LN. (17)

The corresponding evolution equation for is then easily shown to be

 dλdt=325μM2d4λ≡γλ, (18)

whose solution can formally be written as the following integral over orbital frequency

 λ=λ0exp[∫ωω0γ(ω)˙ωdω], (19)

where and are initial values.

## Iv Analytic solutions in absence of radiation reaction

The first step in solving the precession evolution equations is making use of the definition of total angular momentum to eliminate from (14a) and (14b). We then have

 dS1dt = 12d3{[4+3M2M1−3μM1λ]J (20a) −[3MM1−3μM1λ]S2}×S1, dS2dt = 12d3{[4+3M1M2−3μM2λ]J (20b) −[3MM2−3μM2λ]S1}×S2,

where is a constant vector. The usefulness of the conserved quantity is now quite clear: it shows that the coefficients in the coupled evolution system for and are all constants.

In the remainder of this section we solve (20) for both equal mass and unequal mass cases and arbitrary spins, generalizing the analytical analysis of Apostolatos et al. Apostolatos (), valid for equal masses, or when one of the spins is dynamically negligible.

### iv.1 Equal mass case

When , the structure of system (14) simplifies greatly, allowing in fact an exact solution. Defining the phase variable as5

 ψ=∫t012d3[7−32λ]|J|dt′, (21)

and the parameter as

 α=3|J|[4−λ14−3λ], (22)

the evolution of the spins is governed by

 dS1dψ = ^J×S1−αS2×S1, (23a) dS2dψ = ^J×S2−αS1×S2. (23b)

By adding these two equations together, we see immediately that the total spin simply precesses about , i.e.

 dSdψ=^J×S. (24)

The solution to differential equation (24) is

 S=S∥0+S⊥0cosψ+(^J×S⊥0)sinψ, (25)

where and are the projections of the initial total spin [not to be confused with the combination of spins defined in Eq.(10); henceforth will always refer to initial total spin] in directions parallel and perpendicular to respectively. Equation (25) is the first part of the exact solution for the equal mass case. We again emphasize that this solution corresponds to simple precession of the total spin about the total angular momentum axis, with precession frequency .

Consider next the vector . Its evolution is governed by

 dΔdψ = ^J×Δ−2αS2×S1 (26) = ^J×Δ−αS×Δ ≡ Ω×Δ.

Taking one more derivative with respect to yields after some algebra

 d2Δdψ2+(^J−αS)2Δ=(^J⋅Δ)(^J−2αS)+α2(S⋅Δ)S. (27)

Now , and are all constants of motion, which are thus known in terms of input initial conditions. However is not yet known. Hence before (27) can be solved explicitly, we must first solve for . This is done directly by contracting (27) with , giving

 d2(^J⋅Δ)dψ2+α2S2(^J⋅Δ)=α2(S⋅Δ)(S⋅^J). (28)

The solution is immediate

 (^J⋅Δ)=Acosσ+Bsinσ+(^S⋅Δ)(^S⋅^J), (29)

where the phase variable is defined below in Eq.(38). To fix the integration constants and , we use the initial conditions and Eq.(26) to obtain

 (^J⋅Δ) = μJS(^S0⋅Δ0)+[(^J×^S0)⋅(Δ0×^S0)]cosσ (30) +^J⋅(Δ0×^S0)sinσ,

where . Equation (30) shows that the projection of the difference of the spins along the total angular momentum axis oscillates at a frequency slower than the total spin precession frequency, since [cf. Eq.(38) below]. This implies that each individual spin exhibits a nutation motion at frequency . [Here we define nutation motion as non-trivial time evolution of the projection of a given spin vector along .] Since the nutation frequency is proportional to , it is entirely due to spin-spin coupling. Barker, Byrd and O’Connell BBOC () have computed spin nutation frequencies in generic binary systems. However it is difficult to compare their results to ours since they rely on a succession of approximations, along with a somewhat implicit, time-dependent definition of the nutation frequency. In addition it is not entirely clear how their definition of nutation motion maps onto ours. Since our result is exact, up to the same PN order as considered by BBOC (), and provides a clean, constant value for the nutation frequency along with a precise definition of nutation motion, we believe our result is an improvement over the work of BBOC (). However our result is applicable only to binary systems of Kerr black holes.

Substituting (30) into (27), one can then solve (27) using the retarded Green’s function for the simple harmonic oscillator. However it turns out to be more convenient to simply project (27) along the remaining two independent directions, which we choose to be along and , in order to match the form of solution (25) for the total spin. The remaining evolution equations are then

 d2dψ2(S⊥0⋅Δ) = −Ω2(S⊥0⋅Δ)+α2(S⊥0)2(S0⋅Δ0)cosψ−2α(S⊥0)2(^J⋅Δ)cosψ (31) = −Ω2(S⊥0⋅Δ)+α(S⊥0)2(^S0⋅Δ0)(αS−2μJS)cosψ −α(S⊥0)2{[(^J×^S0)⋅(Δ0×^S0)][cos[(1+αS)ψ]+cos[(1−αS)ψ]]} +α(S⊥0)2{[^J⋅(^S0×Δ0)][sin[(1+αS)ψ]−sin[(1−αS)ψ]]},
 d2dψ2[(^J×S⊥0)⋅Δ] = −Ω2[(^J×S⊥0)⋅Δ]+α2(S⊥0)2(S0⋅Δ0)sinψ−2α(S⊥0)2(^J⋅Δ)sinψ (32) = −Ω2[(^J×S⊥0)⋅Δ]+α(S⊥0)2(^S0⋅Δ0)(αS−2μJS)sinψ −α(S⊥0)2{[(^J×^S0)⋅(Δ0×^S0)][sin[(1+αS)ψ]+sin[(1−αS)ψ]]} −α(S⊥0)2{[^J⋅(^S0×Δ0)][cos[(1+αS)ψ]−cos[(1−αS)ψ]]}.

The solutions to these forced harmonic oscillators are

 S⊥0⋅Δ = [S⊥0⋅Δ0]cosϖ+[S⊥0⋅(^Ω0×Δ0)]sinϖ+S(1−μ2JS){(^S0⋅Δ0)(cosψ−cosϖ) (33) +12[(^J×^S0)⋅(Δ0×^S0)][cosϖ−cos(ψ−σ)(1−μJS)−cosϖ−cos(ψ+σ)(1+μJS)]
 (^J×S⊥0)⋅Δ = [(^J×S⊥0)⋅Δ0]cosϖ+[(^J×S⊥0)⋅(^Ω0×Δ0)]sinϖ+S(1−μ2JS){(^S0⋅Δ0)(sinψ−Ω−1sinϖ) (34) +12[(^J×^S0)⋅(Δ0×^S0)][δ−sinϖ−sin(ψ−σ)(1−μJS)−δ+sinϖ−sin(ψ+σ)(1+μJS)] 12[^J⋅(^S0×Δ0)][cos(ψ+σ)−cosϖ(1+μJS)+cos(ψ−σ)−cosϖ(1−μJS)]},

where

 Ω = √Ω2=(1−2αSμJS+α2S2)1/2, (35) δ± = Ω−1(1±αS), (36)

and where the phase variables and are defined as

 ϖ = ∫ψ0Ωdψ′, (37) σ = ∫ψ0αSdψ′. (38)

Here we note that the components of the difference of the spins perpendicular to the total angular momentum axis contain terms which oscillate at four different frequencies. One of these frequencies is again , and the other three are close, but not equal to the total spin precession frequency , namely and . Spin-spin coupling is again responsible for the appearance of these three additional precession frequencies.

The complete vector is then given by combining (30),(33) and (34) to give

 Δ=(^J⋅Δ)^J+(S⊥0⋅Δ)(1−μ2JS)S2S⊥0+[(^J×S⊥0)⋅Δ](1−μ2JS)S2^J×S⊥0. (39)

The individual spins are recovered from . We emphasize once more an important property of this solution: each individual spin exhibits nutation at frequency , which is typically much slower than the precession frequency of the total spin. Since the nutation frequency is directly proportional to , this clearly shows that this effect is entirely induced by spin-spin coupling. The evolution of the components of the spins transverse to is rather complicated, since it contains components oscillating at four different frequencies, which are , and . Again if one neglects spin-spin coupling, all four frequencies of the precession motion reduce to a single precession frequency, namely .

Note also that in the case where, say, , i.e. the spin of body 2 is dynamically negligible, then becomes the total spin and clearly the vector also becomes equal to the total spin. Obviously solution (25) for the total spin is still valid in that situation. It is straightforward to check that by setting in (39), one recovers for all time.

### iv.2 General case

Unlike the equal mass case, the general case of arbitrary spins and unequal masses unfortunately does not allow an exact solution in terms of elementary functions. We first prove this statement. The spin evolution equations in the general case are

 dS1dψ = β1^J×S1−α1S2×S1, (40a) dS2dψ = β2^J×S2−α2S1×S2, (40b)

where

 α1,2 = 1|J|[6M1,2(M−μλ14−3λ)], (41) β1,2 = 4+3M2,1M1,2−3μM1,2λ7−32λ. (42)

Now consider the evolution equations satisfied by and . These are

 d(^J⋅S1,2)dψ = ±α1,2^J⋅(S1×S2), (43) d(S1⋅S2)dψ = β−^J⋅(S1×S2), (44)

where . We see that all right-hand sides in the above equations are proportional to the same function, namely . Defining as

 ρ=∫ψ0^J⋅(S1×S2)dψ′, (45)

we then have

 ^J⋅S1,2 = ±α1,2ρ+(^J⋅S1,2)0, (46) S1⋅S2 = (β−)ρ+(S1⋅S2)0, (47)

where the subscript denotes initial value. By taking a derivative of, say, Eq.(44) and using expressions (46)-(47), one can show that obeys the following differential equation

 d2ρdψ2 = −3β−α1α2ρ2−[β2−+2β−^J⋅ζ−0+(ζ+0)2]ρ (48) −β−[(^J×S1)⋅(^J×S2)]0 +ζ+0⋅[^J×(S1×S2)]0,

where . Thus evolves as an anharmonic oscillator with a cubic term in its potential. In the equal mass case and reduces to a simple harmonic oscillator, thus allowing a solution in terms of trigonometric functions. In the general case one can solve (48) in terms of an elliptic integral6, more specifically in terms of the Weierstrass elliptic function (see e.g. WW ()). However here, for sake of transparency and simplicity, we resort to a perturbative method to solve (48), based on the existence of a small parameter that will be identified shortly. Equation (48) is of the form

 d2ρdψ2=−C2ρ2−C1ρ+C0, (49)

where are all constants, which can be read off directly from (48). Defining and the dimensionless variable , we have

 d2udψ2 = −S1S2C2u2−C1u+C0S1S2 (50) ≡ −B2u2−B1u+B0.

Now the key element to realize is that is a small number, as it scales as

 B2 ∼ (α1S1)(α2S2)∼M2J2S1M1S2M2 (51) ∼ M1M2M2J2χ1χ2 ∼ ηM4J2χ1χ2,

where (41) has been used and where are dimensionless spins and is the symmetric mass ratio.

Let us first consider the case where , i.e. the extreme mass ratio limit. We adopt the convention here that , so . In that case the total angular momentum is dominated by the total spin since . Therefore we have which gives

 B2→ηχ2/χ1→0, (52)

since is finite.

When masses are comparable in the regime where post-Newtonian theory is applicable, the total angular momentum is dominated by the orbital angular momentum, which yields . This then yields

 B2∼χ1χ2ηMr∼χ1χ2ηv2. (53)

Since is typically a number of order unity when masses are comparable, scales as a 1PN correction term compared to the numbers and , which are generically of order unity.

Provided that remains small enough, the anharmonic term in the equation of motion for generically remains a small correction to the harmonic driving term for all time. We shall assume this is the case in what follows, and verify a posteriori that this assumption indeed works quite well. We next define the variable

 x=u−A0, (54)

where

 A0=−B1+√B21+4B0B22B2 (55)

is a constant chosen so that the evolution equation for takes the form

 d2xdψ2=−εx2−w2x, (56)

where and

 w2=B1+2B2A0=√B21+4B0B2. (57)

We now solve (56) by treating the anharmonic term as a small perturbation, with being our expansion parameter. Now it is well-known (see e.g. Goldstein Goldstein ()) that the corrections to the frequency of an oscillator with a cubic perturbation of order in its potential scale as . Thus it is necessary to solve (56) to order accuracy in order to capture this important effect. To that order the solution is

 x(ψ) = −ε[x20+v202w2+εx303w4]+[x0+εx20+2v203w2+ε2x0(29x20−55v20)144w4]cos(Φ) (58) +[v0−ε2x