Statistical studies of Spinning Black-Hole Binaries

2.2 MB·English

by  Carlos O. Lousto

Statistical studies of Spinning Black-Hole Binaries
Carlos O. Lousto, Hiroyuki Nakano, Yosef Zlochower, Manuela Campanelli
Center for Computational Relativity and Gravitation,
and School of Mathematical Sciences, Rochester Institute of Technology,
85 Lomb Memorial Drive, Rochester, New York 14623
Westudythestatisticaldistributionsofthespinsofgenericblack-holebinariesduringtheinspiral
and merger, as well as the distributions of the remnant mass, spin, and recoil velocity. For the
inspiral regime, we start with a random uniform distribution of spin directions S(cid:126) and S(cid:126) over the
1 2
sphereandmagnitudes|S(cid:126) /m2|=|S(cid:126) /m2|=0.97fordifferentmassratios,whereS(cid:126) andm arethe
1 1 2 2 i i
spin-angular momentum and mass of the ith black hole. Starting from a fiducial initial separation
of r = 50M, we perform 3.5-post-Newtonian-order evolutions down to a separation of r = 5M,
i f
where M = m +m , the total mass of the system. At this final fiducial separation, we compute
1 2
the angular distribution of the spins with respect to the final orbital angular momentum, L(cid:126). We
0 perform 164 = 65536 simulations for six mass ratios between q = 1 and q = 1/16 and compute
1 the distribution of the angles L(cid:126)ˆ·∆(cid:126)ˆ and L(cid:126)ˆ·S(cid:126)ˆ, directly related to recoil velocities and total angular
0
momentum. We find a small but statistically significant bias of the distribution towards counter-
2
alignmentofbothscalarproducts. Apost-Newtoniananalysisshowsthatradiation-reaction-driven

n dissipative effects on the orbital angular momentum lead to this bias. To study the merger of
a black-hole binaries, we turn to full numerical techniques. In order to make use of the numerous
J simulations now available in the literature, we introduce empirical formulae to describe the final

