Abstract
Title of Dissertation: Growing Intermediate-Mass Black Holes
with Gravitational Waves
Kayhan G?ultekin, Doctor of Philosophy, 2006
Dissertation directed by: Professor M. Coleman Miller
Department of Astronomy
We present results of numerical simulations of sequences of binary-single scatter-
ing events of black holes in dense stellar environments. The simulations cover a wide
range of mass ratios from equal mass objects to 1000:10:10 Mcircledot and compare purely
Newtonian simulations with a relativistic endpoint, simulations in which Newto-
nian encounters are interspersed with gravitational wave emission from the binary,
and simulations that include the effects of gravitational radiation reaction by using
equations of motion that include the 2.5-order post-Newtonian force terms, which
are the leading-order terms of energy loss from gravitational waves. In all cases,
the sequence is terminated when the binary?s merger time due to gravitational ra-
diation is less than the arrival time of the next interloper. We also examine the
role of gravitational waves during an encounter and show that close approach cross-
sections for three 1 Mcircledot objects are unchanged from the purely Newtonian dynamics
except for close approaches smaller than 10?5 times the initial semimajor axis of
the binary. We also present cross-sections for mergers resulting from gravitational
radiation during three-body encounters for a range of binary semimajor axes and
mass ratios including those of interest for intermediate-mass black holes (IMBHs).
We find that black hole binaries typically merge with a very high eccentricity ?
extremely high when gravitational waves are included during the encounter such
that when the gravitational waves are detectable by LISA, most of the binaries will
have eccentricities e> 0.9 though all will have circularized by the time they are de-
tectable by LIGO. We also investigate the implications for the formation and growth
of IMBHs and find that the inclusion of gravitational waves during the encounter
results in roughly half as many black holes ejected from the host cluster for each
black hole accreted onto the growing IMBH. The simulations show that the Miller
& Hamilton (2002b) model of IMBH formation is a viable method if it is modified
to start with a larger seed mass.
Growing Intermediate-Mass Black Holes
with Gravitational Waves
by
Kayhan G?ultekin
Dissertation submitted to the Faculty of the Graduate School of the
University of Maryland at College Park in partial fulfillment
of the requirements for the degree of
Doctor of Philosophy
2006
Advisory Committee:
Professor M. Coleman Miller, chair
Professor Alessandra Buonanno
Professor Douglas P. Hamilton
Professor Richard F. Mushotzky
Professor Chris S. Reynolds
Professor Derek C. Richardson
c? Kayhan G?ultekin 2006
Preface
Much of the material contained in this dissertation has been published. The work for
Chapter 2 was published in the Astrophysical Journal as ?Growth of Intermediate-
Mass Black Holes in Globular Clusters? (G?ultekin et al. 2004), and the work for
Chapter 3 was published in the Astrophysical Journal as ?Three-Body Dynamics
with Gravitational Wave Emission? (G?ultekin et al. 2006). Both works were parts
of a single study into the dynamics of black holes in dense stellar clusters. This
dissertation brings together those two works with other material for a coherent yet
ongoing study of three-body dynamics with gravitational waves. It examines the
influence of dynamics on gravitational waves, the influence of gravitational waves on
dynamics and their influence on the formation of intermediate-mass black holes.
ii
For Laura and our children.
iii
Acknowledgements
I am grateful to Cole for being my advisor. Cole has been a model
mentor, advisor, and research supervisor throughout my work with him.
He has pushed me intellectually and professionally through his incisive
Socratic questioning, his willingness to pursue a myriad of cerebral cu-
riosities, and his consistent promotion of my work in the astronomical
community. I cannot imagine a better advisor for me.
I am also grateful to Doug for being my collaborator. His advice and
guidance has been helpful and much appreciated.
My time at the University at Maryland has been enriched by all of the
astronomy department faculty, students, and staff.
Finally, IthankLauraforherkindpatienceandwonderfulhelpthroughout
this long ordeal, and I thank my parents for their support and en-
couragement throughout my entire education. I have been lucky to
have had the freedom to follow my heart?s content as a career.
iv
Contents
List of Tables vii
List of Figures viii
1 Introduction 1
1.1 Scientific Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.1.1 Ultraluminous X-ray Sources . . . . . . . . . . . . . . . . . . . 1
1.1.2 Kinematical and Dynamical Evidence for IMBHs . . . . . . . 11
1.2 IMBH Formation Models . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.2.1 Population III Stars . . . . . . . . . . . . . . . . . . . . . . . 14
1.2.2 Dynamics of Stellar Clusters . . . . . . . . . . . . . . . . . . . 14
1.3 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
1.3.1 N-body Simulations . . . . . . . . . . . . . . . . . . . . . . . 18
1.3.2 Gravitational Waves . . . . . . . . . . . . . . . . . . . . . . . 19
1.4 Overview of This Dissertation . . . . . . . . . . . . . . . . . . . . . . 23
2 Growth of Intermediate-Mass Black Holes in Globular Clusters 26
2.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
2.2 Numerical Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
2.3 Simulations and Results . . . . . . . . . . . . . . . . . . . . . . . . . 32
2.3.1 Pure Newtonian Sequences . . . . . . . . . . . . . . . . . . . . 34
2.3.2 General Relativistic Binary Evolution . . . . . . . . . . . . . . 38
2.4 Implications for IMBH Formation and Growth . . . . . . . . . . . . . 42
2.5 Implications for Gravitational Wave Detection . . . . . . . . . . . . . 49
2.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
3 Three-Body Dynamics with Gravitational Wave Emission 54
3.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
3.2 Numerical Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
3.3 Simulations and Results . . . . . . . . . . . . . . . . . . . . . . . . . 60
3.3.1 Individual Binary-Single Encounters . . . . . . . . . . . . . . 60
3.3.2 Sequences of Encounters . . . . . . . . . . . . . . . . . . . . . 68
v
3.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
3.4.1 Implications for IMBH Formation and Growth . . . . . . . . . 71
3.4.2 Implications for Gravitational Wave Detection . . . . . . . . . 77
3.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
4 Summary and Conclusions 83
4.1 Answers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
4.2 Possible Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
4.3 For Posterity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
A Time for IMBH Progenitor to Acquire a Companion 92
B Demonstrationthat2.5PNDragForceReproducesthePetersEqua-
tions for Circular Orbits 95
C Summary of Code 97
C.1 Initial Conditions Generator . . . . . . . . . . . . . . . . . . . . . . . 98
C.2 HNBody-iabl Interface . . . . . . . . . . . . . . . . . . . . . . . . . . 99
C.3 Examination of Integration . . . . . . . . . . . . . . . . . . . . . . . . 100
C.4 Two-Body Approximator . . . . . . . . . . . . . . . . . . . . . . . . . 101
D Other Applications of the Code 102
D.1 Planet in Close Triple System . . . . . . . . . . . . . . . . . . . . . . 102
D.2 Near-Earth Asteroid Binary Disruption . . . . . . . . . . . . . . . . . 105
D.3 IMBH-IMBH Binaries . . . . . . . . . . . . . . . . . . . . . . . . . . 107
Bibliography 109
vi
List of Tables
2.1 Single Encounter Cross-sections for Exchange. . . . . . . . . . . . . . 33
2.2 Sequence Statistics. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
2.3 IMBH Formation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
3.1 Sequence Statistics. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
vii
List of Figures
2.1 Newtonian 1000:10:10 sequence. . . . . . . . . . . . . . . . . . . . . . 35
2.2 Histograms of final semimajor axes for all mass ratios. . . . . . . . . . 37
2.3 Histogram of final eccentricities for 1000:10:10 mass ratio. . . . . . . . 39
2.4 Comparison of sequence eccentricities with thermal distribution. . . . 40
2.5 1000:10:10 gravitational radiation sequence. . . . . . . . . . . . . . . 41
2.6 Total number of black holes ejected as a function of seed mass. . . . . 44
2.7 Probability of IMBH to remain in cluster. . . . . . . . . . . . . . . . 45
2.8 Total time to build up to given mass. . . . . . . . . . . . . . . . . . . 47
2.9 Distribution of eccentricities when binary is in LISA band. . . . . . . 51
3.1 Comparison of HNDrag with Peters equations. . . . . . . . . . . . . . 58
3.2 Two-body capture by means of emission of gravitational waves. . . . 59
3.3 Cross-section for close approach. . . . . . . . . . . . . . . . . . . . . . 63
3.4 Derivative of close approach cross-section curve. . . . . . . . . . . . . 64
3.5 Cross-section for close approach. . . . . . . . . . . . . . . . . . . . . . 65
3.6 Merger cross-section for 10:10:X mass series. . . . . . . . . . . . . . . 67
3.7 Merger cross-section for 1000:1000:X mass series. . . . . . . . . . . . 68
3.8 Merger cross-section for X:10:10 mass series. . . . . . . . . . . . . . . 69
3.9 Merger cross-section for 1000:X:1000 mass series. . . . . . . . . . . . 70
3.10 Growth rate of IMBH progenitor. . . . . . . . . . . . . . . . . . . . . 73
3.11 Number of stellar-mass black holes ejected from cluster. . . . . . . . . 75
3.12 Probability of IMBH-progenitor retention. . . . . . . . . . . . . . . . 76
3.13 Eccentricity of merging binaries while in LISA band. . . . . . . . . . 80
3.14 Eccentricity of merging binaries while in LIGO band. . . . . . . . . . 81
D.1 Methods of forming stellar triple with planet. . . . . . . . . . . . . . 105
D.2 Exchange encounter that leads to a planet in a triple system. . . . . . 106
D.3 NEA binary disruption. . . . . . . . . . . . . . . . . . . . . . . . . . . 107
viii
Chapter 1
Introduction
1.1 Scientific Motivation
1.1.1 Ultraluminous X-ray Sources
Thescientificmotivationforthisdissertationistheevincedexistenceofintermediate-
mass black holes (IMBHs): black holes too massive to be stellar-mass black holes
and not massive enough to be super-massive black holes (20 Mcircledot <? M <? 105 Mcircledot).
The evidence for IMBHs begins with X-ray observations of so-called ultraluminous
X-ray sources (ULXs). ULXs are typically defined as point sources that are not
known to be quasars, galactic nuclei, or supernovae and that have an X-ray lumi-
nosity in the 2 ? 10 keV band of LX > 1039 erg s?1. Others have restricted the
definition essentially to encompass intermediate-mass black hole candidates by re-
quiring a bolometric luminosity LB > 2?1039 erg s?1, an irregular variability in the
X-ray flux, and a resolvable separation from the host galaxy?s nucleus (Kaaret et al.
2004). Here luminosities are equivalent isotropic luminosities, that is, the luminos-
ity assuming that the observed flux is emitted isotropically, and the values given in
the definition of ULXs indicate a mass larger than a stellar-mass black hole by the
1
Eddington argument described below. X-ray observations of galaxies in the early
1980s with the Einstein X-ray satellite revealed X-ray emission not only from active
galactic nuclei (AGNs) but also, unexpectedly, from the centers of otherwise normal
spiral galaxies. With a spatial resolution of 1prime on its Imaging Proportional Counter,
Einstein was not able to determine whether the ULXs were nuclear sources in their
host galaxies or even if they were actually multiple objects (at a distance of 5 Mpc
1prime is 1.5 kpc). In the early 1990s, the successful launch of the ROSAT satellite with
a spatial resolution of 10primeprime on its High Resolution Imager enabled astronomers to
place some of the ULXs outside of the galaxy?s nucleus. With the high resolution
available on the current generation of X-ray telescopes (4primeprime on XMM-Newton and
1primeprime on Chandra), the positions of many more ULXs have narrowed to exclude the
galactic nuclei, and variability observed in multiple observations for some ULXs pre-
cludes the possibility of multiple objects within a single beam since multiple objects
should not be able to coordinate their variability. In addition, these instruments can
provide high spectral resolution observations while excluding contaminating sources.
There are now more than 200 known ULXs with fluxes measured to correspond
to luminosities up to 1041 erg s?1. That these sources are so bright in X-rays and
are variable indicates that the sources are powered by accretion onto a black hole.
The small size of a black hole allows matter to flow deep into its potential well
where friction and tidal forces heat up the material enough to emit X-rays. The
small size also allows for variability on short timescales. The inferred luminosity at
which the source is emitting, however, cannot come from an ordinarily emitting, or-
dinarily accreting, ordinary-sized black hole. Assuming isotropic emission, isotropic
accretion, and an opacity dominated by scattering off free electrons in a completely
ionized medium, it is possible to calculate a given mass?s Eddington luminosity LE,
the maximum luminosity possible without blowing the accreting medium away with
2
radiation pressure. Equating the force due to radiation pressure from Thomson
scattering
Frad = ?TFc = ?TL4picr2 (1.1)
(where Frad is the force from radiation pressure, ?T is the Thomson scattering cross-
section, F is the radiation flux on the medium at a distance r, and L is the intrinsic
radiation luminosity) with the gravitational force on the medium
Fgrav = GMmHr2 (1.2)
(whereFgrav is the gravitational force of the accreting body of massM on a hydrogen
atom of mass mH to get
LE = 4piGcMmH?
T
= 1.26?1038
parenleftBigg M
Mcircledot
parenrightBigg
erg s?1. (1.3)
Taking into consideration that the X-ray luminosity is only a fraction (0.1 - 0.5)
of the bolometric luminosity and given the assumptions, the ULXs are powered by
black holes with masses greater than 20Mcircledot ? greater than 500Mcircledot for the brightest
sources. This is surprising because the maximum mass of black holes that form
from the core collapse of an evolved star of solar metallicity is thought to be about
20 Mcircledot (Fryer & Kalogera 2001). The spatial separation of the ULXs from the
centers of their host galaxies indicates that these objects are not supermassive black
holes (SMBHs) that are in a low spectral state. SMBHs, with masses greater than
105 Mcircledot and less than 1010 Mcircledot, are massive enough to sink to the center of the galaxy
from the effects of dynamical friction on timescales much shorter than the age of
the host galaxies. For example, the brightest ULX in M82, X1, is at a distance of
200 pc from the kinematical center of the galaxy in which the velocity dispersion is
vdisp = 100 km s?1 (Gaffney et al. 1993; Kaaret et al. 2001; Matsumoto et al. 2001).
Assuming that the density is an isothermal sphere, the time for a mass M to sink
3
from a distance r to the center through dynamical friction is
tdf = 1.65ln?r
2vdisp
GM (1.4)
where ln? ? 10 is the Coulomb logarithm (Binney & Tremaine 1987). So a modest
SMBH withM ? 105 Mcircledot orbiting at 500 pc from the nucleus would have sunk to the
center of the galaxy in less than 1010 yr. A SMBH with M ? 106 Mcircledot would sink in
less than 109 yr. Colbert & Mushotzky (1999) found that in spiral galaxies ULXs are
typically separated from the galactic nucleus with an average separation of 390 pc.
ULXs, then, are less massive than SMBHs, and, by the Eddington argument, they
are more massive than stellar-mass black holes. This represents a new category
of black holes: IMBHs. Although there is little doubt left that stellar mass black
holes and SMBHs exist and although the evidence in favor of some IMBHs is strong,
the existence of this new type of black hole is not universally accepted. Universal
acceptance will probably only follow dynamical measurements of the mass of the
black hole that powers a ULX. The fear of writing a dissertation on a possibly
non-existent object notwithstanding, alternative explanations for ULXs must be
examined.
Beaming
The most commonly mentioned alternative to IMBHs as the mechanism to produce
ULXs is the beaming hypothesis of King et al. (2001, see also Reynolds et al. 1997).
In this model, accretion onto stellar mass black holes powers a narrow cone of
emission that is pointed toward us. Because the emission is highly anisotropic, the
inferred luminosities are overestimates of the ULXs? intrinsic luminosities. Thus the
Eddington argument for the lower limit on the mass of the accreting black hole does
not apply to these beamed sources. King et al. (2001) propose that beaming can
be achieved, for example, with a medium with small scattering opacity near the
4
rotational poles of a compact object with a thick disk surrounding it. The beaming
may also be augmented by super-Eddington emission, described below.
Although it would be reckless to say that none of the ULXs is a beamed source,
there are three strong lines of observational evidence against beaming for some
individual ULXs: He II emission from gas surrounding the ULXs, low temperatures
of accretion disks, and observed quasi-periodic oscillations (QPO) in the X-ray fluxes
in two ULXs. He II emission has been found in the nebula surrounding several ULXs.
The most extensively studied such source is the bright ULX in the Holmberg II dwarf
irregular galaxy with inferred X-ray luminosity LX = 5 to 16?1039 erg s?1 (Kaaret
et al. 2004; Pakull & Mirioni 2002). A nebula at the location of the ULX was
detected with emission at ? = 468.6 nm, which comes from the recombination of
fully ionized helium (Pakull & Mirioni 2002). Since the He II emission should be
isotropic, the measurement of the He II flux is a measure of the X-ray luminosity
required to ionize it. The observed ? = 468.6 nm luminosity of 2.5?1036 erg s?1 is
calculated to come from an emitted X-ray luminosity of LX = 3 to 13?1039 erg s?1
(Pakull & Mirioni 2002). This is completely consistent with a roughly isotropic
emission of the observed X-ray flux. Later observations with the Hubble Space
Telescope (HST) were able to resolve the morphology of the optical emission from
this source including He II, [O I], and H? emission (Kaaret et al. 2004). The
morphology shows that the central region of the nebula contains the He II emission,
and the outer regions to the west, where the He II emission ends, contain the [O I]
emission. Both of these regions also show H? emission. This is as expected for an
X-ray photoionized nebula (Kaaret et al. 2004).
The second line of evidence against beamed emission from accreting stellar-mass
black holes is the inference of cool accretion disks. X-ray spectra are typically mod-
eled with either a power-law, a multi-color disk (MCD), or both. The hard power-law
5
component is thought to come from a Comptonized disk or corona. The MCD model
consists of a disk with a temperature profile (for a thin disk T(r) ? r?3/4) so that
each annulus of disk emits as a blackbody. Since the flux emitted from the disk
scales as F ? T4, the bolometric luminosity from the innermost part of the disk
scales as L ? R2inF ? M2T4in for a characteristic radius Rin, which scales with the
Schwarzschild radius and therefore the mass. Then for a black hole that is accreting
at a fixed rate with respect to Eddington, L ? ?M ? M, and the temperature of
the inner disk scales as Tin ? M?1/4. If the accretion rate is allowed to vary for a
fixed M, the temperature scales as Tin ? ?M1/4. A typical disk temperature for a
? 10 Mcircledot stellar-mass black hole binary is ? 107 K. Early spectra of ULXs from
ASCA, which were modeled only with a MCD, revealed disk temperatures that were
too hot to be explained by IMBHs. Chandra and XMM Newton with their higher
spatial resolution, however, were able to eliminate contaminating sources, and with
their higher spectral resolution, they allowed modeling of spectra with combined
power-law and MCD models. With these improvements the models revealed typical
temperatures of T ? 106 K for inferred luminosities of ? 1040 erg s?1(Miller et al.
2004). The combination of such high luminosities with low disk temperatures is
indicative of higher masses. Although this measure of the mass is model-dependent
and simple scalings of commonly accepted properties of X-ray binaries may not ap-
ply, it is difficult to imagine a scenario in which a beamed cone of emission could
produce a spectrum that would look like a disk with a high luminosity and low tem-
perature. One such imagined scenario would be if the beamed emission is optically
thick so that it will cool along the cone of emission to a point in the cone with a large
effective area. In addition, recent work by Winter et al. (2005) shows that ULXs
powered by IMBHs appear to produce X-ray spectra similar to those of Galactic,
stellar-mass black hole systems with high/soft and low/hard states. They also find
6
that some ULXs appear to be very high state stellar-mass black hole systems. So
it appears that some ULXs are consistent with an accreting IMBH, and some are
bright, ordinary X-ray binary systems.
The third line of evidence against beamed emission comes from the detection
of QPOs in two ULXs. The first is an 8.5% rms amplitude oscillation in the X-
ray flux between 2 and 10.0 keV with a centroid frequency of 54 mHz in M82 X-1
(Strohmayer & Mushotzky 2003). Although the QPOs are not completely under-
stood, they are generally thought to be disk phenomena. This means that the disk,
which will be emitting roughly isotropically, is producing a minimum of 8.5% of
the total inferred luminosity LQPO = 0.085 ? 4 ? 1040 = 3.4 ? 1039 erg s?1. That
is, if the disk is emitting only X-rays in the QPO, it is still above the Eddington
luminosity of a 20 Mcircledot object. In addition, the observations show a broad 6.55 keV
Fe K line that is required for all fits to the spectrum. Since this broad line results
from the X-ray continuum illumination of cold, high-density material, it is difficult
to understand how a narrow beam of emission, which would emerge perpendicular
to the disk, could interact with the cold, dense matter. The second report of a QPO
is a 9.3% rms amplitude QPO in the 0.2 to 10 keV band with a centroid frequency
of 203 mHz in the bright ULX in Holmberg IX (Ho IX X-1 Dewangan et al. 2006).
Although not as strong a case as that of M82 X1, combined with the low inferred
disk temperature of T ? 1?4?106 K (Miller et al. 2004), it is difficult to explain
how the spectrum could be created from a beamed source.
7
Super-Eddington Accretion
Another objection to invoking IMBHs to power ULXs is that the Eddington limit
may not apply because of effects from rapidly rotating, thick disks or because of
inhomogeneities in a thin accretion disk?s density. For thin, radiation-pressure-
dominatedaccretiondisks, theso-called?photon-bubbleinstability?canform(Arons
1992; Gammie 1998). Pockets of low density within a surrounding, higher density
medium effectively trap photons by way of radiative diffusion. These photon bub-
bles will grow due to radiation pressure, and the denser regions will grow denser
as radiation leaves them and as magnetic tension prevents the gas from expand-
ing. If the low-density regions are dynamically coupled to the high-density regions,
the disk may remain in an overall state of dynamical equilibrium (Begelman 2002).
This allows the photons to travel along low-density escape routes in the disk without
blowing the disk apart, and thus the system may radiate above the Eddington limit.
Begelman (2006) finds that the maximum luminosity is
Lmax ? 20
parenleftbigg ?
0.03
parenrightbiggparenleftBigg M
Mcircledot
parenrightBigg1/9
LEdd (1.5)
where ? = pm/pr is the ratio between the magnetic pressure and radiation pressure.
For stellar-mass black holes with M = 20 Mcircledot, luminosities may reach up to Lmax ?
30LEdd ? 7 ? 1040 erg s?1. Thus, super-Eddington accretion by means of photon
bubbles alone is just enough to explain the brightest ULXs. One would, however,
still expect disks accreting at super-Eddington to be hotter than is seen in some
ULXs. In addition, these studies are speculative in nature, and it is possible that
photon bubbles may not form at all because of the incoherence and time-dependence
of the magnetic field (Begelman 2006).
8
Frequency of ULXs
The number of identified ULXs is currently around 200 (Colbert & Ptak 2002; Ptak
& Colbert 2004; Swartz et al. 2004). Colbert & Ptak (2002) found 87 ULXs from
ROSAT HRI data, and Swartz et al. (2004) used Chandra data to find 154, of which
they conservatively estimate ? 25% to be background sources. Based on ROSAT
HRI observations, at least 12% of all galaxies contain a ULX withLX > 1039 erg s?1,
and at least 1% of all galaxies contain a ULX with LX > 1040 erg s?1 (Ptak &
Colbert 2004). These frequencies are lower limits because the 10primeprime angular resolution
of HRI limits the number of ULXs found near but not in the nuclei of the galaxies and
because its comparatively soft bandpass of 0.1?2.2 keV will miss absorbed ULXs
behind a column density greater than 1021 cm?2. Consistent with these results is the
possibility that the frequency of ULXs with LX > 2?1039 erg s?1 has been evolving
from ? 36% at a redshift of z ? 0.1 to ? 8% in the nearby universe (Hornschemeier
et al. 2003). ULXs also tend to be found in star-forming galaxies, but of galaxies
with ULXs, ellipticals have the largest number per galaxy, perhaps because of their
larger masses (Colbert & Ptak 2002; Ptak & Colbert 2004). When a spiral galaxy
is the host, the ULXs tend to be found near the nucleus with an average projected
separation of 390 pc (Colbert & Mushotzky 1999). When in ellipticals, however,
ULXs are found in the halo (Irwin et al. 2003).
Luminosity Function of ULXs
The X-ray luminosity function of ULXs shows that ULXs in ellipticals may be
of a different nature than those in spirals. In ellipticals the cumulative luminosity
function is best explained by a single power-law with slope?1.7 while the luminosity
function in spirals is best explained by a broken power-law with a slope ?0.6 below
LX = 1040 erg s?1 and?1.9 above (Swartz et al. 2004). The fact that the population
9
of ULXs in ellipticals can be fit by a single power law suggests that they may be
high-luminosity, stellar-mass X-ray binaries, but with only 7 ULXs with LX >
3 ? 1039 erg s?1 this claim cannot be made with high statistical certainty. The
absence of a break in the luminosity function does not necessarily mean that there
is only a single population of sources. After all, although both accreting neutron
stars and stellar-mass black holes contribute to the X-ray luminosity function, there
is no break at ? 2?1038 erg s?1, the Eddington luminosity for M = 1.4 Mcircledot. Much
has been made about the presence of a ?knee? in the luminosity function of the
population of ULXs in spirals, but there are only ? 15 sources brighter than the
break point, and a power-law with an exponential cutoff also produces a satisfactory
fit (Swartz et al. 2004).
Cluster Associations and Companions
Observations of ULXs and their environments at other wavelengths show an as-
sociation between ULXs and stellar clusters. For example, M82 X-1 is spatially
coincident with the young stellar cluster MGG 11 as determined by near infrared
observations (McCrady et al. 2003). Fabbiano et al. (2001) found a spatial corre-
lation of ULXs with stellar clusters in the merging Antennae system in excess of
that expected from a uniform distribution of ULXs. A comparison of Chandra and
HST images of the CD galaxy NGC 1399 at the center of the Fornax cluster shows a
spatial correlation between many of its X-ray point sources and its globular clusters
(Angelini et al. 2001). These X-ray point sources include two of the three sources
with LX >? 2?1039 erg s?1, and the globular cluster and X-ray positions agree to
within the combined astrometric uncertainties.
In addition to an association with stellar clusters, individual stellar companions
to ULXs have been found in several cases. These companions are thought to be the
10
donor stars in the accreting system. Many of these observations reveal O- or B-giant
or -supergiant stars as companions such as NGC 5204 X-1 (Goad et al. 2002; Roberts
et al. 2001), M81 X-11 (Liu et al. 2002), and M101 ULX1 (Kuntz et al. 2005; Mukai
et al. 2005). In these cases, the ULXs appear to be very luminous X-ray binaries
with black holes that are not significantly more massive than 20 Mcircledot. Intriguingly,
one of the best IMBH candidates, HoII X-1 also contains a hot stellar companion
of spectral type between O4 and B3, but it cannot be determined whether the star
is a dwarf or supergiant. Such searches for companions are of special importance
because the measurement of the period and velocity of the star?s orbit around the
black hole would yield a lower limit to the mass of the system. Knowledge of the
inclination angle gives the mass of the system. The inclination may be found from
evidence of eclipses, which indicate nearly edge-on systems, and from modeling the
light curve of the distorted stellar companion as the orientation of its non-spherical
shape changes with time.
1.1.2 Kinematical and Dynamical Evidence for IMBHs
The presence of non-luminous matter can be inferred from the dynamical and kine-
matical fingerprint it leaves on the luminous matter in the area. While a bound
stellar companion to a black hole is the best way to measure the black hole?s mass,
stars in a dense stellar environment may interact sufficiently to measure the mass of
the system. In this way, evidence from radial velocities of individual stars in M15 as
well as velocity and velocity dispersion measurements in G1, a large stellar cluster
in the Andromeda Galaxy (M31), indicate that these globular clusters may harbor
large dark masses in their cores (? 3.9 [?2.2]?103 Mcircledot and 1.7 [?0.3]?104 Mcircledot,
respectively, Gebhardt et al. 2000, 2002, 2005; Gerssen et al. 2002). For M15 the
data cannot rule out the presence of many compact remnants that contribute to
11
the high mass-to-light ratio, but the most recent observations of G1 can rule out
the absence of dark mass at the 97% confidence level (Gebhardt et al. 2005). One
caveat that should be mentioned is that at least the core of G1 is likely to be relaxed,
which violates the assumption of a collisionless system in the Vlasov equations used
to model the mass distribution of the cluster. In both cases, the observations are
of special interest because they are the only direct dynamical measurements of the
mass of possible IMBHs.
Additionally, the Galactic globular cluster NGC 6752 contains two millisecond
pulsars with high, negative spin derivatives in its core as well as two other millisecond
pulsars well into the halo of the cluster at 3.3 and 1.4 times the half mass radius of
the globular cluster (Colpi et al. 2003, 2002). The pulsars in the cluster core can
be explained by a line-of-sight acceleration by 103 Mcircledot of dark mass in the central
0.08 pc (Ferraro et al. 2003). While the pulsars in the outskirts of the cluster can be
explained by exchange interactions with binary stars, the most likely explanation is
that they were kicked from the core in a close interaction with an IMBH, either a
single IMBH or a binary that contains an IMBH (Colpi et al. 2003, 2002).
Finally, high spatial resolution observations of an infrared source known as
IRS 13E near the galactic center reveal that it is a collection of seven stars, six
of which have similar proper motions (Maillard et al. 2004). The similar proper mo-
tion suggests that these stars may be the remaining core of a cluster that is in the
process of becoming tidally disrupted. If one assumes that the two stars with radial
velocity measurements are gravitationally bound to each other, then one infers the
presence of 1300 Mcircledot of unseen mass. Without further observations of radial ve-
locities of other members of this possible cluster, the IMBH interpretation remains
speculative.
12
1.2 IMBH Formation Models
So there is strong evidence in a number of cases for the existence of IMBHs. One
of the most intriguing questions regarding IMBHs is the method of their formation.
They must form differently than stellar-mass black holes, and they are distinct from
supermassive black holes, whose formation is still not completely understood and
may require an IMBH as a seed. A stellar-mass black hole cannot become an IMBH
through Bondi-Hoyle accretion alone since the accretion rate from the ISM is
?MBH ? 4?10?15
parenleftBigg M
50 Mcircledot
parenrightBigg2parenleftbiggn
ISM
cm?3
parenrightbiggparenleftbigg T
106 K
parenrightbigg?3/2
Mcircledot yr?1 (1.6)
where M is the mass of the accreting black hole, nISM is the number density of the
interstellar medium (ISM), and T is the temperature, assuming thermal velocities
dominate the relative speed between the gas and the black hole. Thus a 50 Mcircledot
black hole accreting from the hot ISM withT = 106 K andnISM = 10?2 cm?3 would
accrete at a rate ?MBH = 4?10?17 Mcircledot yr?1, with growth timescaleM/ ?MBH ? 1018 yr
much longer than a Hubble time. A molecular cloud can have nISM = 103 cm?3
and T = 100 K, but radiation from the accretion will heat the gas to T = 104 K
(see Miller & Colbert 2004). A 50 Mcircledot black hole accreting from such a molecular
cloud has growth timescale M/ ?MBH ? 1.3 ? 1010 yr, which is roughly the age of
the Universe but is much longer than the lifetime of a molecular cloud and the
crossing time of the black hole through the cloud. Thus an IMBH requires a more
exotic method of formation. The three most favored models of IMBH formation
are (1) formation from Population III (Pop III) stars, (2) runaway mergers of stars
during the core collapse of a young stellar cluster, and (3) repeated mergers of
stellar-mass black holes.
13
1.2.1 Population III Stars
Madau & Rees (2001) and Schneider et al. (2002) suggest that IMBHs are the
remnants of the massive (M >?200Mcircledot), first generation of stars. The low metallicity
of these Pop III stars precludes cooling through metal line emission, and since the
Jeans mass scales as MJ ? T3/2, Pop III stars may start their main sequence lives
much more massive than solar-metallicity stars thus giving rise to an extremely top-
heavy initial mass function. These large stars avoid the significant mass loss from
stellar winds, which are driven by metal lines, and nuclear-powered radial pulsations,
which are the result of the hotter core temperatures (Baraffe et al. 2001; Fryer et al.
2001; Kudritzki 2000). Thus a large Pop III star may retain almost all of its mass
through its main sequence life. If the star is greater than ? 250Mcircledot, it can avoid the
electron-positron pair instability that results in a completely destructive supernova
and instead collapses directly to a black hole (Madau & Rees 2001).
Although Pop III stars have not been directly observed, there is little doubt that
a first generation of stars with primordial abundance existed. The high-mass end of
any theoretical initial mass function, however, tends to be the most uncertain, and
this is especially true for zero-metallicity stars. In addition, Pop III stellar remnants
would not be associated with either the globular or young stellar clusters that are
associated with ULXs and other IMBH candidates.
1.2.2 Dynamics of Stellar Clusters
Core-Collapse of Young Stellar Cluster
Many studies have found that IMBHs may form in young stellar clusters where a
core collapse leads to direct collisions of stars (Ebisuzaki et al. 2001; Freitag et al.
2006a,b; G?urkan et al. 2006, 2004; Lee 1993; Portegies Zwart & McMillan 2002).
14
Soon after the first stars in a stellar cluster form, the heaviest stars sink to the
center in a process known as mass segregation (e.g., Fregeau et al. 2002; Sigurdsson
& Hernquist 1993). As two objects undergo an encounter, they tend towards energy
equipartition so that more massive objects tend towards lower velocity and thus sink
in the potential well of the cluster. This central concentration of massive stars can
lead to number densities high enough that direct stellar collisions can occur. The
merger remnant of the first collision has a radius much larger than any other star
and therefore the largest cross-section. Thus it is the most likely to suffer another
collision before any other object. This merger remnant will have again a much
larger cross-section, and the collision process becomes a runaway. The collisions
continue until the massive stars explode in supernovae after which the black hole
remnants will be too small for direct stellar collisions to occur, or the runaway stops
when the supply of stars available for collisions is reduced. Numerical simulations
show that masses up to a few 103 Mcircledot (? 10?3 of the original cluster mass) can be
accumulated through runaway collisions (Freitag et al. 2006b; G?urkan et al. 2004;
Portegies Zwart & McMillan 2002). This merger remnant would then presumably
evolve into an IMBH.
This model successfully explains the connection between IMBHs and young stel-
lar clusters, as well as globular clusters, assuming that the IMBH produced in a
young super star cluster would remain in it until it became a globular. In par-
ticular, this model has been shown to successfully describe why there is a ULX
(M82 X-1) associated with the MGG 11 super star cluster and not MGG 9, which,
unlike MGG 11, has a dynamical friction timescale longer than the ? 3 Myr life-
time of the most massive stars. The long dynamical friction timescale prevents the
central density of the cluster from getting high enough to promote collisions before
the stars have evolved. This model, however, does not fully address the evolution
15
of the merger remnant, which is likely to be nontrivial. For example, the merger
remnant may not fully relax from the effects of a collision before the next. The
merger remnant may remain puffy and allow grazing encounters to have the im-
pacting star?s atmosphere tidally stripped while the star?s core escapes. Another
possibility is that during the collisional period, stellar winds drive significant mass
loss perhaps augmented by metal enrichment from high-energy collisions that mix
material into the merger remnant?s atmosphere. If such things hamper growth, this
process may not be able to create a ? 3000 Mcircledot black hole but, perhaps, a more
modest 100?200 Mcircledot.
Repeated Mergers of Stellar-Mass Black Holes
Miller & Hamilton (2002b) proposed that over a Hubble time stellar-mass black
holes in dense globular clusters may grow by mergers to the inferred IMBH masses.
In their model, the heaviest objects in a globular cluster, the black holes (including
those in binaries), sink to the center of the cluster. If the total mass of black holes
in a cluster is large enough that mass segregation is a runaway process, known as
the ?mass stratification instability,? the black holes dynamically decouple from the
rest of the cluster and interact almost entirely with themselves (e.g., O?Leary et al.
2006). In a typical dense globular cluster with reasonable initial mass function, there
will be roughly one 50 Mcircledot black hole that will act as a seed mass for IMBH growth.
Although a 50 Mcircledot is more massive than is thought to be able to form by evolution
of an isolated star of solar metallicity, if a star is more massive than ? 40 Mcircledot
after mass loss, the supernova will fail to eject the atmosphere of the star (Fryer &
Kalogera 2001). This may be possible with the reduced stellar winds in the lower
metallicity environment of a globular cluster. This would not apply to young stellar
clusters with roughly solar metallicity, but in these clusters, non-runaway collisions
16
of stars could be enough to offset mass loss. Because of the large cross-section of
binaries, this seed black hole will interact with binaries most frequently. Through an
exchange bias in which the most massive objects tend to end up with a companion
after a binary-single encounter, the seed black hole will swap into a binary. By virtue
of being in a binary, the seed mass then interacts mainly with individual black holes.
Heggie?s Law (Heggie 1975) states that encounters with hard (tight) binaries tend
to harden (shrink) the binary further. Thus, these encounters tend to take energy
away from the binary and drive the two black holes in the binary toward each other
until gravitational radiation, which takes energy away from the system, is strong
enough to cause the two black holes to merge. The initial seed black hole has then
increased its mass and may repeat the process of acquiring companions, hardening,
and merging until it runs out of black holes with which to merge.
This method of IMBH formation has a clear connection to globular clusters but
not to very young stellar clusters since it will take a longer time than is available in
the youngest stellar clusters suspected to harbor an IMBH. There are three specific
concerns regarding this model: (1) Is this process fast enough to create IMBHs?
(2) Is this process efficient enough to reach 500 Mcircledot or 1000 Mcircledot with the available
supply of stellar-mass black holes in a globular cluster? (3) Can this mechanism
proceed without ejecting the IMBH progenitor? Additionally, one may ask (see
? 1.3.2 for more): (4) What predictions does this model make about the observable
gravitational wave signals from this process? These four questions are the proximal
scientific motivation for this dissertation.
17
1.3 Background
1.3.1 N-body Simulations
The heart of this work is a series ofN-body simulations to test the Miller & Hamilton
(2002b) model of IMBH formation. While there are many unique aspects to these
simulations, the basic idea has not changed since the first published computational
N-body experiments by von Hoerner (1960, see also Aarseth & Lecar 1975), which
stand in stark contrast to the light-bulb-and-photodetector experiments of Holmberg
(1941). In its simplest form, the gravitationalN-body problem is merely the solution
to the equations
d2
dt2xi = G
Nsummationdisplay
jnegationslash=i
mj xj ?xi|x
j ?xi|
3 (1.7)
where xi are the position vectors of the system?s N particles with mass mi. Though
the problem is simply stated, there is no closed-form analytical solution for N > 2
(Bruns 1887; Poincar?e 1892, see also Whittaker 1989 and Julliard-Tosel 2000).
Therefore, numerical solutions are required. Much work has gone into ingenious
improvements in the efficiency of the calculations as well as the stability of these
techniques, especially for large N.
This work concentrates on binary-single scattering experiments. The three-body
problem has been studied extensively, but with every new generation of computing
power, our understanding of this rich but conceptually simple problem advances
with a wider range of numerical simulations and a changing perspective. Previous
studies of the three-body problem have tended to focus on the case of equal or nearly
equal masses (e.g., Heggie 1975; Hut & Bahcall 1983) though other mass ratios have
been studied (e.g., Fullerton & Hills 1982; Heggie et al. 1996; Sigurdsson & Phinney
1993). The nearly equal mass case does not apply to the case of an IMBH in the
18
core of a stellar cluster. In addition the vast majority of previous work has studied
the effect of a single encounter on a binary though Cruz-Gonz?alez & Poveda (1971)
and Cruz-Gonz?alez & Poveda (1972) studied the dissolution of Oort cloud comets
by simulating the effects of a background of field stars on the Sun and a massless
companion. To determine the ultimate fate of an IMBH, simulations of sequences of
encounters are needed. Furthermore, to our knowledge no previous work has consid-
ered the effects of orbital decay due to gravitational radiation between encounters,
which we expect to be important for very tight or highly eccentric binaries.
1.3.2 Gravitational Waves
Although they would occur in any theory of gravity that obeys special relativ-
ity (Einstein 1905, because information about changes in the gravitational field
propotate at or below the speed of light), one of the most fascinating fallouts of
general relativity (Einstein 1915) is the prediction of gravitational waves (Einstein
1916). Gravitational waves, which carry energy, are ripples of curvature in space-
time. Their generation consistent with general relativity is known to exist from
period measurements of the Hulse-Taylor pulsar (Hulse & Taylor 1975; Taylor et al.
1979). Although the detection of gravitational waves (and thus their propagation)
has not been achieved, the first generation of ground-based gravitational wave obser-
vatories is already operational; the next generation of ground-based observatories
is in development; and the first generation of space-based observatories is being
planned. The latter two are expected to detect gravitational waves from astrophys-
ical sources.
The detection of gravitational waves is difficult because they are so weak, and
they are so weak, in part, because gravitational radiation is quadrupolar. This may
be seen by following Misner et al. (1973, see also Goldstein et al. 2002) in looking at
19
the gravitational equivalent of electric dipole radiation (e.g., Jackson 1999), which
has luminosity
Lelec dipole = 23c3e2?x2 = 23c3 ?d2elec (1.8)
where delec = ex is the electric dipole moment of a charge e with position vector x.
For the gravitational equivalent, start by writing the mass dipole moment
d = summationdisplay
A
mAxA (1.9)
where the sum is carried out over all particles in the system. Since the time rate
change of the mass dipole is the total momentum
?d = summationdisplay
A
mA ?xA = p, (1.10)
then there is no dipole gravitational radiation by dint of conversation of momentum,
?d = ?p = 0. In electrodynamics there is also magnetic dipole radiation, which
comes from the second time derivative of the magnetic moment. The gravitational
equivalent of the magnetic moment is
? = summationdisplay
A
x?mAvA = J, (1.11)
the total angular momentum of the system, which is, of course, also conserved. The
next leading term, quadrupole radiation, has luminosity (or power)
LGW quadrupole = 15Gc5
angbracketleftBig...
I2
angbracketrightBig
(1.12)
where I is the trace-free part of the second moment of the mass distribution. It has
elements
Ijk = summationdisplay
A
mA
parenleftbigg
xAjxAk ? 13?jkx2A
parenrightbigg
(1.13)
where ?jk is the Kronecker delta. The third time derivative of I does not generally
vanish, and thus quadrupole radiation is the leading order of gravitational radiation.
20
The most obvious example of a system that generates gravitational waves is
an orbiting binary pair. Dynamically, gravitational waves carry energy away from a
two-body system: causing bound bodies? orbits to shrink, and causing unbound bod-
ies to become bound or at least ?less unbound.? In this dissertation, the dynamical
effects of gravitational radiation are the same as a dissipative drag force.
With increasing evidence in support of the existence of intermediate-mass black
holes (IMBHs), interest in these objects as gravitational wave sources is also growing
(Hopman & Portegies Zwart 2005; Matsubayashi et al. 2004; Miller 2002; Will 2004).
Orbiting black holes are exciting candidates for detectable gravitational waves. At
a distance d a mass m on a circular orbit of size r around a mass M greatermuchm generates
a gravitational wave amplitude of
h ? Gmrc2 v
2
c2 =
G2
c4
Mm
rd
= 4.7?10?25
parenleftBigg M
Mcircledot
parenrightBiggparenleftBigg m
Mcircledot
parenrightBiggparenleftBigg r
AU
parenrightBigg?1parenleftBigg d
kpc
parenrightBigg?1
, (1.14)
where G is the gravitational constant and c is the speed of light. For comparison,
with one-year integrations both the Laser Interferometer Space Antenna (LISA)
and Advanced LIGO are expected to reach down to sensitivities of 10?23 at fre-
quencies of 10 mHz and 100 Hz, respectively. Thus binaries containing IMBHs with
M >? 100 Mcircledot with small separations at favorable distances are strong individual
sources. During inspiral, the frequency of gravitational waves increases as the orbit
shrinks until it reaches the innermost stable circular orbit (ISCO) where the orbit
plunges nearly radially towards coalescence. Because of the quadrupolar nature of
gravitational waves, the gravitational wave frequency for circular binaries is twice
the orbital frequency. At the ISCO for a non-spinning black hole with M greatermuch m,
where rISCO = 6GM/c2 and h ? Gm/6c2d is independent of the mass of the pri-
21
mary, the gravitational wave frequency is
fGW = 2forb = 2
parenleftBigg GM
4pi2r3ISCO
parenrightBigg1/2
? 4400 Hz
parenleftbiggM
circledot
M
parenrightbigg
. (1.15)
Thus a binary with a 100 Mcircledot black hole will pass through the LISA band (10?4
to 100 Hz, Danzmann 2000) and into the bands of ground-based detectors such as
LIGO, VIRGO, GEO-600, and TAMA (101 to 103 Hz, Ando & the TAMA collab-
oration 2002; Barish 2000; Fidecaro & VIRGO Collaboration 1997; Schilling 1998)
whereas a 1000Mcircledot black hole will be detectable by LISA during inspiral but will not
reach high enough frequencies to be detectable by currently planned ground-based
detectors. After the final inspiral phase, the gravitational wave signal goes through
a merger phase, in which the horizons cross, and a ringdown phase, in which the
spacetime relaxes to a Kerr spacetime (Cutler & Thorne 2002; Flanagan & Hughes
1998a,b). The merger and ringdown phases emit gravitational waves at a higher
frequency with a characteristic ringdown frequency of f ? 104(M/Mcircledot)?1 Hz so
that mergers with more massive IMBHs will still be detectable with ground-based
detectors.
If stellar clusters frequently host IMBHs, then currently planned gravitational
wave detectors may detect mergers within a reasonable amount of time. Optimistic
estimates put the upper limit to the Advanced LIGO detection rate of all black holes
in dense stellar clusters at ? 10 yr?1 (O?Leary et al. 2006). The LISA detection
rate for a 1-yr integration and signal-to-noise ratio of S/N = 10 is (Will 2004)
?det ? 10?6
parenleftBigg H
0
70 km s?1Mpc?1
parenrightBigg3parenleftBiggf
tot
0.1
parenrightBigg
?
parenleftBigg ?
10 Mcircledot
parenrightBigg19/8parenleftBigg M
max
100 Mcircledot
parenrightBigg13/4parenleftbigg
lnMmaxM
min
parenrightbigg?1
yr?1, (1.16)
where H0 is the Hubble constant, ftot is the total fraction of globular clusters
that contain IMBHs, ? is the reduced mass of the merging binary, and Mmin
22
and Mmax denote the range in masses of IMBHs in clusters. If we assume that
Mmin = 102 Mcircledot, Mmax = 103 Mcircledot, ? = 10 Mcircledot, ftot = 0.8 (O?Leary et al. 2006), and
H0 = 70 km s?1Mpc?1, then we get a rate of 0.006 yr?1. This, however, implies
that 103 Mcircledot black holes are continuously accreting 10 Mcircledot black holes, which is
unlikely to be the case. Since the distance out to which LISA can detect a given
gravitational wave luminosity DL scales as the square root of the integration time
T, the volume probed scales asV ?D3L ?T3/2. This means that a 10-yr integration
could yield a rate of 0.2 yr?1, and if IMBHs with mass M = 104 Mcircledot are common,
the rate could be much higher. These rates, however, are optimistic and should be
considered as upper limits. A gravitational wave detection of an IMBH with high
signal-to-noise could also yield the spin parameter and thus shed light on the forma-
tion mechanism of the IMBH (Miller 2002). Because detection of inspiral requires
the comparison of the signal to a pre-computed waveform template that depends on
the orbital properties of the binary, knowing the eccentricity distribution is useful.
For e <? 0.2, circular templates are accurate enough to detect the gravitational wave
signal with LIGO (Martel & Poisson 1999), and this is likely to be the case for LISA
as well. A full understanding of the gravitational wave signals from IMBHs requires
a more detailed study of the complicated dynamics and gravitational radiation of
these systems.
1.4 Overview of This Dissertation
This dissertation was designed to answer the following questions.
1. Is the Miller & Hamilton (2002b) model of IMBH formation fast enough to
make an IMBH before the globular cluster?s supply of stellar-mass black holes
eject themselves from the cluster through dynamical interactions (? 0.4 to
23
1 Gyr)?
2. Can repeated mergers of stellar-mass black holes in a globular cluster produce
an IMBH without ejecting too many stellar-mass black holes via the encounters
that harden the IMBH-progenitor binary?
3. Will a globular cluster retain an IMBH progenitor over the course of its merger
history, or will the binary-single scattering events eject the IMBH progenitor
with great regularity?
4. What is the observable gravitational-wave signature of these merging and in-
spiraling systems?
In the process of answering these questions, several more arose and are also addressed
in this dissertation.
5. What is the reason for the apparent broken power-law in the cross-section for
close approach?
6. How does the cross-section for close approach change when gravitational waves
are included?
7. What is the likelihood of merger because of gravitational radiation during
a binary-single encounter when the leading-order terms of energy loss from
gravitational waves are included?
We investigate these questions with simulations of sequences of binary-single en-
counters as they would occur for black holes in dense stellar clusters. In Chapter 2
we introduce the code (? 2.2) used to integrate the sequences of three-body interac-
tions as well as the major assumptions involved. We then present results of purely
Newtonian simulations and simulations that include gravitational waves between en-
counters, and discuss the implications for IMBH formation and gravitational wave
24
detection. In Chapter 3 we describe the changes to the code necessary to include
gravitational waves during the encounter. We look into close-approach and merger
cross sections in ? 3.3, and present results of simulations that include gravitational
waves as well as their implications for IMBH formation and gravitational wave de-
tection. In Chapter 4 we summarize our main results and conclusions.
25
Chapter 2
Growth of Intermediate-Mass
Black Holes in Globular Clusters
2.1 Overview
In this chapter we present numerical simulations of sequences of high-mass ratio
binary-single encounters. We describe the code used to simulate the encounters in
? 2.2. Next, we present results of the simulations of sequences of encounters on a
range of mass ratios with Newtonian gravity (? 2.3.1) and with gravitational radi-
ation between encounters (? 2.3.2) and show that including gravitational radiation
decreases the duration of the sequence by ? 30% to 40%. In ? 2.4 and ? 2.5 we
discuss the implications of these results for IMBH formation and gravitational wave
detection.
26
2.2 Numerical Method
We perform numerical simulations of the interactions of a massive binary in a stel-
lar cluster. Simulating the full cluster is beyond current N-body techniques, so we
focus instead on a sequence of three-body encounters. Massive cluster objects, such
as IMBHs and tight binary systems, tend to sink to the centers of clusters so that
a single IMBH is very likely to meet a binary (Sigurdsson & Phinney 1995). Ex-
changes in which the IMBH acquires a close companion are common. Such a binary
in a stellar cluster core will experience repeated interactions with additional objects
as long as the recoils from these interactions do not eject the binary. Therefore,
we simulate a sequence of encounters between a hard binary and an interloper. We
perform one interaction and then use the resulting binary for the next encounter.
This is repeated multiple times until the binary finally merges because of gravi-
tational radiation. Because typical velocities involved are non-relativistic and the
black holes are tiny compared to their separations, they are treated as Newtonian
point masses. In Chapter 3 we see how this changes when gravitational radiation
is included throughout. In order to test the influence of the binary?s mass, we use
a range of binary mass ratios. To simplify the problem we study a binary with
mass ratio of N:10 Mcircledot and a 10 Mcircledot interloper, designated as N:10:10, and vary N
between 10 Mcircledot and 103 Mcircledot.
The simulations were done using a binary-single scattering code that was writ-
ten to be as general purpose as possible. See Appendix D for examples of other
applications. Because of the vast parameter space that needs to be covered, the
code uses a Monte Carlo initial condition generator. The orbits are integrated using
HNBody, a hierarchical, direct N-body integrator, with the adaptive fourth-order
27
Runge Kutta integrator option (K. Rauch & D. Hamilton, in preparation)1. Because
we focus on close approaches where a wide range of timescales are important, an
adaptive scheme is often better than symplectic methods, which are time-reversible
and bound the error in energy conservation for a fixed timestep.
In wide hierarchical triples, direct integration can consume a large amount of
computational time. To reduce this, we employ a two-body approximation scheme
that tracks the phase of the inner binary. For a sufficiently large outer orbit, the
orbit is approximately that of an object about the center of mass of the binary. If
both the distance from the outer object to the inner binary?s center of mass and
the semimajor axis of the outer binary are greater than 30 times the inner binary?s
semimajor axis, aib, we calculate this approximate two-body orbit analytically and
keep track of the inner binary?s phase. When the outer object nears the binary
again, we revert to direct numerical integration.
The orbit is integrated until one of three conditions is met: (1) one mass is at a
distance of least 30aib and is departing along a hyperbolic path, (2) the system forms
a hierarchical triple with outer semimajor axis greater than 2000 AU, an orbit so
large that it would likely be perturbed in the high density of a cluster core and not
return, (3) an ionization2 or merger is detected, or (4) the integration is prohibitively
long, in which case the encounter is discarded and restarted with new randomly
generated initial conditions. Roughly 10?4 of all encounters had to be restarted with
most occurring for higher mass ratios where resonant encounters (encounters that
have more than one close approach and are not simple fly-bys) are more common.
1See http://janus.astro.umd.edu/HNBody/.
2The total energy of the binary-single encounters simulated in this dissertation is negative,
which precludes ionization as a possibility, but for the sake of generality this test is included. For
an example in which ionization of the binary is both possible and significant, see ? D.2.
28
It is an unfortunate reality of simulations such as ours that the integration of some
small sample of encounters will require a longer time than is practical. We tested the
effect of our arbitrary limit by extending the allowed integration time by a factor of
102 for one hundred sequences form0 = 100 and 1000Mcircledot. There were no statistically
significant differences. In particular, the average time per sequence, average final
semimajor axis, and average final eccentricity all differed by less than 2%. In half
of the simulations presented in this chapter, we evolve the binary?s orbit due to
gravitational wave emission after each encounter. Since a binary in a cluster spends
most of its time and emits most of its gravitational radiation while waiting for an
encounter rather than during an interaction, we only include gravitational radiation
between encounters. The effects of adding gravitational waves during the encounter
are discussed in Chapter 3. To isolate this effect, we run simulations both with
and without gravitational radiation between encounters. We include gravitational
radiation by utilizing orbit-averaged expressions for the change in semimajor axis a
and eccentricity e with respect to time (Peters 1964):
da
dt = ?
64
5
G3m0m1 (m0 +m1)
c5a3 (1?e2)7/2
parenleftbigg
1 + 7324e2 + 3796e4
parenrightbigg
(2.1)
and
de
dt = ?
304
15
G3m0m1 (m0 +m1)
c5a4 (1?e2)5/2
parenleftbigg
e+ 121304e3
parenrightbigg
, (2.2)
where m0 and m1 (m0 ?m1) are the gravitational masses of the binary pair. Here
G is the gravitational constant, and c is the speed of light. The orbital elements
are evolved until the next encounter takes place, at a time that we choose randomly
from an exponential distribution with a mean encounter time, ??enc? = 1/?nv???,
where n is the number density of objects in the cluster?s core, v? is the relative
velocity, and ? is the cross-section of the binary. If we assume the mass of the
29
binary m0 +m1 greatermuchm2, then
? ?pir2p + 2pirpG(m0 +m1)/v2?, (2.3)
where rp is the maximum considered close approach of m2 to the binary?s center of
mass. For a thermal distribution of stellar speeds, v? = (mavg/m2)1/2vms, where
mavg = 0.4 Mcircledot is the average mass of the main sequence star and vms is the main
sequence velocity dispersion. In our simulations, the second term of Equation 2.3,
gravitational focusing, dominates over the first. Averaging over velocity (assumed
to be Maxwellian) we find
??enc? = 2?107
parenleftbigg v
ms
10 km s?1
parenrightbiggparenleftBigg106 pc?3
n
parenrightBiggparenleftBigg1 AU
rp
parenrightBiggparenleftbigg 1 M
circledot
m0 +m1
parenrightbiggparenleftbigg1 M
circledot
m2
parenrightbigg1/2
yr.
(2.4)
We then subject the binary to another encounter using orbital parameters adjusted
by both the previous encounter and the gravitational radiation emitted between
the encounters. This sequence of encounters continues until the binary merges due
to gravitational wave emission. If orbital decay is not being calculated, then we
determine that the binary has merged when the randomly drawn encounter time is
longer than the timescale to merger, which is approximately
?merge ? 6?1017 (1 Mcircledot)
3
m0m1 (m0 +m1)
parenleftbigg a
1 AU
parenrightbigg4parenleftBig
1?e2
parenrightBig7/2
yr (2.5)
for the high eccentricities of importance in this dissertation.
For binaries with unequal masses, gravitational wave emission before coalescence
imparts a recoil velocity on the binary (Favata et al. 2004). As two objects with
unequal masses or with misaligned spins spiral in towards each other, asymmetric
emission of gravitational radiation produces a recoil velocity on the center of mass of
the binary. Most of the recoil comes from contributions after the masses are inside
30
of the ISCO, where post-Newtonian analysis becomes difficult (Favata et al. 2004).
For non-spinning black holes, the velocity kick from the recoil up to the ISCO is
(Favata et al. 2004)
vr = 15.6 km s?1f(q)f
max
(2.6)
whereq = m1/m0 < 1,f(q) = q2(1?q)/(1+q)5, andfmax ?f(0.38) = 0.018. Favata
et al. (2004) bounded the total recoil to between 20 km s?1 ?vr ? 200 km s?1 for
non-spinning black holes with q = 0.127. Since the recoil velocity scales as q2 for
q lessmuch 1, this may be scaled to other mass ratios. More recently, Blanchet et al. (2005)
argued from high order post-Newtonian expansions that the kick speed for very
small mass ratios q lessmuch 1 was vr/c = 0.043q2, with an uncertainty of roughly 20%.
This is consistent with the results of Favata et al. (2004), but as most of the recoil
originates well inside the ISCO, Blanchet et al. (2005) caution that numerical results
are probably required for definitive answers. Using the effective one-body approach,
Damour & Gopakumar (2006) find a maximum recoil velocity of vr = 74 km s?1.
Most recently, Baker et al. (2006) presented the first accurate, numerical relativity
simulations of the merger of two non-spinning black holes, and for q = 2/3 found
an estimated recoil velocity of vr = 105 km s?1 with an error of less than 10%. All
calculations are consistent with the wide range first reported by Favata et al. (2004).
Because of the disparity among the several calculations, we chose not to include this
effect in our simulations, but we comment on the possible effects it could have in
? 2.4 and ? 3.4.1.
Global energy and angular momentum are monitored to ensure accurate inte-
gration. The code also keeps track of the duration of encounters, the time between
encounters, changes in semimajor axis and eccentricity, and exchanges (events in
which the interloping mass replaces one of the original members of the binary and
the replaced member escapes).
31
As a test of our code, we compared simulations of several individual three-body
encounters with the work of Heggie et al. (1996). As part of a series of works exam-
ining binary-single star scattering events, Heggie et al. (1996) performed numerical
simulations of very hard binaries with a wide range of mass ratios and calculated
their cross-sections for exchange. We ran simulations of one encounter each of a
sample of mass ratios for comparison. To facilitate comparison of encounters with
differing masses, semimajor axes, and relative velocities of hard binaries, Heggie
et al. (1996) use a dimensionless cross-section,
??ex = 2v
2
??
piG(m0 +m1 +m2)a, (2.7)
where v? is the relative velocity of the interloper and the binary?s center of mass
at infinity and ? is the physical cross-section for exchanges. We calculate ? as the
product of the fraction of encounters that result in an exchange (fex) and the total
cross-section of encounters considered: fexpib2max, where bmax is an impact parameter
large enough to encompass all exchange reactions. Our cross-sections are in agree-
ment with those of Heggie et al. (1996) within the combined statistical uncertainty
as seen in Table 2.1.
2.3 Simulations and Results
We used our code to run numerical experiments of three-body encounter sequences
with a variety of mass ratios. The binaries consisted of a dominant body with mass
m0 = 10, 20, 30, 50, 100, 200, 300, 500, or 1000 Mcircledot and a secondary of mass m1 =
10 Mcircledot. Because of mass segregation, the objects that the binary encounters will be
the heaviest objects in the cluster. In order to simplify the problem, we consider
only interactions with interlopers of mass m2 = 10 Mcircledot. The binary starts with a
circular a = 10 AU orbit, and the interloper has a relative speed at infinity of v? =
32
Table 2.1. Single Encounter Cross-sections for Exchange.
m0:m1:m2 Ejected Mass HHM96 This Work
10 : 1 : 1 1 1.086 ? .105 1.054 ? .023
10 ?? ??
10 : 1 : 10 1 7.741 ? .360 7.825 ? .255
10 0.513 ? .087 0.520 ? .043
3 : 1 : 1 1 2.465 ? .170 2.311 ? .073
3 0.072 ? .025 0.059 ? .007
Note. ? This table compares dimensionless cross-sections
for exchange ??ex (see text for details) calculated by Heggie
et al. 1996 and by us. The first column lists the masses, with
binary componentsm0 andm1. Columntwoshowsthemass of
the ejected object. The ejection of the smaller mass is energet-
ically favored so it always has a larger cross-section. There is
general agreement between the two calculations to within the
statistical uncertainty, which we calculate as ??/N1/2ex , where
Nex is the total number of exchanges.
10 km s?1 and an impact parameter, b, relative to the center of mass of the binary
such that the pericenter distance of the hyperbolic encounter would range from rp =
0 to 5a. In order to represent an isotropic distribution of encounters, the distribution
of impact parameters is P(b) ? b. For all binaries, vcirc = [G(m0 +m1)/a]1/2 ?
40km s?1 greatermuchv?, andthus all are considered hard. The Monte Carlo initial condition
generator distributes the orientations and directions of encounters isotropically in
space, and the initial phase of the binary is randomized such that it is distributed
equally in time. In a multi-mass King model for a stellar cluster, the scale height
of bodies with mass m scales as (mavg/m)1/2 (Sigurdsson & Phinney 1995). Thus
black holes with mass m = 10Mcircledot have a scale height 1/5 of the scale height of the
average-mass main sequence star, mavg ? 0.4 Mcircledot. The volume occupied by black
holes, then, is more than 100 times smaller than the volume occupied by average-
mass main sequence stars, and if black holes are at least 10?2 times as numerous,
which would be expected for a Salpeter mass function, then their core number
33
density is at least comparable to that of visible stars (Miller & Hamilton 2002b).
Thus we assume the cluster core has a density ofn = 105 pc?3 and an escape velocity
of vesc = 50 km s?1 for the duration of the simulation. We discuss the consequences
of changing the escape velocity in ? 2.4. For each mass ratio, we simulate 1000
sequences with and without gravitational radiation between encounters.
2.3.1 Pure Newtonian Sequences
Figure 2.1a shows the change of semimajor axis and pericenter distance as a function
of time over the course of a typical Newtonian sequence. The encounters themselves
take much less time then the period between encounters, so a binary spends virtually
all its time waiting for an interloper. Most of the time in this example is spent hard-
ening the orbit from 1 AU to 0.4 AU because as the binary shrinks, its cross-section
decreases and the timescale to the next encounter increases. Figure 2.1b shows the
same sequence plotted as a function of number of encounters. The semimajor axis
decreases by a roughly constant factor with each encounter. This is expected for a
hard binary, which, according to Heggie?s Law (Heggie 1975), tends to harden with
each encounter at a rate independent of its hardness. The eccentricity and therefore
the pericenter distance, rp = a(1?e), however, can change dramatically in a single
encounter (for a discussion on eccentricity change of a binary in a cluster, see Heggie
& Rasio 1996). This sequence ends with a very high eccentricity (e = 0.968), which
reduces the merger time given by Equation 2.5 to less than ?enc.
Table 2.2 summarizes our main results and shows a number of interesting trends.
The average number of encounters per sequence, ?nenc?, increases with increasing
mass ratio since the energy that the interloper can carry away scales as ?E/E ?
m1/(m0 +m1) and since nenc ?E/?E for a constant eccentricity (Quinlan 1996).
Energy conservation assures that every hardening event results in an increased rel-
34
Figure 2.1: Newtonian 1000:10:10 sequence. These panels show the semimajor
axes (upper curves) and pericenter distances (lower curves) as functions of time
(left panel) and number of encounters (right panel) for one sequence of encounters
with no gravitational wave emission. Each change in a and rp is the result of a
three-body encounter. Since the binary is hard, the semimajor axis gradually
tightens by a roughly constant fractional amount per encounter with most of
the time spent hardening the final fraction when close encounters are rare. The
pericenter distance, however, fluctuates greatly due to large changes in eccentricity
during a single encounter. The sequence ends at a very high eccentricity when
the binary would merge due to gravitational radiation before the next encounter.
ative velocity between the binary and the single black hole. If the velocity of the
single black hole relative to the barycenter, and thus the globular cluster, is greater
than the escape velocity of the cluster core (typically vesc = 50 km s?1 for a dense
cluster; see Webbink 1985), then the single mass will be ejected from the cluster.
The average number of ejected masses per sequence, ?nej?, also increases with in-
creasing mass ratio because the higher mass ratio sequences have a larger number of
encounters and because the larger mass at a given semimajor axis has more energy
for the interloper to tap. Conservation of momentum guarantees that if a mass is
ejected from the cluster at high enough velocity, the binary will also be ejected. Ta-
ble 2.2 lists ?fbinej?, the fraction of sequences that result in the ejection of the binary
35
Table 2.2. Sequence Statistics.
m0 Case ?nenc? ?nej? ?fbinej? ?tseq?/106 yr ?af?/AU ?ef?
10 Newt. 51.6 3.9 0.880 82.72 0.164 0.929
GR Evol. 48.7 3.7 0.839 59.89 0.190 0.901
20 Newt. 51.3 6.5 0.835 65.94 0.178 0.924
GR Evol. 47.1 6.1 0.776 43.46 0.230 0.898
30 Newt. 58.9 9.3 0.753 49.11 0.198 0.926
GR Evol. 55.1 8.6 0.676 31.89 0.222 0.892
50 Newt. 73.2 14.6 0.581 33.75 0.230 0.919
GR Evol. 66.7 13.0 0.455 22.73 0.285 0.892
100 Newt. 102.0 24.0 0.229 21.35 0.327 0.936
GR Evol. 93.4 20.1 0.161 14.97 0.357 0.873
200 Newt. 158.4 38.2 0.043 15.13 0.387 0.938
GR Evol. 140.3 31.5 0.026 9.998 0.444 0.872
300 Newt. 208.5 49.1 0.013 11.89 0.468 0.943
GR Evol. 184.0 39.4 0.006 7.822 0.445 0.874
500 Newt. 308.7 71.1 0.001 9.920 0.528 0.944
GR Evol. 269.1 54.9 0 6.225 0.488 0.860
1000 Newt. 562.4 117.3 0 7.363 0.641 0.953
GR Evol. 483.0 88.9 0 4.427 0.556 0.851
Note. ? Table 2.2 summarizes the main results of our simulations of sequences of
three-body encounters. For each dominant mass, m0, we ran 1000 sequences of pure
Newtonian encounters (Newt.) and 1000 sequences of the more realistic Newtonian
encounters with gravitational radiation between encounters (GR Evol.). The columns
list the average number of encounters per sequence ?nenc?, the average number of black
holes ejected from the cluster in each sequence ?nej?, the fraction of sequences in which
the binary is ejected from the cluster, ?fbinej?, the average total time for the sequence
?tseq?, the average final semimajor axis ?af?, and the average final eccentricity ?ef?.
from the cluster. As expected, the fraction decreases sharply with increasing mass
such that virtually none of the binaries with mass greater than 300 Mcircledot escapes the
cluster.
The shape and size of the orbit after its last encounter determine the dominant
gravitational wave emission during the inspiral and are of particular interest to
us. The distribution of pre-merger semimajor axes for all mass ratios is shown in
Figure 2.2. The distributions all have a similar shape that drops off at lowabecause
the binary tends to merge before another encounter can harden it. For large orbits
36
Figure 2.2: Histograms of final semimajor axes for all mass ratios. The solid his-
tograms are pure Newtonian sequences, and the hatched histograms are sequences
with gravitational radiation between encounters. The histograms all have similar
shapes with a sharp drop at low a since the binary tends to merge before another
encounter can harden it, and they have long tails at high a where the binary will
only merge with high eccentricity. The sequences with gravitational radiation
have falloffs at smaller a than those without due to both the circularization and
the extra source of hardening.
the binary will only merge for a high eccentricity, and thus there are long tails in
the histograms towards high a from encounters that resulted in an extremely high
eccentricity. The distributions for lower mass ratios are shifted to smaller a because
for a given orbit, a less massive binary will take longer to merge. This can also be
seen in the mean final semimajor axis, ?af?, in Table 2.2.
Figure 2.3 shows the distribution of binary eccentricities after the final encounter
for one mass ratio. The plot is strongly peaked near e = 1, a property shared by
37
all other mass ratios. This distribution is definitely not thermal, a distribution
in which the probability scales as P(e) = 2e and the mean eccentricity of which
is ?e?th ? 0.7. The high eccentricity before merger results from both the strong
dependence of merger time on eccentricity and the fact that the eccentricity can
change drastically in a single encounter (see Figure 2.1). As the semimajor axis
decreases by roughly the same fractional amount in each encounter, the eccentricity
increases and decreases by potentially large amounts with each strong encounter.
When the eccentricity happens to reach a large value, the binary will merge before
the next encounter. Figure 2.4 shows the eccentricity distribution for all encounters
after the first 10 for all 1000 sequences with a mass ratio of 1000:10:10. The first
ten encounters are excluded so that the binary has faced sufficient encounters to
thermalize. The distribution is roughly thermal up to high eccentricity where the
binaries merge. Thus merger selectively removes high eccentricity binaries from a
thermal distribution.
2.3.2 General Relativistic Binary Evolution
The addition of gravitational radiation between Newtonian encounters is expected
to alter a sequence since it is an extra source of hardening and since it circularizes
the binary. Figure 2.5 shows a typical sequence for the 1000:10:10 mass ratio includ-
ing gravitational radiation. Three-body interactions drive the binary?s eccentricity
up to e = 0.959 and its semimajor axis down to a = 0.713 AU. Then starting at
t = 2.2?106 yr over the course of about ten interactions that only weakly affect the
eccentricity and semimajor axis, gravitational radiation causes the orbit to decay
to a = 0.550 AU and e = 0.946 while the pericenter distance remains roughly con-
stant. The corresponding semimajor axis change in the Newtonian-only sequences
in Figure 2.1 takes 45 encounters and more than twice as long although one must
38
Figure 2.3: Histogram of final eccentricities for 1000:10:10 mass ratio. The solid
histogram is from pure Newtonian sequences, and the hatched histogram is from
sequences with gravitational radiation between encounters. The histogram is cut
at e = 0.8 because ef < 0.8 is rare. The histograms have roughly the same shape
for both cases and for all mass ratios although the gravitational wave sequences
have a consistently lower mean at higher mass ratios because gravitational wave
emission damps eccentricities. The histograms show a decidedly non-thermal
distribution and are strongly peaked near e = 1. Because the timescale to merge
due to gravitational radiation is so strongly dependent on e, the binary will merge
when it happens to reach a high eccentricity.
be careful when comparing two individual sequences. Gravitational waves make the
most difference when the pericenter distance is small, which is guaranteed at the
end of a sequence, but can also happen in the middle as Figure 2.5 shows.
Table 2.2 summarizes the effect of adding gravitational radiation. In general the
effect is greater at higher masses because gravitational radiation is stronger for a
given orbit. Because of the extra energy sink, the binaries merge with fewer encoun-
39
Figure 2.4: Solid line is a histogram of all eccentricities after each encounter
except for the first ten for all pure Newtonian sequences of 1000:10:10. The dashed
line is a thermal distribution of eccentricities. The distribution is roughly thermal
for low eccentricity but deviates for e >? 0.6. The expected thermal distribution
of eccentricities is altered by losses of high eccentricity orbits to merger.
ters, fewer black holes are ejected, and the fraction of sequences in which a binary
is ejected is smaller. The most dramatic change is in the duration of the sequence,
which gravitational radiation reduces by 27% to 40%. The distributions of final
semimajor axes (Figure 2.2) and final eccentricities (Figure 2.3) have similar shapes
to the Newtonian-only distributions. Due to the circularizing effect of gravitational
radiation, binaries of all mass ratios merge with a smaller ?ef? than Newtonian-only
sequences with the largest difference at high mass ratios. Gravitational radiation
also produces a smaller ?af? for m0 >? 300 Mcircledot. This can be seen in Figure 2.2 where
the gravitational radiation simulations display an excess number of sequences with
40
Figure 2.5: 1000:10:10 gravitational radiation sequence. Same as Figure 2.1a
but for a sequence with gravitational radiation between encounters. The effects
of gravitational radiation can be seen between 2.2 and 2.4 ? 106 years. Over
this period, the binary undergoes about ten interactions that do not significantly
affect its orbit. During this time, the semimajor axis decays from a = 0.713 AU to
0.550 AU while the pericenter distance remains small and roughly constant. When
an encounter reduces the eccentricity at 2.4 ? 106 years, gravitational radiation
is strongly reduced. Gravitational radiation becomes important again at the end
of the sequence. The sequence ends with the binary?s merger from gravitational
waves.
low af, which is a consequence of the binaries? lower ef.
41
2.4 ImplicationsforIMBHFormationandGrowth
We can use these simulations to test the Miller & Hamilton (2002b) model of IMBH
formation. We assume that a 50 Mcircledot seed black hole with a 10 Mcircledot companion will
undergo repeated three-body encounters with 10 Mcircledot interloping black holes in a
globular cluster with vesc = 50 km s?1 and n = 105 pc?3. We also assume that
the density of the cluster core remains constant as the IMBH grows. We then test
whether the model of Miller & Hamilton (2002b) can build up to IMBH masses,
which we take to be 103 Mcircledot, (1) without ejecting too many black holes from the
cluster, (2) without ejecting the IMBH from the cluster, and (3) within the lifetime
of the globular cluster. We also test how answers to these questions depend on escape
velocity and seed mass. It should be noted that the 1000 Mcircledot canonical IMBH is not
required to explain the observations of ULXs and is taken as a round number, and
if the kinematical evidence of black holes in G1 and M15 with mass M > 103 Mcircledot
to 104 Mcircledot is flawed as has been suggested, then IMBHs with masses M <? 500 Mcircledot
are sufficient to explain the ULXs. We return to this point in Chapter 4.
If the number of black holes ejected is greater than the total number of black
holes in the cluster core, then the IMBH cannot build up to the required mass by
accreting black holes alone. To calculate the total number of black holes ejected
while building up to large masses, we sum the average number of ejections using a
linear interpolation of the values in Table 2.2. Assuming a cluster escape velocity
of vesc = 50 km s?1, we find that the total number of black holes ejected when
building up to 1000 Mcircledot is approximately 6800 for our Newtonian-only and 5300 for
gravitational radiation simulations. This is far greater than the estimated 102 to 103
black holes available (Portegies Zwart & McMillan 2000). If there were initially one
thousand 10 Mcircledot black holes in the cluster, mergers of the massive black hole with a
42
series of 10Mcircledot black holes would exhaust half of the black holes in? 2.6?108 yr and
would ultimately produce a 240 Mcircledot black hole. Increasing the seed mass increases
the final mass of the IMBH when half of the field black holes run out. If the seed
mass were 100, 200, or 300Mcircledot, then the model would produce a 270, 330, or 410Mcircledot
black hole after exhausting half of the cluster black hole population in 1.9, 1.1, or
0.8 ? 108 yr, respectively. Figure 2.6 shows the number of black holes ejected as
a function of initial black hole mass for a range of escape velocities. Gravitational
radiation recoil velocity, which we did not simulate, however, would increase the
number of ejections at small masses. These growth times are much shorter than
the ? 0.4 ? 1.0 ? 109 yr necessary for stellar-mass black holes to eject each other
from the cluster (O?Leary et al. 2006; Portegies Zwart & McMillan 2000; Sigurdsson
& Hernquist 1993). Therefore, self-depletion of stellar-mass black holes is not a
limiting factor.
Of particular concern is whether the three-body scattering events will eject the
binary from the cluster. The black hole can only merge with other black holes while
it is in a dense stellar environment. The probability of remaining in the cluster after
one sequence is P = 1??fbinej?. As can be seen in Table 2.2, once the black hole
has built up to ? 300 Mcircledot, it is virtually guaranteed to remain in the cluster. When
starting with 50 Mcircledot, we calculate the total probability of building up to 300 Mcircledot to
be 0.0356. Figure 2.7 shows the probability of building up to 300 Mcircledot as a function
of starting mass for different escape velocities for the gravitational radiation case.
Table 2.3 lists probabilities for selected seed masses and escape velocities for the
gravitational radiation case.
In a similar manner, we calculate the total time to build up to 1000 Mcircledot, as-
suming that the supply of stellar-mass black holes and density remain constant,
an assumption which leads to an underestimation of the time. While the time per
43
Figure 2.6: Plot of total number of black holes ejected in building up to 1000 Mcircledot
as a function of seed mass for the gravitational radiation case assuming different
escape velocities. The four curves show different assumed cluster escape velocities
in km s?1. For all but the largest seed masses, the number of black holes ejected
is greater than the estimated ? 103 (indicated by the dashed line) present in a
young globular cluster.
merger is larger for the smaller masses, the total time is dominated at the higher
masses since more mergers are needed for the same fractional increase in mass. For
Newtonian-only simulations the total time is 1.1?109 yr, and for simulations with
gravitational radiation the total time is 7.1?108 yr. These are much less than the
age of the host globular clusters. Figure 2.8 shows the time to reach a specified
mass for both the Newtonian and gravitational radiation cases. If the interactions
of stellar-mass black holes with IMBHs or interactions among themselves eject the
black holes from the cluster core but not from the cluster, however, the number
44
Figure 2.7: Plot of an IMBH?s probability of remaining in the cluster and building
up to 300 Mcircledot as a function of starting mass of the dominant black hole for the
gravitational radiation case assuming different escape velocities labeled in km s?1.
Once the black hole has built up to 300 Mcircledot it is very unlikely that it will be ejected
from the cluster. The lowest mass binaries are much more readily ejected and thus
are very unlikely to survive a sequence of encounters. Miller & Hamilton (2002b)
suggest that IMBHs can be built in this manner with a starting mass ? 50 Mcircledot.
We find that such small initial masses are likely to be ejected from the cluster
core for reasonable escape velocities of dense clusters.
density and therefore the rate of growth decrease. This analysis of the growth time
ignores the time it takes the IMBH progenitor to acquire a companion, but this is
unlikely to cause a large increase in the total time. See Appendix A for an estimate
of the extra time required.
Although there is clearly enough time to build IMBHs as Miller & Hamilton
(2002b) propose, the issues of whether there are enough stellar-mass black holes and
whether the cluster will hold onto the IMBH remain. The combination of an initial
45
Table 2.3. IMBH Formation.
Seed Mass vesc Probability to remain Number of Time
(Mcircledot) (km s?1) in cluster BH ejections (108 yr)
50.0 40.0 0.00264 6414 7.06
50.0 0.0356 5276
60.0 0.129 4038
70.0 0.269 3573
100.0 40.0 0.0821 6312 6.15
50.0 0.290 5188
60.0 0.525 3963
70.0 0.698 3606
200.0 40.0 0.670 5995 4.93
50.0 0.842 4922
60.0 0.932 4077
70.0 0.978 3417
300.0 40.0 1.000 5561 4.05
50.0 1.000 4564
60.0 1.000 3777
70.0 1.000 3164
Note. ? This table lists values for selected seed masses and cluster escape
velocities for the gravitational radiation case. Column 3 lists the probability
for the IMBH to remain in the cluster until it reaches a mass of 300 Mcircledot. The
fourth column lists the total number of black holes ejected in building up to
1000 Mcircledot. Column 5 lists the total time to build up to 1000 Mcircledot. The total
time is not affected by the escape velocity because the density of black holes
in the cluster core is taken to be constant.
mass of 50Mcircledot and an escape velocity of 50 km s?1 is not likely to produce an IMBH
in a globular cluster through three-body interactions with 10Mcircledot black holes, but the
general process could still produce IMBHs. Miller & Hamilton (2002b) argued that
a seed mass of 50 Mcircledot would be retained, but for analytical simplicity they assumed
that every encounter changed the semimajor axis by the same fractional amount
??a/a?. Some encounters, however, can decrease the semimajor axis by several
times the average value and thus impart much larger kicks. The authors therefore
underestimated the minimum initial mass necessary to remain in the cluster. A
hierarchical merging of stellar-mass black holes could, however, still produce an
46
Figure 2.8: Plot of total time to build up to a certain mass when built by mergers
with 10 Mcircledot black holes for Newtonian only results and for runs with gravitational
radiation between encounters. The Newtonian only simulations are slower to build
up, but both cases reach 1000 Mcircledot within about 109 years. The time plotted
assumes a constant density of black holes for the duration of IMBH formation.
IMBH if (1) the initial mass of the black hole were greater than 50Mcircledot, (2) the escape
velocity of the cluster were greater than 50 km s?1, or (3) additional dynamics were
involved. We consider each of these in turn.
If the mass of the initial black hole were, e.g., 250 Mcircledot before the onset of
compact object dynamics, dynamical kicks would not be likely to eject the IMBH,
and it would require fewer mergers to reach 1000 Mcircledot and thus a smaller population
of stellar mass black holes. The initial black hole could start with such a mass if
it evolved from a massive Population III star or from a runaway collision of main
sequence stars (G?urkan et al. 2004; Portegies Zwart & McMillan 2002), or it could
47
reach such a mass by accretion of young massive stars, which would be torn apart
by tidal forces and impart little dynamical kick.
Retention rates could also be increased by the higher escape velocities found in
some globular clusters. NGC 6388, e.g., has an escape velocity of vesc = 78 km s?1
(Webbink 1985). If the escape velocity werevesc = 70 km s?1, the interactions result
in a smaller fraction of ejected binaries and the probability of building from 50 Mcircledot
to 1000 Mcircledot then increases by almost an order of magnitude.
In addition, processes with lower dynamical kicks could prevent ejection. One
promising mechanism is the Kozai resonance (Kozai 1962; Miller & Hamilton 2002a;
Wen 2003). If a stable hierarchical triple is formed, then resonant processes can
pump up the inner binary?s eccentricity high enough so that it would quickly merge
due to gravitational radiation and without any dynamical kick to eject the IMBH
from the cluster. O?Leary et al. (2006) include Kozai-resonance-induced mergers
in their simulations and find that it reduces the number of stellar-mass black holes
ejected, but only by about 10%. Two-body captures (captures in which an interloper
passes close enough to the isolated IMBH that it becomes bound and merges due
to gravitational radiation) would also result in mergers without dynamical kicks.
These kinds of two-body captures that occur during a larger three-body encounter
are included in Chapter 3. Both Kozai-resonance-induced mergers and two-body
captures are devoid of dynamical kicks, but they would suffer a gravitational radia-
tion recoil. A system in which a 10Mcircledot black hole merges into a 130Mcircledot non-rotating
black hole would have a recoil velocity 20 km s?1 ? vr ? 200 km s?1 (Baker et al.
2006; Blanchet et al. 2005; Damour & Gopakumar 2006; Favata et al. 2004). Since
vr ? (m1/m0)2, a merger between a 10 Mcircledot black hole and a seed black hole of mass
of 250 Mcircledot, as discussed above, would experience a recoil velocity <? 50 km s?1.
Mergers with lower mass objects that are torn apart by tidal forces, such as white
48
dwarfs, would impart no gravitational radiation recoil. Finally, a range of interloper
masses instead of the simplified single mass population that we used here may also
affect retention statistics since a smaller interloper would impart smaller kicks while
still contributing to hardening.
Increasing the seed mass and the escape velocity will reduce the number of field
black holes ejected but not by enough. As seen in Figure 2.6, using a seed mass
m0 = 250 Mcircledot and an escape velocity vesc = 70 km s?1 reduces the number of black
holes ejected by 40%, but this is still several times more than are available. The
Kozai-resonance-induced mergers and two-body captures, however, are methods of
merging without possibility of ejecting stellar-mass black holes. In order to reach
our canonical 1000 Mcircledot intermediate mass while ejecting fewer than 103 black holes,
either 70-80% of the mergers must come from these ejectionless methods, or else
there must exist an extra method of hardening.
2.5 Implications for Gravitational Wave Detec-
tion
Our simulations make predictions interesting for gravitational wave detection. After
the last encounter of a sequence, the binary will merge due to gravitational radiation.
As mentioned in ? 1.3.2, as the binary shrinks and circularizes, the frequency of the
gravitational radiation emitted passes through the LISA band and then through
the bands of ground-based detectors. By the time the binaries are detectable by
ground-based instruments, they will have completely circularized, but while in the
LISA band, some will have measurable eccentricities. We calculate the distribution
of eccentricities detectable by LISA by integrating Equation 2.1 and Equation 2.2
until the orbital frequency reaches ?orb = 10?3 Hz at which point the gravitational
49
wave frequency is in LISA?s most sensitive range and is above the expected white
dwarf background. Figure 2.9 shows the distribution of eccentricities for binaries
with gravitational radiation in the LISA band. There are more low eccentricities
at higher mass ratios. This is because at low mass ratios each encounter takes
a fractionally larger amount of energy away from the binary than at high mass
ratios. Thus at low mass ratios, the last encounter will tend to harden the binary
such that it is closer to merger. At high mass ratios, however, encounters take a
smaller fractional amount of energy from the binary, and, thus, the high mass ratio
binaries have more time to circularize during their orbital decay. For the 1000:10:10
mass ratio, a large fraction of the eccentricities are in the range 0.1 <? e <? 0.2
where the binary is eccentric enough to display general relativistic effects such as
pericenter precession, but circular templates may be sufficient for initial detection
of the gravitational wave.
2.6 Conclusions
We present results of numerical simulations of sequences of binary-single black hole
scattering events in a dense stellar environment. We simulate three-body encounters
until the binary will merge due to gravitational radiation before the next encounter.
In half of our simulations, we include the effect of gravitational radiation between
encounters.
1. Sequences of high mass ratio encounters. Our simulations cover a range
of mass ratios including those corresponding to IMBHs interacting with stellar-
mass black holes in stellar clusters. Because the binaries simulated are tightly
bound, the encounters steadily shrink the binary?s semimajor axis until it merges.
The eccentricity, however, jumps chaotically between high and low values over the
50
Figure 2.9: Distribution of eccentricities after integrating the Peters (1964) equa-
tions until in LISA band when the orbital frequency ?orb = 10?3 Hz. The solid
histograms are the Newtonian only sequences, and the hatched histograms are
sequences with gravitational radiation. The sequences with gravitational radia-
tion tend towards lower eccentricity since they have already started to circularize
during the sequence. There is more difference between the two cases at higher
mass ratios since gravitational radiation is stronger. Higher mass ratio binaries
have lower eccentricities than lower mass ratio binaries since the latter start closer
to merger after the final encounter. The 1000:10:10 mass ratio shows that a large
number of detectable binaries would have 0.1 <?e <?0.2 such that they would likely
be detectable by LISA with circular templates yet display measurable pericenter
precession.
51
course of a sequence. Merger usually occurs at high eccentricity since gravitational
radiation is much stronger then.
2. Gravitational wave emission between encounters. The inclusion of gravita-
tional radiation between encounters affects the simulations in several ways. The
extra source of shrinking caused by gravitational wave emission has the effect of
shortening the sequence in terms of both the number of encounters and the total
time, and the circularization from gravitational waves has the effect of decreasing
the final eccentricity of the binary before it merges.
3. IMBH formation. Our simulations directly test the IMBH formation model
of Miller & Hamilton (2002b). We find that there is sufficient time to build up to
1000 Mcircledot when starting from 50 Mcircledot, but our simulations also show that if there
are a thousand 10 Mcircledot black holes in the globular cluster, the seed black hole would
only be able to grow to 240 Mcircledot before exhausting half of the black holes in the
cluster. In addition, the probability of the binary?s remaining in the cluster during
a growth from 50 to 240 Mcircledot is small. In order to avoid ejection from the cluster
with a reasonable probability, either the black hole must have a larger mass at the
onset of dynamical encounters, the cluster?s escape velocity must be larger, or the
black hole must grow by some additional mechanisms such as by Kozai-resonance-
induced mergers, two-body captures, from smaller interlopers, or from interlopers
that are tidally disrupted. Thes results are modified by our inclusion of gravitational
radiation reaction in Chapter 3.
4. Gravitational wave detection. The mergers of binary black hole systems are
strong sources of detectable gravitational waves. We find that the merging binary
will typically start with very high eccentricity. By the time the binary is detectable
by the Advanced LIGO detector, it will have completely circularized, but when
detectable by LISA, it may have moderate eccentricity (0.1 <? e <? 0.2) such that it
52
will display general relativistic effects such as pericenter precession and still possibly
be detectable with circular templates. We find a high rate of mergers in the first few
hundred million years of a globular cluster?s life. This suggests that recently formed,
nearby super star clusters may be promising sources for gravitational radiation from
IMBH coalescence.
53
Chapter 3
Three-Body Dynamics with
Gravitational Wave Emission
3.1 Overview
In this chapter we present a study of the dynamics of black holes in a stellar clus-
ter using numerical simulations that include the effects of gravitational radiation.
We include gravitational radiation reaction by adding a drag force to the Newto-
nian gravitational calculation. Our treatment is similar to that of Lee (1993), but
we focus on individual encounters and sequences of encounters and the resulting
mergers instead of ensemble properties of the host cluster. Chapter 2 incorporated
gravitational radiation by integrating the Peters (1964) orbit-averaged equations for
orbital evolution of a binary that is emitting gravitational waves, but in this chapter
we include the energy loss from gravitational radiation for arbitrary motion of the
masses. Although the vast majority of three-body interactions do not differ greatly
from a purely Newtonian simulation, an important few involve close approaches in
which gravitational waves carry away a dynamically significant amount of energy
54
such that it may cause the black holes to merge quickly in the middle of the en-
counter. This is qualitatively different from the mergers in Chapter 2 which were
caused by gravitational waves emitted by isolated binaries between encounters, and
this new effect is important in considering detectable gravitational waves as well as
IMBH growth.
In ? 3.2 we describe our method of including gravitational waves as a drag force
as well as numerical tests of its accuracy. We present our simulations and major
results in ? 3.3 and discuss the implications for IMBH formation and gravitational
wave detection in ? 3.4.
3.2 Numerical Method
The numerical method we use here is much the same as is described in ? 2.2. In
order to study the dynamics of a massive binary in a dense stellar environment, we
simulate the encounters between the binary and single objects. We include both
individual encounters and sequences of encounters, all of which include gravita-
tional radiation emission. When simulating sequences, we allow the properties of
the binary to evolve from interactions with singles, and we follow the binary until a
merger occurs. A merger is determined to occur when the separation between the
two masses is less than G(m0 + m1)/c2. The simulations are run using the same
code as in Chapter 2 with a few modifications. The integration engine is now HN-
Drag, which is an extension of HNBody (K. Rauch & D. Hamilton, in preparation).
Both HNBody and HNDrag can include the first-order post-Newtonian corrections
responsible for pericenter precession based on the method of Newhall et al. (1983).
HNDrag also has the ability to include pluggable modules that can add extra forces
or perform separate calculations such as finding the minimum separation between
55
all pairs of objects. In this paper we ignore the second-order post-Newtonian terms,
which contribute higher-order corrections to the pericenter precession, but we in-
clude the effects of gravitational radiation on the dynamics of the particles through
the addition of a force that arises from the 2.5-order post-Newtonian equation of
motion for two point masses. We discuss the implications of ignoring lower-order
post-Newtonian corrections in ? 3.3.1. The acceleration on a mass m0 from gravi-
tational waves emitted in orbit around a mass m1 can be written as
dv0
dt =
4G2
5c5
m0m1
r3
parenleftbigg m
1
m0 +m1
parenrightbiggbracketleftBigg
?r(?r?v)
parenleftBigg34
3
G(m0 +m1)
r + 6v
2
parenrightBigg
+v
parenleftBigg
?6G(m0 +m1)r ?2v2
parenrightBiggbracketrightBigg
(3.1)
where r = r1 ? r0 and v = v1 ? v0 are the relative position and velocity vec-
tors between the two masses and ?r and ?v are their unit vectors (Iyer & Will 1993;
Lee 1993). When orbit-averaged, Equation 3.1 gives the Peters (1964) equations
for semimajor axis and eccentricity evolution (Equation 2.1 and Equation 2.2. See
Appendix B for a derivation of Equation 2.1 for the circular case.). We tested the
inclusion of this force in the integrator by comparison with direct, numerical in-
tegration of Equation 2.1 and Equation 2.2 for two different binaries with masses
m0 = m1 = 10 Mcircledot and initial semimajor axis a0 = 1 AU: one with initial eccen-
tricity e0 = 0 and one with initial eccentricity e0 = 0.9 (Figure 3.1). The N-body
integration of these binaries made use of HNDrag?s enhancement factor, which ar-
tificially augments the magnitude of the drag forces for the purposes of testing or
simulating long-term effects. For this test and all numerical integrations with HN-
Drag, we used the fourth-order Runge-Kutta integrator. For both the circular and
the high eccentricity cases, the N-body integrations agree very well with the Peters
(1964) equations. Examination of Equation 3.1 reveals that even though physically
the emission of gravitational radiation can only remove energy from the system, the
56
equation implies ?E > 0 for ?r? ?v > 0 in hyperbolic orbits, becoming worse as the
eccentricity increases (Lee 1993). Integration of Equation 3.1 over an entire orbit,
however, does lead to the expected energy loss. This is because there is an excess of
energy loss at pericenter, which cancels the energy added to the system (Lee 1993).
Thus this formulation does not introduce significant error as long as the integration
is calculated accurately at pericenter, which we achieve by setting HNDrag?s relative
accuracy parameter to 10?13, and the two objects are relatively isolated, which we
discuss below.
We also tested the N-body integration with gravitational radiation for unbound
orbits against the maximum periastron separation for two objects in an initially
unbound orbit to become bound to each other (Quinlan & Shapiro 1989):
rp,max =
?
?85pi
?2G7/2m
0m1 (m0 +m1)
3/2
12c5v2?
?
?
2/7
, (3.2)
where v? is the relative velocity at infinity of the two masses. In Figure 3.2 we plot
the orbits integrated both with and without gravitational radiation for two different
sets of initial conditions that straddle the rp,max threshold. For both sets of initial
conditions, the integrations with gravitational radiation differ from the Newtonian
orbits, and the inner orbit loses enough energy to become bound and ultimately
merge. We used a bisection method of multiple integrations to calculate rp,max, and
our value agrees with that of Quinlan & Shapiro (1989) to a fractional accuracy of
better than 10?5.
For systems of three or more masses, we compute gravitational radiation forces
for each pair of objects and add them linearly. Although this method differs from
the full relativistic treatment, which is nonlinear, the force from the closest pair
almost always dominates. We may estimate the probability of a third object coming
within the same distance by examining the timescales for an example system. A
binary black hole system with m0 = 1000 Mcircledot and m1 = 10 Mcircledot with a separation
57
Figure 3.1: Comparison of HNDrag integration with numerical integration of Pe-
ters (1964) equations for an eccentric binary. Lines are numerical integration of
Equation 2.1 for semimajor axis (solid line) and of Equation 2.2 for eccentricity
(dashed line). The symbols are results from HNDrag integration with gravita-
tional radiation for semimajor axis (diamonds) and eccentricity (squares). The
binary shown has m0 = m1 = 10 Mcircledot with an initial orbit of a0 = 1 AU and
e0 = 0.9. The evolution of the binary?s orbital elements is in very close agreement
for the entire life of the binary.
a = 10?2 AU (? 1000GM/c2) will merge within (Peters 1964)
?merge ? 6?109 (1 Mcircledot)
3
m0m1 (m0 +m1)
parenleftbigg a
10?2 AU
parenrightbigg4parenleftBig
1?e2
parenrightBig7/2
yr ? 600 yr. (3.3)
Though the expression for merger time in Equation 3.3 is valid only for high eccen-
tricities (e ? 1), we calculate the time with e = 0, which serves as an upper limit.
Including a more realistic eccentricity would only decrease the time and help this
argument further. The rate of gravitationally focused encounters with a third mass
58
Figure 3.2: HNDrag-integrated orbits with and without gravitational radiation
inside and outside of two-body capture pericenter. This plot shows orbits of two
10 Mcircledot black holes with relative velocity of 10 km s?1 and pericenter distances of
rp = 0.2rp,max and rp = 1.2rp,max. The lines show the orbits with gravitational
radiation included in the integration, and the diamonds show the Newtonian
orbits for the same initial conditions. The direction of the orbits is indicated
by the arrow. Although it is not apparent for the outer orbit in this plot, both
trajectories differ from their Newtonian counterparts. For the inner orbit, enough
energy is radiated away for the black holes to become bound to each other and
eventually merge.
m2 within a distance r from an isotropic distribution is (Chapter 2)
?enc = 5?10?10
parenleftBigg10 km s?1
v?
parenrightBiggparenleftBigg n
106 pc?3
parenrightBiggparenleftbigg r
10?2 AU
parenrightbiggparenleftBiggm
0 +m1
1 Mcircledot
parenrightBiggparenleftBigg m
2
1 Mcircledot
parenrightBigg1/2
yr?1.
(3.4)
For a number density n = 106 pc?3, a relative velocity v? = 10 km s?1, and
an interloper mass m2 = 10 Mcircledot, the rate of encounters within the same distance
r = a = 10?2 AU is?enc ? 2?10?6 yr?1. Thus the probability of an encounter within
59
the same distance is P ? ?merge?enc ? 10?3 for this mildly relativistic case. For a
separation of 10?3 AU, the probability drops to 10?8. Thus for most astrophysical
scenarios and for all simulations in this chapter, the error incurred from adding the
gravitational radiation force terms linearly is negligible.
3.3 Simulations and Results
3.3.1 Individual Binary-Single Encounters
Close Approach
We begin our study of three-body encounters including gravitational radiation by
calculating the minimum distance between any two objects during a binary-single
scattering event. This quantity has been well studied for the Newtonian case, but
it is still not completely understood (Hut & Inagaki 1985; Sigurdsson & Phinney
1993). FollowingtheinitialconditionchoicesofHut&Inagaki(1985)andSigurdsson
& Phinney (1993), we present 105 simulations of a circular binary with masses
m0 = m1 = 1 Mcircledot and an initial semimajor axis a0 = 1 AU interacting with an
interloper of mass m2 = 1 Mcircledot in a hyperbolic orbit with respect to the center of
mass of the binary. As in Chapter 2, we refer to the mass ratios of three-body
encounters as m0:m1:m2, where m2 is the interloper and the binary consists of m0
and m1 with mbin = m0 +m1 and m0 ?m1. The relative velocity of the binary and
the interloper at infinity is v? = 0.5 km s?1 with an impact parameter randomly
drawn from a distribution such that the probability of an impact parameter between
b and b+db is P(b) ?b and a maximum value bmax = 6.621 AU, which corresponds
to a two-body pericenter distance of rp = 5a0. The encounters are integrated until
finished as determined in Chapter 2 while tracking the minimum distances between
60
all pairs of objects. We follow Hut & Inagaki (1985) and Sigurdsson & Phinney
(1993) in calculating a cumulative, normalized cross-section for close approach less
than r
?(r) = f (r)b
2
max
a20
parenleftbiggv
?
vc
parenrightbigg2
, (3.5)
where
vc ?
radicaltpradicalvertex
radicalvertexradicalbtGm0m1
am2
(m0 +m1 +m2)
(m0 +m1) (3.6)
is the minimum relative velocity required to ionize the system andf(r) is the fraction
of encounters that contain a close approach less than r. We plot ?(r/a0) for the
Newtonian case at several different time intervals within the encounter in Figure 3.3.
Our results for the total cross-section are in almost exact agreement with Sigurdsson
& Phinney (1993) over the domain of overlap, but with the advantage of ten years
of computing advances, we were able to probe down to values of r/a0 that are 102
times smaller. In addition we examine how the total cross-section evolves from
the initial close approach of the binary until the end of the interaction through
subsequent near passes during long-lived resonant encounters. At the time of the
interloper?s initial close approach with the binary, the cross-section is dominated by
gravitational focusing, and thus the bottom two curves in Figure 3.3 are well fit by
power laws with slope of 1. As the interactions continue, resonant encounters with
multiple close approaches are possible, and the cross-section for small values of r/a0
increases. Each successive, intermediate curve approaches the final cross-section by
a smaller amount because there are fewer encounters that last into the next time
bin. A fit of two contiguous power laws to the final curve yields a break at r/a0 =
0.0102 with slopes of 0.85 and 0.35 for the lower and upper portions, respectively.
These values are very close to those obtained by earlier studies (Hut & Inagaki
1985; Sigurdsson & Phinney 1993). There is, however, no reason for a preferred
scale for a Newtonian system, and simple models that assume close approaches
61
are dominated by pericenter passage after an eccentricity kick cannot explain the
lower slope. We numerically calculate d(log?)/d(logr) by fitting multiple lines
to ?(r) in logarithmic space and plot the results in Figure 3.4. The derivative
d(log?)/d(logr) appears to approach unity for very small values of r/a0 where the
close approach can be thought of as a gravitationally focused two-body encounter
within the entire system (Hut & Inagaki 1985). It is surprising that this does not
happen until r/a0 < 10?5.
In order to test the effects of gravitational radiation on the close approach as well
as to test the sensitivity of the results to the phase of the binary, we ran the same
simulations (1) with gravitational radiation, (2) with gravitational radiation and
first-order post-Newtonian corrections, and (3) with just first-order post-Newtonian
corrections. The three new cross-sections are plotted with the Newtonian results
in Figure 3.5. A Kolmogorov-Smirnov test shows the differences between the three
curves to be statistically insignificant (P ? 0.4). Although not statistically signif-
icant, the curves with gravitational radiation appear to drop below the Newtonian
curve for small r/a0 and then climb above for very small r/a0. This is as expected
because gravitational radiation causes this effect by driving objects that become
very close to each other closer still and, in some cases, causing them to merge. A
larger number of simulations could indicate whether the drop is indeed physical
or merely statistical fluctuation. For larger masses, the gravitational radiation is
stronger, and the gravitational radiation curve will differ from the Newtonian curve
at larger r/a0 for a fixed value of a0.
62
Figure 3.3: Cross-section for close approach during binary-single encounters as
a function of rmin/a0. The thick, upper curve is the cross-section for the entire
encounter. The remaining curves are the cross-section at intermediate, equally-
spaced times during the encounter starting from the bottom near the time of
initial close approach. Because we only include 20 intermediate curves, there
is a gap between the last intermediate curve and the final curve. The bottom
two curves have a slope close to 1, consistent with close approach dominated
by gravitational focusing. As the encounters progress, resonant encounters with
multiple passes are more likely to have a close approach at smaller rmin/a0, and
the curves gradually evolve to the total cross-section for the entire encounter.
Merger Cross-Section
The most interesting new consequence from adding the effects of gravitational radi-
ation to the three-body problem is the possibility of a merger between two objects.
Though the two-body cross-section for merger can be calculated from Equation 3.2,
the dynamics of three-body systems increases this cross-section in a nontrivial man-
63
Figure 3.4: Derivative of close approach cross-section curve for the entire en-
counter. Each symbol is the slope for a line segment fit to the top curve from
Figure 3.3 plotted as a function of the midpoint of the range. Because of the small
number of encounters that result in very small close approaches, the multiple line
segments used in the fits cover different ranges in log(r/a0). They were selected
so that each of the 100 line segments covers an additional 1000 encounters that
make up the cumulative cross-section curve. The scatter in the points is indica-
tive of the statistical uncertainty. For smaller close approaches, d(log?)/d(logr)
appears to approach unity. The rise at the right occurs because the cross-section
is formally infinite at rmin/a0 = 1.
ner. We present simulations of individual binary-single encounters for a variety of
masses. As in Chapter 2, the interactions were set up in hyperbolic encounters
with a relative velocity at infinity of v? = 10 km s?1 with an impact parameter
distribution such that the probability of an impact parameter between b and b+db
is P(b) ? b with bmin = 0 and bmax such that the maximum pericenter separation
would be rp = 5a0. The binaries were initially circular with semimajor axes ranging
64
Figure 3.5: Cross-section for close approach like Figure 3.3 including different
orders of Post-Newtonian corrections. The curves are purely Newtonian (solid),
Newtonian plus 2.5-order PN (dotted), Newtonian plus 1-order PN (dashed), and
Newtonian plus 1-order and 2.5-order PN (dash-dotted). The purely Newtonian
and the Newtonian plus 2.5-order PN curves come from 105 encounters each.
The other two curves come from 104 encounters each and show more statistical
fluctuations. The differences among the curves are not statistically significant.
from 10?6 AU to 102 AU, depending on the mass. The masses were picked such that
one of the three mass ratios was unity with all masses ranging from 10Mcircledot to 103Mcircledot
with roughly half-logarithmic steps. For each mass and semimajor axis combina-
tion, we run 104 encounters. We calculate the merger cross-section as?m = fmpib2max
where fm is the fraction of encounters that resulted in a merger while all three ob-
jects were interacting. In Figure 3.6 through Figure 3.9 we plot, as a function of the
semimajor axis scaled to the gravitational radius of the binary ? ? a/(Gmbin/c2),
the cross-section normalized to the physical cross-section of the Schwarzschild radius
65
of the mass of the entire system taking gravitational focusing into account:
??m = ?m
bracketleftBigg
4piGMtotv2
?
GMtot
c2
bracketrightBigg?1
. (3.7)
For all mass ratios ??m increases with ? because hard binaries with wide separations
sweep out larger targets where the interloper can interact with and merge with the
binary components. As ? increases to the point that the binary is no longer hard,
??m will approach the value expected from Equation 3.2. The curves flatten out for
? <? 100 as the cross-section is dominated by the mergers of binary members with
each other because of hardening interactions and eccentricity kicks that bring the
two masses together. For sufficiently small ?, the merger cross-section would be
formally infinite since all binaries would merge quickly. Thus the interesting regime
is where ? is neither too small nor too large; otherwise, the problem is effectively a
two-body problem where one may ignore the interloper (? ? 1) or where one may
ignore one of the binary members (? > 106). For all mass series, as the mass ratios
approach unity, the cross-section increases because complicated resonant encounters,
which produce more numerous and smaller close approaches, are more likely when
all three objects are equally important dynamically.
We note some interesting trends that can be seen in the plots. Note that for
the scalings given, it is only the mass ratios that matter and not the absolute mass
so that the 10:10:10 and 1000:1000:1000 cases only differ because of statistical fluc-
tuations (Figure 3.6 and Figure 3.7). Thus our results can be scaled to others,
e.g., 1000:100:100 would be the same as 100:10:10. For the 10:10:X mass series
(Figure 3.6), the normalized cross-section decreases with increasing interloper mass,
roughly as ??m ? (m2/mbin)?1. This happens because as the interloper dominates
the total mass of the system, complicated resonant interactions with more chances
for close approach are less likely. Thus for the 10:10:1000 case, there are far fewer
chances for a close approach that results in merger. The 1000:1000:X series (Fig-
66
Figure 3.6: Normalized merger cross-sections (Equation 3.7) for individual binary-
single encounters as a function of ? for the 10:10:X mass series. The normalization
is explained in the text. Each symbol represents 104 binary-single encounters.
Error bars given for the top curve are representative for all merger cross-section
curves in Figure 3.6 through Figure 3.9.
ure 3.7) shows a distinct break around ? ? 100. Since the binary mass is the same
for all curves, they all approach the same value for ? <? 100 where the binary mem-
bers merge with each other because of their small separation. For ? >? 100, the
higher mass interlopers are dynamically more important and cause more mergers.
The X:10:10 series curves (Figure 3.8) all approach the 10:10:10 curve for ? >? 105
where the dominant object in the binary has less influence over its companion.
67
Figure 3.7: Normalized merger cross-sections like Figure 3.6 for 1000:1000:X se-
ries. The error bars from Figure 3.6 are representative for the curves in this
figure.
3.3.2 Sequences of Encounters
Because a tight binary in a dense stellar environment will suffer repeated encoun-
ters until it merges from gravitational radiation, we simulate a binary undergoing
repeated interactions through sequences of encounters including gravitational radia-
tion reaction. As in Chapter 2, we start with a circular binary with initial semimajor
axis a0 = 10 AU and a primary of mass m0 = 10, 20, 30, 50, 100, 200, 300, 500,
or 1000 Mcircledot and a secondary of mass m1 = 10 Mcircledot. We simulate encounters with
interloping black holes with mass m2 = 10 Mcircledot. After each encounter, we integrate
Equation 2.1 Equation 2.2 to get the initial semimajor axis and eccentricity for the
68
Figure 3.8: Normalized merger cross-section like Figure 3.6 for X:10:10 series.
The error bars from Figure 3.6 are representative for the curves in this figure.
next encounter. This procedure continues until the binary merges from gravitational
radiation or there is a merger during the encounter. Throughout our simulations
we use an encounter speed of v? = 10 km s?1, an isotropic impact parameter such
that the hyperbolic pericenter would range from rp = 0 to 5a0, and a black hole
number density in the core n = 105 pc?3 (See Chapter 2 for an explanation of
these choices.). For each mass ratio we simulate 1000 sequences of encounters with
gravitational radiation reaction.
Our results are summarized in Table 3.1. The inclusion of gravitational waves
during the encounter makes a significant difference from the results reported in
Chapter 2. The fraction of sequences that result in a merger during an encounterfm
is a good indicator of the importance of gravitational waves. Even form0 = 10Mcircledot, a
69
Figure 3.9: Normalized merger cross-section like Figure 3.6 for 1000:X:1000 series.
The error bars from Figure 3.6 are representative for the curves in this figure.
significant fraction (fm > 0.1) of the sequences merge this way, and form0 > 300Mcircledot
this type of merger is more likely to occur than mergers between encounters, and
thus this effect shortens the sequence significantly. In particular, for m0 = 1000 Mcircledot
compared to the values from Chapter 2, the average number of encounters per
sequence ?nenc? is decreased by 42%; the average number of black holes ejected from
the cluster ?nej? is reduced by 56%; and the average sequence length ?tseq? is 67%
shorter. One caveat for the study of sequences of encounters is that an IMBH in
a cluster of much lower mass objects will gather a large number of companions in
elongated orbits through binary disruptions, and thus the picture of an isolated
binary encountering individual black holes may not hold when the IMBH becomes
very massive (Pfahl 2005a).
70
Table 3.1. Sequence Statistics.
m0/Mcircledot ?nenc? ?nej? fbinej ?tseq?/106 yr ?af?/AU ?ef? fm
10 46.4 3.2 0.652 54.10 0.174 0.904 0.134
20 46.7 5.1 0.515 40.86 0.224 0.900 0.130
30 52.4 7.3 0.457 29.47 0.290 0.898 0.156
50 62.3 10.8 0.329 19.17 0.291 0.897 0.190
100 83.9 16.6 0.103 11.65 0.401 0.893 0.275
200 123.0 24.3 0.011 7.26 0.411 0.885 0.387
300 147.8 26.9 0.002 4.74 0.543 0.881 0.492
500 197.5 33.1 - 3.03 0.611 0.879 0.627
1000 284.2 38.8 - 1.47 0.878 0.914 0.754
Note. ? Main results of simulations of sequences of encounters with gravita-
tional radiation included during the encounter. The columns are: the mass of the
dominant black hole m0, the average number of encounters per sequence ?nenc?,
the average number per sequence of stellar-mass black holes ejected from a stellar
cluster with escape velocity 50 km s?1 ?nej?, the fraction of sequences in which the
binary is ejected fbinej from a stellar cluster with escape velocity 50 km s?1, the
average time per sequence ?tseq?, the average final semimajor axis of the binaries
after the last encounter ?af?, the average final eccentricity of the binaries after the
last encounter ?ef?, and the fraction of sequences that end with a merger during
the encounter fm. Note that ?af? and ?ef? only refer to the binaries that do not
merge during the encounter; these comprise 1?fm of the sequences.
3.4 Discussion
3.4.1 Implications for IMBH Formation and Growth
Our simulations provide a useful look into the merger history of an IMBH or its
progenitor in a dense stellar cluster. As an IMBH grows through mergers with
stellar-mass black holes, it will progress through the different masses that we in-
cluded in our simulations of sequences. We interpolate the results in Table 3.1 to
calculate the time it takes to reach 1000 Mcircledot, the number of cluster black holes
ejected while building up to 1000 Mcircledot, and the probability of retaining the IMBH
progenitor in the cluster for different seed masses and escape velocities of the cluster.
71
The time to build up to 1000Mcircledot is dominated by?tseq?at high masses. Although
each individual sequence is short, far more mergers are required for the same frac-
tional growth in mass. In Figure 3.10 we plot the mass of the IMBH as a function
of time for an initial mass of m0 = 10, 50, and 200 Mcircledot, for which total times to
reach 1000 Mcircledot are 600, 400, and 250 Myr, respectively. Because we assume a con-
stant core density throughout the simulations, the times are unaffected by changing
the cluster?s escape velocity. Without gravitational radiation, the times are roughly
twice as long (Chapter 2) because the length of each sequence is dominated by the
time it spends between encounters at small a when encounters are rarer. With grav-
itational radiation included, mergers that occur during an encounter are more likely
at small separations, and the length of the sequence is shortened. These times are
much shorter than the age of the globular cluster and are smaller than or compara-
ble to timescales for ejection of black holes from the cluster, which we discuss below
(see also O?Leary et al. 2006; Portegies Zwart & McMillan 2000). Thus time is not
a limiting factor in reaching 1000 Mcircledot for an IMBH progenitor that can remain in a
dense cluster with a sufficiently large population of stellar mass black holes.
Each time that an encounter tightens the binary, energy is transfered to the
interloper, which leaves with a higher velocity. If energetic enough, this interaction
will kick the interloper out of the cluster. If the interactions kick all of the interacting
black holes out of the cluster, the IMBH cannot continue to grow. In a massive,
and therefore dense, cluster, there are roughly 103 black holes (Chapter 2). With
gravitational radiation included during the encounters, the number of black holes
ejected is roughly halved compared with the number when gravitational radiation
in included only between encounters (Figure 3.11), but the total number ejected
while building up to 1000 Mcircledot is still a few times the number of black holes available
even for an escape velocity of vesc = 70 km s?1. Thus a black hole smaller than
72
Figure 3.10: Mass of progenitor IMBH as a function of time as it grows through
mergers with 10 Mcircledot black holes in a dense stellar cluster. Solid curves show
results from this work in which gravitational radiation is included, and dashed
curves show results from Chapter 2 in which this effect is only included between
encounters. From bottom to top the curves show the growth of black holes with
initial mass m0 = 10, 50, and 200 Mcircledot. The IMBH progenitors all reach 1000 Mcircledot
in less than 600 Myr, and the inclusion of gravitational radiation significantly
speeds up the growth of the black hole.
m0 <? 600 Mcircledot cannot reach 1000 Mcircledot by this method without additional processes
such as Kozai resonances (G?ultekin et al. 2004; Miller & Hamilton 2002a; Wen 2003).
However, O?Leary et al. 2006 find that Kozai-resonance induced mergers will only
increase the total number of mergers by ? 10%. There is still the potential for
significant growth in a short period of time. If we consider the point at which half
of the black holes have been ejected from the cluster as the end of growth, then a
black hole with initial mass of 50 Mcircledot will grow to 290 Mcircledot in 120 Myr, and a black
73
hole of 200 Mcircledot will grow to 390 Mcircledot in less than 100 Myr (compared to 240 Mcircledot
and 330 Mcircledot without gravitational radiation during the encounter). In addition,
this ejection of stellar-mass black holes by a binary with a large black hole is faster
than by self-ejection from interactions among stellar-mass black holes calculated by
Portegies Zwart & McMillan (2000), who find that ? 90% of black holes are ejected
in a few Gyr. O?Leary et al. (2006), however, find that the inclusion of a mass
spectrum of black holes further speeds up the ejection of stellar-mass black holes by
a small amount.
For every kick imparted on an interloper, conservation of momentum ensures a
kick on the binary. Even with a large black hole, extremely large kicks can eject the
binary from the cluster, at which point the IMBH progenitor can no longer grow.
We can calculate the probability of IMBH retention for an individual sequence as
1 ? fbinej, from which we interpolate the probability of remaining in the cluster
while growing to 300 Mcircledot when the binary is essentially guaranteed to remain in
the cluster. We plot this probability as a function of seed mass for several different
escape velocities in Figure 3.12. The inclusion of gravitational waves during the
encounter increases the retention probability for small masses. For m0 = 50 Mcircledot
the cluster retains the binary more than 12% of the time, and 49% of the time for
m0 = 100 Mcircledot. Because the energy that an interloper can carry away from the
system scales as
?E ? m1m
0 +m1
|EB| = m1m
0 +m1
Gm0m1
2a , (3.8)
the encounters at the end of the sequence, when a is smallest, are the most likely
to impart a kick large enough to eject the binary from the cluster. This is also
the point at which effects from gravitational radiation are strongest and at which
close encounters are most likely to cause a merger. When the encounter ends in a
merger, there can be no more ejections. The mergers from gravitational radiation
74
Figure 3.11: Number of black holes ejected in building up to 1000 Mcircledot (solid
curves) and to 500 Mcircledot (dashed curves) as a function of seed mass for different
cluster core escape velocities, given in units of km s?1. The dotted line indicates
the expected number of black holes in a large globular cluster. The dot-dashed
curve from Chapter 2 shows the number of black holes ejected from the cluster
in building up to 1000 Mcircledot for a cluster escape velocity of 50 km s?1 without the
effects of gravitational waves during the encounter. For all but the largest seed
masses, the globular cluster does not contain enough black holes for the IMBH to
reach 1000 Mcircledot. There are, however, a sufficient number of black holes to build up
to 500 Mcircledot for a seed mass greater than 225 Mcircledot or an escape velocity of at least
60 km s?1. The inclusion of gravitational radiation during the encounter roughly
halves the number of ejections.
75
Figure 3.12: Probability for a binary with an IMBH to remain in the cluster until
building up to 300 Mcircledot as a function of seed mass for different cluster core escape
velocities given in units of km s?1. Solid curves are results from this work, and the
dashed curve is from Chapter 2 for an escape velocity of 50 km s?1. The inclusion
of gravitational radiation significantly increases the retention probability.
decrease the number of ejections by decreasing the number of encounters and thus
the number of possible ejections as well as cutting off what would otherwise be the
end of the sequence, in which ejections are more likely to occur.
As mentioned in ? 2.2, our analysis of the ejection of stellar-mass black holes as
well as of IMBH progenitors does not include the effects of gravitational radiation
recoil. Based on the four most recent calculations (Baker et al. 2006; Blanchet et al.
2005; Damour & Gopakumar 2006; Favata et al. 2004), a seed mass ofm0 = 150Mcircledot
merging withm1 = 10Mcircledot companions will produce a recoil velocityvr <?50 km s?1,
and if the numerical relativity simulations are the most accurate predictions, a seed
76
mass of m0 = 100 Mcircledot is sufficient. Therefore, a seed mass greater than 100 Mcircledotis
sure to avoid ejection from both dynamical interactions and gravitational radiation.
3.4.2 Implications for Gravitational Wave Detection
In addition to the likelihoods and rates of growth of black holes in dense stellar
systems, our simulations shed light on the gravitational wave signals that come
from the mergers of these black holes. Making optimistic assumptions, O?Leary
et al. (2006) calculate upper limits for Advanced LIGO detection rates for black hole
mergers in stellar clusters that formed at a redshift z = 7.8. For their wide range of
cluster properties, they find detection rates ranging from ?AdLIGO ? 0.6 to 10 yr?1.
For cluster parameters that most closely resemble those used in Chapter 2 and in
this work (GMH model series), they find ?AdLIGO ? 2 to 4 yr?1. Our simulations
show that when gravitational radiation is included in the integration the number of
black holes ejected per merger decreases for all mass ratios. With fewer black holes
ejected from the cluster, the overall rate of black hole mergers increases. For the
10:10:10 case, the number of ejections per merger decreases by ? 10%, and for the
1000:10:10 case the number decreases by more than a factor of 2, thus increasing
the rates found by O?Leary et al. (2006). The exact increase in rate is difficult
to estimate because the total number of mergers is dominated by mergers between
stellar-mass black holes, yet the most easily detected mergers involve black holes
with larger masses.
Because dynamical interactions strongly affect the eccentricity of a binary and
because the timescale for merger is a such a strong function of eccentricity, binaries
in a cluster tend to have very high eccentricities after their last encounter (Chap-
ter 2, O?Leary et al. 2006). With the addition of gravitational radiation during the
encounter, we find that the merging binaries become more eccentric because a sig-
77
nificant fraction of the mergers (fm in Table 3.1) occur during the encounter. These
mergers typically happen between two black holes that are not bound to each other
until they come close to each other and emit a significant amount of gravitational
radiation, after which the two black holes are in an extremely high eccentricity orbit
(1?e <? 10?3).
To see how these high eccentricities affect the detectability of the gravitational
wave signal, we integrate Equation 2.1 and Equation 2.2 until the binaries are de-
tectable by LISA and then Advanced LIGO. For circular orbits, the frequency of
gravitational wave emission is twice the orbital frequency, but masses in eccentric
orbits emit at all harmonics: fGW = n?/2pi, where n is the harmonic number and
? =
bracketleftBiggG(m
0 +m1)
a3
bracketrightBigg1/2
(3.9)
with peak harmonic for e > 0.5 at approximately n = 2.16(1?e)?3/2 (Farmer &
Phinney 2003). We consider the binary to be detectable by LISA when the peak
harmonic frequency is between 2 mHz and 10 mHz. We plot the distribution of ec-
centricities in Figure 3.13. The distributions are essentially a combination of those
from Figure 2.9 and a sharp peak near e = 1, which comes from the mergers during
the encounter. The number in the sharp peak increases with mass as fm increases
such that for 1000:10:10 more than 75% of the merging binaries detectable by LISA
have an eccentricity greater than 0.9. Between 15% and 25% of all of the merging
binaries have eccentricities so high that the peak harmonic frequency is above the
most sensitive region of the LISA band, but they should still be emitting strongly
enough at other harmonics to be detectable. Such high eccentricity presents chal-
lenges for the detection of these signals from the data of space-based gravitational
wave detectors because (1) it requires a more computationally expensive template
matching that includes non-circular binaries and (2) the binaries only emit a strong
amount of gravitational radiation during the short time near periapse as they merge.
78
For a given semimajor axis, these extremely high eccentricities will also increase the
gravitational wave flux emitted and thus increase the distance out to which LISA
can detect them, but the detection rate may be compensated by the fact that more
parameters are required (Will 2004). We also integrate the orbital elements of the
binaries until they are in the Advanced LIGO band (40 Hz < fGW < fISCO) or
within a factor of 2 of their ISCO frequency for m0 > 100 Mcircledot. We find that they
have almost completely circularized (Figure 3.14). A tiny fraction (< 0.5%) of the
runs with m0 = 500 and 1000 Mcircledot have merging binaries with eccentricities such
that 1?e <? 10?6.
3.5 Conclusions
1. Gravitational radiation in N-body. We present results of numerical simulations
of binary-single scattering events including the effects of gravitational radiation dur-
ing the encounter. We include gravitational radiation by adding the 2.5-order post-
Newtonian force term (Equation 3.1) to the equation of motion within the HNDrag
framework. The code reproduces the expected semimajor axis and eccentricity evo-
lution, and it gives the expected two-body capture radius.
2. Close approach and merger cross-sections. We use the new code to test the
effects of gravitational radiation on a standard numerical experiment of binary-single
encounters. We probe the close approach cross-section to smaller separations than
has been simulated previously and find that the inclusion of gravitational radiation
makes little difference except for extremely close encounters (rp < 10?5a), at which
point gravitational radiation drives the objects closer together. We also present the
cross-section for merger during binary-single scattering events for a variety of mass
ratios and semimajor axes.
79
Figure 3.13: Histogram of eccentricities of merging binary while in the LISA band
(fGW = 2 mHz to 10 mHz) out of a total of 1000 sequences. The histograms
show a combination of the binaries that merged after the last encounter with
eccentricities concentrated around 0 < e <? 0.3 and the black holes that merged
quickly during the encounter with eccentricities very close to unity. The peaks in
the rightmost bin in all plots lie above the range of the plots.
3. IMBH growth in dense stellar clusters. We simulate sequences of binary-single
black hole encounters to test for the effects of gravitational radiation and to test
formation and growth models for intermediate-mass black holes in stellar clusters.
We find that the inclusion of gravitational radiation speeds up the growth of black
holes by a factor of 2, increases the retention of IMBH progenitors by a factor of 2,
and decreases the ejection of stellar-mass black holes by a factor of 2. All of these
effects act to enhance the prospects for IMBH growth.
4. Detectability of gravitational waves. We analyzed the merging binaries from
80
Figure 3.14: Histogram of eccentricities of merging binary while the gravitational
wave frequency is detectable from current and future ground-based detectors. The
upper limit of the frequency range is the ISCO frequency. We used a lower limit
for frequency range of 100 Hz for m0 = 10, and 20 Mcircledot; 35 Hz for m0 = 30, 50, and
100 Mcircledot; and half the ISCO frequency for the higher mass binaries. The binaries
are very close to circular once they are in the frequency range of ground-based
detectors. The peaks in the leftmost bin in all plots lie above the range of the
plots.
the simulations of black holes in dense stellar clusters to look at the detectability
of the gravitational wave signals from these sources. We find that the mergers that
occur rapidly during the encounter as opposed to those that occur after the final
encounter are an important source of black hole mergers, becoming the dominant
source of mergers at the higher mass ratios. The mergers that do occur during the
encounter tend to have extremely high eccentricity (e > 0.9) while in the LISA
band, presenting challenges for their detection. When the gravitational wave signal
81
from the merging black holes is in the Advanced LIGO band, the orbit will have
completely circularized.
82
Chapter 4
Summary and Conclusions
4.1 Answers
Let us return to the questions raised in ? 1.4, which may be simplistically para-
phrased as, ?Does the Miller & Hamilton (2002b) model of IMBH formation work;
can we see any unique gravitational waves from this model; and what else have
you learned?? These may be equally simplistically answered as, ?Yes, with slight
modifications; yes and no; and quite a bit.? More specifically, the first question was:
1. Is the Miller & Hamilton (2002b) model of IMBH formation fast enough to
make an IMBH before the globular cluster?s supply of stellar-mass black holes
eject themselves from the cluster through dynamical interactions (? 0.4 to
1 Gyr)?
The simulations show that the process of repeatedly merging stellar-mass black holes
with a seed mass can happen quickly enough. The simulations in this dissertation
can be seen as building on each other with an additional layer of gravitational
physics that the previous layer did not include. The first set of simulations in
Chapter 2 used Newtonian integrations with a relativistic endpoint. The second
83
set of simulations in Chapter 2 added orbital decay between Newtonian encounters
in a sequence with a relativistic endpoint, and the simulations in Chapter 3 added
gravitational wave emission to the encounters. Each additional layer decreased the
amount of time needed to reach 1000 Mcircledot: from 1.1 Gyr to 0.71 Gyr to 0.41 Gyr,
assuming a seed mass of 50Mcircledot. This is comparable to the 0.4?1.0 Gyr that it takes
for the stellar-mass black holes to eject themselves (O?Leary et al. 2006; Portegies
Zwart & McMillan 2000; Sigurdsson & Hernquist 1993). This process happens
more quickly than originally anticipated by Miller & Hamilton (2002b) because the
authors assumed that the binaries would merge with eccentricity e ? 0.7. Instead,
the binaries merge with a high eccentricity, which allows the binaries to merge with
a much larger semimajor axis and, therefore, in a much shorter time. The times
listed above assume that the density of stellar-mass black holes remains constant
throughout the process. In reality, the density will be decreasing as the stellar-
mass black holes eject themselves and as the IMBH progenitor ejects them from the
cluster core. This leads to the next question.
2. Can repeated mergers of stellar-mass black holes in a globular cluster produce
an IMBH without ejecting too many stellar-mass black holes via the encounters
that harden the IMBH-progenitor binary?
There are roughly 103 stellar-mass black holes in a typical large globular cluster, and
the process of binary-hardening and merging will eject more than this when growing
to 1000Mcircledot. As before, when each additional layer of gravitational physics is added,
the number decreases: from 6800 to 5300 to 2900. This assumes that the interactions
continue regardless of the number of black holes in the cluster. Another way of
looking at this is that an IMBH (or its progenitor) in a binary in a stellar cluster
is the dominant source of stellar-mass black hole ejections. Before all of the black
holes are ejected, however, all scenarios allow 15 to 25 mergers before exhausting
84
roughly half of the black holes, and thus all scenarios allow for significant growth.
These numbers also assume an escape velocity from the cluster of vesc = 50 km s?1,
which is a typical value. The escape velocity, however, has been found to be much
higher for some Galactic globular clusters. For escape velocities of vesc = 60 and
70 km s?1, the most realistic simulations find 2300 and 1900 black holes ejected,
respectively. The total number of black holes in a very massive globular cluster
could reach ? 2000, and such a cluster would have a higher-than-average escape
velocity. Thus it is entirely plausible that the Miller & Hamilton (2002b) model
could reach 1000 Mcircledot in such a cluster. Finally, as seen in Figure 3.11, an escape
velocity of vesc ? 60 km s?1 will eject fewer than 1000 black holes when building up
to 500Mcircledot. Such a black hole is certainly an IMBH and is massive enough to explain
all but the brightest ULXs. Even the brightest ULX, M82 X-1, can be explained if
its brightest fluxes come from short periods of accreting at ? 2 times the Eddington
rate.
3. Will a globular cluster retain an IMBH progenitor over the course of its merger
history, or will the binary-single scattering events eject the IMBH progenitor
with great regularity?
The answer to this question depends on the seed mass of the IMBH progenitor. The
classical Miller & Hamilton (2002b) model calls for a 50 Mcircledot black hole as the seed
mass. For the three layers of gravitational physics the probability of a seed mass?s
retention is 0.00766, 0.0356, and 0.124, respectively. The retention probability is the
biggest quantitative change among the three classes of simulations: from negligible
to respectable. If roughly 10?1 of all globular clusters can hold on to their IMBH
progenitor, then only 10?2 of the IMBHs produced need to be currently accreting
in order to explain the ? 10% of nearby galaxies with ULXs, assuming 100 globular
clusters per galaxy. So the classical model can withstand ejection by dynamical
85
kicks from three-body encounters. Withstanding ejection by gravitational radiation
recoil is a different matter. As discussed in ? 2.2, gravitational radiation carries
momentum with it, and the asymmetric radiation from unequal masses causes the
center of mass of the binary to spiral outwards from its original position as seen by
an observer far from the binary. Using the fitting formula of Fitchett (1983) and
the numerical relativity simulations of Baker et al. (2006), the recoil velocity is
vrecoil = 9?103q
2 (1?q)
(1 +q)5 km s
?1 (4.1)
forq ?m1/m0 ? 1, where the factor in front comes from settingvrecoil = 105 km s?1
for q = 0.67. For a 50 Mcircledot black hole that only merges with 10 Mcircledot black holes,
q = 0.2, and Equation 4.1 implies vrecoil ? 115 km s?1. Thus the classical Miller
& Hamilton (2002b) model cannot keep the seed mass in the cluster after the first
merger. If, however, the seed mass were m0 = 100 Mcircledot, then the recoil velocity
when merging with m1 = 10 Mcircledot black holes would be vrecoil ? 50 km s?1, or just
at the escape velocity of a typical globular cluster. It is reasonable, however, for
black holes to have mass m1 = 20 Mcircledot. To avoid ejection from mergers with such
black holes, the seed mass would need to be m0 > 200 Mcircledot. If we only consider
clusters with escape velocity vesc > 60 km s?1, then a seed mass with m0 > 170 Mcircledot
would avoid ejection. Such 100Mcircledot to 200Mcircledot black holes are already exotic objects
(IMBHs by our own definition!), and they require their own formation mechanism.
One possibility is if two 50 Mcircledot black holes exist in a single cluster. They would
likely find each other through interactions and end up in a binary. As equal-mass
objects, they would suffer no gravitational radiation recoil when they merged, and
then the resulting 100 Mcircledot black hole might grow to ? 170 Mcircledot by merging with
10 Mcircledot black holes. Similarly, a 200 Mcircledot black hole could be created by hierarchical
merging of four 50 Mcircledot black holes (or eight 25 Mcircledot black holes, etc.). While such
scenarios are possible, they require finesse to prevent the wrong black holes from
86
meeting up at the wrong time. In addition, slightly unequal masses can still lead
to recoil velocities greater than the escape velocity of the cluster, and misaligned
spins would contribute as well. The most obvious method of providing a massive
seed is from the core collapse phase of a stellar cluster. If the runaway collisions
of stars during this phase produce a merger remnant that loses 80% to 90% of its
mass in its evolution to a black hole, then it will produce a black hole with mass
200 to 300 Mcircledot. If the merger remnant always retains its mass, a 200 to 300 Mcircledot
seed can be formed if the collisions do not become a true runaway process. If the
core collapse phase only results in 1 or 2 collisions of ? 100 Mcircledot stars, the resulting
200 to 300 Mcircledot black hole could act as a seed in the host cluster.
It has been said that the Miller & Hamilton (2002b) model is an inefficient
process for making IMBHs (e.g., O?Leary et al. 2006), but the most conservative
conclusion one can draw from the simulations presented in this dissertation is that
a 200 Mcircledot seed mass will grow to 500 Mcircledot in less than 200 Myr, within the available
supply of black holes with a retention probability greater than 0.9. Put differently,
it is unclear how such a seed mass could avoid such growth with any regularity. The
model requires another mechanism to create the seed mass, but the core-collapse
collisions appear to be a promising source. Finally, the Miller & Hamilton (2002b)
model can ?save? the runaway collisions method of forming IMBHs if the merger
remnants cannot turn much of their mass into a black hole.
4. What is the observable gravitational-wave signature of these merging and in-
spiraling systems?
Though the expected detection rate of these events with LISA is low, detectable
inspirals are likely to be from highly eccentric binaries. Again, the inclusion of
gravitational radiation makes a difference for the expected eccentricity. From Fig-
ure 2.9, one can see the difference between the histograms for high masses, and in
87
Figure 3.13 one can see that gravitational radiation included during dynamical en-
counters completely changes the histogram. For LIGO events, the binaries will have
almost completely circularized, but the increased rate of mergers for simulations
that include gravitational radiation reaction forces indicates that detection rates
may be higher than previously estimated.
5. What is the reason for the apparent broken power-law in the cross-section for
closest approach in a three-body encounter?
The curves in Figure 3.3 show that as time progresses within a three-body encounter
the cross-section for closest approach evolves from that expected for a single pass
in a gravitationally focused encounter into the apparent broken power-law curve.
This suggests that the reason is that multiple passes by resonant encounters with
additional close approaches increase the number of encounters with small rmin/a0.
The slope of the curve approaches unity, as expected, but it does not happen until
very small values of rmin/a0. This is surprising because one would expect that as
two objects get within a distance much smaller than the semimajor axis of the
binary (the next smallest scale distance), one could approximate the encounter as
a two-body encounter. These selected two-body encounters, however, probably do
not come from the isotropic parameter space that would produce a curve with slope
unity.
6. How does the cross-section for close approach change when gravitational waves
are included?
The curves in Figure 3.5 show that gravitational radiation reaction forces do not
significantly alter the cross-section except for extremely close encounters that are
driven all the way to merger.
88
7. What is the likelihood of merger because of gravitational radiation during
a binary-single encounter when the leading-order terms of energy loss from
gravitational waves are included?
Figure 3.6 through Figure 3.9 show the shape of the merger cross-section curves.
There appears to be no simple scaling that takes into account the different regimes
of mass ratios possible.
4.2 Possible Future Work
Though it covers a lot, this dissertation does not represent all that could be done
with this technique to answer these questions. The method of simulating an IMBH
progenitor?s evolution through dynamics and mergers with isolated encounters may
be extended with additional considerations. One such consideration is to include
simulations of the acquisition of the IMBH progenitor?s binary companion. Although
Appendix A shows that this will probably not siginificantly increase the time to grow
an IMBH, it is possible that including the acquisition of the companion will decrease
the time slightly if the typical initial orbit of the binary is much tighter or extremely
eccentric. Including a mass spectrum instead of the single population of 10Mcircledot black
holes will also change the number of black hole ejections and retention probability.
If the population contains a large number of black holes smaller than 10 Mcircledot, the
number of stellar-mass black holes ejected increases while the probability of the
IMBH?s ejection decreases, and vice versa for a large number of black holes more
massive than 10 Mcircledot. It is unclear what will happen if the population contains large
numbers of both.
Binary-binary interactions are also important. Table 3.1 shows that ? 2?104
encounters are required to grow from 50 Mcircledot to 1000 Mcircledot. This means that even
89
for a binary fraction fbin ? 10?2, there will be ? 200 encounters with a binary. In
binary-binary encounters of stars, collisions between stars are very frequent (Fregeau
et al. 2004); so it is likely that mergers are important as well. Another intriguing
aspect of binary-binary interactions is the possibility of long-lived triple systems as
an outcome. Such triple systems may be subject to the Kozai (1962) mechanism
that causes the inner binary?s eccentricity to oscillate. Under the influence of gravi-
tational radiation, the binary would merge when the eccentricity reached a very high
value. The wide outer binary also increases the interaction cross-section, which may
result in resonant single-triple encounters that result in mergers.
Finally, it would be possible to model the changing number density of the black
hole population as it interacts with itself and an IMBH. This could be done with
loss-cone theory (e.g., Yu & Tremaine 2003), and it may allow a more accurate
calculation of the time and merger rate.
4.3 For Posterity
If one wishes to reduce this dissertation to a general, minimalist sound bite, it would
be: In order to understand the dynamics of compact objects in dense clusters or their
gravitational waves, their influence on each other must be studied because Newtonian
dynamics informs the generation of gravitational waves, and gravitational radiation
alters the dynamics.
Although the scientific motivation for this dissertation is sound science, it is
possible that there is some fundamental assumption that is incorrect. There are,
after all, dissertations concerning the ether (e.g., Gilbert 1901; Webster 1913). What
is the relevance of this dissertation if the assumptions are wrong? If ULXs or
cluster kinematics are shown not to be evidence of IMBHs, the process of repeated
90
mergers that increase the mass of black holes must still proceed to some extent
if there is a seed mass. Even if runaway-collisions and the initial mass function
of the cluster conspire to preclude massive black holes as seeds to this process,
there are other possible channels to a seed, for example, binary black holes or non-
runaway collisions such as blue stragglers. The overall efficiency of making big
black holes will be reduced without a large seed, though. If very solid evidence that
IMBHs definitely do not exist in clusters is found, then the question becomes what
prevents repeated mergers of black holes in clusters? The merger cross-sections
and close approach cross-sections should be relevant unless the solid underlying
physics of general relativity or Newtonian dynamics does not hold. Even then, the
mathematical aspects of the three-body problem or three-body plus drag problem
exist.
91
Appendix A
Time for IMBH Progenitor to
Acquire a Companion
Because the simulations of sequences of encounters do not include the interactions
that result in the initial acquisition of a companion, the time needed to do so is
not included in the analysis of IMBH growth. Here we estimate the time needed.
Because of the exchange bias that leads to more-massive objects in binaries and less-
massive objects ejected from the system, an IMBH progenitor is favored to pick up a
companion in a strong encounter with a binary. In order to reach the starting point
of the simulations in Chapter 2 and Chapter 3, the IMBH progenitor needs to acquire
a black hole companion. Analogous to the derivation of Equation 2.4, the average
time to interact with a black hole in a binary system is ??comp? = 1/?nbinv??(r)?,
where nbin is the number density of black holes with a companion of any mass. If
the black holes in binaries are distributed evenly throughout the entire population
of black holes, then nbin = nfbin, where fbin is the fraction of black holes with
a companion. In reality, black holes with companions will be more massive than
average and therefore have a smaller scale height and larger density, but we will
92
ignore the mass of the companion for now: mbin = mBH +mcomp ?mBH = 10 Mcircledot.
As in Equation 2.3, ?(r) is the gravitationally focussed cross section for interac-
tion
?(r) ?pir2 + 2pirGmIMBH/v2? ? 2pirGmIMBH/v2?, (A.1)
for gravitationally-focusing-dominated interactions for an IMBH progenitor of mass
mIMBH. For a strong interaction to result in an exchange, the distance must be
of order the resulting semimajor axis of the binary after the encounter af. If the
binding energy of the binary does not change with the encounter, the distance of
the encounter is r ? af ? a0mIMBH/mcomp for an initial semimajor axis a0. The
distribution of semimajor axes among black holes with companions in a dense stellar
cluster is difficult to determine, but the largest binaries would be near the hard/soft
boundary, which is
a0 = GmBHv2
?
? 100 AU, (A.2)
where we have assumed an average v? = 10 km s?1. If semimajor axes are dis-
tributed with equal number per decade as they are with main sequence binaries,
then a conservative average semimajor axis would be ?a0? = 3 AU. The average
time is then
??comp (mIMBH)? =
angbracketleftBigg v
?mcomp
2pinBHfbinGm2IMBHa0
angbracketrightBigg
= 0.04 Myrf
bin (mIMBH/50 Mcircledot)
2 (?a
0?/3 AU)
,
(A.3)
where we have used nBH = 105 pc?3, mcomp = 0.4 Mcircledot. The total time spent
acquiring binaries while building from 50 Mcircledot to 1000 Mcircledot is
T =
96summationdisplay
i=1
??comp (mIMBH,i)?f?1ex , (A.4)
wheremIMBH,i = 50Mcircledot+(i?1)mBH is the mass of the growing IMBH progenitor and
fex is the fraction of such encounters that results in the desired exchange. For the
mass ratios considered here, fex >? 0.05 (Heggie et al. 1996). Because of the strong
93
dependence on the mass of the IMBH progenitor, the total time is dominated at
small masses. Even for fbin = 0.05 and ?a0? = 3 AU, T ? 80 Myr, which is less
than 25% of the smallest total growth time when starting at 50 Mcircledot. In Figure 3.10,
a seed mass of 200 Mcircledot grows to 1000 Mcircledot in 250 Myr, but because Equation A.4
is dominated by the time when the progenitor mass is small, the total time spent
acquiring companions is less than 10 Myr.
If anything, these estimated times are over-estimates. For example, the IMBH
progenitor has a ? 2 ? 5% probability of acquiring the lower-mass object as a
companion, which may then be exchanged for a stellar-mass black hole during one
of the many hardening encounters that will take energy out of the binary. A similar
process may occur for a tight binary with two ? 5 Mcircledot stars, which would also have
sunk to the center of the stellar cluster with the other ? 10 Mcircledot objects on the
order of ? 106 yr. Therefore, ignoring the time to acquire a companion introduces
only a small systematic error to the total IMBH growth time calculated in ? 2.4 and
? 3.4.1.
For comparison, we may also calculate the time for the IMBH progenitor to
acquire a black hole companion by two-body capture via gravitational radiation.
This time is ?2?body = 1/?nv??(rp,max)?, where rp,max is the maximum periastron
separation for two objects to become bound to each other, given by Equation 3.2.
For m0 = mIMBH = 50 Mcircledot, m1 = mBH = 10 Mcircledot, and the values given above, the
maximum periastron separation is rp,max = 3.3?10?4 AU, and the average time is
?2?body ? 40 Gyr. Thus, three-body encounters with exchanges are the dominant
source of binary acquisition.
94
Appendix B
Demonstration that 2.5PN Drag
Force Reproduces the Peters
Equations for Circular Orbits
In this appendix we show that Equation 3.1 produces the Peters (1964) orbit-
averaged equations for semimajor axis evolution (Equation 2.1) for circular orbits.
Peters (1964) has the derivation for the general case, but the circular case is elegant.
In the barycentric frame, the velocity vectors of the two masses m0 and m1 are
always anti-aligned, and m0v0 +m1v1 = 0, so that the relative velocity between the
two masses is
v ? v1 ?v0 = ?(1 +m0/m1)v0. (B.1)
The semimajor axis of a binary may be expressed as a function of the relative
velocity between the two masses:
a = GMv2 , (B.2)
where M = m0 +m1 is the total mass of the binary. Differentiating Equation B.2
95
with respect to time yields
da
dt = ?
2GM
v3
dv
dt. (B.3)
Substituting Equation B.1 into Equation B.3, we get
da
dt =
2GM
v3
parenleftbiggM
m1
parenrightbiggdv
0
dt . (B.4)
Since ?r?v = 0 for a circular orbit, Equation 3.1 becomes
dv0
dt =
4G2
5c5
m0m1
a3
parenleftbiggm
1
M
parenrightbiggbracketleftbigg
v
parenleftbigg
?6GMa ?2v2
parenrightbiggbracketrightbigg
, (B.5)
where we have taken r = a. Substituting Equation B.5 into Equation B.4 gives
da
dt = ?
2GM
v3
bracketleftBigg4G2
5c2
parenleftbiggm
0m1
a3
parenrightbigg
v
parenleftbigg
6GMa + 2v2
parenrightbiggbracketrightBigg
. (B.6)
Using Equation B.2 and rewriting, this becomes
da
dt = ?
64
5
G3
c5
m0m1M
a3 , (B.7)
which is equivalent to Equation 2.1 for e = 0.
96
Appendix C
Summary of Code
The code used in this dissertation is called iabl, which stands for ?It?s a binary?s
life.? It was written to be as general purpose as possible though it is optimized
for the simulations run for this dissertation. The source for iabl consists of 12
files, which total over 3800 lines of C code. The code uses HNBody or HNDrag
(HNBody hereafter) for the actual integration of the orbits with its classical Runge-
Kutta integrator, but the rest of the tasks are done by iabl itself. The code was
developed to be modular and reusable in other contexts. Because the computation
is dominated by the actual N-body integration, the rest of the routines could be
designed for readability and simplicity as long as the integrations were efficient. The
code can be broken down into four primary sections: an initial conditions generator,
the interface between HNBody and iabl, examination of integration, and the two-
body approximator. We now describe these primary sections.
97
C.1 Initial Conditions Generator
The initial conditions generator uses drand48, a standard random number generator
to produce a Monte Carlo sampling of the multi-dimensional parameter space of the
three-body problem. The initial conditions generator is supplied with the masses
of the three objects, m0, m1, and m2, the initial semimajor axis of the binary,
a0, the initial eccentricity of the binary, e0, the relative speed of the encounter
between the center of mass of the binary and the interloper at infinity, v?. These
are either supplied as command-line arguments or derived from the final conditions
of the previous encounter in the sequence. The parameters generated are the true
anomaly of the binary, f, the impact parameter of the interloper with respect to
the center of mass of the binary as measured from infinity, b?, the azimuth of the
interloper?s initial velocity vector with respect to the pericenter of the binary, ?, the
zenith angle of the interloper?s initial velocity vector with respect to the plane of
the binary?s orbit, ?, an angle to determine in which direction the impact parameter
takes effect, ?.1 A limitation of the code is that the distribution of encounters is
assumed to be either isotropic in three dimensions or isotropic and coplanar. This
could in principle be changed to accommodate any distribution of encounters. The
above parameters are converted to Cartesian coordinates in the barycentric frame.
The above parameters are stored in a structure within the main function, and are
output with other output dumps for analysis.
1Note that the use of some of the symbols in this Appendix is inconsistent with their use
elsewhere in this dissertation.
98
C.2 HNBody-iabl Interface
Without full access to the HNBody?s source code, there were two choices available
for the interface: (1) write a new driver for HNBody using the documented API or
(2) use system calls to HNBody with its mature, distributed driver. We chose the
latter because it limited the initial work required to produce a functioning program
despite the loss of efficiency. As mentioned above, because the integration of the
N-body equations of motion dominates the computational load, the efficiency lost
from using system calls instead of linked function calls does not significantly affect
the overall efficiency.
Before calling HNBody, iabl writes to files the input data from the initial condi-
tions generator as well as the amount of time for HNBody to integrate. The initial
guess for how long the encounter will take is that of a simple hyperbolic encounter
of the interloper about the center of mass of the binary from an initial distance
ri = 30(a20 +b2)1/2 to a final distance rf = 30a0. The time is given by
t0 =
radicalBigg
a3
GM [e(sinhFi + sinhFf)?(Fi +Ff)] (C.1)
where a and e are the semimajor axis and eccentricity of the hyperbolic orbit and
F is the hyperbolic analog of the eccentric anomaly given by
F = cosh?1
bracketleftbigg1
e
parenleftbigg
1 + ra
parenrightbiggbracketrightbigg
(C.2)
(Goldstein et al. 2002; Richardson et al. 1998). If the encounter is determined not
to have finished in this time, the integration time is doubled. This time doubling
repeats 10 times after which the time for the next integration grows as ti = 1.2ti?1.
After 25 time extensions, iabl resets the encounter with new initial conditions.
99
C.3 Examination of Integration
One of the most sophisticated parts of iabl is its series of routines used to determine
whether the integration has completed. The first test that iabl runs is to check
for merger by gravitational waves. Though iabl currently uses the HNBody module
ExitCondition to determine merger status, the simulations in this dissertation looked
for errors from HNBody either from the exit value or the output, which indicates
that objects came too close to each other for HNBody to continue.
If a merger has not occurred, iabl checks for an ionization. In this check iabl (1)
looks to see that the total energy of the system is positive; (2) checks to see that no
pair of objects are bound to each other; (3) checks to see that no object is bound
to the center of mass of the other two objects; (4) checks to see that each particle?s
velocity vector is pointed away from the center of mass of the system (r?v < 0 in
barycentric coordinates). If all four tests are passed, a static variable flag is set, and
the integration is extended. If iabl finds that the systems passes all tests twice in a
row, an ionization is determined to have occurred. The twice-in-a-row test is used
because of pathological cases in which no ionization occurred but it appears to be
so. It is very unlikely for such already unlikely cases to occur at the second point
of integration.
If no merger and no ionization has occurred, iabl checks to see if the encounter
has resolved itself. First iabl determines the closest pair of objects and hypothesizes
that they are the binary. The encounter is determined not to have finished, and
the integration is extended if (1) the third object is within a distance of 30a0, (2)
the closest pair of objects is not bound to each other, (3) the third object is bound
to the center of mass of the binary with semimajor axis a3 < 2000 AU, or (4) the
third object is unbound but approaching the center of mass of the system. If the
100
third object (1) is bound to the center of mass of the binary with a semimajor axis
30aib < a3 < 2000 AU, where aib is the semimajor axis of the inner binary, (2) is
separated from the center of mass of the binary by a distance r > 30aib, and (3)
has its velocity vector pointed away from the center of mass of the system, then
the two-body approximation is run. If all of these tests fail, then the encounter is
determined to have finished.
C.4 Two-Body Approximator
The two-body approximator in iabl is conceptually simple yet speeds up typical
simulations by a factor between 7 and 15. Empirically, we have found that roughly
10% of the encounters in a simulation take 90% of the total time. These encounters
are long, resonant encounters in which all three objects become temporarily bound
to each other. When the result of a binary-single encounter, these systems are
unstable and eventually disrupt themselves. The two body approximator speeds
this process up by treating hierarchical triples as two sets of binaries: the inner and
the outer. The outer binary approximates the inner binary as a single particle at
the location of its center of mass. The approximator advances the outer binary until
it is returning toward the center of mass at a distance r = 30aib. The approximator
calculates the time elapsed and advances the inner binary?s orbit. Then iabl calls
HNBody to integrate the orbit again.
101
Appendix D
Other Applications of the Code
Though the primary purpose of iabl was for the research in this dissertation, it has
been used for several other studies, three of which we outline here.
D.1 Planet in Close Triple System
Therecentlyfoundplanetinthehierarchicaltriple-stellarsystemHD188753(Konacki
2005) motivated the study of dynamical interactions that can create observed plan-
ets in multiple systems. The planet has a minimum mass of 1.14 times that of
Jupiter, 1.14 MJ, and is in an a = 0.045 AU orbit around a 1.06 Mcircledot star. This star,
the optical primary of the system, is in an a = 12.3 AU eccentric orbit (e = 0.5)
around a tight, single-line spectroscopic binary. The tight binary is composed of a
0.96 Mcircledot star with a 0.67 Mcircledot companion in an a = 0.67 AU orbit with an eccentric-
ity e = 0.1. The ? 6 AU minimum separation between the primary and the tight
binary poses a problem for the in situ formation of the planet, which is expected to
have formed at a distance of ? 3 AU before it migrated close to the star. Because
of its high eccentricity, the tight binary will truncate the protoplanetary disk at a
distance of 1.3 AU, preventing it from forming giant planets (Jang-Condell 2005;
102
Pichardo et al. 2005).
If the system could not have formed as we currently observe it, then it must
have evolved to this state. Since most stars, especially stars in multiple systems,
form in clusters, dynamical interactions are likely to play an important role in the
evolution of this system. In addition, while far from definitive, the eccentricity of
the primary with respect to the spectroscopic binary (e = 0.5) is suggestive of a
dynamical interaction. Binaries that have undergone a significant interaction have
their eccentricities distributed such that after a single encounter the probability of
having an eccentricity e is P(e) = 2e. For a binary with an initially low eccentricity,
dynamical interactions are the most likely method to boost the eccentricity. Jang-
Condell (2005) suggests that the tight binary may have captured the planet-bearing
star after the planet had formed, but purely dynamical captures involving three
stars cannot lead to a stable triple-system. A dynamical interaction that ends in a
stable triple system must use a fourth star that does not end up in the system. The
three most likely interactions involving four stars are outlined in Figure D.1. In all
cases the planet forms before the dynamical interaction and, due to its small mass
and small separation, plays only a minor role in the essential dynamics of creating
a triple system.
In method E, ?Formation by Exchange,? the star with the planet (A) encounters
a close, hierarchical triple system composed of a tight binary (B) and another star
(X), with masses MA, MB, and MX, respectively. Star A interacts with the triple
and exchanges into the binary thus ejecting star X. Although triple-single exchange
interactions are not well studied, we may approximate the tight binary as a single
object with its combined mass. In binary-single interactions, the most massive
objects tend to end up in the binary (e.g., Heggie et al. 1996). Thus if MX <
MA,MB, the exchange reaction is more likely to happen.
103
In method H, ?Formation by Hardening,? star A is in a hierarchical triple system
with binary B with a separation wide enough for a planet to form (> 50?100 AU).
The triple system then faces repeated encounters with several different stars X. The
interactions extract energy from the wide binary and carry it away from the system.
This hardens the binary in a method analogous to that of Miller & Hamilton (2002b).
The circular speed of A about B in HD 188753 is vcirc = (GM/a)1/2 = 13.8 km s?1,
and the typical velocity dispersion of an open cluster isvclust ? 0.3 to 3 km s?1 (e.g.,
Binney & Tremaine 1987). For the less-studied triple-single interactions, it is not
clear whether the wide or tight binary is more likely to harden.
Method EH, ?Formation by Exchange and Hardening,? involves both an ex-
change and a hardening in a single interaction. Star A starts in a wide binary
system with star X and interacts with binary B. Binary B exchanges into the sys-
tem and ejects star X. Although binary-binary interactions are more likely to occur
in a cluster than triple-single interactions, this process is the least likely to occur.
The reason is that in order for the exchange to be favored, MX must be the smallest
mass, but for an equal-energy orbit that increases in mass, the semimajor axis will
increase. Each of the above three methods may be complicated further by replacing
any star with a binary with the same total mass and with separation much smaller
than the next closest object.
As a proof of concept, iabl was used to integrate a triple-single-plus-planet en-
counter that displays the desired exchange and the observed configuration (Fig-
ure D.2). This is done by starting with the observed triple system and integrating
the encounter backwards. The final results become the initial conditions for the
desired reaction. For simplicity, the ejected star is assumed to have the same mass
as star A MX = MA and a final velocity at infinity vf = vclust = 1 km s?1. The
initial configuration that leads to the triple system has a wide binary semimajor axis
104
Figure D.1: Three methods of forming a stellar triple system with a planet. The
methods are labeled, from top to bottom, E, H, and EH.
of 10.1 AU, eccentricity e = 0.61 and a relative velocity of 9.8 km s?1 with respect
to star A. Although integrating backwards from a known configuration is useful for
testing the crude feasibility of such scenarios, it is not appropriate for determining
the probability of these exchanges because the discovered initial conditions will not
be representative of the environment in the host cluster, such as the high relative
velocity in this example. A broad parameter search of likely initial conditions would
be necessary. The conclusions found in this study were generally the same as those
found by Portegies Zwart & McMillan (2005) and Pfahl (2005b)
D.2 Near-Earth Asteroid Binary Disruption
As part of a large project studying the tidal disruption of near-Earth asteroids
(NEAs), iabl was used to compute the fraction of binaries disrupted on close ap-
proach to Earth for a variety of parameters, including impact parameter, semimajor
105
Figure D.2: Orbits of triple-single-plus planet encounter that results in an ex-
change into the configuration observed for HD 188753 by method E. The magenta
curve is star A (and planet), which comes in from the upper-left. The initial triple
system comes in from the bottom-right, and contains the close binary system (blue
and red) and the star that is eventually ejected (black).
axis, relative velocity, and mass ratio of the binary. Figure D.3 shows an example
curve compared to simulations by Bottke & Melosh (1996) for a relative velocity
v? = 12 km s?1, binary asteroid mass ratio q = 0.125 with an a = 3 km circular
orbit for asteroids of mass density ? = 2.6 g cm?3.
106
FigureD.3: Plotoffractionofbinariesdisruptedasafunctionofimpactparameter
in units of the Earth?s radius. Symbols are from simulations by Bottke & Melosh
(1996), and the line is from simulations run by iabl. See text for details. The
agreement between the two sets of simulations is very good. Data provided by
K. Walsh.
D.3 IMBH-IMBH Binaries
Recent simulations of young stellar clusters have shown that collisional runaways
can lead to two or more IMBHs (G?urkan et al. 2006). These IMBHs would sink to
the center and eventually form a binary with each other and, through subsequent
interactions with stars in the cluster, harden and ultimately merge. During inspiral,
these massive binaries are extremely loud LISA signals and may possibly be used
as star-formation tracers to the young stellar clusters in which they formed. To test
the lag time between the formation of an IMBH-IMBH binary and merger, we ran
107
simulations similar to those in this dissertation but with initial binary m0 = m1 =
1000Mcircledot, a = 20000 AU with encounters withm2 = 10Mcircledot stars that were assumed
to be tidally disrupted if they got too close to either IMBH. We found that the time
to merge was roughly 2 Myr and that the merging binaries were not likely to be
very eccentric when detectable by LISA.
108
Bibliography
Aarseth, S. J. & Lecar, M. 1975, ARA&A, 13, 1
Ando, M. & the TAMA collaboration. 2002, Classical and Quantum Gravity, 19,
1409
Angelini, L., Loewenstein, M., & Mushotzky, R. F. 2001, ApJ, 557, L35
Arons, J. 1992, ApJ, 388, 561
Baker, J. G., Centrella, J., Choi, D.-I., Koppitz, M., van Meter, J. R., & Miller,
M. C. 2006, ArXiv Astrophysics e-prints, arXiv:astro-ph/0603204
Baraffe, I., Heger, A., & Woosley, S. E. 2001, ApJ, 550, 890
Barish, B. C. 2000, Advances in Space Research, 25, 1165
Begelman, M. C. 2002, ApJ, 568, L97
?. 2006, ArXiv Astrophysics e-prints, arXiv:astro-ph/0602022
Binney, J. & Tremaine, S. 1987, Galactic dynamics (Princeton, NJ, Princeton Uni-
versity Press, 1987, 747 p.)
Blanchet, L., Qusailah, M. S. S., & Will, C. M. 2005, ApJ, 635, 508
Bottke, W. F. & Melosh, H. J. 1996, Icarus, 124, 372
Bruns, H. 1887, Acta Mathematica, 11, 25
Colbert, E. J. M. & Mushotzky, R. F. 1999, ApJ, 519, 89
Colbert, E. J. M. & Ptak, A. F. 2002, ApJS, 143, 25
Colpi, M., Mapelli, M., & Possenti, A. 2003, ApJ, 599, 1260
109
Colpi, M., Possenti, A., & Gualandris, A. 2002, ApJ, 570, L85
Cruz-Gonz?alez, C. & Poveda, A. 1971, Ap&SS, 13, 335
Cruz-Gonz?alez, C. & Poveda, A. 1972, in Gravitational N-Body Problem, Proceed-
ings of IAU Colloq. 10, held in Cambridge, England, 12-15 August, 1970. Edited
by Myron Lecar. D. Reidel Publishing Company, 1971., p.99, ed. M. Lecar, 99?
Cutler, C. & Thorne, K. S. 2002, ArXiv General Relativity and Quantum Cosmology
e-prints, arXiv:gr-qc/0204090
Damour, T. & Gopakumar, A. 2006, ArXiv General Relativity and Quantum Cos-
mology e-prints, arXiv:gr-qc/0602117
Danzmann, K. 2000, Advances in Space Research, 25, 1129
Dewangan, G. C., Griffiths, R. E., & Rao, A. R. 2006, ArXiv Astrophysics e-prints,
arXiv:astro-ph/0602472
Ebisuzaki, T., Makino, J., Tsuru, T. G., Funato, Y., Portegies Zwart, S., Hut, P.,
McMillan, S., Matsushita, S., Matsumoto, H., & Kawabe, R. 2001, ApJ, 562, L19
Einstein, A. 1905, Annalen der Physik, 17, 891
?. 1915, Sitzungsberichte der K?oniglich Preu?ischen Akademie der Wissenschaften,
844
?. 1916, Sitzungsberichte der K?oniglich Preu?ischen Akademie der Wissenschaften,
688
Fabbiano, G., Zezas, A., & Murray, S. S. 2001, ApJ, 554, 1035
Farmer, A. J. & Phinney, E. S. 2003, MNRAS, 346, 1197
Favata, M., Hughes, S. A., & Holz, D. E. 2004, ApJ, 607, L5
Ferraro, F. R., Possenti, A., Sabbi, E., Lagani, P., Rood, R. T., D?Amico, N., &
Origlia, L. 2003, ApJ, 595, 179
Fidecaro, F. & VIRGO Collaboration. 1997, in General Relativity and Gravitational
110
Physics; Proceedings of the 12th Italian Conference, edited by M. Bassan, V.
Ferrari, M. Francaviglia, F. Fucito, and I. Modena. World Scientific Press, 1997.,
p.163, 163
Fitchett, M. J. 1983, MNRAS, 203, 1049
Flanagan, ?E. ?E. & Hughes, S. A. 1998a, Phys. Rev. D, 57, 4535
?. 1998b, Phys. Rev. D, 57, 4566
Fregeau, J. M., Cheung, P., Portegies Zwart, S. F., & Rasio, F. A. 2004, MNRAS,
352, 1
Fregeau, J. M., Joshi, K. J., Portegies Zwart, S. F., & Rasio, F. A. 2002, ApJ, 570,
171
Freitag, M., G?urkan, M. A., & Rasio, F. A. 2006a, MNRAS, 335
Freitag, M., Rasio, F. A., & Baumgardt, H. 2006b, MNRAS, 324
Fryer, C. L. & Kalogera, V. 2001, ApJ, 554, 548
Fryer, C. L., Woosley, S. E., & Heger, A. 2001, ApJ, 550, 372
Fullerton, L. W. & Hills, J. G. 1982, AJ, 87, 175
Gaffney, N. I., Lester, D. F., & Telesco, C. M. 1993, ApJ, 407, L57
Gammie, C. F. 1998, MNRAS, 297, 929
Gebhardt, K., Kormendy, J., Ho, L. C., Bender, R., Bower, G., Dressler, A., Faber,
S. M., Filippenko, A. V., Green, R., Grillmair, C., Lauer, T. R., Magorrian, J.,
Pinkney, J., Richstone, D., & Tremaine, S. 2000, ApJ, 543, L5
Gebhardt, K., Rich, R. M., & Ho, L. C. 2002, ApJ, 578, L41
?. 2005, ApJ, 999, L999
Gerssen, J., van der Marel, R. P., Gebhardt, K., Guhathakurta, P., Peterson, R. C.,
& Pryor, C. 2002, AJ, 124, 3270
Gilbert, N. E. 1901, Ph.D. Thesis
Goad, M. R., Roberts, T. P., Knigge, C., & Lira, P. 2002, MNRAS, 335, L67
111
Goldstein, H., Poole, C., & Safko, J. 2002, Classical Mechanics (San Francisco:
Addison-Wesley, 3rd ed.)
G?ultekin, K., Miller, M. C., & Hamilton, D. P. 2004, ApJ, 616, 221
?. 2006, ApJ, 640, 156
G?urkan, M. A., Fregeau, J. M., & Rasio, F. A. 2006, ApJ, 640, L39
G?urkan, M. A., Freitag, M., & Rasio, F. A. 2004, ApJ, 604, 632
Heggie, D. C. 1975, MNRAS, 173, 729
Heggie, D. C., Hut, P., & McMillan, S. L. W. 1996, ApJ, 467, 359
Heggie, D. C. & Rasio, F. A. 1996, MNRAS, 282, 1064
Holmberg, E. 1941, ApJ, 94, 385
Hopman, C. & Portegies Zwart, S. 2005, ArXiv Astrophysics e-prints, arXiv:astro-
ph/0506181
Hornschemeier, A. E., Bauer, F. E., Alexander, D. M., Brandt, W. N., Sargent,
W. L. W., Bautz, M. W., Conselice, C., Garmire, G. P., Schneider, D. P., &
Wilson, G. 2003, AJ, 126, 575
Hulse, R. A. & Taylor, J. H. 1975, ApJ, 195, L51
Hut, P. & Bahcall, J. N. 1983, ApJ, 268, 319
Hut, P. & Inagaki, S. 1985, ApJ, 298, 502
Irwin, J. A., Athey, A. E., & Bregman, J. N. 2003, ApJ, 587, 356
Iyer, B. R. & Will, C. M. 1993, Physical Review Letters, 70, 113
Jackson, J. D. 1999, Classical Electrodynamics (New York: Wiley, 3rd ed.)
Jang-Condell, H. 2005, ArXiv Astrophysics e-prints, arXiv:astro-ph/0507356
Julliard-Tosel, E. 2000, Celestial Mechanics and Dynamical Astronomy, 76, 241
Kaaret, P., Prestwich, A. H., Zezas, A., Murray, S. S., Kim, D.-W., Kilgard, R. E.,
Schlegel, E. M., & Ward, M. J. 2001, MNRAS, 321, L29
Kaaret, P., Ward, M. J., & Zezas, A. 2004, MNRAS, 351, L83
112
King, A. R., Davies, M. B., Ward, M. J., Fabbiano, G., & Elvis, M. 2001, ApJ, 552,
L109
Konacki, M. 2005, Nature, 436, 230
Kozai, Y. 1962, AJ, 67, 591
Kudritzki, R.-P. The First Stars. Proceedings of the MPA/ESO Workshop held at
Garching, Germany, ed. , A. WeissT. G. Abel & V. Hill, 127
Kuntz, K. D., Gruendl, R. A., Chu, Y.-H., Chen, C.-H. R., Still, M., Mukai, K., &
Mushotzky, R. F. 2005, ApJ, 620, L31
Lee, M. H. 1993, ApJ, 418, 147
Liu, J.-F., Bregman, J. N., Irwin, J., & Seitzer, P. 2002, ApJ, 581, L93
Madau, P. & Rees, M. J. 2001, ApJ, 551, L27
Maillard, J. P., Paumard, T., Stolovy, S. R., & Rigaut, F. 2004, A&A, 423, 155
Martel, K. & Poisson, E. 1999, Phys. Rev. D, 60, 124008
Matsubayashi, T., Shinkai, H., & Ebisuzaki, T. 2004, ApJ, 614, 864
Matsumoto, H., Tsuru, T. G., Koyama, K., Awaki, H., Canizares, C. R., Kawai, N.,
Matsushita, S., & Kawabe, R. 2001, ApJ, 547, L25
McCrady, N., Gilbert, A. M., & Graham, J. R. 2003, ApJ, 596, 240
Miller, J. M., Fabian, A. C., & Miller, M. C. 2004, ApJ, 614, L117
Miller, M. C. 1984, Gifted Children Newsletter, 10, 15
?. 2002, ApJ, 581, 438
Miller, M. C. & Colbert, E. J. M. 2004, International Journal of Modern Physics D,
13, 1
Miller, M. C. & Hamilton, D. P. 2002a, ApJ, 576, 894
?. 2002b, MNRAS, 330, 232
Misner, C. W., Thorne, K. S., & Wheeler, J. A. 1973, Gravitation (San Francisco:
W.H. Freeman and Co., 1973)
113
Mukai, K., Still, M., Corbet, R. H. D., Kuntz, K. D., & Barnard, R. 2005, ApJ, 634,
1085
Newhall, X. X., Standish, E. M., & Williams, J. G. 1983, A&A, 125, 150
O?Leary, R. M., Rasio, F. A., Fregeau, J. M., Ivanova, N., & O?Shaughnessy, R.
2006, ApJ, 637, 937
Pakull, M. W. & Mirioni, L. 2002, ArXiv Astrophysics e-prints, arXiv:astro-
ph/0202488
Peters, P. C. 1964, Physical Review, 136, 1224
Pfahl, E. 2005a, ApJ, 626, 849
?. 2005b, ApJ, 635, L89
Pichardo, B., Sparke, L. S., & Aguilar, L. A. 2005, MNRAS, 359, 521
Poincar?e, H. 1892, Les methodes nouvelles de la mecanique celeste (Paris, Gauthier-
Villars et fils, 1892-99.)
Portegies Zwart, S. F. & McMillan, S. L. W. 2000, ApJ, 528, L17
?. 2002, ApJ, 576, 899
?. 2005, ApJ, 633, L141
Ptak, A. & Colbert, E. 2004, ApJ, 606, 291
Quinlan, G. D. 1996, New Astronomy, 1, 35
Quinlan, G. D. & Shapiro, S. L. 1989, ApJ, 343, 725
Rauch, K. P. & Hamilton, D. P. in preparation, AJ
Reynolds, C. S., Loan, A. J., Fabian, A. C., Makishima, K., Brandt, W. N., &
Mizuno, T. 1997, MNRAS, 286, 349
Richardson, D. C., Bottke, W. F., & Love, S. G. 1998, Icarus, 134, 47
Roberts, T. P., Goad, M. R., Ward, M. J., Warwick, R. S., O?Brien, P. T., Lira, P.,
& Hands, A. D. P. 2001, MNRAS, 325, L7
Schilling, R. 1998, in AIP Conf. Proc. 456: Laser Interferometer Space Antenna,
114
Second International LISA Symposium on the Detection and Observation of Grav-
itational Waves in Space, 217?221
Schneider, R., Ferrara, A., Natarajan, P., & Omukai, K. 2002, ApJ, 571, 30
Sigurdsson, S. & Hernquist, L. 1993, Nature, 364, 423
Sigurdsson, S. & Phinney, E. S. 1993, ApJ, 415, 631
?. 1995, ApJS, 99, 609
Strohmayer, T. E. & Mushotzky, R. F. 2003, ApJ, 586, L61
Swartz, D. A., Ghosh, K. K., Tennant, A. F., & Wu, K. 2004, ApJS, 154, 519
Taylor, J. H., Fowler, L. A., & McCulloch, P. M. 1979, Nature, 277, 437
von Hoerner, S. 1960, Zeitschrift fur Astrophysik, 50, 184
Webbink, R. F. 1985, in IAU Symp. 113: Dynamics of Star Clusters, 541?577
Webster, D. L. 1913, Ph.D. Thesis
Wen, L. 2003, ApJ, 598, 419
Whittaker, E. T. 1989, A Treatise on the Analytical Dynamics of Particles and Rigid
Bodies (Cambridge, UK: Cambridge University Press)
Will, C. M. 2004, ApJ, 611, 1080
Winter, L. M., Mushotzky, R. F., & Reynolds, C. S. 2005, ArXiv Astrophysics
e-prints, arXiv:astro-ph/0512480
Yu, Q. & Tremaine, S. 2003, ApJ, 599, 1129
115