1 remnant black hole mass, spin, and recoil velocity for merging black-hole binaries with arbitrary
2 massratiosandspins. Ourformulaearebasedonthepost-Newtonianscaling,tomodeltheplunge
phase, with amplitude parameters chosen by a least-squares fit of recently available fully nonlinear

] numericalsimulations,supplementedbyinspirallossesfrominfinitytotheISCO.Wethenevaluate
c
thoseformulaeforrandomlychosendirectionsoftheindividualspinsandmagnitudesaswellasthe
q
binary’smassratio. Thenumberofevaluationshasbeenchosensuchthatthereare10configurations
-
r per each dimension of this parameter space, i.e. 107. We found that the magnitude of the recoil
g velocity distribution decays exponentially as P(v) ∼ exp(−v/2500 km s−1) with mean velocity
√
[ < v >= 630 km s−1 and standard deviation <v2 >−<v>2 = 534 km s−1, leading to a 23%
4 probability of recoils larger than 1000 km s−1, and a highly peaked angular distribution along the
v final orbital axis. The studies of the distribution of the final black-hole spin magnitude show a
7 universal distribution highly peaked at Sf/m2f = 0.73 and a 25◦ misalignment with respect to
9 the final orbital angular momentum, just prior to full merger of the holes. We also compute the
1 statistical dependence of the magnitude of the recoil velocity with respect to the ejection angle.
3 The spin and recoil velocity distributions are also displayed as a function of the mass ratio. We
. finallyalsocomputetheeffectsoftheobserverorientationwithrespecttotherecoilvelocityvector
0
to take into account the probabilities to measure a given redshifted (or blueshifted) radial velocity
1
of accretion disks with respect to host galaxies.
9
0
: PACSnumbers: 04.25.Dm,04.25.Nx,04.30.Db,04.70.Bw
v
i
X
I. INTRODUCTION effect [4], a repulsive spin-orbit interaction, that delays
r the merger of black-hole binaries (BHB) when the spins
a
are aligned with the orbital angular momentum, and
Astrophysical black-hole (BH) binaries are character-
simultaneously causes the system to radiate excess an-
ized by the mass ratio q = m /m ≤ 1 of the smaller to
1 2
gular momentum, leading to a remnant BH with sub-
larger BH, where m is the mass of BH i, the total mass
i
maximal spin. The same mechanism produces an addi-
M =m +m ,eccentricitye(assumedtobeverysmall),
1 2
tionalattractiveeffectwhenthespinsarecounter-aligned
and spins S(cid:126) (where S /m2 <1). In addition, it is often
i i i withtheorbitalangularmomentum,leadingtoaprompt
convenient to parameterize the binary with the symmet-
merger. Thus the radiation of angular momentum and
ric mass ratio η =q/(1+q)2, specific spins α(cid:126) =S(cid:126) /m2,
i i i energy is asymmetric with respect to the relative orien-
total spin S(cid:126) = S(cid:126)1 +S(cid:126)2, ∆(cid:126) = M(S(cid:126)2/m2 −S(cid:126)1/m1), and tations of the total spin angular momentum vector and
orbital angular momentum L(cid:126). the orbital angular momentum.
The relevance of the spins to the dynamics of black- When the spins are not exactly aligned or counter
hole (BH) mergers was recognized soon after the break- aligned, new effects appear (the hangup effect is still
throughinnumericalrelativity[1–3]thatallowedforthe present). Precession of the spins is important dynami-
long term, stable numerical evolution of such systems. cally because it cause the orbital plane to strongly pre-
Notable examples of the early findings are the ‘hangup’ cessjustpriortomerger[5]. Thefinalspinofthemerged
2
holecanflipwithrespecttothedirectionsoftheindivid- available,butthereremainssignificantuncertaintiescon-
ual ones, mainly due to the addition of the orbital angu- cerning their accuracy outside this range of parameters.
lar momentum [5]. While the spin-orbit coupling leads In this paper we propose a set of formulae that incorpo-
to strong precessional effects near merger, the magni- rate the benefits of both approaches in a unified way.
tudes of the spins are not affected to the same degree. Due to the large astrophysical interest of computing
In particular spin-orbit interactions are too weak to in- remnantrecoilvelocities,themodelingofrecoilvelocities
duce the binary to corotate (or maintain corotation of followed an independent path, particularly since the dis-
an initially corotating binary) at the last stages of the covery[7,24]thatthespinsoftheblackholesplayacru-
merger because the timescale for the radiation driven cial role in producing recoils of up to 4000 km s−1. The
inspiral is much smaller than the spin-orbit interaction importanceofmodelingtherecoilvelocitiesasafunction
timescale [6]. of the astrophysical parameters of the progenitor binary
Numerous other papers have studied different spin ef- was quickly realized [7, 8, 24].
fects, such as the large recoil velocities acquired by the The news that the merger of binary black holes can
remnant of the merger of two spinning black holes [7–12] produce recoil velocities up to 4000 km s−1, and hence
and long term evolutions of generic BHBs (i.e. unequal allow the remnant to escape from major galaxies, led
mass and unequal, randomly-oriented spins) [13, 14] to to numerous theoretical and observational efforts to find
citeafewofthenearlyonehundredpaperspublishedon traces of this phenomenon. Several studies made pre-
the subject since 2006. dictions of specific observational features of recoiling su-
The characterization of the remnant black hole (BH) permassive black holes in the cores of galaxies in the
as the by-product of a generic BH binary (BHB) merger electromagnetic spectrum [25–31] from infrared [32] to
isofgreatastrophysicalinterestasitallowsonetomodel X-rays [33–35] and morphological aspects of the galaxy
the growth of BHs during the evolution of the universe cores [36–38]. Notably, there began to appear observa-
and their effect on the dynamical evolution of galactic tions indicating the possibility of detection of such ef-
cores and globular clusters, as well as the collisions of fects [39–41], and although alternative explanations are
galaxies and stellar size binary systems. Thanks to the possible[42–45],thereisstilltheexcitingpossibilitythat
recent breakthroughs in Numerical Relativity [1–3] one these observations can lead to the first confirmation of a
can now precisely compute the masses, spins and recoil predictionofGeneralRelativityinthehighly-dynamical,
velocitiesofthesemergedBHBs from fullynonlinearnu- strong-field regime.
merical simulations. In our approach to the recoil problem [7, 24] we chose
The modeling of the remnant black hole using fully- to use post-Newtonian theory as a guide to model the
numerical techniques was pioneered by the ‘Lazarus recoil dependence on the physical parameters of the pro-
method’ [15] for spinning black holes followed by the genitorBHB(SeeEqs. (3.31)in[46]),whilearguingthat
breakthrough ‘moving puncture’ approach. In Refs. [4– only full numerical simulations can produce the correct
6] the authors studied BHBs characterized by equal- amplitude of the effect. Bearing this in mind, we pro-
mass, equal-spin individual BHs, with the spins aligned posed an empirical formula for the total recoil veloci-
or counter-aligned with the orbital angular momentum, ties (see Eq. (22) below.) Our heuristic formula describ-
using fully nonlinear numerical calculations and found ing the recoil velocity of a black-hole binary remnant
a simple ad hoc expression relating the final mass and as a function of the parameters of the individual holes
spinoftheremnantwiththespinsoftheindividualBHs. has been theoretically verified in several ways. In [24]
This scenario was later revisited in [16, 17] and the for- thecosΘdependencewasestablishedandwasconfirmed
mula for the remnant spin was generalized (by assuming in [47] for binaries with different initial separations. In
that the angular momentum is only radiated along the [48] the decomposition into spin components perpendic-
orbital axis, and neglecting the energy loss) in [18] for ular and parallel to the orbital plane was verified, and
arbitrary BH configurations (although in the latest pa- in[49]itwasfoundthatthequadratic-in-spincorrections
per of this sequel this condition was removed [19].) In to the in-plane recoil velocity are less than 20 km s−1.
[20] a more general ad hoc fitting function was proposed. Recently in [50] we confirmed the leading η2 (where
A more comprehensive approach was proposed in [21]; η = (mm1+1mm22)2 is the symmetric mass ratio) dependence
where a generic Taylor expansion, reduced by the phys- of the large recoils out of the orbital plane.
ical symmetries of the problem, was used to fit the ex- Since the magnitude (and direction) of the recoil ve-
isting full numerical simulations. A different approach locity of the remnant black holes depend so sensitively
waspresentedin[22]wheretheparticlelimitapproxima- on the spin orientation just around the time of the for-
tion was extended to the equal-mass case and the effects mation of a common event horizon, it is important to
of post-ISCO (Innermost Stable Circular Orbit) gravi- establish that random oriented spins of individual black
tational radiation were neglected. This approach was holes at large separations (as a plausible initial astro-
further improved in [23] by taking binding energies into physical scenario) lead to randomly oriented black holes
account. Alloftheseapproachesshowacertaindegreeof near merger, or if there is some bias in their orientations
agreement with the remnant masses and spins obtained by the time they get very close together (i.e. at typical
in the few dozen fully nonlinear numerical simulations numerical simulations separations of a few M, where M
3
is the total mass). It has been argued recently [51–53] butions of mass ratios and individual spins (magnitudes
that the presence of gas and accretion of the individ- andspins). Thoseevaluationsleadtoalargerecoilveloc-
ual black holes during the inspiral phase for long time ity tail in the distribution with non negligible probabil-
scales can lead to a preferential alignment of spins with ities for v > 1000 km s−1, and highly peaked about the
theorbitalangularmomentum,andhencetoaconfigura- direction of orbital angular momentum at merger. Like-
tionthatleadstomodestrecoilvelocities(afewhundred wise,evaluationsofthefinalspinformulaeleadtoawide
km s−1). The latest results by Dotti et al [53] indicate distribution peaked at magnitudes of S /M2 ≈0.73 and
f f
that spin alignment occur in the scale of a few million orientationspeakedatanangle∼25o withrespecttothe
years within 10 degrees of the orbital angular momen- orbitalangularmomentum. Wecompletethepaperwith
tum for cold disks and 30 degrees for warmer disks. The a discussion of the astrophysical consequences of these
proportion of wet to dry mergers in the universe still results and include an appendix with the computation
needs to be established. Current rough estimates give of the innermost stable circular orbit radius, energy and
comparable percentages for both kinds of mergers. With angular momentum around Kerr black holes (needed for
3to5mergerspergalaxyduringtheirlifetimeforcurrent the remnant formulae), with an analytic solution for the
spiral and elliptic galaxies respectively [102]. While the equatorial and polar orbits.
verificationoftheseclaimsfor‘wet’mergersisunderway,
inthispaperwewouldliketoexplorethepossibilitythat
such alignment (or counter-alignment) mechanism exists
for purely gravitational interactions (‘dry mergers’). In II. INSPIRAL PHASE OF BHBS
general we will seek to find if there is any bias in the
individual spin distributions of black holes at close sep-
A. PN techniques
arations (5−8M for starting typical full numerical evo-
lutions) when starting evolutions with post-Newtonian
methods at large radii with random spin orientations. We construct the PN equations of motion using the
These distributions, in turn, will then help in choosing formulae provided in Refs. [54–57]. To obtain the con-
configurationforfullnumericalsimulationsofclosebina- servative part of the PN equations of motion, we use the
ries. following Hamiltonian,
Thepaperisorganizedasfollows. InSectionIIwede-
scribe the post-Newtonian formalism to analyze the in- H = H +H +H +H
O,Newt O,1PN SO,1.5PN O,2PN
spiralstageofthebinaryevolutions. WeusetheHamilto-
+H +H +H
nianformulation(upto3.5PNorder)toderivetheequa- SS,2PN SO,2.5PN O,3PN
+H +H , (1)
tions of motion in the ADM-TT gauge. Conservative S1S2,3PN S1S1(S2S2),3PN
and radiative effects of the spins are included up to the
next leading PN order. We also include a purely ana- where H contains the terms associated with the orbital
O
lytic analysis of the projection of the quantity ∆(cid:126) along motion up to 3PN order, HSO contains the spin-orbit
the orbital angular momentum L(cid:126), which has a strong coupling terms up to 2.5PN order, and HSS contains the
effect on the recoil velocity, to qualitatively predict a spin-spin coupling term up to 3PN order. Note that
slight bias towards counter-alignment of these two vec- Porto and Rothstein has discussed the spin-spin inter-
tors. Theresultsofthestatisticsofnumericalintegration action by using effective field theory techniques [58–61].
of the post-Newtonian equations of motion (EOM) fol- Theseareaverypowerfulapproachtosystematicallydis-
lows. We performed integrations from initial separations cuss the dynamics of finite size objects.
of r = 50M with 164 spin orientation chosen at ran- The equations of motion are then obtained via,
dom and magnitudes fixed at large astrophysical values,
i.e. Si/m2i = 0.97 for different mass ratios in the range dXi ∂H
1/16 ≤ q = m /m ≤ 1. The results quantitatively con- = {Xi,H}= , (2)
1 2 dt ∂P
firm the bias towards counter-alignment of ∆(cid:126) and total i
dP ∂H
spin S(cid:126) with respect to the orbital angular momentum L(cid:126). i = {P ,H}+F =− +F , (3)
dt i i ∂Xi i
Section III deals with the merger phase, when the black
holes are much closer to each other and in a few orbits dS(cid:126)1 = {S(cid:126) ,H}= ∂H ×S(cid:126) , (4)
willmergeintoasinglelargerone. Thisisthetypicalsce- dt 1 ∂S(cid:126) 1
1
nario that full numerical simulations assume. The bulk dS(cid:126) ∂H
properties of the remnant black hole can be summarized 2 = {S(cid:126) ,H}= ×S(cid:126) , (5)
dt 2 ∂S(cid:126) 2
in terms of empirical remnant formulae that describe its 2
total mass, spin and recoil velocity. We proposed formu-
lae for these quantities based on post-Newtonian scaling where{···,···}denotesthePoissonbrackets,Xi =xi−xi
1 2
with amplitudes fixed by the full numerical simulations. andPiarerelativecoordinatesandlinearmomentaofthe
Withtheseformulaeathand,weperformstatisticalstud- binary, S(cid:126) and S(cid:126) are the spins of each body, and F is
1 2 i
ies by evaluation of these expression for random distri- theradiationreactionforce. Theradiationreactionforce
4
F(cid:126) is given by [54], Here since we focus only on the dissipative effect, we ig-
nore ∆(cid:126)˙ and |∆(cid:126)|· This is because there is no radiation
1 dE
F(cid:126) = P(cid:126) reaction term in Eqs. (4) and (5). The radiation reac-
ω|L(cid:126)| dt tion effect are introduced by the evolution equation of
8 v8 (cid:26)(cid:18) m (cid:19) the linear momentum given in Eq. (3). Furthermore, we
+ η2 ω 61+48 2 P(cid:126) ·S(cid:126)
15 |L(cid:126)|2R m1 1 expect that the time evolution of the spin directions due
to the conservative force will cancel out in a statistical
(cid:18) (cid:19) (cid:27)
m
+ 61+48 1 P(cid:126) ·S(cid:126) L(cid:126) , (6) treatment. Hence, we have
m 2
2
where L(cid:126) =X(cid:126) ×P(cid:126), R=|X(cid:126)|, v =(Mω)1/3 and ω is the
ω
orbital frequency. We use the following notation:
M = m1+m2, (7) (L(cid:126)ˆ·∆(cid:126)ˆ)·dis = L(cid:126)˙|L(cid:126)di|s|∆(cid:126)·∆(cid:126)| − L(cid:126)|·L(cid:126)∆(cid:126)|2||L(cid:126)∆(cid:126)||·dis . (14)
δM = m −m , (8)
1 2
m m
η = 1 2 , (9)
M2
S(cid:126) = S(cid:126) +S(cid:126) , (10)
1 2
(cid:32) (cid:33)
S(cid:126) S(cid:126)
∆(cid:126) = M 2 − 1 , (11) The dissipative effect on the angular momentum is
m m
2 1 given by
δm
S(cid:126) = 2S(cid:126) + ∆(cid:126)
0 M
(cid:18) (cid:19) (cid:18) (cid:19)
m m
= 1+ 2 S(cid:126) + 1+ 1 S(cid:126) . (12)
m1 1 m2 2 L(cid:126)˙ = X(cid:126) ×F(cid:126)
dis
To calculate dE/dt, the instantaneous loss in energy, we = 1 dEL(cid:126)ˆ
use the formulae given in Refs. [62] [103]. ω dt
8 v8 R(cid:26)(cid:18) m (cid:19)
− η2 ω 61+48 2 P(cid:126) ·S(cid:126)
15 |L(cid:126)|2 m1 1
1. PN prediction of distribution of L(cid:126)ˆ·∆(cid:126)ˆ (cid:18) m (cid:19) (cid:27)
+ 61+48 1 P(cid:126) ·S(cid:126) P(cid:126) (15)
m 2
2
The time derivative of the inner product L(cid:126)ˆ ·∆(cid:126)ˆ where
L(cid:126)ˆ and ∆(cid:126)ˆ are the unit vector corresponding to L(cid:126) and ∆(cid:126),
respectively, is given by
where we have used the quasi-circular assumption and
(L(cid:126)ˆ ·˙∆(cid:126)ˆ)= L(cid:126)˙ ·∆(cid:126) + L(cid:126) ·∆(cid:126)˙ Eq. (6).
|L(cid:126)||∆(cid:126)| |L(cid:126)||∆(cid:126)|
L(cid:126) ·∆(cid:126) |L(cid:126)|· L(cid:126) ·∆(cid:126) |∆(cid:126)|· Using this dissipation of the angular momentum, we
− − . (13) obtain
|L(cid:126)|2|∆(cid:126)| |L(cid:126)||∆(cid:126)|2
(L(cid:126)ˆ·∆(cid:126)ˆ)· = − 8 vω11 q 1 (cid:26)−q2 (61q+48)(P(cid:126)ˆ·α(cid:126) )2+(61+48q)(P(cid:126)ˆ·α(cid:126) )2
dis 15 M (1+q)4 |α(cid:126) −qα(cid:126) | 1 2
2 1
(cid:27)
+q [(61q+48)−(61+48q)](P(cid:126)ˆ·α(cid:126) )(P(cid:126)ˆ·α(cid:126) ) , (16)
1 2
in the leading PN order calculation. Here α(cid:126)1 = S(cid:126)1/m21, in L(cid:126)˙dis cancels out, and we have
α(cid:126) =S(cid:126) /m2 and q =m /m . Note that the dE/dt term
2 2 2 1 2
5
(L(cid:126)ˆ·S(cid:126)ˆ)· = − 8 vω11 q 1 (cid:26)q3 (61q+48)(P(cid:126)ˆ·α(cid:126) )2+(61+48q)(P(cid:126)ˆ·α(cid:126) )2
dis 15 M (1+q)4 |q2α(cid:126) +α(cid:126) | 1 2
1 2
(cid:27)
+q [(61q+48)+q (61+48q)](P(cid:126)ˆ·α(cid:126) )(P(cid:126)ˆ·α(cid:126) ) , (17)
1 2
Next, we consider the time integration from t = t to t=t .
i f
(cid:90) tf(L(cid:126)ˆ·∆(cid:126)ˆ)· dt = − 5 (1+q)2 (cid:90) Rf(L(cid:126)ˆ·∆(cid:126)ˆ)· (cid:18)M(cid:19)−3dR
dis 64 q dis R
ti Ri
(cid:26)
= − 1 1 1 −q2 (61q+48)(P(cid:126)ˆ·α(cid:126) )2+(61+48q)(P(cid:126)ˆ·α(cid:126) )2
36(1+q)2 |α(cid:126) −qα(cid:126) | 1 2
2 1
+q [(61q+48)−(61+48q)](P(cid:126)ˆ·α(cid:126) )(P(cid:126)ˆ·α(cid:126) )(cid:27)(cid:34)(cid:18)M (cid:19)3/2−(cid:18)M(cid:19)3/2(cid:35) ,
1 2 R R
f i
(cid:90) tf(L(cid:126)ˆ·S(cid:126)ˆ)· dt = − 5 (1+q)2 (cid:90) Rf(L(cid:126)ˆ·S(cid:126)ˆ)· (cid:18)M(cid:19)−3dR
dis 64 q dis R
ti Ri
(cid:26)
= − 1 1 1 q3 (61q+48)(P(cid:126)ˆ·α(cid:126) )2+(61+48q)(P(cid:126)ˆ·α(cid:126) )2
36(1+q)2 |q2α(cid:126) +α(cid:126) | 1 2
1 2
+q [(61q+48)+q (61+48q)](P(cid:126)ˆ·α(cid:126) )(P(cid:126)ˆ·α(cid:126) )(cid:27)(cid:34)(cid:18)M (cid:19)3/2−(cid:18)M(cid:19)3/2(cid:35) , (18)
1 2 R R
f i
where we considered only the evolution of vω, i.e., the consider the evolution of P(cid:126)ˆ·α(cid:126)2 in Eq. (16),
inspiral,andhaveusedtheleadingradiationreactionand
the Newtonian velocity, (cid:16)P(cid:126)ˆ ·˙(cid:126)α (cid:17)= −vω (cid:16)X(cid:126)ˆ ·α(cid:126) (cid:17) . (21)
2 R 2
dR 64 q (cid:18)M(cid:19)3
= − ,
dt 5 (1+q)2 R This equation means that the direction of α(cid:126)2 does not
(cid:114) change,i.e.,thereisnoprecessionofthespin. Therefore,
M
vω = R . (19) we may replace (P(cid:126)ˆ·α(cid:126)2)2 in Eq. (18) by the one-orbit av-
erage<(P(cid:126)ˆ·α(cid:126) )2 > of(P(cid:126)ˆ·α(cid:126) )2. Althoughtheadiabatic
In the above integration, we derived the formula as- 2 t 2
suming P(cid:126)ˆ·α(cid:126) = constant (i = 1, 2). However, since the evolution of <(P(cid:126)ˆ·α(cid:126)2)2 >t is present, its effect comes in
i at higher PN order in Eq. (18). In this case, it should
spins precess, we need to check the evolution of P(cid:126)ˆ ·α(cid:126)i. be noted that we may consider a test particle orbiting
FromtheevolutionequationsforspinsintheleadingPN around a Kerr black hole with the spin S(cid:126) . According
2
order, the evolution equations for P(cid:126)ˆ·α(cid:126) are given by to [63] in the black hole perturbation approach, the par-
i
ticle’s angular momentum and the black hole’s spin tend
(cid:16)P(cid:126)ˆ·α(cid:126) (cid:17)· = (cid:20)−vω + Mvω 1 (cid:18)2q+ 3(cid:19)(cid:21) to be anti-parallel.
1 R R2 (1+q)2 2 On the other hand, in the case of comparable mass
×(cid:16)X(cid:126)ˆ ·α(cid:126) (cid:17) , binaries, thedirectionofα(cid:126)i changesonatimescalemuch
1 shorter than the integration time. Hence, Eq. (18) is not
(cid:16)P(cid:126)ˆ·α(cid:126)2(cid:17)· = (cid:20)−vRω + MRv2ω (1+qq)2 (cid:18)2+ 32q(cid:19)(cid:21) expInecTteadblteoI,bweeaschcuowratteheinqqde→pe1ndliemnicte.ofEq.(18)when
we ignore the spin precession. Here, we take the average
×(cid:16)X(cid:126)ˆ ·α(cid:126) (cid:17) , (20) with respect to the direction of two spins to represent
2
the randomly oriented spins. We also present the spin
We note that in the limit q → 0, we only need to amplitude dependence in Table II.
6
TABLEI:Theqdependenceintheevolutionof(L(cid:126)ˆ·∆(cid:126)ˆ) and 350
dis
(L(cid:126)ˆ·S(cid:126)ˆ) fromr=50M tor=5M. Weset|α(cid:126) |=|α(cid:126) |=0.97. 300
dis 1 2
q (L(cid:126)ˆ·∆(cid:126)ˆ) (L(cid:126)ˆ·S(cid:126)ˆ) 250
dis dis
1.00 0.0000 −0.0283
200
0.75 −0.0111 −0.0287
0.50 −0.0224 −0.0310 150
0.25 −0.0343 −0.0366
0.125 −0.0406 −0.0412 100
0.0625 −0.0440 −0.0441
0.00 −0.0475 −0.0475 50
Μ
(cid:45)0.5 0 0.5 1
TABLEII:Theamplitudedependenceofthespinintheevo- FIG. 1: The P(µ = L(cid:126)ˆ·∆(cid:126)ˆ) distribution for q = 1 at r = 5M
lution of (L(cid:126)ˆ·∆(cid:126)ˆ)dis and (L(cid:126)ˆ·S(cid:126)ˆ)dis from r = 50M to r = 5M. starting from a uniform distribution at r = 50M. Here we
We set q=0.25 and |α(cid:126)1|=|α(cid:126)2|=α. plot the number of events in the given range of µ out of 164
total events.
α (L(cid:126)ˆ·∆(cid:126)ˆ) (L(cid:126)ˆ·S(cid:126)ˆ)
dis dis
0.97 −0.0343 −0.0366
√
0.97/ 2 −0.0242 −0.0259
0.97/2 −0.0171 −0.0183
0.97/4 −0.0086 −0.0092
0.97/8 −0.0043 −0.0046
0.97/16 −0.0021 −0.0023
B. Statistical Results
For our PN evolutions with used an adaptive fourth-
order Runge-Kutta time-integration scheme with a rela-
tive tolerance of 10−13. The initial data for the simula- mentumforthegivenmassratios. Toanalyzethesedata
tions were generated using the 3PN conservative equa- quantitatively, we bin the data from µ = −1 to µ = 1
tions for quasi-circular orbits with orbital frequency with bin widths of δµ = 0.01. We fit the resulting data
MΩ = 0.00275, which corresponds to an orbital radius P(µ) to a linear function P(µ) = P(0)+ dP| µ for each
dµ 0
of 50±2M. In most cases we stopped the PN simula- mass ratio. The results are summarized in Table III and
tions at a fixed orbital radius of 5M, but also performed plotsofthefitsaregiveninFigs.9-14. Weperformasim-
a set of simulations that terminated at r = 8M in or- ilar analysis for the angle that S(cid:126) makes with the orbital
der to see the effect of the final orbital radius on the angularmomentum(seeFig.17). Wealsoperformasim-
distributions. To obtain the initial PN orbital param- ilaranalysis,butwithqfixedtoq =1/4andα =α =α
√ 1 2
eters, we used uniform distributions of α(cid:126) and α(cid:126) over
1 2 reducedbyfactorsof 2and2,respectively(SeeFigs.7-8
the sphere (by choosing uniform random distributions
and 15-16), and fit the resulting slope dP/dµ as a func-
in µ = cosθ and φ) with fix amplitude α = 0.97. We
tion of α. Here the fit favors a leading-order linear de-
produced 65536 random spin configurations for each fix
pendenceinαoveraleading-orderquadraticdependence
mass ratio q = 1,3/4,1/2,1/4,1/8,1/16. Each run took
(where the constant term is assumed to be zero) (See
approximately10minutes. Inadditionweperformedsets
√ Fig. 18). If we set the constant in the fit to zero, then a
of 65536 run for q = 1/4 and α = 0.97/ 2, α = 0.97/2,
√ lineardependenceis dP/dµ=−(0.02491±0.00098)α for
aswellasα=0.97/ 2butterminatingatr =8M rather thedistributionofL(cid:126)ˆ·∆(cid:126)ˆ anddP/dµ=−(0.0301±0.0041)α
than r =5M. We denote these three latter distributions
in Table III by 0.25S1, 0.25S2, and 0.25F, respectively. forthedistributionoftheL(cid:126)ˆ·S(cid:126)ˆ. Notethattheskewingof
In the following section we examine the distribution the distributions takes place at smaller radii, as can be
of the angle µ = L(cid:126)ˆ ·∆(cid:126)ˆ that ∆(cid:126) = M(S(cid:126) /m −S(cid:126) /m ) seenbydifferencesinthe0.25S1and0.25Fdistributions,
2 2 1 1 whichdifferonlyintheorbitalradius(5M for0.25S1and
makes with the orbital angular momentum (at r =5M).
8M for 0.25F) where the distributions are measured. In
Atr =50M thisdistributionisuniform(sinceS(cid:126) andS(cid:126)
1 2 Table IV we show fits for the distributions of the angles
arechosenfromauniformdistributiononthesphere). In Sˆ ·Lˆ and Sˆ ·Lˆ for the same set of runs. Note that
1 2
Figs. 1-6, we show histograms of the distribution of the the distribution of Sˆ ·Lˆ (the smaller BH’s spin) become
angle L(cid:126)ˆ ·∆(cid:126)ˆ that ∆(cid:126) makes with the orbital angular mo- essentially uniform fo1r q <1/4.
7
350 350
300 300
250 250
200 200
150 150
100 100
50 50
Μ Μ
(cid:45)0.5 0 0.5 1 (cid:45)0.5 0 0.5 1
FIG.2: TheP(µ=L(cid:126)ˆ·∆(cid:126)ˆ)distributionforq=3/4atr=5M FIG.5: TheP(µ=L(cid:126)ˆ·∆(cid:126)ˆ)distributionforq=1/8atr=5M
starting from a uniform distribution at r = 50M. Here we starting from a uniform distribution at r = 50M. Here we
plot the number of events in the given range of µ out of 164 plot the number of events in the given range of µ out of 164
total events. total events.
350 350
300 300
250 250
200 200
150 150
100 100
50 50
Μ Μ
(cid:45)0.5 0 0.5 1 (cid:45)0.5 0 0.5 1
FIG.3: TheP(µ=L(cid:126)ˆ·∆(cid:126)ˆ)distributionforq=1/2atr=5M FIG.6: TheP(µ=L(cid:126)ˆ·∆(cid:126)ˆ)distributionforq=1/16atr=5M
starting from a uniform distribution at r = 50M. Here we starting from a uniform distribution at r = 50M. Here we
plot the number of events in the given range of µ out of 164 plot the number of events in the given range of µ out of 164
total events. total events.
400
350
300
300
250
200
200
150
100
100
50
Μ Μ
(cid:45)0.5 0 0.5 1 (cid:45)0.5 0 0.5 1
FIG.4: TheP(µ=L(cid:126)ˆ·∆(cid:126)ˆ)distributionforq=1/4atr=5M FIG.7: TheP(µ=L(cid:126)ˆ·∆(cid:126)ˆ)distributionforq=1/4atr=5M
starting from a uniform distribution at r = 50M. Here we starting fro√m a uniform distribution at r = 50M and α1 =
plot the number of events in the given range of µ out of 164 α2 =0.97/ 2. Hereweplotthenumberofeventsinthegiven
total events. range of µ out of 164 total events.
8
P(cid:72)Μ(cid:76)
350 0.60
300
250 0.55
200
150 0.50
100
0.45
50
Μ
(cid:45)0.5 0 0.5 1
Μ
FIG.8: TheP(µ=L(cid:126)ˆ·∆(cid:126)ˆ)distributionforq=1/4atr=5M (cid:45)1.0 (cid:45)0.5 0.5 1.0
starting from a uniform distribution at r = 50M and α1 = FIG.11: ThefittothenormalizedP(µ=L(cid:126)ˆ·∆(cid:126)ˆ)distribution
α =0.97/2. Here we plot the number of events in the given
2 at r = 5M for q = 1/2. The data have been binned with a
range of µ out of 164 total events.
bin width of δµ=0.01
P(cid:72)Μ(cid:76) P(cid:72)Μ(cid:76)
0.60 0.60
0.55 0.55
0.50 0.50
Μ
(cid:45)1.0 (cid:45)0.5 0.45 0.5 1.0 Μ
(cid:45)1.0 (cid:45)0.5 0.5 1.0
0.40 0.40
FIG. 9: The fit to the normalized P(µ = L(cid:126)ˆ·∆(cid:126)ˆ) distribution FIG.12: ThefittothenormalizedP(µ=L(cid:126)ˆ·∆(cid:126)ˆ)distribution
at r =5M for q =1. The data have been binned with a bin at r = 5M for q = 1/4. The data have been binned with a
widthofδµ=0.01andnormalizedtoatotalprobabilityof1. bin width of δµ=0.01 and normalized to a total probability
of 1.
P(cid:72)Μ(cid:76) P(cid:72)Μ(cid:76)
0.60 0.60
0.55 0.55
0.50 0.50
Μ
(cid:45)1.0 (cid:45)0.5 0.5 1.0 0.45 Μ
(cid:45)1.0 (cid:45)0.5 0.5 1.0
0.40 0.40
FIG.10: ThefittothenormalizedP(µ=L(cid:126)ˆ·∆(cid:126)ˆ)distribution FIG.13: ThefittothenormalizedP(µ=L(cid:126)ˆ·∆(cid:126)ˆ)distribution
at r = 5M for q = 3/4. The data have been binned with a at r = 5M for q = 1/8. The data have been binned with a
bin width of δµ=0.01 and normalized to a total probability bin width of δµ=0.01 and normalized to a total probability
of 1. of 1.
9
P(cid:72)Μ(cid:76)
0.60 TABLE III: The distribution P(µ) of the angle µ = cosθ
between∆(cid:126) andtheL(cid:126) atr=5M startingfromauniformdis-
tribution at r = 50M (top), and the similar distribution for
the angle between S(cid:126) and L(cid:126) (bottom). The 0.25S1 configura-
0.55 √
tionshadα =α =0.97/ 2andthe0.25S2hasα=0.97/2,
1 2 √
whilethe0.25Fconfigurationshaveα=0.97/ 2andprovide
thedistributionsatr=8M (ratherthanr=5M),allothers
0.50 had α =α =0.97.
1 2
q P(µ)
1.00 0.5000±0.0018+(0.0009±0.0031)µ
0.45
0.75 0.5000±0.0019−(0.0138±0.0034)µ
0.50 0.5000±0.0019−(0.0180±0.0033)µ
Μ
(cid:45)1.0 (cid:45)0.5 0.5 1.0 0.25 0.5000±0.0018−(0.0251±0.0031)µ
0.40 0.125 0.5000±0.0020−(0.0248±0.0035)µ
0.0625 0.5000±0.0019−(0.0226±0.0033)µ
FIG.14: ThefittothenormalizedP(µ=L(cid:126)ˆ·∆(cid:126)ˆ)distribution 0.25S1 0.5000±0.0020−(0.0156±0.0035)µ
at r =5M for q =1/16. The data have been binned with a 0.25S2 0.5000±0.0012−(0.0123±0.0031)µ
bin width of δµ=0.01 and normalized to a total probability 0.25F 0.5000±0.0021−(0.0108±0.0037)µ
of 1. 1.00 0.5000±0.0021−(0.0345±0.0037)µ
0.75 0.5000±0.0020−(0.0284±0.0035)µ
P(cid:72)Μ(cid:76)
0.50 0.5000±0.0019−(0.0286±0.0031)µ
0.60 0.25 0.5000±0.0019−(0.0261±0.0034)µ
0.125 0.5000±0.0018−(0.0249±0.0034)µ
0.0625 0.5000±0.0019−(0.0225±0.0033)µ
0.55 0.25S1 0.5000±0.0020−(0.0162±0.0034)µ
0.25S2 0.5000±0.0019−(0.0125±0.0034)µ
0.25F 0.5000±0.0020−(0.0103±0.0034)µ
0.50
0.01
D
Μ S
(cid:45)1.0 (cid:45)0.5 0.5 1.0
0.00
0.40
−0.01
FIG.15: ThefittothenormalizedP(µ=L(cid:126)ˆ√·∆(cid:126)ˆ)distributionat m
d
r=5M for q =1/4 and α1 =α2 =0.97/ 2. The data have P/
d
beenbinnedwithabinwidthofδµ=0.01andnormalizedto
−0.02
a total probability of 1.
P(cid:72)Μ(cid:76)
−0.03
0.60
−0.04
0.55 0 0.2 0.4 0.6 0.8 1
q
0.50 FIG. 17: The dependence of the slope in the distribution of
the angle between ∆(cid:126) and the orbital angular momentum, as
Μ well as the angle between S(cid:126) and the orbital angular momen-
(cid:45)1.0 (cid:45)0.5 0.45 0.5 1.0 tum as a function of mass ratio.
An important consequence of choosing uniform distri-
0.40
butions for the directions of S(cid:126) and S(cid:126) (with magnitude
1 2
FIG.16: ThefittothenormalizedP(µ=L(cid:126)ˆ·∆(cid:126)ˆ)distribution |S(cid:126)i|=0.97)isthattheinitialdistributionsforthesquares
atr=5M forq=1/4andα1 =α2 =0.97/2. Thedatahave of the magnitudes of S(cid:126) and ∆(cid:126), P(S2) and P(∆2), are
beenbinnedwithabinwidthofδµ=0.01andnormalizedto uniform in the range [α(m2 −m2)]2 to [α(m2 +m2)]2,
2 1 2 1
a total probability of 1. and zero outside this range (i.e. there is an equal prob-
10
0.000
TABLE IV: The distribution P(µ) of the angle µ = cosθ
(D) Fit of dP(m)/dm = Ba
betweenS(cid:126)1 andtheL(cid:126) atr=5M startingfromauniformdis- (D) Fit of dP(m)/dm = Ba2
tribution at r = 50M (top), and the similar distribution for −0.005 (S) Fit of dP(m)/dm = Ba
theanglebetweenS(cid:126) andL(cid:126) (bottom). The0.25S1configura- (S) Fit of dP(m)/dm = Ba2
2 √ D
tionshadα =α =0.97/ 2andthe0.25S2hasα=0.97/2,
1 2 √ −0.010 S
whilethe0.25Fconfigurationshaveα=0.97/ 2andprovide 4)
thedistributionsatr=8M (ratherthanr=5M),allothers =1/
q
had α1 =α2 =0.97. Note that the distribution of angles for at −0.015
the smaller component S(cid:126)1 becomes uniform as q→0. d (m
P/
q P(µ) d −0.020
1.00 0.5000±0.0019−(0.0278±0.0033)µ
0.75 0.5000±0.0020−(0.0129±0.0034)µ
0.50 0.5000±0.0019−(0.0189±0.0033)µ −0.025
0.25 0.5000±0.0019−(0.0044±0.0034)µ
0.125 0.5000±0.0019−(0.0000±0.0033)µ
−0.030
0.0625 0.5000±0.0019−(0.0026±0.0033)µ 0.4 0.5 0.6 0.7 0.8 0.9 1
0.25S1 0.5000±0.0020−(0.0019±0.0034)µ a
0.25S2 0.5000±0.0018−(0.0008±0.0031)µ
0.25F 0.5000±0.0021−(0.0007±0.0036)µ
FIG. 18: The dependence of the slope in the distribution of
1.00 0.5000±0.0020−(0.0237±0.0035)µ theanglebetween∆(cid:126) andtheorbitalangularmomentumL(cid:126) for
0.75 0.5000±0.0018−(0.0252±0.0032)µ
q=1/4 as a function of |α(cid:126) |=|α(cid:126) |=α, as well as the angle
1 2
0.50 0.5000±0.0020−(0.0259±0.0034)µ between S(cid:126) and the orbital angular momentum as a function
0.25 0.5000±0.0020−(0.0261±0.0034)µ
of α. In all fits the constant term is taken to be zero. The
0.125 0.5000±0.0019−(0.0249±0.0034)µ
data here favor a linear dependence in α.
0.0625 0.5000±0.0019−(0.0225±0.0033)µ
0.25S1 0.5000±0.0021−(0.0162±0.0037)µ
0.25S2 0.5000±0.0020−(0.0125±0.0034)µ
ing quasi-circular orbits. This equation is coupled with
0.25F 0.5000±0.0020−(0.0105±0.0034)µ
the spin and angular momentum precession equations,
which include the leading order spin-orbit and spin-spin
couplings. On the other hand, in our calculation, the
ability of finding any given value of S2 or ∆2 in this
PN equations of motion are derived from the Hamilto-
range). HoweverthedistributionsP(∆)andP(S)there-
nian and include radiation reaction effects. These have
fore contain an additional linear factor in ∆ and S (i.e.
higherPNorderspin-orbitandspin-spincouplingterms.
P(x) = 2xP(x2) for any variable x), respectively. One
Furthermore, the second term of the right hand side of
immediate consequence is that the distributions P(∆)
Eq. (6) has a significant effect in the PN evolutions. Al-
and P(S) are both maximized for the largest allowed though the evolution of L(cid:126)ˆ in [65] is determined only by
values of S and ∆. Given the observation that large ∆
the conservative dynamics, we have also considered the
in the orbital plane [7] leads to very large recoils, this
dissipativeeffectduetotheradiationreaction. Wefindin
bias, if present in nature, would favor observations of
the PN prediction that this dissipative effect creates the
largerecoils. SeeSec.IVforfurtheranalysisoftherecoil
statistically significant counter-alignment of the spins.
distribution.
Schnittman in Ref. [64] has studied the evolution of
spins in binary systems using orbit-averaged PN equa-
tions of motion what allowed longer term evolutions III. MERGER PHASE OF BHBS
(from separations up to 1000M). The results indicate
strong correlations of the late angle among spins when Unlike in the earlier inspiral phase, during the plunge
one starts fixing the initial direction of the spin of the and merger the PN equations of motion do not provide
primaryobjectandchoosethesecondary’sspindirection a quantitatively accurate description of the merger dy-
at random (See Figs. 6 and 7 in [64].) Bogdanovic et namics, and therefore do not provide robust estimates of
al revisit this scenario in Ref. [51] and find that if one thefinalremnantmass,spin,andrecoil. However,analy-
isallowedtochooseinitialrandomdistributionsforboth sis of the recoil in particular shows that PN analysis can
spins the resulting evolution leads to close to isotropic be used to derive heuristic formulae (based on how PN
distributionsofthelatedirectionsofthespins(Seetheirs predictions scale with spins and masses) that give quan-
Fig. 1). In our paper we find an small but statistically titatively correct predictions [7, 24, 66] and incorporate
significant bias towards counteralignment of the spins the symmetries of the problem. We will use this model-
with the orbital angular momentum (See Figs. 9-16.) ing in the case of the total radiated energy and angular
More recently, Herrmann et al. presented numerical momentum. Inparticularwewillsupplementtheinspiral
studies of the PN equations on GPUs [65]. They used losses, modeledbytheenergyandangularmomentumof
the evolution equation for the orbital frequency assum- the ISCO in the particle limit (extended to the compa-