GRAVITATIONAL-WAVE
ASTRONOMY WITH
COALESCING COMPACT
BINARIES: DETECTION AND
PARAMETER ESTIMATIONWITH ADVANCED DETECTORSRORY JAMES EDWIN SMITH

A thesis submitted to the
University of Birmingham
for the degree of
DOCTOR OF PHILOSOPHY

Astrophysics and Space Research
School of Physics and Astronomy
College of Engineering and Physical Sciences
University of Birmingham
January 2013

Acknowledgements
It is my sincere pleasure to thank the following people: My supervisors Ilya
Mandel and Alberto Vecchio for their mentorship, support and patience; Chad
Hanna for his mentorship, and for always keeping an open door despite his busy
schedule; The members of the Bayesian and CBC groups for all their help; All
my friends in the ASR group at Birmingham, especially Charlotte Bond, Mark
Burke, Ludovico Carbone, Chris Collins, Paul Fulda, Kat Grover, Carl-Johan
Haster, Konstantinos Kolokythas, Ilya Mandel, Chiara Mingarelli, Trevor Sidery
and Will Vousden for years of fun and adventures; Andy Pearsall, Brigid Jones,
Amelia Perry, Charlotte Kimbly, Stacey McGowan, Jim Phillips, Peter Rothery,
Patrick Scott, Jim Smyth and Mat Turnbull for their friendship over the years;
The Dog & Dolphin, refuge from real world problems; And my family, especiallymy parents, for their love and support.
Several of the chapters in this thesis were the result of collaborations andbenefited from discussions with several people. Chapter 2 was based on work
done in collaboration with Ilya Mandel and Alberto Vecchio, and benefited fromdiscussions with Eliu Huerta. Chapter 3 was based on work done in collaboration
with Kip Cannon, Chad Hanna, Drew Keppel and Ilya Mandel, and benefited
from discussions with Tyson Littenberg. Chapter 4 was based on work done incollaboration with Chad Hanna, Ilya Mandel and Alberto Vecchio, and benefited
from discussions with several people, in particular: Kipp Cannon, Scott Field,
Drew Keppel, Vivien Raymond, Manuel Tiglio and Alan Weinstein. Finally, I
am grateful for the many helpful comments on this thesis from Stephen Fairhurst
and Andreas Freise.
My work has been supported by the STFC and the Perimeter Institute for
Theoretical Physics.

iAbstract

The current generation of interferometric gravitational-wave detectors, LIGO and Virgo, are undergoing upgrade to their so-called advanced phase. These instruments, together with new instruments inJapan and India, KAGRA and LIGO India, will form a network ofadvanced gravitational-wave detectors with which detections are expected to become routine. Amongst the prime sources for gravitationalwave astronomy are coalescing compact binaries consisting of neutron stars and/or black holes. Filtering detector data to detect thesesources relies on precise templates of the expected gravitational-wave
signals. In addition, estimating the parameters encoded in the signals(masses, spins etc...) requires sophisticated Bayesian inference techniques. Templates are typically computationally expensive to generate and can be a bottle-neck in data analysis.
Here we focus on two aspects of gravitational-wave astronomy using coalescing compact binaries. The first part of this thesis focuses on studying the requirements of template waveforms to detect
intermediate-mass black holes through the coalescence of a stellarmass companion into an intermediate-mass black hole in Advanced
LIGO. The second part of this thesis focuses on numerical and analytic techniques to improve the efficiency of (Bayesian) parameter
estimation on coalescing binaries when parameter estimation is dominated by template waveform generation. Such efficiency improvements to parameter estimation are crucial for gravitational-wave astronomy using advanced detectors.

Contents
AcknowledgementsiContents

iii

List of Figuresv

List of Tablesvii

1 INTRODUCTION

12 DETECTING INTERMEDIATE MASS-RATIO BINARIES WITHADVANCED DETECTORS
21
2.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

21

2.2

Matched filtering of gravitational-wave signals . . . . . . . . . . .

25

2.3

Waveforms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

28

2.3.1

EOBNR . . . . . . . . . . . . . . . . . . . . . . . . . . . .

29

2.3.1.1

Inspiral-plunge waveform . . . . . . . . . . . . .

29

2.3.1.2

Merger-ringdown waveform . . . . . . . . . . . .

30

2.3.2

HG waveforms . . . . . . . . . . . . . . . . . . . . . . . . .

30

2.3.3

TaylorT4 . . . . . . . . . . . . . . . . . . . . . . . . . . . .

32

2.4

SNR from inspiral, merger and ringdown . . . . . . . . . . . . . .

32

2.5

E↵ectiveness of inspiral-only templates for IMRAC searches . . .

38

2.6

Comparison of inspiral-only waveforms . . . . . . . . . . . . . . .

44

2.7

Discussion and conclusion . . . . . . . . . . . . . . . . . . . . . .

47

iii

CONTENTS

3 TOWARDS RAPID COMPACT BINARY PARAMETER ESTIMATION USING INTERPOLATED WAVEFORMS
49
3.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

49

3.2

CBC parameter estimation . . . . . . . . . . . . . . . . . . . . . .

51

3.2.1

MCMC . . . . . . . . . . . . . . . . . . . . . . . . . . . .

54Singular value decomposition of gravitational-wave template banks
for CBC sources . . . . . . . . . . . . . . . . . . . . . . . . . . . .

56

3.4

Interpolating template waveforms using the SVD . . . . . . . . .

59

3.5

Parameter estimation using interpolated waveforms . . . . . . . .

62

3.6

Results: Comparison of parameter estimates using interpolated
and non-interpolated waveforms . . . . . . . . . . . . . . . . . . .

67Discussion and conclusion . . . . . . . . . . . . . . . . . . . . . .

71

3.3

3.7

4 PARAMETER ESTIMATION USING INTERPOLATED TIMEDOMAIN TEMPLATE WAVEFORMS
73
4.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

73

4.2

Interpolating inspiral-only time-domain waveforms . . . . . . . . .

75

4.2.1

Parameter estimation using interpolated TaylorT4 waveforms: 5-dimensional case . . . . . . . . . . . . . . . . . .

78

4.2.1.1

Prior distributions . . . . . . . . . . . . . . . . .

78

4.2.1.2

SVD set up . . . . . . . . . . . . . . . . . . . . .

79

4.2.1.3

Results . . . . . . . . . . . . . . . . . . . . . . .

79

4.3

Interpolating inspiral-merger-ringdown waveforms . . . . . . . . .

81

4.4

Future work and Discussion . . . . . . . . . . . . . . . . . . . . .

84

5 RAPIDLY EVALUATING THE COMPACT BINARY LIKELIHOOD FUNCTION VIA INTERPOLATION
85
5.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

85

5.1.1

Directly interpolating the likelihood function . . . . . . . .

88

5.1.2

Likelihood interpolation: Examples . . . . . . . . . . . . .

89

5.1.2.1

Example 1: High-mass binary black holes . . . .

90

5.1.2.2

Example 2: Binary neutron stars . . . . . . . . .

92

5.1.2.3

Practical considerations . . . . . . . . . . . . . .

92

iv

CONTENTS

5.1.2.4

Discussion and conclusion . . . . . . . . . . . . .

6 SUMMARY

94
97

A FREQUENTLY USED QUANTITIES
101
A.1 Definitions of discrete representations of continuous functions . . . 101
A.2 Post-Newtonian phasing formulae . . . . . . . . . . . . . . . . . . 103
A.2.1 TaylorF2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
A.2.2 TaylorT4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
B SINGULAR VALUE DECOMPOSITION OF M ⇥ N MATRICES105
B.1 Definitions of the singular value decomposition of M ⇥ N matrices 105
B.1.1 Step-by-step example . . . . . . . . . . . . . . . . . . . . . 106
B.1.2 Geometric picture of the SVD . . . . . . . . . . . . . . . . 108
REFERENCES

111v

CONTENTS

vi

List of Figures
1.1

Inspiral phase of a gravitational wave from compact binary coalescence. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

1.2

Anticipated Advanced LIGO noise sources. . . . . . . . . . . . . .

8

1.3

Noise amplitude spectral densities of various Earth based gravitationalwave detectors. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
92.1

Advanced LIGO high-power, zero-detuned noise PSD . . . . . . .

33

2.2

RMS strain of IMRAC sources in Advanced LIGO. . . . . . . . .

35

2.3

SNR of IMRAC sources in Advanced LIGO as a function of total
mass and mass ratio. . . . . . . . . . . . . . . . . . . . . . . . . .

37

E↵ectiveness of inspiral-only templates to filter inspiral, merger
and ringdown signals. . . . . . . . . . . . . . . . . . . . . . . . . .

41

E↵ectiveness of inspiral-only EOBNR templates at filtering HG signals (left) and TaylorT4 signals (right) as a function of the sourcecomponent masses encoded in the signal. . . . . . . . . . . . . . .

45

2.4
2.5

3.1

Coalescing compact binary system relative to a ground based gravitationalwave detector. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

3.2

Initial LIGO noise power spectral density. . . . . . . . . . . . . .

64

3.3

Patching of parameter space. . . . . . . . . . . . . . . . . . . . . .

67

3.4

Marginalized one-dimensional posterior probability density functions produced using interpolated and non-interpolated template
waveforms. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

68

vii

LIST OF FIGURES

4.1
4.2

4.3

5.1

Histogram of mismatches between (normalized) interpolated andnon-interpolated TaylorT4 waveforms. . . . . . . . . . . . . . . .
Marginalized one-dimensional posterior probability density functions produced using interpolated and non-interpolated TaylorT4
waveforms. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Histogram of mismatches between (normalized) interpolated andnon-interpolated EOBNR waveforms. . . . . . . . . . . . . . . . .
Interpolated and non-interpolated log likelihoods (top), and percentage error (bottom) for a data set containing a gravitationalwave signal from the coalescence of binary black holes. . . . . . .

77

80
83

91

B.1 Action of a matrix on the unit circle. . . . . . . . . . . . . . . . . 108

viii

List of Tables
2.1
2.2

3.1

3.2

4.1

E↵ectiveness of inspiral-only searches. . . . . . . . . . . . . . . . .
E↵ectiveness of inspiral-only EOBNR template waveforms in recovering signals modelled using di↵erent waveform families. . . . .
Mean posterior parameter estimates (and standard deviations) of
the marginalized PDFs using interpolated and non-interpolated
SPA waveforms (Fig. 3.4). . . . . . . . . . . . . . . . . . . . . . .
Computational time of template waveform generation in units ofcomputational time of interpolated waveforms. . . . . . . . . . . .

43
46

68
70

Mean posterior value parameter estimates (and standard deviations) of the marginalized PDFs using interpolated and non-interpolated
TaylorT4 waveforms (Fig. 4.2). . . . . . . . . . . . . . . . . . . . 81ix

LIST OF TABLES

xChapter 1
INTRODUCTION
Gravitational-wave astronomy will open an entirely new observational window to
the Universe with the direct detections of gravitational waves from high energy astrophysical phenomena. A network of ground-based interferometric gravitationalwave observatories (LIGO-Hanford and LIGO-Livingston (USA) and Virgo (Italy)[1]) are undergoing upgrades to their so-called “advanced” configurations which
will be an order of magnitude more sensitive than their initial configurations
and are scheduled to begin collecting data around 2015 [2]. In addition, two
new ground-based interferometric detectors in Japan and India, KAGRA [3] and
LIGO India [4] respectively, are expected to be online around 2020. This willusher in the “advanced detector era” in which gravitational-wave detections are
expected to become routine [5].
The advanced detector network is expected to detect gravitational waves froma variety of fascinating astrophysical and cosmological events; from the late stagesof the coalescence of compact binaries [5; 6; 7; 8] to a stochastic background ofgravitational waves of astrophysical and cosmological origin [9; 10]. Of all thetarget sources, coalescing compact binaries are perhaps the the most likely candidates for detection [5]. Compact binary systems, consisting of neutron stars
and/or black holes are believed to be in abundance throughout the Universe. Indeed, the first indirect evidence for the existence of gravitational radiation camefrom the observation of the binary pulsar system PSR B1913+16, known as the“Hulse-Taylor” binary pulsar after its discoverers, which is undergoing orbital
decay through the emission of gravitational waves [11]. For their discovery, Hulse

11. INTRODUCTION

and Taylor were awarded the Nobel Prize in 1993. The Hulse-Taylor binary pulsar is undergoing orbital decay very slowly and will not merge for around threehundred million years. In contrast the advanced-detector network will directly
observe the gravitational waves from the final minutes to seconds of the coalescence of compact binaries before they merge to form a single object. This makes
compact binary coalescence (CBC) an ideal laboratory for probing dynamical,
strong-field gravity.
With observations of gravitational waves from coalescing compact binaries
we will be able to make precise measurements of their masses and spins as well
as of their demographics [12; 13]. These measurements are essential for astrophysics/astronomy and fundamental physics using gravitational waves from compact binaries. Determination of the masses and spins will distinguish neutronstars and black holes and can be used to test the null-hypothesis that black holes
are described by the Kerr metric [14; 15]. In addition, determination of source
demographics, together with merger rates, places constraints on the formation
mechanisms of binary systems [16]. While there are around ten known binary
neutron star (NS-NS) systems [17] there are currently no known binary black hole
(BH-BH) or black-hole-neutron-star (BH-NS) systems [18]. Despite this, BH-BH
and BH-NS systems are believed to be fairly common throughout the Universe:
population synthesis models suggest that realistic detection rates of BH-BH and
BH-NS systems in Advanced LIGO/Virgo (aLIGO/Virgo) could be as high as
40 yr 1 [5]. Unsurprisingly, binary black holes are difficult to observe in the electromagnetic spectrum and gravitational waves are likely the only means by which
they can be detected.
There are numerous open issues relating to black hole systems which gravitational -wave detections will be able to make important contributions to. In
particular gravitational-wave detectors will be able to search for binary systems
with total masses up to around 500 M [19] and will therefore be able to provide
the first unambiguous observations of “intermediate-mass” black holes (IMBHs).
Intermediate-mass is a blanket term which refers to black holes with masses(roughly) in the range 50 103 M . Believed to form in globular clusters [e.g.,20; 21; 22], IMBHs are extremely difficult to observe in the electromagnetic spectrum and there are currently only tentative results suggesting their existence.

2

Detecting IMBHs will thus have important implications for globular cluster dynamics. In addition, detections of compact binary coalescence will be importantin addressing the so called black hole “mass-gap” problem [23]: The stellar massdistribution is continuous in the range 0.1 100 M and one would expect thesame to be true of ensuing compact remnants. However, while neutron stars
have been observed to have masses up to around 2 M , the lightest black holes
are observed to have masses around 5 M , with no compact objects with masses
in-between [24]. Gravitational-wave observations of BH-BH or BH-NS systems
will therefore be able to shed light on whether there is a genuine gap between theheaviest neutron stars and the lightest black holes.
The central theme of the work presented here is the detection of, and parameter estimation on gravitational waves from CBC sources. Gravitational-wave
astronomy using compact binaries will ultimately be limited by our ability to robustly detect and estimate the parameters of the sources. Here we focus on two
challenges for detection and parameter estimation in aLIGO. Firstly, we study thepotential to detect IMBHs using aLIGO through the coalescence of a stellar-masscompanion into an IMBH. We refer to this kind of binary as an “intermediate
mass-ratio coalescence” (IMRAC) source. Secondly, we address the computational efficiency of performing parameter estimation on gravitational waves fromCBC sources. Parameter estimation on gravitational waves from CBC sources
is currently highly computationally intensive [25], and the computational costof performing parameter estimation on aLIGO data will be exacerbated if notdirectly addressed.

Detecting gravitational-wave signals from coalescing compact binaries and estimating the source parameters
Despite the scientific potential of gravitational-wave detectors, detection and parameter estimation are extremely difficult in practice. A passing gravitational
wave creates a transverse strain in spacetime which we aim to detect using interferometers. An interferometric gravitational-wave detector in the path of such agravitational-wave thus experiences a relative change in the length of its arms.
This strain is the observable quantity upon which gravitational-wave astronomy

3

1. INTRODUCTION

relies. The exact form of the strain from a gravitational-wave from binary coalescence depends sensitively on the parameters of the binary, e.g. its masses andspins. However, in general a gravitational-wave signal from a coalescing binary
has a characteristic signature. The coalescence of a binary system can be approximately split into three stages; (i) the so-called “inspiral” phase during which thecomponents of the binary evolve along a quasi-circular orbit, “spiralling” into one
another, (ii) a rapid merger where the two objects form a single object, and (iii) a
quasi-normal ringing phase where the (merged) object behaves like a damped oscillator [26]. While it is possible to have inspiralling binaries on elliptical orbits,
binaries observable over the bandwidth of aLIGO/Virgo are expected to havecircularized, see e.g. [27]. These three stages of coalescence imprint a distinct
signature to the gravitational-wave signal. Schematically it is useful to consider
a gravitational-wave signal h(t), as seen in a detector, to have the formh(t) = hinspiral (t)⇥(⌧t)+ hmerger RD (t)⇥(t

⌧) ,

(1.1)

where hinspiral (t) is the part of the waveform from the inspiral phase of coalescenceand hmerger RD (t) is the part of the waveform from the merger and ringdown (RD)phases of coalescence. Here ⇥(.) is the usual Heaviside step function. The time
t = 0 is taken to be the time at a fiducial gravitational-wave frequency, fgw . The
time ⌧ corresponds to the time from fgw until the binary reaches a point close to
the “inner-most stable circular orbit” (ISCO), after which the binary undergoes aplunge before the components merge [28]. Typically, ⌧ is referred to as the “chirptime”. To leading order the chirp time is given by [29]⌧ = 2.18 s!5/3

1.12 MMc

100 Hz
fgw

!8/3

,

(1.2)

where Mc = (m1 m2 )3/5 /(m1 + m2 )1/5 is known as the “chirp mass” of the binary.
For a binary consisting of two neutron stars, each with mass m1 = m2 = 1.4 M ,at 10 Hz, roughly the low frequency cut o↵ of aLIGO, one has access to aroundt ⇡ 20 min of the gravitational-wave signal from the inspiral phase of coalescence.

4

At leading order, the gravitational-wave frequency at which the binary reaches
the ISCO is✓◆
M
fISCO = 4.4 kHz
,
(1.3)
M
where M is the total mass of the binary. Thus a binary consisting of two neutronstars with total mass of 2.8 M will reach ISCO at around 1.5 kHz after which
the components begin to merge.

The inspiral phase of coalescenceTo leading order a coalescing binary is a mass-quadrupole radiator and thegravitational-wave strain from the inspiral phase (at leading order) is controlled
by the second time-derivative of the binary’s quadrupole moment, with a characteristic 1/r fall-o↵ where r is the distance from the source to an observer. Fora binary with chirp mass Mc , angular frequency ⌦, and orbital phase (t), theinspiral phase of the strain signal of a face-on, overhead binary is to leading order
given by [29]hinspiral (t) = hinspiral(t) cos (t) ,
0
!5/3
1Mpc
M
c
hinspiral(t) ⇡ 10 23⌦2/3 ,
0
r
M

(1.4)

It is simple to estimate the typical strain from the inspiral of a compact binary
source using Eq. (1.4). The gravitational-wave frequency is twice the orbital
frequency, hence ⌦ = ⇡ fgw where fgw is the gravitational-wave frequency. Thus
for a binary neutron star system with component masses m1 = m2 = 1.4 M
located at a fiducial distance of 15 Mpc (roughly the distance to the Virgo galaxy
cluster) at a gravitational-wave frequency of 10 Hz, around the low frequency
cut o↵ of Advanced LIGO, the strain amplitude is approximately 10 23 . For aninterferometric gravitational-wave detector with arms of length l = 4 Km, thestrain is l/l ⇡ 10 23 where l is the change in the arm length. Thus the change
in length of the arms which one aims to measure is approximately l ⇡ 10 20 m.
As the binary evolves in time the frequency evolution is driven by the energy
radiated through gravitational waves. The gravitational-wave frequency evolution5

1. INTRODUCTION

Figure 1.1: Inspiral phase of a gravitational-wave signal emitted from the coalescence of a face-on, overhead binary black hole system. The component massesof the binary are m1 = m2 = 5 M . The waveform’s amplitude and frequency
increase up to the merger. This waveform was generated using the LIGO Scientific Collaboration Algorithm Library’s (LAL) implementation of the TaylorT4
waveform approximant [26].

(to leading order) scales like [29]f˙gw /

Mc
M

!5/3

11/3
fgw
.

(1.5)

As the frequency of the gravitational-wave signal increases, both hinspiral(t) and0
f˙gw increase. This endows the inspiral phase of the coalescence signal with acharacteristic “chirp” signature, i.e., a signal which increases in both amplitudeand frequency. An example of the inspiral phase of a gravitational-wave signal
from binary coalescence is shown in Fig. 1.1.

6

The merger and ringdown phases of coalescence
After the inspiral phase the components of the binary merge and undergo quasinormal ringing, with the gravitational-wave strain amplitude decreasing exponentially over a short time scale. The characteristic time scale is given by [30]
T =

2 1
,
⇡ f0
(1.6)

f0 = 10 kHz

M
.
M

(1.7)

To leading order the ringdown waveform (for a binary with non-spinning components) has the form [30]hringdown (t) ⇡ 10 20

1 Mpc Me t/T cos(2⇡f0 t) .
r
M

(1.8)

For a double neutron star binary with component masses m1 = m2 = 1.4 M theamplitude decreases by an e-fold after time t ⇠ 5 ⇥ 10 5 s.
Detecting the strain signal from a coalescing binary
Detecting gravitational waves thus requires measuring a signal whose amplitudeis physically very small. Noise sources in the detectors typically cause strainsin the instrument at the frequency of the gravitational-wave strain signal, whichmakes the task of detecting the strain from a gravitational wave challenging. The
dominant noise sources for aLIGO include thermal (or Brownian) noise, quantum
noise from the uncertainty in the photon count in the lasers (also called shot noise)and seismic noise from the Earth’s seismic activity [2]. The strain spectrum
p(strain/ Hz) of the various noise sources of aLIGO is shown in Fig. 1.2. Eachnoise source typically acts over a certain frequency range which is illustrated inFig. 1.2. For example, thermal noise in the suspensions and mirrors is present
over a frequency range between O(10 Hz) and O(1 kHz). Conversely, seismic noise
typically dominates the noise below frequencies of O(10 Hz) which e↵ectively
sets the low frequency cut o↵ of the instruments. The combined noise (shownas the black curve) is known as the detector noise amplitude spectral density7

1. INTRODUCTION

Figure 1.2: Anticipated Advanced LIGO noise sources generated using the online
GWINC tool [31]. The combined noise (shown as the black curve) is known
as the detector noise amplitude spectral density (ASD) and defines the overallsensitivity of the instrument. The noise ASD corresponds to the instrument in
its high-powered zero-detuned configuration.

(ASD) and defines the overall sensitivity of the instrument. The peak strainsensitivity of aLIGO can be read o↵ from Fig. 1.2 and is around 10 23 Hz 1/2
over a frequency range 10 Hz . f . 103 Hz. In Fig. 1.3 we show the noise
amplitude spectral densities for various Earth based gravitational-wave detectors.
From Fig. 1.3 one can see the improvement in sensitivity between LIGO/Virgo,
and their advanced configurations aLIGO/Virgo, which are around an order ofmagnitude more sensitive.
The problem of detecting signals in noise is a classic of signal analysis. The
noise properties are well understood and over time scales around the typicallength of a gravitational-wave signal it is e↵ectively stationary [1; 29]. In addition, for many practical purposes it is often assumed that the noise will not becorrelated with a gravitational-wave signal, for example, by considering it to be
Gaussian. However, in practice real noise is not Gaussian and a large e↵ort isspent on characterizing the non-Gaussian features of the detector noise (see e.g.,[33; 34; 35]). Secondly, the functional form of gravitational-wave signals from8Figure 1.3: Noise amplitude spectral densities of various Earth based
gravitational-wave detectors [32]. Note the increase in sensitivity by around anorder of magnitude between LIGO and Virgo, and Advanced LIGO and Advanced
Virgo.

9

1. INTRODUCTION

CBC sources can be modelled accurately [26], and one can filter the data for a
signal of known functional form. This is known as “matched filtering” [5; 6; 7; 8]
and is e↵ective for detecting gravitational waves from CBC sources precisely because the gravitational-wave signal can be accurately modelled. In the presenceof a gravitational-wave signal h(t), the data d(t) has the formd(t) = n(t) + h(t) ,(1.9)

where n(t) is the noise. To detect gravitational waves using matched filtering
one compares the data with a set of theoretical gravitational waveforms {T (t)}(known as templates), normalized by the square of the detector’s ASD, which act
as filters of the data. The filter which maximizes the detection-statistic is thenchosen to be the candidate signal. A commonly used detection statistic is thesignal-to-noise ratio (SNR) which is a ratio of the power of the gravitational-wave
signal to that of the noise. The SNR of the data set d in Eq. (1.9) is [29]SNR = 4<1

Z 1

˜ )/Sn (f ) ,df T̃ ⇤ (f ) d(f(1.10)

0where T̃ (f ) is the Fourier-transform of the template T (t) (and similarly for thedata) and Sn (f ) is the detector power spectral density (PSD), which is the square
hR
i1/2
1
of the ASD. The normalization is given by = 0 df T̃ ⇤ (f ) T̃ (f )/Sn (f ).
The optimal SNR is achieved when the template waveform is same as the signal
in the data set and hence one requires the templates to be a faithful description of real gravitational waveforms. Using templates which are not accurate
descriptions of gravitational waves contained in the data degrades the SNR, as
the optimal SNR is recovered only when the template corresponds exactly to thesignal. Using inaccurate templates also leads to losses in detection rates which
can be understood as follows. The gravitational-wave signal from binary coalescence has a 1/r fall-o↵, c.f. Eq. (1.4), where r is the distance from the Earth
to the source. Using a template which degrades the SNR by a fixed fractionalammount ✏ can be compensated for by increasing the overall amplitude of thetemplate. This is achieved by reducing the distance to the source by the sameamount as the SNR loss. Reducing the distance to the source by a factor of ✏10thus reduces the observable volume to a factor of ✏3 of its original size. Assuming
an isotropic distribution of sources, the fractional loss in detection rates should
thus scale like (1 ✏3 )One of the key considerations for a gravitational-wave search is the templates
one uses to filter the data. For example, searches which target so called “lowmass” systems with total masses less than 25 M [7] typically only require waveforms which describe the inspiral phase of coalescence. This can be understoodby considering the frequency at which the binaries merge; a system with total
mass of 5 M will merge at around 1 kHz, c.f. Eq. (2.1). For initial LIGO thisinspiral phase sits in the frequency band over which the detector is most sensitiveand hence the inspiral phase can dominate the SNR, c.f. Fig. 1.3. Conversely,
heavier systems are liable to merge at frequencies at which the detector is mostsensitive and so searches for these systems need to employ waveforms which also
describe the merger and ringdown. The so-called “high mass” searches in initial
LIGO [8] which aimed to detect compact binaries with masses up to 100 M employed such waveforms. Systems with even higher masses may have a very short
detectable inspiral phase, or may have even merged at the point at which they
are detectable. For such systems, the ringdown can be the dominant source of the
SNR and hence one may require waveforms which describe solely the ringdown.
The “ringdown search” during initial LIGO [30] employed such ringdown-only
waveforms to search for binaries with total masses up to 390 M .

Extracting the binary’s parameters
The gravitational waves encode the parameters of the source and for a binary
system with non-spinning components there are nine such parameters: the twocomponent masses, the distance to the Earth, the two angles which specify thesky position, the binary’s inclination, the “polarization phase”, the time of coalescence and the orbital phase at coalescence. Spinning binaries also have thetwo spin-vectors of each component which increases the number of parameters to
fifteen. Because the filters are parameterized, the filter which maximizes the SNR
can be used to estimate the parameters of the source. This is thes maximum likelihood estimator of the source parameters [29] and in this fashion one can extract11

1. INTRODUCTION

information about the source of a candidate gravitational-wave signal.
A limitation of these kinds of parameter estimates is that one only gets a“point estimate” of the parameters of the gravitational-wave source, the maximum likelihood estimator. A more informative approach would also provide aquantitative picture of the uncertainty of parameter estimates over the full domain of the parameter space. For this one subjects a stretch of data to a full
Bayesian analysis and computes the posterior probability distributions of thesource parameters of a candidate gravitational-wave signal, see e.g. [12; 36; 37].
This approach allows one to assign meaningful probability distributions to themeasured properties of the candidate source. The Bayesian approach also allows
for “model selection”, by comparing the null hypothesis that the data contains
only noise, to the hypothesis that the data also contains a signal. This is theratio of the probability that a data set is described by just noise, to the probability that the data contains a gravitational-wave signal. This is commonly known
as the “Bayes factor” [12] and provides a useful way of quantifying the relative
probability that the data contains a gravitational-wave signal.

Challenges of waveform modelling for detection and parameter estimation
There are a number of theoretical challenges associated with gravitational-wave
detection and (Bayesian) parameter estimation. Firstly, in order to robustly
detect and make parameter estimates of gravitational-wave sources, one requires
templates which are accurate descriptions of real gravitational waves. While there
are several schemes in which gravitational waveforms can be modelled, such as thepost-Newtonian (PN) perturbative expansion and numerical solutions to the Einstein field equations (numerical relativity) [26], there is not yet a complete analytic
description of a gravitational waveform for binary systems with arbitrary mass
ratios and spins. For example, the PN expansion is valid when the binary system
has components with comparable masses and the expansion is itself only formallyvalid up to a point before the binary’s components merge, when the expansion
parameter, v/c becomes large. For orbital velocities v/c = (M ⇡f )1/3 & 0.2 thePN energy flux deviates significantly from more accurate numerical expressions,12see [e.g., 38; 39]. Hence waveforms computed using the PN expansion only describe the “inspiral” of the coalescence, ignoring the merger and the subsequentquasi-normal ringing of the merged objects. In contrast, numerical relativity is anapproach which solves the Einstein field equations numerically and can describe
the full coalescence. However, numerical relativity is extremely computationallyexpensive, with a single simulation taking between several hundred, to thousandsof CPU-hours [40].
The limitations of both approaches owe to the difficulty of the two-body problem in general relativity. Because of this there has been a strong e↵ort to produce
waveforms which are a hybrid of analytic and numerical approaches, stitching
together analytic fits and numerical expressions to produce an approximation
scheme which can describe the full coalescence. An example of such a scheme
which is commonly employed for gravitational-wave detection is the e↵ective onebody numerical relativity (EOBNR) [26; 41] approach. Despite its success, it is
only known to be accurate for binaries with mass ratios q ⌘ m2 /m1 1/6 though
it is also formally valid in the test-particle limit q ! 0. For binaries with extreme
mass-ratios (q ⌧ 1), consisting of a stellar-mass black hole and a super-massive
black hole, several authors ([e.g., 42; 43]) have modelled the two body motion
as that of a point particle in Kerr/Schwarzschild space-time, with corrections toaccount for the gravitational radiation-reaction from solutions to the Teukolsky
equation. Like the PN approach, this also only describes the inspiral phase of the
coalescence. More recently, the authors in [44] have added high-order-in-massratio corrections to extreme mass-ratio inspiral waveforms so that they describe
binaries with so called “intermediate mass-ratios”. For the work here we consider
an intermediate mass-ratio q to be in the range 1/200  q  1/10. However, todate only one numerical simulation of an intermediate mass-ratio binary exists[45] and analytic approximations to describe intermediate mass-ratio inspirals areat an early stage of development. This renders the problem of robustly detecting
sources which consist of an intermediate-mass black hole and a small stellar masscompanion difficult. For matched filtering to be successful, the template must
be able to track the phasing of the gravitational wave in the data to within awave-cycle, over many wave cycles. Errors in the phasing of the filters can renderthem ine↵ective and can thus lead to losses in detection rates.

13

1. INTRODUCTION

Template waveform generation can be a computational burden to performing
matched-filtered searches or Bayesian parameter estimation. In order to perform
a matched-filtered search for gravitational waves, or a full Bayesian parameter
estimation, one requires a set of template waveforms which have to be computed
numerically, known as a “template bank”. Templates generated using the PN orEOBNR approach are computed by first solving the underlying two-body equations of motion from which the gravitational waves are extracted. This is somewhat undesirable for gravitational-wave data analysis because of the sheer numberof template waveforms required to construct a template bank for detection. In
order to e↵ectively search for gravitational waves, one typically constructs a dense
template bank on the intrinsic parameter-space (the two masses and spins which
control the phase evolution) to ensure that one does not lose a significant amountof SNR through under-dense gridding [46]. The typical number of template waveforms required for detection can be on the order of hundreds of thousands [47]. In
addition, the accuracy requirements of Bayesian parameter estimation are somewhat stricter than for detection, and a standard analysis requires computing of
millions of template waveforms [25]. Whereas in a gravitational-wave search one
performs a maximization procedure to find the maximum-likelihood estimator, in
Bayesian parameter estimation one is interested in computing the full posterior~ H)
probability density function (PDF) of the source parameters. The PDF p(✓|d,of a set of parameters ✓~ which parameterize a model, H, assumed to describe~ H), the a priori probability disa data set d, is related to the likelihood p(d|✓,~tribution of the source parameters p(✓|H)
and the “evidence” p(d|H) via Bayes’
theorem [48]:~~~ H) = p(✓|H) p(d|✓, H) .
p(✓|d,(1.11)
p(d|H)As discussed earlier, the parameter vector ✓~ for coalescing binaries is nine dimensional for binaries with non-spinning components and fifteen dimensional forbinaries with spinning components. It is immediately clear that a fine griddingin each dimension is not feasible and so one must make use of algorithms topreferentially sample regions of the parameter space where the bulk of the posterior probability is contained. There are currently several stochastic samplingalgorithms which are dedicated to parameter estimation on gravitational waves

14

from CBC sources. These include Markov-chain Monte Carlos [13] and nested
sampling [12] algorithms. A typical stochastic Bayesian analysis can take up to
several weeks to run on a single stretch of data, around ten seconds in duration.
The dominant cost of parameter estimation is waveform generation, and the costof the analysis scales linearly with the cost of waveform generation.

Methods to improve the efficiency of detection and parameter estimation are
therefore crucial for gravitational-wave astronomy. On the side of detection it has
recently been shown that there is a large redundancy in the number of filters used
for gravitational-wave searches [49] and the authors employed a singular value decomposition (SVD) of a template bank of waveforms used for a gravitational-wave
search to show that the total number can be reduced by around an order of magnitude without incurring significant losses in SNR. In addition, the decompositionof the (discrete) template bank can be used to generate a continuous descriptionof the bank (e↵ectively an interpolation technique) [50]. The interpolation bypasses the need to perform ODE solving to generate template waveforms andcould o↵er improvements to the efficiency of detection and parameter estimation
when the computational cost is dominated by template waveform generation. It
is all the more crucial to address the issue of the computational cost of template
generation in the run-up to the advanced detector era as the in-band signal duration will be longer than in the initial configurations of the instruments. This willrender computing templates more computationally expensive if not directly dealt
with. For example, a binary neutron star system with total mass of 2.8 M at a
gravitational-wave frequency of 40 Hz, roughly the low frequency cut o↵ of initial
LIGO, would have an in band signal duration of around 25 s (c.f. Eq. (1.2)).
Conversely, the same system at a gravitational-wave frequency of 10 Hz, roughly
the low frequency cut o↵ of aLIGO, has an in band signal duration of around 20
min. Computing hundreds of thousands, or millions of these waveforms will render detection and parameter estimation an even more computationally expensive
procedure.

15

1. INTRODUCTION

Overview of the thesis
Chapter II
Chapter II focuses on the detection of intermediate-mass black holes in aLIGO.
As a detection channel we consider the coalescence of a stellar-mass black hole
into an intermediate-mass black hole (IMBH) which we refer to as an intermediate
mass-ratio coalescence (IMRAC). We consider binaries in which the IMBH has
mass in the range 24 M 200 M with a stellar-mass companion having masses in
the range 1.4 M
18.5 M , and with mass ratios 1/140  q  1/10. This range
is well within the detectable mass-range of aLIGO. We investigate the relative
contribution to the signal-to-noise ratio (SNR) of the three di↵erent phases of
the coalescence - inspiral, merger and ringdown - of IMRACs as would be seen
by aLIGO. We find that merger and ringdown play the dominant role over thebulk of the mass space though there is a small portion in which the inspiralportion of the coalescence signal dominates. We identify three regions in the
mass-space in which (i) inspiral-only searches could be performed with losses indetection rates, L, in the range 10% . L . 27%, (ii) searches based on inspiralonly templates lead to a loss in detection rates in the range 27% . L . 50%, and(iii) templates that include merger and ringdown are essential to prevent losses indetection rates greater than 50%. In addition, we provide a comparison of existing
inspiral-only waveforms over the mass range considered in the study. We consider
waveforms based on the EOBNR and PN approximation schemes, and a waveformfamily developed specifically to describe intermediate mass-ratio inspirals [44]which we refer to as the “Huerta-Gair” waveform family after its authors. We
measure the consistency of waveforms in di↵erent families by computing theirnormalized inner-products, maximized over the parameters which control thephase evolution of the waveforms. We find that the di↵erent families are 90%
consistent with each other which reinforces the importance of several numerical
relativity simulations in this mass-range with which semi-analytical models of
waveforms could be constructed in order to produce robust templates for detecting
IMRACs. This study is restricted to binaries on (quasi-) circular orbits with nonspinning components.
The text of Chapter II is based on a paper by Rory Smith, Ilya Mandel and16

Alberto Vecchio [19]. My contribution to this work was: (i) the implementation of
the Huerta-Gair waveform family, (ii) implementation of all the calculations (iii)organized and led the project (iv) wrote the paper. The paper benefited fromdiscussions with Eliu Huerta. This paper has arXiv number arXiv:1302.6049[astro-ph.HE].

Chapter III
Chapter III is the first of three chapters on techniques to improve the efficiencyof Bayesian parameter estimation. In Chapter III we provide a proof of principle of a technique which could be used to improve the efficiency of Bayesianparameter estimation on CBC sources when the computational cost of the analysis is dominated by template waveform generation. The work employs a singular value decomposition (SVD) of a template bank of waveforms which we use
to interpolate templates for parameter estimation. The interpolated templates
have already been employed in gravitational-wave searches [49; 50]. Interpolating template waveforms can be a less intensive computational procedure than
generating templates via di↵erential equation solving, for example. Because theaccuracy requirements of gravitational-wave searches are typically less stringent
than for parameter estimation we (a) study the bias in posterior parameter estimates incurred through using interpolated template waveforms, and (b) computethe relative computational cost of computing interpolated templates and noninterpolated templates. We find that the bias in mean posterior parameter estimates is negligible on the scale of statistical measurement uncertainty and thatinterpolating templates is around an order of magnitude faster than generating
standard time-domain waveforms, e.g., TaylorT4 and EOBNR [26] waveforms.
Our study uses a simple class of interpolated inspiral-only template waveforms
and is restricted to a small patch of the parameter space. We consider future workrequired in order for interpolated template waveforms to be integrated within parameter estimation pipelines.
Chapter III grew out of a collaboration between Rory Smith, Kipp Cannon,
Chad Hanna, Drew Keppel and Ilya Mandel while on a fellowship at the Perimeter
Institute for Theoretical Physics. The text of Chapter III is based on a paper by17

1. INTRODUCTION

Rory Smith, Kipp Cannon, Chad Hanna, Drew Keppel and Ilya Mandel [51]. My
contribution to this work was: (i) writing the SVD-interpolation software, (ii) integrating the SVD-interpolation software into the LSC Algorithm Library (LAL),(iii) the comparison of PDFs generated using interpolated and non-interpolated
waveforms in Sec. 3.5, (iv) the comparison of computational time of waveformfamilies in Sec. 3.6, (v) organized and led the project, (vi) wrote the paper. This
paper has arXiv number arXiv:1211.1254 [astro-ph.HE].
Chapter IV
Chapter IV builds on the previous chapter and represents continuing work on
improving the efficiency of Bayesian parameter estimation. Here we extend theinterpolation used in Chapter III to more general classes of template waveforms;
inspiral-only time-domain waveforms, and time-domain waveforms which describe
the full coalescence signal. We find that our implementation of the interpolation
is successful for inspiral-only time-domain waveforms but fails for time-domain
waveforms which describe the full coalescence.
All the work in Chapter IV is my own. It grew out of the work in Chapter III
which was a collaboration with Kipp Cannon, Chad Hanna, Drew Keppel andIlya Mandel, and benefited from discussions with the afore mentioned people.
Chapter V
Chapter V builds on the previous two chapters on improving the efficiency ofBayesian parameter estimation. While the previous two chapters focussed on interpolating template waveforms, here we show how to efficiently compute the continuous likelihood function on the three-dimensional subspace of parameters on
which it has a non-trivial dependence – the chirp mass, symmetric mass ratio andcoalescence time – via interpolation. We show that subsequently, sampling thisinterpolated likelihood function is a significantly cheaper computational process
than directly evaluating the likelihood; we report improvements in computationaltime of between two and three orders of magnitude while keeping likelihoods accurate to . 0.025 %. Generating the interpolant of the likelihood function overa significant portion of the CBC mass space is computationally expensive but

18highly parallelizable. This suggests that the overall wall time of generating theinterpolant of the likelihood function can in principle be very small compared tothe time of a full parameter estimation analysis.
The work in Chapter V is based on a paper by Rory Smith, Chad Hanna, Ilya
Mandel and Alberto Vecchio. This work benefited from discussions with several
people, in particular Kipp Cannon, Scott Field, Drew Keppel, Vivien Raymond,
Manuel Tiglio and Alan Weinstein. This paper has arXiv number arXiv:1305.3798[astro-ph.HE].

19

1. INTRODUCTION

20Chapter 2
DETECTING INTERMEDIATE
MASS-RATIO BINARIESWITH ADVANCEDDETECTORS
2.1

Introduction

Observations of ultra-luminous X-ray sources and simulations of globular cluster
dynamics suggest the existence of intermediate-mass black holes (IMBHs) [5;
20; 52; 53; 54; 55]. However, observational evidence for their existence is still
under debate, see e.g. [56; 57]. Gravitational waves from binary coalescences
involving IMBHs with masses 50 M . M . 500 M are potentially detectableby advanced detectors – including Advanced LIGO [58], Advanced Virgo [59], and
KAGRA [60] – with a low frequency cuto↵ of around 10 Hz. If IMBHs do exist,
one likely contribution to gravitational-wave detections is believed to be throughthe coalescence of a compact stellar-mass companion (black hole or neutron star)with an IMBH, at a possible rate of up to ⇠10 yr 1 [5; 14; 61], with masses of
the IMBH in the range 50 M . M . 350 M [14]. We will denote these signals

21

2. DETECTING INTERMEDIATE MASS-RATIO BINARIESWITH ADVANCED DETECTORS

as intermediate mass-ratio coalescences (IMRACs)1 .
Given that IMBHs in this mass range have proved extremely difficult to observe in the electromagnetic spectrum, gravitational-wave detections may provide
the first unambiguous observations of such objects. Such observations would form
an important channel for probing the dynamical history of globular clusters. Furthermore, Advanced LIGO/Virgo (aLIGO/Virgo) may be able to provide measurements of the quadrupole moment of a black hole [14; 15], which would allow
a null-hypothesis test of the Kerr metric for IMBHs.
The gravitational waveform from the coalescence of two compact objects canbe divided into three phases: a gradual inspiral, a rapid merger, and the quasinormal ringdown of the resulting black hole. The relative contribution to the
expected coalescence signal from inspiral, merger and ringdown is an importantconsideration for gravitational-wave searches. To leading Newtonian order thegravitational wave frequency at the inner-most stable circular orbit (ISCO) is
fISCO = 4.4 kHz

✓

M
M

◆

,

(2.1)

where M is the total mass of the binary. For advanced detectors with a low
frequency cut-o↵ of ⇠ 10 Hz, we may only have access to either the very latestages of the inspiral, or solely merger and ringdown for the heaviest IMRAC
systems. While the power in the merger and ringdown is suppressed by a factorof the mass ratio relative to the power in the inspiral [62], the fact that IMRAC systems are liable to merge either in-band, or at the low frequency limitof the bandwidth, means that merger and ringdown may be significant over alarge portion of the detectable mass-space. Additionally, for cases where IMRAC waveforms are inspiral-dominated, it is useful to know where inspiral-only
searches could be targeted.
Detecting IMRACs through gravitational waves will require template gravitationalwaveform families adapted to highly asymmetrical mass-ratio systems. However
the development of numerical relativity simulations and perturbative techniques
1

In the literature, the term frequently used for this class of objects is intermediate massratio inspirals or IMRIs, see e.g. [14; 61]. However, in the context of ground-based observations,in particular with second-generation instruments, we will show that the full coalescence isimportant for these systems, and it therefore seems to be more appropriate to call them IMRAC.

22

in this regime is at an early stage which is potentially problematic. The issue
of appropriate template waveform families is thus central to the detection of IMRACs through gravitational waves.
The e↵ective-one-body approach, calibrated to numerical relativity, has led
to template waveforms, known as EOBNR [41], that describe the full inspiral, merger and ringdown coalescence-signal for comparable mass-ratio binaries;
EOBNR waveforms should also be accurate at extreme mass ratios. However, todate only one full numerical simulation exists for mass-ratio q = 1/100 binaries[45]. EOBNR waveforms have not yet been compared to NR simulations at suchmass ratios, so their validity in the IMRAC regime remains to be demonstrated.
Meanwhile, in the context of extreme mass-ratio binaries, several authors havemodelled the two-body motion by computing radiative and conservative self-force
corrections to Kerr geodesic motion [42; 43]. Waveforms computed within thisscheme are inspiral-only and are only developed to lowest order in the mass ratio.
These waveforms have been adapted to describe intermediate mass-ratio inspirals
by including higher-order-in-mass-ratio corrections in [63] and have been used to
study the detection of intermediate mass-ratio inspirals in the Einstein Telescope[44]. We refer to these intermediate mass-ratio inspiral waveforms as the “Huerta
Gair” (HG) waveform family after its authors. This waveform family should bephysically well motivated to describe the inspiral of IMRACs.
Typically one does not have an exact representation of “true” gravitationalwave signals but requires templates which are sufficiently e↵ective at filteringsuch signals. A common metric for quantifying how well approximate waveformfamilies are at filtering gravitational-wave signals is known the “e↵ectiveness”,
or fitting factor [26]. This measures the fraction of the theoretical maximum
signal-to-noise ratio (SNR) that could be recovered by using non-exact template
waveforms.
The work in this Chapter focuses on the requirements of template waveforms
for the detection of IMRACs in aLIGO. Firstly, by computing the e↵ectivenessof inspiral-only template waveforms at filtering the full coalescence signal, we
determine the relative importance of the inspiral and merger-ringdown phases.
We identify three regions in the m1 -m2 plane in which: (a) inspiral-only searches
are feasible with losses in detection rates L in the range 10% . L . 27%, (b)232. DETECTING INTERMEDIATE MASS-RATIO BINARIESWITH ADVANCED DETECTORS

searches are limited by the lack of merger and ringdown in template waveforms
and are liable to incur losses in detection rates in the range 27% . L . 50%,and (c) where merger and ringdown are essential for searches in order to prevent
losses in detection rates greater than 50%.
Secondly, to gain insight into the accuracy of the inspiral portion of IMRAC
waveforms we compute the e↵ectiveness of the inspiral-only portion of EOBNR
waveforms at filtering gravitational-wave signals as described by the HG waveform family. We find that there is a non-negligible discrepancy between EOBNR
and HG inspirals in the regime where inspiral-only searches could be consideredsufficient. For reference we also compare EOBNR inspirals to a post-Newtonian(PN) [64; 65] waveform family known as TaylorT4 [26]. The PN expansion is liable to be a poor choice of approximant for IMRACs because of the large numberof cycles spent at small radii. We find that EOBNR and HG are in better agreement with each other than to TaylorT4, as might be expected from the previous
observation.
Our approach does not directly address the accuracy of template waveforms,
because none of the waveforms considered have been matched to full numerical
waveforms. However, assuming that the waveform families we consider “bracket”
the correct gravitational waveforms in the IMRAC regime, this approach provides a useful heuristic for quantifying the e↵ectiveness of existing gravitational
waveforms for IMRAC searches. Further numerical relativity simulations will
be important in the continuing development of accurate template waveforms for
IMRACs.
Our analysis improves upon previous work to determine the detectability of
IMRAC sources [66] which only considered the faithfulness of template waveforms,
i.e., the e↵ectiveness of template waveforms evaluated at the signal parameters.
Additionally, that study only considered inspiral-only waveforms and focused on
the LISA regime.
This chapter is organized as follows. In Sec. 2.2 we review the basics ofmatched filtering as applied to gravitational-wave detection. In Sec. 2.3 we describe the waveform families used in our study. In Sec. 2.4 we compute thecontributions to signal-to-noise ratio (SNR) from the inspiral and merger andringdown phases of EOBNR waveforms in the intermediate mass-ratio regime. In

24

Sec. 2.5 we study the e↵ectiveness of inspiral-only waveforms to filter inspiralmerger-ringdown gravitational wave signals from IMRAC sources and identify
the three regions in which di↵erent searches could be conducted. In Sec. 2.6 wecompare the inspiral portion of EOBNR waveforms to HG and TaylorT4 waveforms. In Sec 2.7 we consider the implications of our results for future searches
in advanced detectors.

2.2

Matched filtering of gravitational-wave signals

In order to quantify the requirements of gravitational waveforms for IMRAC
searches, we summarize the key concepts of matched filtering which is the core
principle behind gravitational-wave detection. The fundamental object is the
optimal filter of the data which is used to define the maximum signal-to-noise
ratio achievable from a stretch of data which contains a gravitational-wave signal.
The output of a gravitational-wave detector which contains a gravitationalwave signal will generally be of the form s(t) = h(t) + n(t). Here s(t) is adata-stream consisting of a gravitational-wave signal h(t) and noise n(t). We
take the noise to be Gaussian and stationary for simplicity. We note that ingeneral the noise will not be Gaussian and stationary in the detector network,
though will not consider the e↵ects of non-Gaussianity in the noise here.
One can filter the data for a signal h(t) if we know the form of h(t). Consider
the following filtering procedure. We take a signal containing a rapidly oscillating
function of time, h(t), to be the signal buried in Gaussian and stationary noise.
The convolution of the detector output, s(t), with the signal, averaged over anobservation time T , is
1
T

Z T
0

1
dt s(t)h(t) =
T

Z T
0

1
dt h (t) +
T
2

Z T

dt n(t)h(t).

(2.2)

0

Both h and n are separately oscillating functions of time. The first term on theright hand side of Eq. (2.2) is positive definite and grows, for large values of T , as
T . Averaged over the observation time T , its value will be the average amplitude25

2. DETECTING INTERMEDIATE MASS-RATIO BINARIESWITH ADVANCED DETECTORS
RTof the signal, 1/T 0 dt h2 (t) ⇠ h20 , where h0 is the average amplitude of h(t).
Conversely, the noise n(t) and the signal h(t) are uncorrelated. The second termon the right hand side of Eq. (2.2) grows at most, for large T , as T 1/2 so thatRT
1/T 0 dt n(t)h(t) ⇠ (⌧ /T )1/2 n0 h0 , where n0 is the characteristic amplitudeof the noise and ⌧ is its characteristic period of the noise. Hence in the limit
T ! 1, the second term on the right hand side of Eq. (2.2) averages to zero andthe signal is filtered out of the data-stream.

More formally, we work with the filtered detector output ŝ given by
ŝ =

Z 1

dt s(t)K(t),

(2.3)

1

where s(t) is a data stream and K(t) is the filter function. A commonly used
detection statistic is the “signal-to-noise” ratio (SNR), defined as the ratio of the
expected value of ŝ when a signal is present, S, to the rms value of ŝ when thesignal is absent, N. We thus seek the filter K(t) which maximises the SNR. The
expectation value of the filtered signal, S, is [29]ZS =dt hs(t)i K(t) ,Z=

dt h(t) K(t) ,Z=

df h̃⇤ (f )K̃(f ) ,(2.4)

where we have used the fact that hn(t)i = 0 on the second line. The Fouriertransform, h̃(f ) of h(t), and its inverse are defined in Eq. (A.1). The rms valueof ŝ when the signal is absent, N can be computed from [29]
N

2

=

h

2

hs (t)ihs(t)i

2

i

h=0

,

= hs2 (t)ih=0 ,Z
1=
df Sn (f )|K̃(f )|2 ,
2

(2.5)

where the quantity in the Sn (f ) is known as the detector “power spectral density”

26

(PSD), and is defined by the ensemble average of the noise [67]:11
Sn (f ) = h|ñ(f )|2 i ,
2
T

(2.6)

where T is the duration of an observation time over which the noise would typically be collected in an experiment. The SNR is thus given by [67]
S
= hR
N
1

R1

1

1

df h̃(f )K̃ ⇤ (f )df (1/2)Sn(f )|K̃(f )|2

i1/2 .

(2.7)

To find the filter K which maximizes Eq. (2.7) we define the “noise-weighted
inner-product”:
Z 1ã⇤ (f )b̃(f )(a|b) = 4 <
df
,
(2.8)
Sn (f )
0in terms of which the SNR, Eq. (2.7), is given by
S
(u|h)
=
,
N
(u|u)1/2(2.9)

where u(t) is a function with Fourier transform ũ(f ) = (1/2) Sn (f )K̃(f ). Geometrically, we have a vector of unit norm û = u/(u|u)1/2 with which we are taking
the inner-product with h. Thus, the SNR is maximized when û is parallel to h.
The optimal filter, K, which maximizes the SNR is therefore given byK̃(f ) = const.

h̃(f ).
Sn (f )

(2.10)

The constant is irrelevant and the best filter is thus given by the signal itself. Thisis a crucial point as it states that we must know the form of the signal which wewant to detect in order to obtain the optimal SNR. In terms of the noise-weighted
inner-product the optimal SNR is expressed succinctly as
⇣S⌘
N

max

= (h|h)1/2 .

(2.11)

Central to the task of detecting gravitational waves from IMRACs is the choice27

2. DETECTING INTERMEDIATE MASS-RATIO BINARIESWITH ADVANCED DETECTORS

of model waveform h which we require to filter the data. In the section below weoutline three model waveform families which we use in our study.

2.3

Waveforms

In this section we summarise the key concepts entering the construction of the
waveforms used in this study. Throughout the chapter, for a binary system withindividual component masses m1 and m2 (with m2 < m1 ) we define the total
mass as M ⌘ m1 + m2 , and mass ratio and symmetric mass ratio as q ⌘ m2 /m1
and ⌘ ⌘ m1 m2 /(m1 + m2 )2 , respectively.
We first consider the family of waveforms constructed by calibrating thee↵ective-one-body approach to numerical relativity (EOBNR) [41; 68; 69]. The
EOBNR family describes the full inspiral-merger-ringdown signal; it is currently
used in searches that reach the IMBH mass range, so far up to 100 M [8]. The
free coefficients in the family have been fitted to comparable mass ratio numerical
relativity simulations, and by construction this family is deemed to be faithful inthe test particle limit. For this work, we use the implementation provided by the
LIGO Scientific Collaboration Algorithm Library (LAL) that corresponds to theapproximant EOBNRv2 [41].

We also consider a waveform family based on test particle motion in Kerr
/ Schwarzschild space-time with radiative and conservative self-force corrections
which we refer to as the Huerta-Gair (HG) family [44]. The approximation scheme
is constructed specifically to handle highly-asymmetrical mass-ratio binaries andis therefore a physically well motivated approximation scheme for intermediate
mass-ratio inspirals. These waveforms have been compared against Teukolsky
based waveforms and the match exceeds 95% over a large portion of the parameter space [63]. These waveforms have been used to study detection of intermediate mass-ratio inspirals in Einstein telescope in [44]. The Huerta-Gairwaveforms describe only the inspiral portion of the coalescence signal. There isno corresponding LAL approximant for the HG waveform family.
Finally, as a reference we use the standard inspiral-only post-Newtonian approximation corresponding to the LAL approximant TaylorT4 at three-and-a-halfpost-Newtonian order in phase [26].

28The TaylorT4 waveforms used here are computed in the so-called “restricted”
amplitude approximation, which assumes the waveform amplitude to be zerothpost-Newtonian order and only includes the leading second harmonic of the orbital phase. We only include the leading second harmonic of the orbital phase in
the EOBNR waveforms. We do not consider the e↵ects of spin or eccentricity in
any of the waveform families, as we restrict to circular orbits and non-spinningblack holes.

2.3.1

EOBNR

EOBNR uses an e↵ective one body (EOB) Hamiltonian to evolve the system’s
trajectory through the inspiral and plunge phase, while numerical simulations andblack hole perturbation theory are used to construct the merger and ringdown
phases respectively. The PN Hamiltonian generates a waveform hinsp plunge (t)which is matched at some time, tmatch , to a merger and ringdown waveform,
hmerger RD . The full EOBNR waveform is given by (see [26] and references
therein),h(t) = hinsp plunge (t)⇥(tmatch+ hmerger RD (t)⇥(tt)

tmatch ),(2.12)

where ⇥(.) is the Heaviside step function.
2.3.1.1

Inspiral-plunge waveform

To calculate the inspiral trajectory one solves Hamilton’s equations given the
Hamiltonian of the EOB metric and the inspiral-plunge waveform at leading
Newtonian-order is [26],h(t)insp plunge = A !
ˆ 1/3 sin [2 (t)] ,(2.13)

where the “angular frequency” !
ˆ is given by Eq. (3.23b) in [26]. The overallscaling A depends on the source location with respect to the detector. In this29

2. DETECTING INTERMEDIATE MASS-RATIO BINARIESWITH ADVANCED DETECTORS

chapter we will use normalized template waveforms and so we leave this factor
unspecified.
2.3.1.2

Merger-ringdown waveform

The merger-ringdown waveform is made up of a superposition of quasi-normal
modes (QNMs) [41; 70; 71]. We only consider the dominant (2, 2) mode and themerger-ringdown waveform is given by
hmerger RD (t) = A

N
X1

22

A22n e i 22n (t tmatch ) ,(2.14)

n=0

where n is the overtone number of the Kerr QNM, N is the number of overtones,t22
match is the time at which the (2, 2) mode is matched to the inspiral-plunge
waveform and the coefficient A22n is found by a matching procedure. We takei↵n where !n > 0 are the oscillation frequencies and ↵n > 0 are asso22n ⌘ !nciated with the inverse of the decay times of the QNM. The full inspiral-mergerringdown waveform is then given by matching the inspiral- plunge waveform, Eq.
(2.13), to the merger-ringdown waveform, Eq. (2.14), waveform close to the EOB
light-ring [72].

2.3.2

HG waveforms

HG waveforms aim to capture the main features of the binary inspiral when the
mass-ratio of the system is small. Physically, this means we work in a test particle
limit. For non-spinning bodies, this limit guarantees that the inspiralling particle
instantaneously follows Schwarzschild geodesics while the flux of the orbital parameters is adiabatic. The orbital parameters (E, Lz ) are evolved slowly using aprescription for the flux of energy, E, and angular momentum projected along thecentral black holes spin axis, Lz . Combining these fluxes with the geodesic equations for a test particle in Schwarzschild space-time yields the particles trajectoryin Schwarzschild coordinates, (t, p, ✓, ). Here we specialize to (quasi-)circularorbits around a non-spinning central black hole.
The angular momentum per unit mass of the test particle are

30

⇣p⌘
Lz
M
p
=
,
m2
M
1 (3M/p)(2.15)and the inspiral trajectory is generated by the flux of the above quantity which
is given by [63],L̇z⇢⇣ M ⌘3/2
32 m22 ⇣ M ⌘7/2
1247 M
=
1
+ 4⇡
5 M r
336 p
p
!)
⇣ M ⌘1/2⇣M ⌘44711 ⇣ M ⌘2
+ ca1 + cb1+ cc1
9072 pp
p

The values of the coefficients are ca1 =(2.16)

28.1517, cb1 = 60.9607, cc1 = 40.9998.

Conservative self force corrections modify the orbital frequency and are givenby, [63]d1:= ⌦ =
dt
M
1
=
M
+ m2

✓
✓

M
p
M
p

◆3/2
◆3/2

d0 + d1

+ ⌦,
(
1+

✓

M
p

◆

+ d1.5

✓

M
p

◆3/2

+ d2

✓

M
p

◆2 ! ). (2.17)

The values of the coefficients are d0 = 1/8, d1 = 1975/896, d1.5 = 1027 ⇡, d2 =1152343/451584. To compute the inspiral trajectory at each instant in time weevolve Lz and according to Eq. (2.16) and Eq. (2.17) and compute p from Lzfrom Eq. (2.15). The gravitational-wave strain, h(t), is given by,h(t) = A ⌦(t) p(t)2sin [2 (t)] .

(2.18)The overall scaling A depends on the source location with respect to thedetector. In this chapter we will use normalized template waveforms and so weleave this factor unspecified.

31

2. DETECTING INTERMEDIATE MASS-RATIO BINARIESWITH ADVANCED DETECTORS

2.3.3

TaylorT4

The TaylorT4 approximant is based on the post-Newtonian (PN) expansion (seesection III B. of Ref. [26] and references therein). The orbital phase is given bysolving the coupled di↵erential equations of the orbital phase and velocity, (t)and v(t), respectively,dv 3 (t)
= 0,
dt
M
dv
= K(t),
dt(2.19)(2.20)

where the acceleration K(t) is given Eq. (A.14).

The gravitational wave strain h(t) can then be expressed directly ash(t) = Av 2 sin [2 (t)] .

(2.21)The overall amplitude A depends on the source location with respect to thedetector. The amplitude is given in [73]. In this chapter we will use normalized
template waveforms and so we leave this factor unspecified.

2.4

SNR from inspiral, merger and ringdown

In this section we consider the relative contributions of the di↵erent portions of
the gravitational-wave coalescence signal to the SNR as a function of the IMRAC’smass. We take the instrument noise power spectral density (PSD) Sn (f ) to be
the aLIGO high-power, zero-detuned noise PSD [74] which is shown in Fig. 2.1.
The expectation value of the optimal matched filtering SNR, in the case when
the signal and template waveforms are identical, is given by [67]
⇣S⌘
N

max

= (h|h)1/2 ,
" Z

fmax

=

4<

fmin

f |h̃(f )|p
f Sn (f )!2d ln f#1/2,(2.22)

where limits of integration correspond to the bandwidth of the detector. Writ-

32

Figure 2.1: Advanced LIGO high-power, zero-detuned noise PSD [74]33

2. DETECTING INTERMEDIATE MASS-RATIO BINARIESWITH ADVANCED DETECTORS

ing the maximum SNR in the form above clearly separates it into contributions
from the signal strain, f |h̃(f )|, and the root-mean-squared (rms) noise spectralpamplitude, f Sn (f ), which is the strain signal associated with the detector noise.

One can gain insight into the relative contributions to the SNR from inspiral,merger and ringdown by comparing the gravitational-wave strain to the noiserms value. In Fig. 2.2 we show the strain for selected overhead and face-on (i.e.,
optimally-located and orientated) IMRAC sources at a fiducial distance of 1Gpc
as described by the EOBNR waveform family, and noise rms amplitude. The
ISCO frequency of each signal is shown as a vertical line. The strain from merger
and ringdown is thus the portion after the ISCO frequency. The contribution
to the strain from merger and ringdown from binaries with component masses(m1 , m2 ) = [(200, 20) M , (200, 2) M ], is greater than that of the noise rms
amplitude (black curve with triangles). In general, systems with ISCO frequencies
between 30 100 Hz, merge in the “bucket” of the noise curve, i.e., where thedetector is most sensitive. For example, for the (m1 , m2 ) = (200, 20) M system(red dotted curve in Fig. 2.2), the merger and ringdown contribute the bulk of
the SNR. Conversely, the inspiral contribution to the SNR is strongly suppressed
for such massive systems.
In Fig. 2.3 we show the maximum SNR, Eq. (2.22), as a function of the binary’s
total mass produced by inspiral-only and full EOBNR waveforms at four di↵erent
mass-ratios in the range 1/200  q  1/10. We construct inspiral-only EOBNR
waveforms by Fourier transforming the full waveform into the frequency domain(see Eq. (A.1) ) and truncating it at the ISCO frequency. We have considered
the SNR for optimally-located and orientated sources at a fiducial distance of 1
Gpc. The lower bound mass of the smaller body is set to m2 = 1.4 M which isthe canonical neutron-star mass.
We note that the IMRAC systems we are considering are liable to produce anappreciable SNR at distances greater than 1 Gpc. Such cosmological distances
correspond to red shifts z
0.2 [75]. Gravitational waves from binaries which
propagate in an expanding FRW universe experience a Doppler shifting of theirfrequency such that the observed gravitational-wave frequency fobs is related to
the gravitational-wave frequency at the source, fs , by fobs = (1/(1 + z)) fs [29].
Such a red shifted frequency corresponds to a binary with a chirp mass which

34

Figure 2.2: Strain of optimally-located and orientated IMRAC sources at afiducial distance of 1 Gpc as described by EOBNR waveforms, and aLIGO
noise. The corresponding ISCO frequency of each signal is shown as a vertical dashed line. The strain from the merger and ringdown from each source
contributes after the ISCO frequency. For the sources with component masses(m1 , m2 ) = [(200, 20) M , (200, 2) M ], the strain from merger and ringdown sits
above the noise spectrum. The SNR from the full EOBNR waveform and fromits inspiral-only portion are shown in Fig. 2.3.

35

2. DETECTING INTERMEDIATE MASS-RATIO BINARIESWITH ADVANCED DETECTORS

is larger by a factor of (1 + z) which is typically referred to as the “red shiftedchirp mass” of the binary. Thus, systems at red shift z = 0.2 appear to be 1.2
times heavier. Hence, when we refer to binary masses in the following pages, we
strictly mean the red shifted masses of the binary.
The lowest total mass for the q = 1/50, 1/100 and 1/200 subplots is set by
fixing the mass of the smaller body to m2 = 1.4 M . For the q = 1/10 subplot
in Fig. 2.3 the smallest total mass is set to M = 35 M as the inspiral phaseaccounts for ⇠ 100% of the SNR below this value. The lower limit of integration of
Eq. (2.22) is 10 Hz and the upper limit is 2048 Hz, which is the Nyquist frequency
of discretely sampled EOBNR waveforms generated at a sampling rate of t =4096 s in the time-domain. We only consider systems with total masses suchthat the ISCO frequency is greater than 10 Hz (our low frequency cut-o↵). The
highest total mass for each of the subplots in Fig 2.3 is set to M = 300 M which
ensures the ISCO frequency is greater than 10 Hz.
As anticipated from Fig. 2.2 there is a significant di↵erence in the SNR between inspiral-only and full EOBNR waveforms which can be seen at all fourmass-ratios. We also note that for systems with mass-ratios of q = 1/10 withtotal masses below around M = 35 M the inspiral phase is the dominant sourceof SNR. If we consider 3% as a fiducial value of the di↵erence between the full
SNR and the one associated to the inspiral-only waveform – which leads to a loss
in detection rates of 10% – this happens at M ⇡ 35 M for q = 1/10. For binaries with q = 1/50, 1/100 and 1/200, the minimum di↵erence in SNR between
inspiral-only and full waveforms is ⇡ 6%, 15% and 40% respectively for the mass
ranges considered in Fig 2.3.
In summary, we have shown that inspiral-only templates will miss a significant portion of the total SNR of IMRAC signals over the bulk of the detectablemass-space. Future searches will therefore require templates that can match thefull inspiral-merger-ringdown. However, there is a small region of the parameter
space for which inspiral-only templates may suffice for searches, without inducing drastic losses in detection rates. In the following section we quantify thee↵ectiveness of inspiral-only templates for searching for full coalescence signals in
aLIGO.

36

Figure 2.3: SNR of optimally-located and orientated IMRAC sources at afiducial distance of 1 Gpc vs total mass for four di↵erent mass ratios; q =(1/10, 1/50, 1/100, 1/200). The solid line is the SNR from full EOBNR waveforms and the dashed line from inspiral-only EOBNR waveforms truncated at
the ISCO frequency in the frequency domain. The lower bound mass of thesmaller body is set to m2 = 1.4 M which is the canonical neutron-star mass.
The lowest total mass for the q = 1/50, 1/100 and 1/200 subplots in Fig. 2.3 isset by fixing the mass of the smaller body to m2 = 1.4 M . For the q = 1/10
subplot the smallest total mass is set to M = 35 M as the inspiral accounts
for ⇠ 100% of the SNR below this value. We only consider systems with total
masses such that the ISCO frequency is greater than 10 Hz (our low frequency
cut-o↵). The highest total mass for each of the subplots in Fig 2.3 is set to M =300 M which ensures the ISCO frequency is greater than 10 Hz. We find that
there is a non-negligible contribution to the SNR from merger and ringdown in
IMRAC signals above a total mass of around M = 35 M . The di↵erence in
SNR between inspiral-only and full waveforms is at the 3% level at around M =35 M at q = 1/10. For binaries with q = 1/50, 1/100 and 1/200 , the minimumdi↵erence in SNR between inspiral-only and full waveforms are at the 6%, 15%
and 40% levels, respectively, over our mass range of interest. For IMRACs of astrophysical interest, more extreme mass ratios correspond to a greater total mass,
which can place the merger and ringdown at a frequency where the detector has
the greatest sensitivity.

37

2. DETECTING INTERMEDIATE MASS-RATIO BINARIESWITH ADVANCED DETECTORS

2.5

E↵ectiveness of inspiral-only templates for
IMRAC searches

We have shown in Sec. 2.4 that the SNR from merger and ringdown will provide
a significant contribution to the total SNR over a broad portion of the IMRAC
mass-space, c.f. Fig. 2.3. There is however a small portion of the parameter
space where the SNR is dominated by the inspiral phase which can be seen,e.g. for q = 1/10 binaries in Fig. 2.3 with total masses M  35 M . Thus
it is important to quantify the e↵ect of using inspiral-only templates to search
for IMRAC signals which contain inspiral, merger and ringdown phases. The
use of template waveforms that are not exact representations of the signals theyfilter degrades the SNR, as the optimal SNR can be recovered only when thetemplate waveform corresponds exactly to the signal h, cf. Eq. (2.22). In practice
however, we do not have access to an exact representation for h. Using a nonexact template waveform T to filter h caps the maximum recoverable SNR to
[76]
⇣S⌘
N

~
(h(~ )|T (✓)),~ (✓))~ 1/2✓~ (T (✓)|T
⇣S⌘
= ✏
,
N max
= max

(2.23)

where ~ and ✓~ represent the parameter vector of the signal and template, respectively. We define ✏ as the e↵ectiveness of a template waveform family T at
recovering the maximum SNR from a gravitational-wave signal h; by definition0  ✏  1. This quantity is also referred to as the “fitting factor” in the literature [76]. It is convenient to define waveforms normalized to unit norm asâ(f ) = ã(f )/(a|a)1/2 so that (ĥ|ĥ) = (T̂|T̂) = 1 and the e↵ectiveness can be
written succinctly as~ .
✏ = max(ĥ(~ )|T̂(✓))(2.24)✓~

Using normalized waveforms also has the advantage of eliminating the dependenceof the waveforms on the source orientation and distance, which enter as an overallscaling.

38

Calculating the e↵ectiveness, Eq. (2.24), requires maximizing over the component masses (m1 , m2 ) and the time and phase at coalescence. We can efficiently
maximize over the time and phase by Fourier transforming the integrand of thenoise-weighted inner-product [77],z(tc ) = 4

Z fmaxfmin

dfã(f )b̃⇤ (f ) 2⇡if tce,
Sn (f )

(2.25)

which yields a complex time series whose elements correspond to the innerproduct of a and b as one of the signals is time-shifted with respect to the other.
We can efficiently find the time at coalescence, tc , by finding the time at which thenorm of this time series is a maximum. The phase at coalescence c is then automatically given by finding the argument of the time-series at its peak amplitude.
We thus modify the inner product (a|b):(a|b) ! (a|b)0 = max |z(tc )| ,tc

(2.26)

which we will adopt as the definition of the inner-product for the remainder ofthis chapter.
To compute the e↵ectiveness of an inspiral-only IMRAC search we evaluate
Eq. (2.24) for signals covering the IMRAC mass space. We take as our signal
waveform, h, the full inspiral-merger-ringdown EOBNR waveform. We take thetemplate, T , to be an inspiral-only EOBNR waveform, formed by truncating thefull EOBNR waveform at the ISCO frequency in the frequency domain. With
such signals and templates the e↵ectiveness provides a measure of the maximum
SNR which could be achieved through using an inspiral-only template to filter
full coalescence-signals. To get a broad coverage of the IMRAC mass space wecompute Eq. (2.24) for signals whose source masses cover the ranges 1.4 M
 m2  18.5 M and 24 M  m1  200 M , with mass ratios spanning therange q := m2 /m1 2 [1/140, 1/10]. For each signal we evaluate the e↵ectiveness,
Eq. (2.24), where the template T describes the inspiral-only portion of an EOBNR
waveform. We maximize over time and phase by maximizing the inner product of
the signal with an inspiral-only EOBNR template, Eq. (2.26). The maximization
over the masses is performed by finding the largest inner product between the39

2. DETECTING INTERMEDIATE MASS-RATIO BINARIESWITH ADVANCED DETECTORS

signal and a bank of template waveforms. The template bank is characterised byintrinsic parameters (M , ⌘) which are the chirp mass and symmetric mass-ratioof the system: M = (m1 m2 )3/5 /(m1 + m2 )1/5 , ⌘ = (m1 m2 )/(m1 + m2 )2 , wherem1 and m2 are the component masses of the binary. The bank spans an extended
mass range 3 M  M  30 M and 0.0065  ⌘  0.082. The results of thee↵ectiveness of an inspiral-only IMRAC search are shown in Fig. 2.5.
Inspiral-only templates are ⇠ 98% e↵ective at filtering full coalescence signals
for total masses M . 30 M . Such systems have an ISCO frequency 150 Hz .
fISCO which is well within the peak sensitivity of the noise curve. However for thebulk of the mass space the e↵ectiveness is below 75%. This is unsurprising given
the SNR curves in Fig. 2.3 which clearly show the importance of the contributionof merger and ringdown to the SNR.
The loss in SNR incurred through using inspiral-only templates directly affects detection rates. Because the SNR scales inversely with the distance, theobservable volume will scale with the cube of the e↵ectiveness. Assuming that
GW sources are isotropically distributed in the sky, the fractional loss in detection
rates will be 1 ✏3 . The percentage loss in detection rates through using inspiralonly EOBNR templates to recover the full coalescence signal is also shown in
Fig. 2.5. Over a broad portion of the mass-space inspiral-only templates incur
losses in detection rates between 60 85%. As the total mass of the binary approaches 440 M the Schwarzschild ISCO frequency, Eq. (2.1), approaches 10 Hz
which is near the low-frequency cut-o↵ of the detectors. Hence the relative contribution of the inspiral phase to the coalescence signal of heavier systems diminishes
until the only contribution is from the merger and ringdown. This is a striking
indication of the need of merger and ringdown in IMRAC template waveforms.
This suggests the importance of full numerical simulations in this regime in order to construct such a reliable waveform family including inspiral, merger, andringdown phases.
We identify three regions in the m1 -m2 plane in which various searches could
be constructed. The regions are defined by contours of constant e↵ectivenesspwhich are approximately given by C = (m1 /M ) m2 /M , which are foundpurely empirically, with 1.4 M  m2  18.5 M and mass-ratios q 2 [1/140, 1/10].
The e↵ectiveness is related to C by ✏ ⇡ 1/100 ⇥ (1.6 C 7.3 ⇥ 10 3 C 2 ).

40

Figure 2.4: E↵ectiveness of inspiral-only EOBNR templates to filter full inspiral,merger and ringdown EOBNR signals as a function of the source component
masses, and corresponding losses in detection rates. The diagonal corresponds
to mass-ratio q = 1/10. Inspiral-only EOBNR templates are constructed bytruncating the full waveform at the ISCO frequency in the frequency-domain.
For the bulk of the parameter space inspiral-only templates are . 75% e↵ective
at filtering inspiral, merger and ringdown signals. Inspiral-only templates are⇠ 97 98% e↵ective for total masses M . 30 M . Inspiral-only templates
within the 90%-e↵ectiveness contour should be sufficient for IMRAC searches
without incurring greater than 30% losses in detection rates.

41

2. DETECTING INTERMEDIATE MASS-RATIO BINARIESWITH ADVANCED DETECTORS

Between the 97%- and 90%-e↵ectiveness contours, the losses in detection rates
are between 10% . L . 27% and so an inspiral-only search could be sufficient
without incurring drastic losses in detections. The region bound from below ine↵ectiveness by the 90%-e↵ectiveness contour is defined by C  100, with thee↵ectiveness increasing with decreasing C. Between the 90%- and 80%- e↵ectiveness contours the losses in detection rates are around 27% . L . 50%. Thus,
within this region searches will be limited by the lack of merger and ringdown
in template waveforms, though an inspiral-only search would be feasible in principle. This contour is defined by 100 . C . 150. Below the 80%-e↵ectiveness
contour inspiral-only searches will incur losses in detection 50% < L and somerger and ringdown will be crucial for searches. The region bound from abovein e↵ectiveness by the 80%-e↵ectiveness contour which is defined by 150 . C,with e↵ectiveness decreasing with increasing C. The results are summarized inTable 2.1.

42

43

27% . L . 50%

80% . ✏ . 90%

Contours of constant ef- Implication for searches
fectiveness in mp
1 -m2 plane(C = (m1 /M ) m2 /M )within mass range of interest
C . 100
Inspiral-only search sufficient but with nonnegligible loss in detection
rates.
100 . C . 150
Inspiral-only search possible but limited by lackof merger and ringdown.
Could potentially lose halfof signals with inspiral-only
templates.
150 . C
Merger and ringdown crucial for searches. Will missover half of signals withinspiral-only templates.

Table 2.1: E↵ectiveness of inspiral-only searches, the corresponding loss in detection rates and the region in the m1 m2 plane bounded by constant-e↵ectiveness contours. For a given region in the m1 -m2 plane bounded by constante↵ectiveness contours we summarize the implications for IMRAC searches.

50% . L

10% . L . 27%

90% . ✏ . 97%✏ . 80%Loss in detection
rates, L(%)

E↵ectivenessof inspiral-only
search, ✏(%)

2. DETECTING INTERMEDIATE MASS-RATIO BINARIESWITH ADVANCED DETECTORS

2.6

Comparison of inspiral-only waveforms

We have shown that merger and ringdown are crucial for e↵ective searches overa large portion of the IMRAC mass space, though there is a small region inwhich an inspiral-only search could be constructed without incurring losses indetection rates greater than around 27%. For this region, it is therefore important
to study whether currently available waveforms are sufficiently accurate. The
inspiral phase can be computed using perturbative expansions and thus it isinteresting to quantify the consistency of di↵erent expansions.
To assess the e↵ectiveness of the EOBNR inspiral, we employ a waveform family designed to approximate intermediate mass-ratio inspirals which we refer toas “Huerta-Gair” (HG) waveforms [44]. HG waveforms describe only the inspiralportion of the coalescence signal.
We repeat the study done in the previous Section using now the HG waveformfamily as the signal h and inspiral-only EOBNR as the template T . The results
are reported over the whole parameter space in Fig. 2.5. For completeness, inTable 2.2 we also show the values of the e↵ectiveness, Eq. 2.24, for selected mass
combinations of EOBNR inspiral-only templates for filtering full EOBNR and
HG signals respectively.

44

45

5

10

15

50

m1

100

150

200

0.92

0.94

0.96

0.98

5

10

15

50

m1

100

150

200

0.92

0.94

0.96

0.98

Effectiveness
1.00

Figure 2.5: E↵ectiveness of inspiral-only EOBNR templates at filtering HG signals (left) and TaylorT4 signals (right)
as a function of the source component masses encoded in the signal.

m2

Effectiveness
1.00

m2

2. DETECTING INTERMEDIATE MASS-RATIO BINARIESWITH ADVANCED DETECTORS

m1
(M )
50
100
200
50
100
200
5
20
100

m2
(M )
5
5
5
1.4
1.4
1.4
5
20
100

Signal waveforms
full EOBNR Huerta-Gair TaylorT4
0.90
0.95
0.96
0.76
0.97
0.97
0.53
0.99
0.99
0.96
0.96
0.90
0.86
0.98
0.94
0.67
0.98
0.98
0.99
0.89
0.99
0.99
0.92
0.99
0.52
0.51
0.50

Table 2.2: Summary of e↵ectiveness of inspiral-only EOBNR template waveforms
in recovering signals modelled using di↵erent waveform families – full EOBNR,
Huerta-Gair and TaylorT4 – for selected component masses. Merger and ringdown become more prominent in the coalescence signal as the total-mass of thesystem is increased. The EOBNR inspiral is typically better at matching HG
signals in the IMRAC regime than TaylorT4 signals. Results for equal mass-ratio
systems are shown for reference below the horizontal line.

For the high-mass part of the mass-space the e↵ectiveness of the EOBNR
inspirals with respect to the HG waveforms is close to 100%. This is perhaps
unsurprising because very high mass systems will have short inspirals and possibledi↵erences in the waveforms will not produce a significant degradation of SNR
when matched over a small number of wave cycles. However, for lighter systems
the e↵ectiveness can be as low as 90 %, which occurs in the region of mass spacein which inspiral-only searches would be most feasible (see Table 2.1).
For reference we also compare inspiral-only EOBNR templates to TaylorT4
signal waveforms (which are inspiral-only). We construct signal waveforms on thesame grid in m1 m2 as for HG waveforms and use the same template bank ofinspiral-only EOBNR waveforms. The results are summarized in the right panelof Fig. 2.5, and in Table 2.2 for selected masses. We find that the EOBNR inspiralhas good filtering efficiency for the TaylorT4 waveform family. However, EOBNR
is clearly a better match to the HG waveform family over a larger range of masses
and mass ratios than to TaylorT4. This can be seen more clearly by comparing
the subplots in Fig. 2.5. This is unsurprising given that the PN expansion is46unreliable at high velocities and highly asymmetrical mass-ratios. For orbital
velocities v/c = (M ⇡f )1/3 & 0.2 the PN energy flux deviates significantly fromnumerical results, see [e.g., 38; 39]. A binary at its ISCO frequency has v/c ⇠ 0.4
which is well beyond the region of validity.

2.7

Discussion and conclusion

We have shown that over the bulk of the IMRAC mass space, merger and ringdown contribute significantly to the gravitational-wave coalescence signal. This
happens despite the suppression of the power in the merger and ringdown in signals from binaries with very asymmetric mass ratios. The importance of merger
and ringdown is due to the greater sensitivity to these waveform portions forhigh-mass signals, for which most of the inspiral may fall at frequencies below
the detector’s sensitive band.
However there is a relatively large patch in mass space in which the inspiralonly waveforms are more than 90% e↵ective. We identified three regions in whichdi↵erent searches could be considered appropriate based on thresholds of acceptable losses in detection rates. The mass space splits into a region in which
inspiral-only searches could be feasible, incurring losses in detection rates of up
to ⇠ 27%; a region in which searches would be limited by lack of merger and ringdown in template waveforms, incurring losses in detection rates up to 50%; anda region in which merger and ringdown are critical to prevent losses in detection
rates over 50%. The search regions are summarized in Table 2.1.
We have further shown that in the region of the IMRAC mass space in which
inspiral-only searches are feasible, approximants adapted to asymmetric massratio binaries are important, as here the binary is liable to have highly relativistic velocities v & 0.2. We considered a waveform family designed to describeintermediate mass-ratio binaries which we referred to as the “Huerta-Gair” (HG)
waveform family By computing the e↵ectiveness of inspiral-only EOBNR waveforms to filter signals described by the HG waveform family, we showed losses
in recovered SNR could be as great as 10%. In Table 2.2 we summarize thee↵ectiveness of the signal–template combinations used in the chapter.
We believe that template waveforms for IMRAC searches will benefit from47

2. DETECTING INTERMEDIATE MASS-RATIO BINARIESWITH ADVANCED DETECTORS

calibration to several numerical simulations. We note that there already exists
one very short numerical waveform of a q = 1/100 binary which we have not usedin our study, and which EOBNR is not currently calibrated to [45].

48Chapter 3
TOWARDS RAPID COMPACT
BINARY PARAMETERESTIMATION USINGINTERPOLATEDWAVEFORMS
3.1

Introduction

Astronomy and tests of fundamental physics with gravitational waves from compact binary coalescence (CBC) will ultimately be limited by our ability to estimate the binary’s source parameters from the gravitational-wave signal [e.g.,78; 79; 80]. CBC sources with total masses in the range 2 M . M . 500 M
will be amongst the prime sources for Advanced LIGO [58] and Advanced Virgo[59] when they begin operating around 2015 [5].
In a Bayesian treatment of parameter estimation, one is interested in theposterior probability distribution of the set of source parameters of the underlying
model given an observed stretch of data. Waveform computation represents the
majority of the computation cost in the Bayesian analysis of CBC sources, so thetotal computational cost scales roughly linearly with waveform generation. This49

3. TOWARDS RAPID COMPACT BINARY PARAMETERESTIMATION USING INTERPOLATED WAVEFORMS

becomes burdensome when one needs to explore a large dimensional parameter
space as the number of waveform computations is large, e.g., O(106 ) [25].

Recently, Cannon et al. [49] showed that a truncated singular value decomposition (SVD) can be applied to gravitational-wave template banks which span
the two parameters describing the masses of the coalescing binary. The SVD decomposes the bank into a set of “basis templates” and projection coefficients. In
general, the number of basis templates is much less than the total number of templates in the bank. Furthermore, the projection coefficients can be interpolated
across the domain of the bank [50]. Template waveforms can thus be interpolated. Because the intrinsic parameter space of CBC sources with non-spinningcomponents is two-dimensional (the two mass parameters), it is possible to set
up the waveform computation for parameter estimation such that the waveform
calculations are done by interpolation alone. However, the errors incurred frominterpolation could, in principle, a↵ect parameter-estimation accuracy.
In this chapter, we describe the application of SVD-interpolated waveforms to
CBC parameter estimation. For a simulated data set containing a gravitational
wave signal we provide a proof of principle that SVD-interpolated waveforms canbe used for parameter estimation without significantly a↵ecting the accuracy of
the inferred probability distributions of the source parameters. We further showthat the computational cost of using interpolated waveforms is around an orderof magnitude less than that of commonly-used time-domain waveform families.
This technique has the potential to increase the computational efficiency of CBC
parameter estimation when the computational cost is dominated by waveformcomputation. Our application of the SVD is limited to a small patch of parameter
space about the injected signal value.
This chapter is organized as follows. In Sections 3.2 and 3.4 we outline theprinciples of parameter estimation for CBCs and interpolating template waveforms based on the SVD, respectively. In Sec. 3.5 we describe the applicationof SVD-interpolated waveforms to parameter estimation. In Sec. 3.6 we compare the results of using interpolated and non-interpolated waveforms for parameter estimation and compare the computational cost of interpolation to using
non-interpolated waveform families. In Sec. 3.7 we consider the future of using
SVD-interpolated waveforms for parameter estimation and discuss the technical50requirements of implementing these waveforms in parameter-estimation pipelines.

3.2

CBC parameter estimation

The central quantity of interest in Bayesian parameter estimation is the posterior
probability density function (PDF) of a set of parameters ✓~ which parameterize amodel, H, assumed to describe a data set d. The PDF is related to the likelihood
function and prior probability via Bayes’ theorem and is given by~ H) =p(✓|d,~~ H)
P(✓|H)
L(d|✓,,
p(d|H)(3.1)~ H) is the likelihood function and P(✓|H)~
where L(d|✓,is the prior probability~ The quantity in thewhich encodes our a priori belief in the distribution of ✓.
denominator, p(d|H), is known as the evidence and is an overall normalization
factor which we will not deal with here.
The CBC parameter vector ✓~ is high-dimensional. The phasing and amplitudeof a waveform from a non-spinning coalescing compact binary source is controlled
by two mass parameters, the chirp mass M = (m1 m2 )3/5 /(m1 + m2 )1/5 and symmetric mass ratio ⌘ = (m1 m2 )/(m1 + m2 )2 , where m1 and m2 are the component
masses of the binary. In addition, a gravitational wave source with respect to the
Earth is specified by location dependent parameters. These are the distance fromthe Earth D, inclination ◆, right ascension ↵, declination , polarization phase 'and time and phase at coalescence, tc and c . In general, the CBC parameter
vector ✓~ is nine-dimensional for circular binaries with non-spinning components.
An illustration of a coalescing binary system, with the set of source parameters isshown in Fig. 3.1. One of the goals of gravitational-wave astronomy is to estimate
the PDF of the parameters of a candidate gravitational wave source in order toassign a meaningful probability to our measurements of the source properties anddemographics. To compute the right hand side of Eq. (5.1), we directly evalu~ H). Under the hypothesis that the data, d,
ate the likelihood function, L(d|✓,~ theconsists of Gaussian, stationary noise n and a gravitational-wave signal h(✓),

51

3. TOWARDS RAPID COMPACT BINARY PARAMETERESTIMATION USING INTERPOLATED WAVEFORMS

Figure 3.1: Cartesian coordinate system (x, y, z) attached to a ground based
interferometer and a coalescing binary system in that frame [81]. The binary’s
components are labelled by two mass parameters m1 , m2 . The sky position isdetermined by the two polar angles (✓, ) which fixes the position vector N̂ . The
inclination is given by cos ◆ = N̂ .L̂ where L̂ is the orbital angular momentum
vector of the binary.

52

likelihood is a Gaussian [12]:~~~ H) = e (d h(✓)|d h(✓))/2 ,L(d|✓,(3.2)

In the above, (a|b) is the usual noise-weighted inner-product,(a|b) = 4<

Z fmaxfmin

dfã(f )b̃⇤ (f ),Sn (f )

(3.3)

and Sn (f ) is the detector’s noise power spectral density (PSD). The limits of integration correspond to the bandwidth of the detector. A significant computational
cost of evaluating the likelihood comes from computing the template waveform~ at each point in the parameter space.
h(✓)
The full PDF is multi-dimensional and to get estimates on an individual pa~ Writing
rameter we work with the marginalized PDF of a single parmeter ✓A 2 ✓.
~ the marginalized one-dimensional PDF of ✓A is thus:✓~ = (✓A , ⇥),p(✓A |d, H) =

Z

~
⇥

~ H).
~ p(✓|d,d⇥(3.4)

To efficiently evaluate the likelihood function we typically use a stochastic sampling algorithm. Here we employ Markov-chain Monte Carlo (MCMC), whose
application to gravitational-wave parameter estimation is described in detail in[13], and in Sec. 3.2.1.
We use the stationary phase approximation (SPA) inspiral waveforms for both
our simulated signal and template model. This allows us to directly “inject”
the signal waveform into simulated frequency-domain noise without performing
an additional Fourier transformation, which could introduce spurious artifacts
related to the abrupt in-band termination of a time-domain waveform. The postNewtonian frequency-domain waveform for a face-on, overhead binary has theform
A(M; f ) 2⇡if tc i c +i (M,⌘;f )h̃(f ) =
e
,
(3.5)
D
where expressions for the phase (M, ⌘; f ) is given by Eq. (A.11) for TaylorF2
post-Newtonian waveforms at 3.5 post-Newtonian order and the amplitude is53

3. TOWARDS RAPID COMPACT BINARY PARAMETERESTIMATION USING INTERPOLATED WAVEFORMS

given by [29]1 ⇣ 5 ⌘1/2 5/6 7/6
M f
,
(3.6)
⇡ 2/3 24in geometric units, G = c = 1. In general, the template model will be determinedby the needs of a particular analysis. Below, we describe the MCMC algorithm
which we use to compute the posterior PDF, Eq. 5.1.
A(M; f ) =

3.2.1

MCMC

Our goal is to compute the posterior PDF of the set of astrophysical parameters✓~ which describe a compact binary, given a data set containing a gravitationalwave signal from such a source. The posterior PDF is given by Eq. (5.1). In
principle, one could simply grid up the parameter space of interest, and computethe posterior PDF at each grid point. However, a number of problems become
immediately apparent with this approach. Firstly, the compact binary parameter
space is high dimensional (9 dimensional for binaries with non-spinning components, and 15 dimensional for binaries with spinning components). Assuming anequally spaced grid of N points in each dimension, one would require between
N 9 N 15 computations. Even a modest grid with N ⇠ O(100) points in eachof the nine dimensions would require a minimum of 1009 ⇠ 1018 computations.
This coupled with the fact that evaluating the posterior PDF is an expensive
computation renders this approach hopeless. Furthermore, the problems withthe approach are exacerbated when the dimensionality is increased to 15, forobvious reasons. Hence, it is useful to consider algorithms which stochasticallysample the parameter space preferentially, in regions of high posterior probability.
Ultimately we want the output of such an algorithm to be a set of sample
points {P} which are drawn from the posterior PDF. In the context of MCMC,
this set is known as a “Markov chain” [82], and the goal of the MCMC is to
produce such a Markov chain via a stochastic process. The stochastic process
we will consider is the classical Markovian random walk through the parameter
space. By “Markovian” we mean that the random walk is such that given the ith
state in the Markov chain, the next state (i + 1) is determined solely by i, andno others. In the context of parameter estimation this means that the (i + 1)th
sample of the posterior PDF in the Markov chain, P (✓~i+1 ), depends only on the54ith sample, P (✓~i ). Thus, by Bayes’ theorem we haveP (✓~i+1 ) =P (✓~i ) P (✓~i |✓~i+1 ).
P (✓~i+1 |✓~i )(3.7)

The above equation provides a criterion known as “detailed balance” [82] which
the Markovian random walk must obey:
P (✓~i+1 ) P (✓~i+1 |✓~i ) = P (✓~i ) P (✓~i |✓~i+1 ) .

(3.8)

In order to utilize the random walk, encoded in the detailed balance condition
Eq. (3.8), to form the Markov chain of samples, we need to be able to specify
a means of computing the conditional probability of state (i + 1) given state i,
P (✓~i+1 |✓~i ), and vice versa. For practical purposes one ought to be able to proposea transition from state i to state (i + 1) and accept the proposal provided thatdetailed balance is satisfied. The proposal should be arbitrary in the sense thatit should reflect the requirements of a particular analysis. For example, whenestimating the parameters of compact binaries we may want to explore known
correlations or degeneracies between parameters, or pairs of parameters. We
thus require that we have an additional distribution g(✓~i+1 |✓~i ) which proposes
which point in parameter space should be explored next. This is often referred
to as a “jump proposal” in the literature [13; 25]. However, this alone is not
sufficient to satisfy detailed balance, Eq. (3.8), because not all transitions willbe allowed. Hence we can only accept transitions with a particular probabilityA(✓~i+1 |✓~i ) which is known as the “acceptance probability”. Thus we can write
the conditional probability of state (i + 1) given state i as
P (✓~i+1 |✓~i ) = g(✓~i+1 |✓~i ) A(✓~i+1 |✓~i )(3.9)

and the detailed balance condition, Eq. (3.8), is given by
P (✓~i+1 ) g(✓~i+1 |✓~i ) A(✓~i+1 |✓~i ) = P (✓~i ) g(✓~i |✓~i+1 ) A(✓~i |✓~i+1 )55(3.10)

3. TOWARDS RAPID COMPACT BINARY PARAMETERESTIMATION USING INTERPOLATED WAVEFORMS

Defining the ratio of acceptance probabilities,r :=A(✓~i+1 |✓~i )P (✓~i+1 ) g(✓~i |✓~i+1 )=,A(✓~i |✓~i+1 )P (✓~i ) g(✓~i+1 |✓~i )(3.11)

we accept a transition with probability⇣⌘
acceptance probability = min 1, r .

(3.12)We denote the ratio r as the “acceptance ratio”. Together with a (user specified)jump proposal, the above condition provides us with all the information required
to draw samples from the posterior distribution. We can compute the ratio
P (✓~i+1 )/P (✓~i ) as this is just a ratio of likelihoods, which can be computed fromEq (3.2), and once we specify a jump proposal g(.), we can evaluate the acceptanceratio.
The method for drawing samples from the posterior PDF which we have explained above was first proposed by Metropolis and Hastings [83; 84] and hence
is referred to as the Metropolis-Hastings algorithm for MCMC. In practice, specifying good jump proposals is extremely difficult and we do not discuss this topic
here. For the purposes of using MCMC for compact binary parameter estimation,
we use the implementation within the LAL library [85], known as lalinference.

3.3

Singular value decomposition of gravitationalwave template banks for CBC sources

The SVD which we will apply to gravitational waveforms for parameter estimationis adapted from a technique used for gravitational-wave detection [49; 50]. To detect gravitational-waves from binaries with non-spinning components, one filters
the data with a large number of templates which span the space of the CBC chirpmass and symmetric mass-ratio. The component masses of the gravitational-wave
signal source are not known a priori and so one must scrutinize the continuous
mass space of CBC signals to high fidelity. Templates in the space are distributed
such that there is a fixed loss in SNR between nearest-neighbour templates which56

is small. This is the so called “minimal match” between neighbouring templates,
which is typically set to be greater than 97% [6; 7; 8; 49; 50] and corresponds to aloss in SNR of less than 3%. With such a template bank one might expect somecorrelation between the templates as the SNR between neighbouring templatesdi↵ers by only 3 % [49]. These correlations can be exploited by performing a SVD
of a template bank which is used to generate an orthogonal set of basis templates
and coefficients from which the original bank can be recovered. In particular,
we will demonstrate that many of the basis vectors are redundant and hence theoriginal bank can be recovered using a relatively small number of bases. Below
we illustrate the SVD as applied to template waveforms. In Appendix B we provide supplementary material on the singular value decomposition, including an
example of the SVD of a simple matrix.
The optimal filter for a matched filtered search is a gravitational-wave template, weighted by the inverse of the amplitude spectral density, a process known
as “whitening”. We define the whitened frequency series asã(f )ãw (f ) = p,
Sn (f )

(3.13)

where a(f ) is the un-whitened frequency series. The study in this chapter deals
exclusively with discretely sampled data and template waveforms. Because ofthis we find it useful to adopt the following notation for the data and waveforms.
We denote the discretely sampled whitened data-stream, and the ↵th template
waveform in a discrete template bank, as d~ and ~h↵ respectively, where it is understood that these correspond to whitened versions of the raw, un-whitened data
and template respectively. Both the data-stream and template are discretely sampled at a fixed sampling rate, t, in the time-domain. The elements of the data
and waveforms at time i t are di and h↵i respectively. We adopt the conventionof [49] and define the SNR of the ↵th template as the vector inner product⇢↵ = < ~h↵ . d~ ,(3.14)

where the definition of the vector inner product is given by Eq. (A.10). Note that
this definition of the SNR implicitly assumes that the templates ~h↵ are normalized57

3. TOWARDS RAPID COMPACT BINARY PARAMETERESTIMATION USING INTERPOLATED WAVEFORMS

to unit norm.
For non-spinning CBC sources the gravitational-wave phasing is controlled by
the component masses of the binary. The template bank for a CBC search is
generated by constructing a grid in the M ⌘ plane such that for neighbouring
templates indexed by ↵ and ↵0 we have ~h↵ . ~h↵0 > 0.95. Instead of filtering thedata using N real-valued filters we wish to remove the redundancy in the set of
templates by reproducing ⇢↵ to the desired accuracy ⇢0↵ by linearly combining theoutput of a smaller set of real-valued filters, ~uµ . Hence, we want to truncate thenumber of unique filters such that [49]
⇢0↵ =

0
N
X1

~ ,
A↵µ (~uµ . d)

(3.15)µ=0

where A is a real-valued “reconstruction matrix” and the total number of innerproducts required for the matched-filtered search is reduced from N to N 0 . To
find the basis vectors ~uµ we compute the SVD of a real-valued N ⇥ L matrix H,with N  L, whose columns are template waveforms:
H = {~h0 , ~h1 , . . . , ~hN 1 } .

(3.16)The SVD factors H as [49]
H⌫j =

N
X1

v⌫µ µ uµj ,(3.17)µ=0

where v is a unitary matrix of reconstruction coefficients, ~ is a vector of singularvalues and u is a matrix of orthonormal basis vectors ~uµ . The matrix H can beapproximately reconstructed such that there are fewer basis vectors than templates in the original template bank. We define the approximate reconstructionof H, H0 , as0N
X10H⌫j :=v⌫µ µ uµj ,(3.18)µ=0

where N 0 < N which reduces the number of columns of u. The reconstructed58

SNR, Eq. (3.15), is then given by
⇣

⇢0↵ =

⌘
~ 0 . d~
H
↵

0
N
X1

=

(3.19)~ .
v↵µ µ (~uµ . d)

(3.20)µ=0

Comparing the above expression with Eq. (3.15), it follows that the SVD allows
the ↵th template waveform, ~h↵ , to be written as
~h↵ =

0
N
X1

Mµ↵ ~uµ ,(3.21)µ=0

where we have defined the “projection coefficient” M⌫↵ asMµ↵ := v↵µ µ .

(3.22)The number of bases N 0 can be set by fixing the fractional loss in SNR which isincurred through truncating the sum in Eq. (3.18). The expected fractional loss
in recovered SNR of any filter is given by [49]
D ⇢E
⇢

N 1

=

1 X
2N µ=N 02µ ,(3.23)and the required number of filters N 0 is thus fixed by Eq. (3.23).

3.4

Interpolating template waveforms using the
SVD

The accuracy requirements for gravitational-wave parameter estimation are typically higher than for searches. While a (dense) discrete template bank is sufficientfor searches, parameter estimation requires that we are able to stochastically sample from a continuous parameter space in order to map the full likelihood-surface
where the bulk of the posterior support is contained. Constructing a discrete tem-59

3. TOWARDS RAPID COMPACT BINARY PARAMETERESTIMATION USING INTERPOLATED WAVEFORMS

plate bank is impractical for parameter estimation because one does not know apriori the shape of the surface. Conversely, for detection one is only interested
in a single number, the maximum SNR, and the fractional loss in SNR of any
stretch of data can be estimated from the error associated with the finite griddingof the bank.
We can construct a continuous representation of any template waveform ina given domain using the SVD of a template bank of gravitational waveforms.
Recall that the SVD allows any template waveform to be written as a linear
combination of basis waveforms and projection coefficients, c.f. Eq. (3.21). While
this representation is discrete, the set of basis waveforms ought to be complete
over the entire domain of the (discrete) template bank provided that neighbouring
templates are sufficiently similar. It therefore follows that a template waveform~h↵0 , which is not in given template bank, should also be expressible as a sumof the same basis waveforms and its own projection coefficients. The projection
coefficients can be found as follows. By the orthonormality of the bases, the innerproduct of ~h↵0 with a basis template ~uµ produces the µth projection coefficientassociated with ~h↵0 :~h↵0 . ~uµ = Mµ↵0 .
(3.24)
This suggests that there is a continuous description of the projection coefficients
across the domain of the template bank, and hence they ought to be amenable
to interpolation. Because the template bank spans the M ⌘ plane a single
index ↵ is somewhat cumbersome for enumerating the projection coefficients,
which are function of two variables, M and ⌘. Without loss of generality we alter
the notation of the projection coefficients and re-define the SVD of a template
waveform as
0N
X1~h(Mk , ⌘l ) =Mµ (Mk , ⌘l ) ~uµ ,(3.25)µ=0

where k and l enumerate points on a grid in the M ⌘ plane. An “interpolated
waveform”, ~h0 , can thus be constructed according to Eq. (3.25) from a linear
combination of interpolated coefficients and basis waveforms:~h0 (M, ⌘) =X

Mµ0 (M, ⌘) ~uµ ,µ60(3.26)

where Mµ0 (M, ⌘) is an interpolated projection coefficient evaluated at an arbitrary
point (M, ⌘) within the domain of the template bank. In order to interpolate
the projection coefficients, we follow [50] and use Chebyshev polynomials of thefirst kind, although in principle other interpolation schemes could be used. The
Chebyshev polynomials of the first kind are defined as [50]TJ (x) =

(x

p

x2

1)J + (x +
2w

p

x2

1)J

,

(3.27)

pwhere w = (1 + J0 )(Jmax + 1)/2 is a normalization constant and J0 is the
Kroenecker delta. The polynomials TJ (x) satisfy the discrete orthogonality condition
JX
max
TI (xj )TJ (xj ) = IJ .
(3.28)
j=0

For interpolation, we use a net of points, scaled such that each dimension coversthroot of the Chebyshev polynomials, comthe interval [ 1, 1], located at the Jmax
monly referred to as interpolation nodes. This net of points is used in order tomitigate oscillatory behaviour at the boundary of the interpolation interval. Fora single dimension, the nodes occur at the locations
j + 12
xj = cos ⇡Jmax + 1!

.

(3.29)In order to obtain the reconstruction coefficients at these locations, we project
waveforms from these locations (scaled to the appropriate points in M and ⌘)onto the basis vectors. We generate waveforms at the (scaled) Chebyshev nodes
using the following maps!
xj + 1
M [xj ] = Mmin +M,2!
yl + 1⌘ [yl ] = ⌘min +⌘,2(3.30)(3.31)

which project Chebyshev nodes on the interval [ 1, 1] onto the M and ⌘ dimensions. Here Mmin is the lower boundary in the M-dimension and M =61

3. TOWARDS RAPID COMPACT BINARY PARAMETERESTIMATION USING INTERPOLATED WAVEFORMS

Mmax Mmin where Mmax is the upper boundary of the M-dimension, andvice versa for the ⌘-dimension. From the values on the net, we interpolate toother positions in the M ⌘ plane using 2D-Chebyshev interpolation for each
set of reconstruction coefficients Mµ (M, ⌘). These reconstruction coefficients areprojected onto the Chebyshev polynomials [50]CKLµ =

K
max L
max
X
X
k=0

TK (xk )TL (yl )Mµ (M [xk ] , ⌘ [yl ]) ,(3.32)

l=0

which provide weights for the 2D-Chebyshev interpolation. The values of Kmax
and Lmax correspond to the number of templates in the M and ⌘ directions
respectively. The interpolated reconstruction coefficients are then given by a
combination of the weights together with the Chebyshev polynomials [50]Mµ0 (M, ⌘) =

K
max L
max
X
X

CKLµ TK (x [M])TL (y [⌘]) ,(3.33)

K=0 L=0

where the inverse maps x [M] and y [⌘] project points on M and ⌘ onto theinterval [ 1, 1]. The inverse maps are given byx [M] =y [⌘] =

M
⌘

Mmin
M
!

⌘min
⌘

!

1,

1.

(3.34)(3.35)

Below we illustrate the application of the SVD to parameter estimation.

3.5

Parameter estimation using interpolated waveforms

We will compare the marginalized one-dimensional PDFs, Eq. (3.4), obtainedby using an interpolated template waveform family to those generated using thestandard, non-interpolated template waveform family for the same data set. For
illustration we consider a toy example with five free parameters. In addition,62

we consider data from a single detector, as opposed to a network, for simplicity.
~ andWe generate a single-detector data set d~ containing a signal waveform h(✓)
Gaussian stationary noise n with a power spectral density typical of initial LIGO
[1]. By only having five free parameters we e↵ectively set the prior on the otherfour to be delta functions centered on the signal values. We chose to fix the sky
position (↵, ) and the inclination and polarization angle (◆, ') of the template
waveforms such that they are not searched over.
Because interpolation is carried out in the mass space only, we focus on exploring the e↵ects of interpolation on mass parameters and parameters that are
known to be very strongly correlated with masses (time and phase of coalescenceand distance). If the accuracy of the recovery of these parameters is una↵ectedby interpolation, we do not expect the angular parameters to be a↵ected, either.
However, it is important to realize that the absolute accuracy with which someparameters, particularly distance, are estimated is improved by fixing sky location and orientation parameters and lifting corresponding degeneracies. Thus,
the measurement uncertainties inferred below should not be considered typical.
Since we demand that systematic biases from using interpolated templates are
smaller than statistical measurement uncertainties, the improvement in the accuracy of distance measurement means that we are being conservative in evaluating
the quality of SVD-interpolated parameter estimation.
The signal contained in the data set has source parameters (M, ⌘, D, tc , c ) =(7.45 M , 0.247, 33 Mpc, 0 s, 2.16) and we have omitted the others which are
not searched over. The signal has a signal-to-noise ratio SNR = 14.8. The
frequency-domain data is sampled at f = 1/32 Hz with a maximum frequencyof 512 Hz.
The choice of prior distribution is an important factor in Bayesian inference.
The prior incorporates the ranges of the model parameters, and the a priori
probability distribution on those parameters before the data are analysed. We
use a uniform prior on log D and ⌘ with ranges D 2 [1 Mpc, 100 Mpc] as we haveno a priori preferred distance to the source. We similarly use uniform priors on
0.1 s  tc  0.1 s, respectively.
c and tc over the range 0  0  2⇡ and
We use a prior on M of the form P(M|H) / M 11/6 in the range M 2[7.20 M , 7.60 M ] and a uniform prior on ⌘ in the range ⌘ 2 [0.175, 0.250]. The

63

3. TOWARDS RAPID COMPACT BINARY PARAMETERESTIMATION USING INTERPOLATED WAVEFORMS

prior on M corresponds to the Je↵reys prior for the waveform family describedby Eq. (3.5) [12]. The Je↵reys prior is used to assign equal prior volume to points
in the parameter space, which is necessary when dealing with a space which has
non-trivial curvature.
For the mock data set we ran a MCMC comprising five parallel Markov chainsin order to compute the PDF p(✓~ = (M, ⌘, D, tc , c )|d, H). The limits ofintegration of the likelihood function, Eq. (3.2), are fixed to fmin = 40 Hz, fmax =512 Hz. To extract the posterior samples from the raw MCMC output we discard
the first 10, 000 samples as burn-in.
We measure the convergence of the parallel chains using the Gelman-Rubin
R-statistic [86]. For well converged chains this should be close to R = 1 and we
regard the MCMC to be complete once R  1.001 for all parameters.

SVD Setup

Figure 3.2: Noise PSD used to whiten template waveforms, Eq. (3.13), roughlymatching that of initial LIGO [1].

64The input to the SVD is a set of discretely sampled, whitened time-domain
waveforms [50]. The frequency-domain SPA waveforms are first whitened in thefrequency domain with the ASD and transformed into the time-domain (seeEq. (A.4) for our definition of the discrete Fourier transform and its inverse).
We use a noise PSD roughly matching that of initial LIGO [1] shown in Fig. 3.2.
To improve the interpolation we chose to align the waveforms such that they all
go through frequency at which the noise PSD is a minimum at the same time.
This is achieved by shifting the waveforms by a time tshift in the frequency domain
via the transformationh̃(f ) ) h̃(f )e 2⇡if tshift .
(3.36)The value of tshift is computed from the time is takes for the binary to reachcoalescence from the frequency at which the PSD is a minimum. This time is
known as the “chirp time” and for the SPA waveform family it is given in [87].
The frequency at which the noise PSD is a minimum is fmin = 152.84 Hz.
By carrying out the interpolation in the time domain, we show that the technique can be applied to time-domain waveform families. Time-domain waveforms
are typically computationally expensive for parameter estimation (see Sec. 3.6), so
this approach allows us to assess the computational savings associated with interpolating them. It is also consistent with the work in [49; 50], where time-domain
waveforms generated by Fourier transforming frequency-domain waveforms wereinterpolated.
We ensure that all templates are of the same length, equal to the next highest
power of two of the longest time-domain waveform in the set, which in our caseis 2 s. For the proof of principle we apply the SVD to a small patch in M ⌘space bounding the signal parameters. This region is set by the prior range on
M and ⌘ described above, chosen to be broad enough so that, for our data set,
there is no likelihood support near the boundaries.
The number of computations for the SVD of a N ⇥ L matrix with N  L
scales like O(LN 2 ). For the purposes of constructing the SVD we have found it
efficient to split the M ⌘ space into four equally sized patches, with a separate
SVD applied to each patch. The waveforms in each bank lie on an evenly spaced(15 + 1) ⇥ (15 + 1) grid. We add one to the grid in each dimension because the65

3. TOWARDS RAPID COMPACT BINARY PARAMETERESTIMATION USING INTERPOLATED WAVEFORMS

interpolation can be unpredictable at the boundaries of the space. The grid ischosen by requiring that the normalized inner product between non-interpolated
waveforms and interpolated waveforms generated within each patch is at least
99.9%. Such normalized mismatches of < 0.001 between interpolated and noninterpolated waveforms should ensure that parameter-estimation accuracy is notcompromised as long as the the signal-to-noise ratio does not exceed ⇠ 20 (sothat twice the mismatch times the square of the SNR is less than unity [88; 89]),
although parameter estimation could remain accurate at even higher SNRs. The
patching is shown in Fig. 3.3. Each waveform in the template bank corresponds toa face-on, overhead binary generated at a fiducial distance of 1 Mpc. We choose
to truncate the SVD so that the norm of any reconstructed template is conserved
to a level of ⇠ 10 5 (c.f. Eq. (3.23)). The truncated-SVD of the template bankin each patch uses 20 basis waveforms.

One subtlety of the implementation is that the interpolated waveforms arewhitened time-domain filters. The likelihood function Eq. (3.2) is computed using
un-whitened frequency series which in our example correspond to the Fouriertransform of a time series. To recover the appropriate frequency series we firsthave to Fourier transform the interpolated waveform back into the frequency
domain.

We further note that we do not use waveforms nomalized to unit norm as
input to the SVD, as was done in [49; 50]. Over large areas of the mass-space it isnecessary to normalize the waveforms because there may be significant di↵erencesin the power of waveforms across the space and so the SVD can disproportionallyweight waveforms with more power if they are not normalized to unit norm. We
are not a↵ected by this issue because all the waveforms in our example haveroughly the same power due to the limited extent of the template bank in massspace. In general one would use waveforms normalized to unit norm for SVD
over a larger parameter space; however, to recover the overall amplitude of theinterpolated waveforms the normalization coefficients would themselves have to
be interpolated.

66

Figure 3.3: Patching of the parameter space. The space is patched into four
regions. The SVD is applied to each patch (shown in red) which is padded at theboundary by one waveform and so contains a (15+1)⇥(15+1) grid of template
waveforms. The SVD of the template banks consists of 20 basis templates in eachpatch.

3.6

Results: Comparison of parameter estimatesusing interpolated and non-interpolated waveforms

The marginalized PDFs for complete MCMC runs using non-interpolated andinterpolated waveforms are shown in Fig. 3.4. We have omitted the marginalized
one-dimensional PDF of the coalescence phase c as it of little physical interest.
Each run required around 1.5 ⇥ 106 waveform computations. The mean posteriorvalues of the distributions along with the signal values are shown in Table 3.1.
The chirp-mass distribution computed using interpolated waveforms is clearlybiased. This is corroborated by a two-sample KolmogorovSmirnov (K-S) test [90]which reveals that the two sets of samples are not consistent with arising fromthe same distribution with overwhelming odds. Nevertheless, the systematic bias
in the mean posterior parameter estimate of chirp mass is a factor of four smaller
than the statistical measurement uncertainty. Thus, we pass a commonly-usedthreshold for sufficient waveform-model accuracy [e.g., 89]. We note that theaccuracy could be increased by, for example, using a higher density template bank67

3. TOWARDS RAPID COMPACT BINARY PARAMETERESTIMATION USING INTERPOLATED WAVEFORMS

Figure 3.4: Marginalized PDFs, Eq. (3.4), produced using non-interpolated SPA
waveforms (blue or dark grey) and interpolated SPA waveforms (red or lightgrey). The signal value is shown as the dashed red vertical line.
Param

M (M )⌘
D (Mpc)tc (s)

Mean
posterior value(interpolated
SPA)
7.472
(2.5⇥10 2 )0.2457(7.1⇥10 3 )32.39 (2.11)
1.0 ⇥ 10 3(4⇥10 4 )Meanterior(SPA)

posvalue

7.467(2.5⇥10 2 )0.2457(7.2⇥10 3 )32.40 (2.13)
1.0 ⇥ 10 3(4⇥10 4 )Signal Value

7.450
0.247333.00
0Table 3.1: Mean posterior parameter estimates (and standard deviations) of
the marginalized PDFs using interpolated and non-interpolated SPA waveforms(Fig. 3.4).

68

or using waveforms normalized to unit norm as input to the SVD. In general therequired accuracy can be estimated from the detection trigger SNR [89].
The two-sample K-S test marginally fails for the coalescence-time distribution,
but there is no evidence of a systematic bias on the scale of statistical measurement errors. We find that the sets of posterior samples for the other two PDFs
in Fig. 3.4, symmetric mass ratio and distance, are consistent with arising fromsame distribution as quantified by the K-S test.

Computational cost of template waveform generation
Two commonly-used time-domain waveform families which are relevant for parameter estimation are the inspiral-only post-Newtonian approximant TaylorT4[26] and the e↵ective-one-body family calibrated to numerical relativity (EOBNR,[26]) that includes inspiral, merger, and ringdown phases. The latter are typicallymore computationally intensive.
For technical reasons, our comparison uses interpolated waveforms which are
based on the SPA approximation rather than time-domain waveform families.
This is inconsequential because we Fourier transformed the SPA waveforms into
the time-domain for the SVD procedure. We have observed that applying thetruncated SVD to TaylorT4 time-domain templates in the mass space used toconstruct the SVD in Sec. 3.5 yields the same number of basis vectors as when
applying it to inverse-FFT’d SPA templates. The computational cost of interpolation will therefore be identical.
Our measure to compare the computational costs of interpolated TaylorT4
and EOBNR waveforms is the time it takes to compute a single interpolated
waveform. While this does not estimate the theoretical minimum number of
FLOPs of the process, and is also hardware dependent, it does provide a usefulheuristic for comparing the relative speed of each waveform family. Recall thatthe interpolated waveforms are a time-domain approximant and hence the comparison is to determine the computational savings for time-domain waveforms.
We restrict our comparison to waveforms generated in the mass space in Fig. 3.3.
The length of TaylorT4 and EOBNR waveforms will in general depend on thespecific source masses. For a fair comparison we compare waveforms which have69

3. TOWARDS RAPID COMPACT BINARY PARAMETERESTIMATION USING INTERPOLATED WAVEFORMS

approximately the same number of data points. Because EOBNR must be generated at a sampling rate of 4096 Hz, we ensure that the interpolated and TaylorT4
waveforms are sampled at this frequency. Each waveform is approximately 2 s induration.
The results of the comparison are shown in Table 3.2. For reference we also
show the computational time of standard SPA waveforms. We find that on average, the interpolated waveforms are ten times faster to generate than TaylorT4
and fifteen times faster than EOBNR, a significant increase in computational
efficiency. However, for the waveform parameters considered here, inspiral-only
waveforms could be generated at lower sampling rates than the 4096 Hz requiredfor EOBNR waveforms; therefore, the cost of constructing interpolated or TaylorT4 waveforms can be around four times smaller relative to EOBNR than thevalues quoted in Table 3.2.
Waveform Family
SPA
Interpolated
TaylorT4
EOBNR

Computational Time (T )
0.2
1
10
15

Table 3.2: Computational time of template waveform generation in units of computational time of interpolated waveforms, T. EOBNR, TaylorT4 and interpolated waveform families are generated at a sampling rate of 4096 Hz and havea duration of 2 s. The interpolated waveforms consist of 20 pre-computed basis
vectors. SPA waveforms are generated in the frequency domain; to ensure the
SPA waveforms contain the same number of sample points they are generated ata sampling frequency f = 1/2 Hz and have a maximum frequency of 2048 Hz.

We also estimate the cost of pre-computing the SVD interpolation. We havepreviously noted that the computational cost of an SVD of an N ⇥ L matrix with(N  L) scales like O(N 2 L). One also needs to compute the N ⇥ L matrix oftemplate waveforms as input to the SVD. The cost of computing a waveform oflength L is typically O(L), possibly with a very large pre-factor. Thus, the total
cost of pre-computing interpolation coefficients will be less than O(N 2 ) times
the cost of an individual waveform computation. For instance, in our example,
N = 16 ⇥ 16 = 256, so interpolation can reduce overall MCMC costs for any

70

time-domain waveform templates by an order of magnitude or more when the
typical MCMC chain length of & 106 samples is taken into account.

3.7

Discussion and conclusion

We have provided a proof of principle that interpolated waveforms can be used
for parameter estimation without unacceptable loss of accuracy. Our example
was restricted to a five-dimensional search over the source chirp mass M andsymmetric mass ratio ⌘, the distance to the source D and the time and phase atcoalescence tc and c . We further restricted the prior ranges on M and ⌘ to M 2[7.20 M , 7.60 M ] and ⌘ 2 [0.175, 0.250], respectively. The systematic biases
observed when using interpolated waveforms are demonstrated to be smaller thanstatistical measurement uncertainties. Thus, SVD-interpolated waveforms satisfy
the stringent waveform-model accuracy criteria imposed by parameter-estimation
requirements.
The relative computational times of generating interpolated waveforms andtime-domain TaylorT4 and EOBNR waveforms are shown in Table 3.2. Interpolated waveforms can be generated at around an order of magnitude more cheaply
than TaylorT4 or EOBNR. This suggests that the computational cost of parameter estimation can be significantly reduced by employing SVD-interpolated
waveforms for likelihood computation when the latter is dominated by the costof waveform generation.
In order for interpolated templates to be viable for parameter estimation
pipelines we need to apply the SVD-interpolation technique to a significantly
larger region of the CBC mass space than in the example considered here. Searchesof gravitational waves from low-mass systems look for binaries with a maximum
total mass of 25 M and a minimum component mass of 1 M [7] and high mass
searches target binaries with total masses between 25 M and 100 M [8]. To beable to apply our parameter estimation technique to triggers from such searches
in a single step, without first determining the more limited mass region where
there is significant likelihood support, we will need to efficiently patch the parameter space over a large mass range so that the computational cost of generating
the SVD can be minimized. A necessary condition for setting up the SVD in all71

3. TOWARDS RAPID COMPACT BINARY PARAMETERESTIMATION USING INTERPOLATED WAVEFORMS

patches is that its computational cost, plus the cost of generating interpolated
waveforms, must be less than the cost of performing the parameter estimation
with non-interpolated waveforms. This will be the subject of future work.
For simplicity, we have considered analysing a single-detector data set. If
one were to extend the use of interpolated waveforms to analysing data from a
network of interferometers, one would require independent sets of interpolated
waveforms for each data set of the network. This is because our technique relies
on a SVD of template waveforms, whitened by the detector noise PSD. Because
each detector in general has a di↵erent PSD, one would thus require separate setsof whitened filters, each whitened with the appropriate noise PSD. For a typicaldetector network of three interferometers, one would thus need to generate three
template banks and perform three SVDs. While this increases the computationalcost, we note that computing independent sets of template banks and SVDs ishighly parallelizable, and hence the wall-time of computing three banks could
easily be made the same as computing a single template bank.
Furthermore, we have to be able to extend the SVD to generic waveform families. Particularly interesting is the potential to extend the technique to EOBNR
waveforms, which are currently expensive to generate, and waveform families
which describe binaries with spinning components. The latter class of waveforms
have an intrinsic parameter space with up to six more independent parameters(two spin vectors) and the current technique of interpolation within the intrinsic parameter space of waveforms may become costly in large-dimensional spaces.
However, it is interesting to consider the potential to apply the technique to spinaligned/anti-aligned waveforms [e.g., 87; 91] as this class of waveforms have only
one extra parameter, the reduced spin of the binary.The analysis of data fromadvanced LIGO and Virgo detectors, which may have lower-frequency cuto↵s
close to 10 Hz [58], will require template waveforms that are several minutes induration. This technique is likely to be highly relevant to parameter estimationin that context.

72

Chapter 4
PARAMETER ESTIMATIONUSING INTERPOLATEDTIME-DOMAIN TEMPLATE
WAVEFORMS
4.1

Introduction

Techniques to improve the efficiency of Bayesian parameter estimation are crucial
in the run-up to the advanced-detector era. Interpolated gravitational waveforms
based on the singular value decomposition have the potential to significantly increase the efficiency when parameter estimation is dominated by waveform computation.
For interpolated waveforms to be viable for parameter estimation on real data,
a number of extensions of the technique demonstrated in Chapter III are required.
While it has been demonstrated that interpolated waveforms can be used for parameter estimation without incurring significant bias in parameter estimates, theproof-of-principle study was limited to one waveform family and a small patch ofparameter space. We begin this chapter by investigating the extension to genericwaveform families which describe CBC sources with non-spinning components. In
Chapter III we interpolated a waveform family constructed by Fourier transform-73

4. PARAMETER ESTIMATION USING INTERPOLATEDTIME-DOMAIN TEMPLATE WAVEFORMS

ing frequency-domain SPA waveforms into the time-domain, which we used as
templates for parameter estimation. We demonstrated that because the interpolation technique is the same regardless of the waveform family being interpolated
one could expect to reduce the computational cost of generating time-domain
waveforms, e.g. TaylorT4 and EOBNR [26], via interpolation.
It is occasionally useful for parameter estimation to utilize time-domain waveforms. For example, the NINJA project [92], and its sequel (NINJA2 [93]) currently underway, is an exercise in detecting and performing parameter estimationon (time-domain) waveforms generated within the framework of numerical relativity. Here EOBNR waveforms should be useful as a template waveform family
because they have been matched to similar numerical waveforms. One of the reasons we did not initially apply the interpolation technique in Chapter III directly
to time-domain waveforms was due to a number of potential issues which could
complicate the interpolation. For example, because the waveforms are whitenedin the frequency domain one first has to Fourier transform a time-domain waveform into the frequency domain to perform the whitening process. This will
introduce artifacts into the Fourier series which can compromise the accuracy of
the SVD, and hence interpolated waveforms. In addition, EOBNR waveforms
may be difficult to interpolate because of the relative time-scales of the inspiralphase and the ringdown.
We first demonstrate that parameter estimation can be performed using interpolated TaylorT4 waveforms, which are inspiral-only time-domain waveforms.
We find that while there are clear biases in the recovered posterior probability density functions, the bias of the mean posterior parameter estimates aremuch less than the scale of the statistical uncertainty and so we pass a common
threshold for model waveform accuracy [89]. Furthermore, the improvement incomputational time of generating interpolated TaylorT4 waveforms is around afactor of 10 with respect to non-interpolated TaylorT4 waveforms. In general,
we find that applying the SVD to time-domain waveforms which describe thefull inspiral, merger and ringdown phases of coalescence is challenging and thequality of interpolated EOBNR waveforms is too low for parameter estimation,
but could suffice for gravitational-wave searches.

74

4.2

Interpolating inspiral-only time-domain waveforms

To demonstrate the subtleties of interpolating time-domain waveforms we firstconsider a simple inspiral-only time-domain waveform family known as TaylorT4.
We refer the reader to Sec. 2.3.3 for details on the construction of the TaylorT4
waveform family. Spectral leakage is a general feature of Fourier transforming atime-series which has discontinuities from a hard beginning and ending. (Such atime-series can be thought of as the convolution of an infinite, smooth time-series
with a box window). Thus, the Fourier transform of the convoluted time-series
will contain additional high- and low-frequency components from the box window,
which is a linear superposition of oscillating functions. This can be ameliorated
by generating the time series at a lower sampling frequency which can move manyof the artifacts to a region in which one is not interested.
There is also another problem associated with using time-domain waveforms
which is symptomatic of our implementation. Inspiral-only time-domain waveforms are integrated from an initial frequency up to the ISCO frequency. Timedomain waveforms generated within the LIGO Scientific Collaboration Algorithm
Library (LAL) are aligned such that the phase at last time-sample of di↵erent
waveforms is the same. However, because the time series are discretely sampled,the frequency at the last time sample is in general not the same and is a reflectionof inaccurate time-frequency evolution of the waveforms. This may be an issue forinterpolating time-domain waveforms because the waveforms have to be aligned
in frequency domain with the minimum of the noise PSD and so there may be
inconsistencies in the alignment of the waveforms for the SVD. This can be ameliorated by noting that in the limit that the sampling frequency is infinitely fine,the frequency at the last sample will tend to the same value. Hence at higher
sampling frequencies we may expect the time-frequency evolution of time-domain
waveforms to be more accurate.
To study the e↵ects of discrete sampling on interpolated waveforms, we compute the SVD of the waveforms over a range of sampling rates. As a measureof the quality of the interpolated waveforms we compute the inner products between interpolated and non-interpolated TaylorT4 waveforms on a uniform grid.

75

4. PARAMETER ESTIMATION USING INTERPOLATEDTIME-DOMAIN TEMPLATE WAVEFORMS

We chose to work with the “mismatch” between (normalized) interpolated andnon-interpolated waveforms which is simply
MM = 1(hinterp |hnon interp ) ,(4.1)

and is zero when hinterp = hnon interp .
Recall that the input to the SVD is a set of whitened waveforms, in this casewaveforms described by the TaylorT4 family. The waveforms are first whitened in
the frequency domain with the square root of the noise PSD and transformed backinto the time-domain for interpolation. We use the same noise PSD as in Chapter
III which is typical of initial LIGO [1] and is shown in Fig. 3.2. The waveforms
are generated from a starting frequency of 30 Hz in the time-domain which we use
to mitigate the e↵ects of spectral leakage in the Fourier-domain waveform when
whitening, and are truncated at a low frequency cut-o↵ of 40 Hz before backtransforming into the time-domain for the SVD. To improve the interpolation wechose to align the waveforms such that they all go through frequency at which
the noise PSD is a minimum at the same time. The frequency at which the noise
PSD is a minimum is fmin = 152.844 Hz. This is achieved by time-shifting the
waveforms in the frequency domain by an amount tshift which can be calculated
from the time is takes for the binary to reach coalescence from fmin . This time iscalculated by computing Eq. (A.14) to find the frequency evolution of the binary
as a function of time.
To compute the SVD, we work within the region of parameter space usedin Chapter 3, which is defined by 7.2 M  M  7.6 M , 0.175  ⌘  0.250.
We divide the parameter space into four equally sized patches as in Fig. 3.3.
In each patch we construct a template bank of (25 + 1) ⇥ (25 + 1) whitened
waveforms on a uniform grid. In Fig. 4.1 we histogram the set of mismatches
between interpolated and non-interpolated waveforms, for waveforms generated
at three di↵erent sampling rates: 512 Hz, 1024 Hz, 2048 Hz and 4096 Hz. The
SVD in each case uses 20 basis templates per patch of parameter space.
To compute the mismatch between interpolated and non-interpolated TaylorT4 waveforms we generate the waveforms on a uniform 50 ⇥ 50 grid in theregion of parameter space defined by 7.2 M  M  7.6 M , 0.175  ⌘  0.250.

76

Figure 4.1: Histogram of mismatches between (normalized) interpolated and noninterpolated TaylorT4 waveforms generated on a uniform (50 ⇥ 50) grid in apatch of (M, ⌘) defined by 7.2 M  M  7.6 M , 0.175  ⌘  0.250. Eachpanel corresponds to waveforms generated at di↵erent sampling rates: 512 Hz,
1024 Hz, 2048 Hz and 4096 Hz. We set a threshold mismatch of 10 3 between
interpolated and non-interpolated waveforms, shown as a red vertical line, above
which we deem interpolated waveforms to be insufficiently accurate for parameter
estimation. The best case occurs for a sampling rate of 4096 Hz for which the
majority of mismatches are less than the threshold.

77

4. PARAMETER ESTIMATION USING INTERPOLATEDTIME-DOMAIN TEMPLATE WAVEFORMS

We set a threshold mismatch of 10 3 above which the accuracy of interpolated
waveforms may seriously compromise parameter estimation accuracy [89]. We
find that the most accurate interpolated waveforms are those generated at asampling frequency of 4096 Hz for which the vast majority of waveforms are
more accurate than the threshold value. We note that there are around ⇠ 20
waveforms which are less accurate than the threshold, however this is unlikely tohave a large impact over the bulk of the parameter space. We further note that
increasing the number of templates in the template bank on which we performthe SVD has almost no e↵ect in improving the quality of interpolated templates.
Below we demonstrate the use of interpolated TaylorT4 waveforms as templates
for parameter estimation.

4.2.1

Parameter estimation using interpolated TaylorT4
waveforms: 5-dimensional case

Our goal is to estimate the parameters of a gravitational-wave signal containedin Gaussian stationary noise using both the interpolated and non-interpolated
TaylorT4 waveform family as templates. The gravitational-wave signal containedin the data is described by the non-interpolated TaylorT4 waveform family. We
generate waveforms in the time-domain at a sampling frequency of t = 1/4096 s.
The signal contained in the data set has source parameters (M, ⌘, D, tc , c ) =(7.45 M , 0.247, 33 Mpc, 0 s, 2.16) and we have omitted the others which are
not searched over. The signal has a signal-to-noise ratio SNR = 14.8. We use thesame prior distributions as in Sec. 3.5. To recap, the prior distributions are set
as follows.
4.2.1.1

Prior distributions

We use a uniform prior on log D and ⌘ with ranges D 2 [1 Mpc, 100 Mpc] as wehave no a priori preferred distance to the source. We similarly use uniform priors
on c and tc over the range 0  0  2⇡ and 0.1 s  tc  0.1 s, respectively.
We use a prior on M of the form P(M|H) / M 11/6 in the range M 2[7.20 M , 7.60 M ] and a uniform prior on ⌘ in the range ⌘ 2 [0.175, 0.250]. The

78

prior on M corresponds to the Je↵reys prior for the waveform family describedby Eq. (3.5) [12].
For the mock data set we ran a MCMC comprising five parallel Markov chainsin order to compute the PDF p(✓~ = (M, ⌘, D, tc , c )|d, H). The limits ofintegration of the likelihood function, Eq. (3.2), are fixed to fmin = 40 Hz, fmax =2048 Hz. The upper limit of integration is the nyquist frequency of waveforms
sampled in the time-domain at t = 1/4096 s. To extract the posterior samples
from the raw MCMC output we discard the first 10, 000 samples as burn-in.
We measure the convergence of the parallel chains using the Gelman-Rubin
R-statistic [86]. For well converged chains this should be close to R = 1 and we
regard the MCMC to be complete once R  1.001 for all parameters.
4.2.1.2

SVD set up

To generate interpolated template waveforms we take the SVD of a templatebank of whitened TaylorT4 waveforms over a region of parameter space definedby our prior range in M and ⌘: M 2 [7.20 M , 7.60 M ], ⌘ 2 [0.175, 0.250]. We
patch the parameter space into four equally sized regions as in Fig. 3.3 in which
we generate a separate template bank. The waveforms in the template banks are
generated at a sampling frequency of t = 1/4096 s. Each template bank uses
(25 + 1) ⇥ (25 + 1) waveforms.
4.2.1.3

Results

The marginalized PDFs for complete MCMC runs using non-interpolated andinterpolated TaylorT4 waveforms are shown in Fig. 4.2. We have again omitted
the marginalized one-dimensional PDF of the coalescence phase c as it of little
physical interest. Each run required around 1 ⇥ 106 waveform computations to
be fully converged. The mean posterior values of the distributions along with thesignal values are shown in Table 4.1.
As with interpolated stationary phase approximation (SPA) waveforms, thechirp-mass distribution computed using interpolated TaylorT4 waveforms is clearlybiased. However it is encouraging that the systematic bias in the mean posterior
chirp mass estimate is around a factor of forty smaller than the statistical mea-

79

4. PARAMETER ESTIMATION USING INTERPOLATEDTIME-DOMAIN TEMPLATE WAVEFORMS

Figure 4.2: Marginalized PDFs, Eq. (3.4), produced using non-interpolated TaylorT4 waveforms (blue or dark grey) and interpolated TaylorT4 waveforms (redor light grey). The signal value is shown as the dashed red vertical line.

80ParamM (M )⌘
D (Mpc)tc (s)

Mean
posterior value(interpolated
TaylorT4)
7.424(3.7⇥10 2 )0.236(1.3⇥10 2 )31.03 (2.4)
1.0 ⇥ 10 3(4⇥10 4 )Meanposterior value(TaylorT4)

Signal Value

7.425(3.3⇥10 2 )0.236(1.2⇥10 2 )31.09 (2.6)
1.0 ⇥ 10 3(4⇥10 4 )

7.450
0.247
33.00
0

Table 4.1: Mean posterior value parameter estimates (and standard deviations) of
the marginalized PDFs using interpolated and non-interpolated TaylorT4 waveforms (Fig. 4.2).

surement uncertainty. Again we pass a commonly-used threshold for sufficient
waveform-model accuracy [e.g., 89] and at a higher threshold than with interpolated SPA waveforms. In addition, the time distribution is also clearly biased,
however as with the chirp mass, there is no bias of the mean posterior time parameter estimate on the scale of statistical measurements errors. This is also thecase with the other two parameters, distance and symmetric mass ratio. We findthat none of the sets of posterior samples for the PDFs in Fig. 3.4 are consistent
with arising from same distribution as quantified by the K-S test.
We find that the reduction in computational time in using interpolated TaylorT4 waveforms is around a factor of 10, which is consistent with Table 3.2.

4.3

Interpolating inspiral-merger-ringdown waveforms

We now extend the interpolation technique to waveforms which describe the fullinspiral, merger and ringdown phases of coalescence. We take as our model timedomain waveform EOBNR which was used in Sec. 2.3 and is described in detailin [26]. Unlike TaylorT4 the EOBNR family of waveforms must be generated

81

4. PARAMETER ESTIMATION USING INTERPOLATEDTIME-DOMAIN TEMPLATE WAVEFORMS

at a minimum sampling rate of 4096Hz in order to accommodate the ringdown
phase. We use the implementation provided by the LIGO Scientific Collaboration
Algorithm Library (LAL) that corresponds to the approximant EOBNRv2 [41].
An clear subtlety in applying the SVD to EOBNR is the alignment of the
waveforms. Unlike SPA or TaylorT4 waveforms, there is no well defined end of
EOBNR waveforms. Because of this we chose to align the peak amplitude of
the waveforms with the minimum of the noise PSD in the frequency-domain.
We use the same noise PSD as in Chapter III which is typical of initial LIGO
[1] and is shown in Fig. 3.2. In order to perform the alignment we first haveto compute the time-frequency evolution associated with EOBNR waveforms.
This is computed within the waveform generation routine within LAL for theapproximant EOBNRv2 [41].We again take the reference frequency to be fmin =152.84375 Hz, which is the frequency at which the PSD is a minimum.
To determine the quality of interpolated EOBNR waveforms we compute themismatch , Eq. (4.1), between interpolated and non-interpolated EOBNR waveforms. We apply the SVD to a set of whitened EOBNR waveforms generated in
the same region of parameter space shown in Fig. 3.3 defined by 7.2 M  M 
7.6 M , 0.175  ⌘  0.250. We divide the parameter space into four equally sized
patches as in Fig. 3.3. As input to the SVD, we use (35 + 1) ⇥ (35 + 1) whitened
waveforms in each patch’s template bank. The waveforms are generated froma starting frequency of 30 Hz in the time-domain which we use to mitigate thee↵ects of spectral leakage in the Fourier-domain waveform when whitening, andare truncated at a low frequency cut-o↵ of 40Hz before back-transforming into thetime-domain for the SVD. The SVD of each template bank uses 22 basis waveforms. To compute the mismatches we use a uniform grid of 50⇥50 waveforms
over the parameter space defined by 7.2 M  M  7.6 M , 0.175  ⌘  0.250.
We compute the mismatch for waveforms generated at the following sampling
frequencies: 4096 Hz, 8192 Hz, 16384 Hz and 32768 Hz.
We observe that for the range of sampling frequencies the mismatches have awide distribution about the threshold mismatch of 10 3 as is shown in Fig. 4.3.
These mismatches suggest that the interpolated EOBNR waveforms are not suitable for parameter estimation as they are not sufficiently accurate. The waveforms we have considered have a relatively long (⇠ 2 s) inspiral followed by a82

Figure 4.3: Histogram of mismatches between (normalized) interpolated and noninterpolated EOBNR waveforms generated on a uniform (50 ⇥ 50) grid in a patchof (M, ⌘) defined by 7.2 M  M  7.6 M , 0.175  ⌘  0.250. Each panel
corresponds to waveforms generated at di↵erent sampling frequencies: 4096 Hz,
8192 Hz, 16384 Hz and 32768 Hz. We set a threshold mismatch of 10 3 between
interpolated and non-interpolated waveforms, shown as a red vertical line, above
which we deem interpolated waveforms to be insufficiently accurate for parameter
estimation. The best case occurs for sampling rate of 32768 Hz, though it is nearly
identical to the 16384 Hz case. In all cases there is a wide distribution about
the threshold mismatch of 10 3 which indicates that these interpolated EOBNR
waveforms are not suitable for parameter estimation. We note that while thesemismatches indicate that the interpolated EOBNR waveforms are not suitablefor parameter estimation, they would be suitable to form a template bank forgravitational-wave searches.

83

4. PARAMETER ESTIMATION USING INTERPOLATEDTIME-DOMAIN TEMPLATE WAVEFORMS

short (⇠ 10 2 ) ringdown. We note that increasing the number of basis vectors
used for interpolation, and increasing the number of waveforms in the templatebank on which we perform the SVD yields almost no improvement in the qualityof interpolated waveforms.
While this is a set back for parameter estimation, there are alternative approaches to parameter estimation in which arbitrary waveforms could be usedwhile still significantly decreasing its computational cost. We discuss these in thefollowing chapter.
We note that while the interpolated EOBNR waveforms are insufficient forparameter estimation, they are sufficient to form a template bank for detection
where the accuracy requirements are typically lower. For example, gravitational
wave template banks are constructed such that neighbouring templates have amismatch of around 3%.

4.4

Future work and Discussion

We have considered in detail the potential to use interpolated waveforms as one
means to do this. For inspiral only waveforms we have found that the using interpolated waveforms leads to biases in the posterior probability distributions of
the parameters. However, the bias in the mean posterior parameter estimates are
negligable on the scale of statistical measurement uncertainties. In general, parameter estimation on real gravitational-wave signals from CBC sources requiresthat we are also able to use waveforms which describe the full inspiral merger andringdown phases of coalescence. Our attempt to interpolate such waveforms has
failed in this sense because the quality of the interpolated waveforms is insufficient for parameter estimation, though in principle could be used for searches. In
the following chapter we consider how to circumvent the problems encounteredin interpolating time-domain template waveforms by directly interpolating the
likelihood function, without first computing an intermediate template waveform.

84

Chapter 5
RAPIDLY EVALUATING THE
COMPACT BINARY
LIKELIHOOD FUNCTION VIAINTERPOLATION
5.1

Introduction

In the previous two chapters we considered using interpolated template waveforms as a means to reduce the computational cost of parameter estimation.
While the technique is promising in principle, we found that our implementationof interpolated waveforms is somewhat limited as we were unable to successfullyinterpolate time-domain waveforms which describe the full inspiral, merger andringdown phases of compact binary coalescence. While there may be other, morerobust, interpolation techniques which would allow generic template waveforms
to be interpolated, we will not pursue them here. Instead, we will consider another means with which one could significantly reduce the computational cost ofparameter estimation, again based on interpolation using the SVD.
We demonstrate that the SVD can be applied directly to the likelihood function itself which in turn allows the likelihood function to be directly interpolated.
With this formulation of the likelihood function the only waveform computations85

5. RAPIDLY EVALUATING THE COMPACT BINARY
LIKELIHOOD FUNCTION VIA INTERPOLATION

required for parameter estimation are those needed to construct the SVD of thelikelihood function. We find that for a test-case, we are able to speed up theevaluation of the likelihood function by two orders of magnitude, whilst keeping
likelihoods accurate to . 0.025%. In addition, we will demonstrate that this improvement in computational time is conservative, and we discuss how additional
computational savings would be achieved in practice. Crucially, the method wewill present is totally transparent to the form of the template waveforms used forthe likelihood computation.
Recall that the goal of Bayesian parameter estimation is to compute the pos~ which underlie
terior probability density function (PDF) of a set of parameters, ✓,a model assumed to describe a data set d. The PDF is related to the likelihood
function and prior probability via Bayes’ theorem:~ =p(✓|d)~ L(d|✓)
~
P(✓)
,
p(d)

(5.1)

~ is the likelihood and P(✓)
~ is the prior probability which encodes ourwhere L(d|✓)~ The quantity in the denominator, p(d), isa priori belief in the distribution of ✓.
known as the “evidence”. Computing (5.1) requires evaluating the likelihood.
For binaries with non-spinning components ✓~ is nine dimensional. Exploringsuch a high dimensional space requires sophisticated stochastic Bayesian inferencetechniques [12; 13; 25] which preferentially sample the parameter space in regionsof high posterior probability. The bulk of the computational cost of evaluating the
likelihood function comes from computing template waveforms. Analyses on firstgeneration interferometer data require computing O(106 ) such waveforms [25; 51].
Sampling techniques such as Markov chain Monte Carlo (MCMC) [13; 25] andnested sampling [12; 94] evaluate likelihoods, and hence compute template waveforms, serially. Thus the total computational time to fully sample the parameter
space scales linearly with the total time spent generating template waveforms. It
can take hours to weeks to analyse a single stretch of data of a few seconds induration, depending on the choice of the template waveform family. This problem
will be exacerbated when analysing second-generation interferometer data as the
waveforms will be forty times longer in duration if the starting frequency fmin ischanged from 40 Hz to 10 Hz.

86

For binaries with non-spinning components, the frequency-domain waveformh̃(✓~ ; f ) has the schematic formh̃(✓~ ; f ) =XAµ (✓~L )h̃0 (M , ⌘; f )e2⇡if tc ,(5.2)µ=+ ,⇥

where A+ ,⇥ denotes the (scalar) amplitudes of the “plus-” and “cross-” polarization states of the waveform. In general h̃0 depends on the waveform family being
used and can be computed by Fourier transforming the associated time-domain
representation of the waveform family. The parameters which describe the binary
are the chirp mass and symmetric mass-ratio, M and ⌘, the time at coalescence tc
and a set of parameters which describe the location and orientation of the binary✓~L .

Evaluating the likelihood function on the three-dimensional subspace of parameters (M , ⌘ , tc ) represents the largest computational burden to parameter estimation on gravitational waves from CBC sources with non-spinning components
because the likelihood function depends non-trivially on these parameters, and sorequires a new waveform evaluation. In [51], we considered interpolation between
waveforms over the mass parameter space as a way to reduce computational cost.
Here, we demonstrate that the evaluation of an interpolated likelihood function
over the (M , ⌘ , tc ) subspace is a much faster computational procedure than thestandard calculation of the likelihood (5.3) by using either full or interpolated
waveforms. For the purposes of parameter estimation, one is not interested intemplate waveforms per se, but rather in the posterior probability distributionsof the underlying parameters of the template waveforms that are assumed todescribe the data. By directly using interpolated likelihood functions, one e↵ectively bypasses template waveform generation during the sequential steps of an
MCMC. This likelihood-interpolation technique is robust and could, in principle,
be generalized to arbitrary template waveform families, in particular those that
describe CBCs with spinning components.

87

5. RAPIDLY EVALUATING THE COMPACT BINARY
LIKELIHOOD FUNCTION VIA INTERPOLATION

5.1.1

Directly interpolating the likelihood function

We wish to generate a representation of the likelihood function over the continuous M, ⌘ and tc subspace. To achieve this we will interpolate the likelihood
function over M, ⌘ and tc . The likelihood function that describes the probabilityof observing a data stream d = h + n containing a given gravitational-wave signal~ t) and Gaussian and stationary noise n(t) is [12]
h(✓;
i
1h ~
~ + (d|d) ,h(✓)|h(✓))
2~ = (d|h(✓))~
log L(d|✓)(5.3)

where (a|b) is the usual noise-weighted inner product [46]. We define the complexvalued time-series corresponding to the inner product between two time series a(t)and b(t) as one is shifted by an amount tc with respect to the other:z[a , b](tc ) := 4Z fmaxdf

fminã(f )b̃⇤ (f ) 2⇡if tce,
Sn (f )

(5.4)

In the above, ã(f ) is the Fourier transform of a(t) and Sn (f ) is the detector noise
power spectral density (PSD). The limits of integration are in general specified
by the bandwidth over which an analysis is being conducted. In terms of z(tc )the inner products in (5.3) are succinctly expressed as~~ = <A(✓~L )z[h0 (M, ⌘), h0 (M, ⌘)](0),(h(✓)|h(✓))~ = <B(✓~L )z[d , h0 (M, ⌘)](tc ) ,(d|h(✓))

(5.5)
(5.6)

and A(✓~L ) and B(✓~L ) are known projections which contain the ✓~L dependence inthe likelihood function.
We have previously interpolated template waveforms over the mass parameters[50; 51], and here we show that the same technique can be applied to interpolating
the time series z[d , h0 ](tc ). The interpolation of z[d , h0 ](tc ) is based on the SVD
of a set of (discretely sampled) time series distributed on a two-dimensional grid.
In this case the two-dimensional grid spans M and ⌘ and the time parameter is
tc . We use the notation ~z [d , h0 ] to describe the discretely sampled z[d , h0 ](tc ).
Recall that the SVD of the discretely sampled time series ~z [d , h0 ] allows it to be
written as a linear superposition of orthonormal basis vectors ~uµ and projection88coefficients Mµ [49]:~z [d , h0 (M, ⌘)] =X

Mµ (M , ⌘) ~uµ .

(5.7)

µ

The coefficients Mµ can be interpolated over M and ⌘ and we follow the method
in [50] which uses Chebyshev polynomials of the first kind.
Interpolation of z[h0 , h0 ](0) over M and ⌘ is simple as it is scalar valued
and we again use Chebyshev polynomials of the first kind. Below we provide an
example of the interpolation technique outlined here.

5.1.2

Likelihood interpolation: Examples

We compare interpolated likelihood functions to those generated by direct evaluation of waveforms and inner products. We consider two test cases: (i) the
coalescence of binary black holes, and (ii) the coalescence of binary neutron stars.
We generate a discretely sampled, simulated data set d~ for a single interferometer consisting of Gaussian and stationary noise ~n and a gravitational-wave
signal ~h. The data set is 32 s in duration and has a sampling rate in the time
domain of 4096 Hz. For binary black holes we model the gravitational-wave signal ~h using the e↵ective one-body approach calibrated to numerical relativity
simulations (EOBNR) [41]. Such a gravitational-wave signal describes the fullinspiral, merger and ringdown phases of coalescence. For binary neutron stars wemodel the gravitational waveform using a post-Newtonian (PN) model computed
to 3.5 PN order in phase [26], which describes the inspiral phase of the coalescence only. We use an implementation of EOBNR and post-Newtonian waveforms
from the LSC Algorithms Library (LAL) [85] corresponding to the approximants
EOBNRv2 and TaylorT4 respectively.
Generating the interpolant of the likelihood function requires the following
stages: (i) patch the mass space into smaller domains, (ii) generate a set ofwaveforms over a dense grid in each patch, (iii) filter the data with the template
waveforms to compute the likelihoods, (iv) pack the likelihoods into a matrix andperform the SVD, (v) build the interpolant in each patch. Only after these stages
have been completed can the interpolated likelihood function be sampled.

89

5. RAPIDLY EVALUATING THE COMPACT BINARY
LIKELIHOOD FUNCTION VIA INTERPOLATION

We first construct a discrete, uniform grid of template waveforms in M ⌘parameter space. We will use a small region around the parameters of the signal,
as M and ⌘ are typically constrained to . 1% and . 10%, respectively, dependingon the signal parameters and signal-to-noise ratio (SNR) [51; 95]. The region in
M ⌘ where the posterior has significant support can be found quickly during theburn-in phase of the MCMC, which requires a small fraction of the total numberof samples necessary to evaluate the posterior probability distribution function.
We use the Chebyshev interpolation described in [50] to interpolate z[h0 , h0 ](0)for waveforms across the grid. To interpolate ~z [d , h0 ] we first find the basis vectors ~uµ by constructing a matrix from the set of {~z [d , h0 ]}, the columns of whichcorrespond to a unique ~z [d , h0 ] on the grid of waveforms, which we factor using
the SVD. After performing the SVD, we truncate the number of basis vectors
such that on average the norm of each ~z is conserved to one part in 105 [50].This
can significantly reduce the number of basis vectors.We then apply the Chebyshev
interpolation [50] to interpolate projection coefficients across the M ⌘ grid.
5.1.2.1

Example 1: High-mass binary black holes

The signal is parameterized by ✓~s = (M = 15.01 M , ⌘ = 0.205 , D = 100 Mpc , ◆ =0 , = 0 , ↵ = 0 , = 0 , tc = 0.1 s , c = 0). We use a noise PSD typical of initial
LIGO [1]. The signal has an SNR of ⇡ 15. In order to interpolate the likelihood
function across M , ⌘ and tc , we work within a small region of M ⌘ space whoseboundaries are given by 14.56 M  M  15.46 M and 0.143  ⌘  0.25.
Assuming a statistical measurement uncertainty on M and ⌘ of 1% and 10%,respectively, the parameter ranges correspond to a ⇠ 3 range about the signal
value. Note that we cannot go about ⌘ = 0.25 in the ⌘ interval. We further
restrict our range in tc to be in a ±0.2 s window about the trigger time, which is
a common time prior in Bayesian parameter estimation [12].
In Fig. (5.1) we compare a likelihood function generated via direct evaluationof inner products, to one which we have generated via SVD-interpolation. We
find that we are able to reconstruct the log likelihood function by interpolation
to within a fractional percentage error of at most 0.025 %. While we have only
plotted an interpolated likelihood function at the signal values of M and ⌘, the90

Figure 5.1: Interpolated and non-interpolated log likelihoods (top), and percentage error (bottom) for a data set containing a gravitational-wave signal from thecoalescence of binary black holes.

91

5. RAPIDLY EVALUATING THE COMPACT BINARY
LIKELIHOOD FUNCTION VIA INTERPOLATION

errors quoted here are typical across the mass range we have considered. Meanwhile, for this waveform model and parameters, computing the likelihood via theinterpolation procedure is around two orders of magnitude faster than generatinga template waveform and directly evaluating the inner products in (5.3).
5.1.2.2

Example 2: Binary neutron stars

The signal is parameterized by ✓~s = (M = 1.217 M , ⌘ = 0.2497 , D = 20 Mpc , ◆ =0 , = 0 , ↵ = 0 , = 0 , tc = 0.1 s , c = 0). We again use a noise PSD typicalof initial LIGO [1], and the signal has SNR ⇡ 15. We interpolate the likelihood function over a small region of M ⌘ space whose boundaries are given
by 1.199 M  M  1.235 M and 0.212  ⌘  0.25. Assuming a statisticalmeasurement uncertainty of 0.5% on M and 5% on ⌘, these parameter ranges
correspond to a ⇠ 3 range about the signal value. Note that we cannot go about⌘ = 0.25 in the ⌘ interval. We restrict our range in tc to be in a ±0.2 s window
about the trigger time.
Again, we find that we are able to reconstruct the log likelihood function
to within a fractional percentage error of at most 0.025 %. For the binary neutron star case, we find that computing the likelihood via interpolation is aroundthree orders of magnitude faster than direct evaluation. This di↵erence is larger
than for the higher-mass binary black hole case because the waveform duration
is significantly longer for binary neutron stars, whereas the cost of computinginterpolated likelihoods remains fixed.
Below we discuss practical issues pertaining to incorporating interpolated likelihoods into real gravitational-wave parameter-estimation pipelines.
5.1.2.3

Practical considerations

For our interpolation technique to be viable for real data analyses, the total
computational time of first constructing the interpolated likelihood function, andthen sequentially sampling the interpolated likelihood function, must be less than
the time for sequentially sampling the likelihood function directly.
Parameter estimation requires sampling the parameter space until the sampler
has met its convergence criterion. The total number of likelihood evaluations

92

for convergence is typically ⇠ O(106 ) [25; 51]. When directly evaluating the
likelihood function, the number of likelihood evaluations is a reasonable proxy
for the number of waveform evaluations, which dominate the computational cost.
To sample the interpolated likelihood function there is an additional upfront
cost of constructing the interpolant of the likelihood function. This cost will
depend on the region of the parameter space over which the likelihood functionneeds to be interpolated and template waveforms must be computed. However,building the interpolant is highly parallelizable and computing it over an extended
region of parameter space could be split into multiple independent subsets. This
could greatly reduce the wall time of computing the interpolant. We have notedthat one can restrict the range in parameter space over which the interpolant isbuilt by using an MCMC to sparsely explore the parameter space in regions ofhigh posterior probability. In practice, the number of samples for this “burn-in”
is often ⇠ O(104 ) [51], and the likelihood has significant support in a relatively
small patch in parameter space. The likelihoods computed during the burn-in
evaluation could thus be stored for future interpolation.
One could also interpolate ~z (tc ) over many patches covering the parameter
space in parallel. We have not investigated optimal patching, nor the required
denseness of likelihood template calculation in order to generate a good basis for~z (tc ); this will be the subject of future work.
Once the interpolant is constructed, the cost of sampling the parameter space
will depend on that of computing the interpolated likelihood function. In our
example we found that computing the interpolated likelihood function is between
two and three orders of magnitude cheaper than directly evaluating the likelihood
function, depending on the region in parameter space in which the likelihood
function is being computed. The actual improvement will depend on the typicalcost of waveform computation, which is a function of both the template waveformfamily used and the waveform parameters.
Here we had used the SVD to find a basis for the set of ~z (tc ). The SVD is
not a unique technique for finding a basis set, and we note that the authors in[96] and [97] employ a greedy algorithm to efficiently generate a set of bases forgravitational waveforms which could in principle be applied to a set of ~z (tc ).
We have so far discussed interpolation in the mass parameters. It may also be

93

5. RAPIDLY EVALUATING THE COMPACT BINARY
LIKELIHOOD FUNCTION VIA INTERPOLATION

necessary to interpolate the quantities z[d , h0 ](tc ) in the tc direction, because thecoalescence time in a particular interferometer may lie in between discretely sampled time points. Second-order interpolation provides sufficient accuracy when
the waveform is sampled at 4 kHz.

5.1.2.4

Discussion and conclusion

We have demonstrated a method to sample the CBC likelihood function via
interpolation, with improvements of two to three orders of magnitude in efficiency.
Our method utilizes a SVD of the likelihood function on a three-dimensional
subspace of parameters: the chirp mass M, symmetric mass-ratio ⌘ and timeat coalescence tc . The SVD factors the likelihood function into a set of scalarcoefficients which describe a surface in M and ⌘, and a set of orthonormal basis
vectors which describe how the surface is translated along tc . The projection
coefficients can be interpolated on the M ⌘ plane and then trivially scaledby elements of the basis vectors to generate the likelihood at (M , ⌘ , tc ). This
provides an efficient means to interpolate in three dimensions.
We note that while we have chosen an interpolation technique based on the
SVD, it is by no means unique and other interpolation techniques have been applied to gravitational-wave data analysis [e.g., 97]. Notably, Mitra et al. [98] considered interpolating the matched-filtered output of gravitational-wave searches.
They interpolated the signal-to-noise ratio, which is e↵ectively a component of
the likelihood function, and so their method could, in principle, be extended tointerpolate likelihood functions. The key di↵erence with our approach is that we
use a decomposition of the likelihood as a function of time, while [98] treat it asa scalar quantity. This provides us with an efficient means of reducing the total
data needed for interpolation, exploiting correlations along the tc direction by
ranking the basis vectors in order of importance in reconstructing the likelihood
function. Hence, we can e↵ectively exorcise redundant information based on our
accuracy requirements. The number of bases needed to approximately reconstruct the likelihood to high accuracy using the SVD is generally small comparedto the number of raw likelihood vectors which we decompose.
Likelihood interpolation appears to be more robust than waveform interpola-94

tion [51], and so utilizing interpolated likelihood functions may also be a steppingstone to tackling the more difficult issue of rapidly estimating the parameters ofbinaries with spinning components.

95

5. RAPIDLY EVALUATING THE COMPACT BINARY
LIKELIHOOD FUNCTION VIA INTERPOLATION

96

Chapter 6
SUMMARY
We have considered two aspects of gravitational-wave astronomy; the ability
to detect intermediate-mass ratio coalescences (IMRACs) with Advanced LIGO(aLIGO), and methods to improve the efficiency of parameter estimation on gravitational waves from coalescing compact binaries using the singular value decomposition (SVD).
In Chapter II we studied the requirements of model waveforms for detecting IMRACs using aLIGO. We considered binaries in which the intermediatemass black hole has mass in the range 24 M
200 M with a stellar-masscompanion having masses in the range 1.4 M
18.5 M and with mass ratios
1/140  q  1/10. In this mass space we studied the contributions to the signalto-noise ratio (SNR) from the three stages of binary coalescence; inspiral, merger
and ringdown. We have shown that over the bulk of the mass space, merger andringdown contribute significantly to the coalescence signal, despite the suppression of the power in the merger and ringdown in the coalescence signal. Thisis because the systems we have considered can merge in the frequency range inwhich aLIGO is most sensitive. For the mass space we considered, the inspiralportion of the coalescence signal may fall at frequencies below the detectors sensitive band and hence the bulk contribution to SNR comes from the merger andringdown. However, we have also shown that despite the importance of merger
and ringdown, there is a small portion of the mass space in which the inspiralportion dominates the coalescence signal. We identified three regions in the massspace in which di↵erent searches could be considered appropriate based on thresh-97

6. SUMMARY

olds of acceptable losses in detection rates. The mass pace splits into a region inwhich inspiral-only searches could be feasible, incurring losses in detection rates
of up to ⇠ 27%; a region in which searches would be limited by lack of merger
and ringdown in template waveforms, incurring losses in detection rates up to50%; and a region in which merger and ringdown are critical to prevent losses indetection rates over 50%. The search regions are summarized in Table 2.1.
We also considered the accuracy of model waveforms which are currently
used for gravitational-wave searches for high-mass binaries in the IMRAC regime.
Model waveform families based on the e↵ective one-body approach, matched tonumerical relativity - known as EOBNR [26] - are only known to be accurate
at describing binaries with mass ratios q > 1/6 and so to detect IMRACs one
needs waveforms which are able to describe binaries which have smaller mass
ratios. We thus compared the EOBNR waveforms to a waveform family designed
for such small mass-ratio binaries which we refer to as the “Huerta-Gair” (HG)
[44; 63] waveform family. Because HG waveforms only describe the inspiral portion of the coalescence signal, we only compared the EOBNR inspiral to HG
waveforms. By computing the fitting factor between these waveform families, wefound that EOBNR-inspirals are
90% e↵ective at describing the HG inspiral.
This implies that EOBNR waveforms may be suitable for conducting a search forgravitational waves from IMRACs. However, this result reinforces the need todevelop approximation schemes which faithfully describe IMRACs. In particular,
due to the importance of merger and ringdown in the coalescence signal, numerical relativity simulations in this regime will be important in order to produce
template waveforms which describe the full coalescence signal.
In Chapters III, IV and V we considered methods to improve the efficiencyof parameter estimation of gravitational waves from coalescing binaries using
the singuar value decomposition (SVD). We applied the SVD in two ways: (i)to interpolate template waveforms and (ii) to directly interpolate the likelihood
function. The SVD has already been used to interpolate template waveforms forgravitational-wave searches [49] which reduces the computational cost of waveform generation in searches. Because the waveform accuracy requirements ofparameter estimation are typically higher than for gravitational-wave searches,in Chapter III we investigated the potential to use interpolated waveforms in pa-

98

rameter estimation as a means to reduce the computational cost associated withwaveform generation. We provided a proof-of-principle that interpolated template waveforms can be used for parameter estimation. While using interpolated
waveforms leads to biases in the posterior probability density functions (PDFs),
the bias in the mean posterior parameter estimates are smaller than the scale ofstatistical measurement uncertainty and thus we pass a common accuracy threshold for model waveform accuracy [89]. Furthermore, we provided a comparison of
the computational times of various waveform families. We found that in general,interpolating template waveforms is at least an order of magnitude faster thanusing non-interpolated waveforms and so using interpolated template waveforms
for parameter estimation is a promising technique to improve its efficiency.
Despite the potential of using interpolated waveforms, the study in Chapter
III was limited to one waveform family. We used the SVD to interpolate a class ofwaveforms known as TaylorF2 [26] which are already efficient to generate. Thiswas because of the simplicity of implementing TaylorF2 waveforms. In Chapter
IV we considered interpolating generic classes of template waveform. This extended the proof-of-principle in Chapter III to include inspiral-only time-domain
waveforms and time-domain waveforms which describe the full coalescence signal(inspiral, merger and ringdown). We found that we were able to use an interpolated inspiral-only time-domain waveform family known as TaylorT4 [26] as
templates for parameter estimation. We reported negligible biases in mean parameter estimates on the scale of statistical measurement uncertainty. We foundthat interpolated TaylorT4 waveforms are around an order of magnitude fasterto generate than non-interpolated TayorT4 waveforms, which is consistent with
our results in Chapter III. We also were able to interpolate time-domain waveforms which describe the full coalescence signal. We used a waveform family
which uses the “e↵ective one-body” approach, matched to numerical simulationsof the merger and ringdown, known as EOBNR [26]. However, we found that
we were not able to generate the interpolated templates at sufficient accuracy
for parameter estimation, though they could be considered accurate enough forgravitational-wave searches. There were a number of reasons for the loss in accuracy. Firstly, the relative time scales of the inspiral and ringdown are significantlydi↵erent. For the mass range we considered the inspiral phase is around 2 s in99

6. SUMMARY

duration while the ringdown is O(10 2 ) s in duration. In addition, we encountered a difficulty in our implementation of interpolated waveforms which makes
it difficult to extend to time-domain waveforms. When performing the SVD of
a set of template waveforms we found it useful to first align the waveforms intime at the frequency at which the noise PSD is a minimum. The time-frequency
evolution of time-domain waveforms can be inaccurate unless one uses a very
high sample rate, and hence the alignment of the waveforms we used can be o↵by a small amount. This finding was reinforced by the fact that at high sampling
rates, the accuracy of interpolated time-domain waveforms becomes higher.
In order to address the issue of efficient parameter estimation with genericwaveforms, in Chapter V we showed that the SVD enables one to directly interpolate the likelihood function. We analytically formulated the interpolated
likelihood function and showed that parameter estimation based on an interpolated likelihood function solves numerous computational issues. For example,
waveform generation is only required in order to generate the SVD of the likelihood function and evaluating the interpolated likelihood function requires no new
waveform computations. In addition likelihood interpolation should be totallytransparent to the template waveform family being used. We reported (conservative) improvements in the efficiency of evaluating the interpolated likelihood
function between two and three orders of magnitude, whilst keeping likelihoods
accurate to . 0.025%. We believe that this technique is highly relevant to parameter estimation in aLIGO because the in-band signals are liable to be around
two orders of magnitude longer than initial LIGO which will exacerbate the costof template waveform generation, and hence parameter estimation. Further work
on directly interpolating the likelihood function is well merited and is a promising
approach to rapidly estimating the parameters of compact binaries.

100Appendix A
FREQUENTLY USEDQUANTITIES
A.1

Definitions of discrete representations of continuous functions

Here we define the conventions that we adopt for discrete representations of continuous functions, which are frequently used.
For a real-valued function a(t), the convention for the Fourier Transform isã(f ) =Z +1

dt a(t)e 2⇡if t ,(A.1)

df ã(f )e+2⇡if t .

(A.2)

1

with the inverse transforma(t) =Z +1
1

For the case of discrete time/frequency series, consider a (discrete) real-valued
time series ~a sampled at discrete time intervals t = 1/(2fNyq ), where fNyq is
the Nyquist frequency. The discrete time series has a duration T in the timedomain. The total number of data samples of the discrete, real-valued timeseries is therefore N = T / t = 2T fNy . The data points at discrete time tj and101

A. FREQUENTLY USED QUANTITIES

frequency fk of ~a, and its (discrete) Fourier transform ã, are therefore:a(tj ) = a(j t) = ajã(fk ) = ã(k f ) =(A.3)
t ⇥ ãk ,(A.4)

where we have defined the Fourier series as:
ãk =

X

aj e 2⇡ijk/N ,

(A.5)

1 X
ãk e+2⇡ijk/N .
N k

(A.6)

j

aj =

We approximate the inner product (a|b) between two (real) functions a and b
in the discrete case as:(a|b) = 4<

Z 1
0

⇡ 4<

f

ã⇤ (f )b̃(f )df ,Sn (f )
N/2 ⇤Xã (fk )b̃(fk )k=Sn (fk )

(A.7)
,
(A.8)

where Sn is the detector noise power spectral density (PSD). The inner productcan further simplified by using “whitened” frequency series. We define a whitened
frequency series as the raw frequency series weighted by the square root of thenoise PSD:ã(fk )ãw (fk ) = p,(A.9)
Sn (fk )Using whitened waveforms, we define the vector inner product between two discretely sampled, whitened time-series ~aw and ~bw , as~aw . ~bw := 4 f

N/2Xã⇤w (fk )b̃w (fk ) ,(A.10)

k=0

which in general is complex. The standard definition of the inner product,
Eq. A.8, is just the real part of the vector inner product.

102

A.2

Post-Newtonian phasing formulae

We provide useful formulae for post-Newtonian gravitational waveforms. We use
v as our expansion parameter which is given by v = (⇡M f )1/3 in geometric units,
G = c = 1, where M is the total mass of the binary.

A.2.1

TaylorF2

For the frequency-domain representation of the gravitational-wave phasing, we
use the formula corresponding to the approximant TaylorF2. The phasing up to
3.5PN order is given by [26]

✓
◆
⇡
3
20 743 11+1++ ⌘ v 2 16⇡v 3c
4 128 ⌘ v 5
9 336
4
✓
◆
✓
◆
3058673 5429
617 2 438645 65
+ 10+⌘+
⌘ v +⇡⌘
1016064 1008
144
756
9
⇢
✓
◆
⇢
v
11583231236531 640 2 6848
6848
5
1 + 3 log
v +
⇡
vlso
4694215680
3
21
21
✓
◆
2
15737765635 2255 ⇡
76055 2 127825 3 6log (4 v)
+
⌘+
⌘
⌘ v
3048192
12
1728
1296
✓
◆
77096675 378515
74045 2 7
+⇡
+
⌘
⌘ v ,
254016
1512
756

(F2)
3.5 (f ) = 2⇡f tc(A.11)
where ⌘ = m1 m2 /(m1 + m2 )2 is the symmetric mass ratio, m1 and m2 are thecomponent masses, and = 0.577216 . . . is the Euler constant.

A.2.2

TaylorT4

The TaylorT4 orbital phase is given by solving the coupled di↵erential equationsof the orbital phase and velocity, (t) and v(t), respectively,
d
v3
= 0,
dt
M
dv
= K(t),
dt

103(A.12)(A.13)

A. FREQUENTLY USED QUANTITIES

where the acceleration K(t) is given by [26]


✓
◆
✓
◆
34103 13661
32 ⌘ 9743 11
59 2 42
3
K(t) =v 1+ ⌘ v + 4⇡v ++⌘+ ⌘ v
5 M
336
4
18144
2016
18
✓
◆
✓
4159 189
16447322263 16 2 1712+
⌘ ⇡v 5 ++ ⇡672
8139708800
3
105
✓
◆
◆
451 2 56198689
541 2 5605 3 856
2
+
⇡
⌘
⌘
⌘
log(16v ) v 6
48217728
896
2592
105
✓
◆
4415 358675
91495 2⌘⌘ ⇡v 7 .
(A.14)
4032
6048
1512

104

Appendix B
SINGULAR VALUEDECOMPOSITION OF M ⇥ N
MATRICES
This appendix provides supplementary material to Sec. 3.3. We give examples of
how the singular value decomposition (SVD) can be performed in practice for simple M ⇥ N matrices. In the examples, we will provide a geometric interpretationof the singular value decomposition.

B.1

Definitions of the singular value decomposition of M ⇥ N matrices

Consider a M ⇥ N matrix A, where M is the number of rows of A and N is
the number of columns. For simplicity we assume M
N . The singular value
decomposition of A is a factorization of the form [99]A = U⌃VT ,(B.1)

where U is an orthonormal M ⇥ N matrix, V is an orthonormal square N ⇥ N
matrix and ⌃ is a diagonal square N ⇥ N matrix. The matrices U, V and ⌃ aredefined as follows [100]:105

B. SINGULAR VALUE DECOMPOSITION OF M ⇥ N MATRICES• The diagonal elements of ⌃ are known as the “singular values” of A and are
denoted by i . The singular values correspond to the positive square rootsof the eigenvalues of the associated matrix K = AT A. For convenience, thesingular values are ordered such that 1...
2
N.
• The N columns of V and U are known as the right- and left-singular
vectors of A respectively. The right-singular vectors, ~vi are the normalized
eigenvectors of AT A and the left-singular vectors, ~ui , are the normalized
eigenvectors of AAT and satisfy
A ~vi =

ui ,
i~

(B.2)

AT ~ui =

vi .
i~

(B.3)

Below we provide a step-by-step example of computing the SVD of a simple
matrix A.

B.1.1

Step-by-step example

Compute the SVD of the matrix A given by
A=

"

2 2
1 2

#

.

(B.4)

Step 1: Compute the singular values of A
To compute the singular values we first need to compute AT A:
K = AT A =The eigenvalues,

"

5
2

2
8

#

.

(B.5)

i , are given by solving for the characteristic polynomial,det(KI2 ) = 0 ,(B.6)

where I2 is the 2 ⇥ 2 identity matrix. It is easy to show that the eigenvalues are
1 = 9,

106

2 = 4,(B.7)

and hence the singular values

i are

1 = 3,

2 = 2

(B.8)

Step 2: Compute the right-singular vectors of A
The normalized eigenvectors, ~vi , of K are the right-singular vectors of A. It
is simple to show that these are given by1~v1 = p
5

"

1
2

#

,

1
~v2 = p
5

"

2
1

#

.

(B.9)

Step 3: Compute the left-singular vectors of A
Having computed the singular values and right-singular vectors of A, we caneasily compute its left-singular vectors using Eq. (B.2). Since 1 , 2 6= 0 we canimmediately calculate ~u1 and ~u2 from Eq. (B.2):
11~u1 = A~v1 = p
5
1

"

2
1

#

,

1

1
~u2 = A~v2 = p
5
2

"

1
2

#

.

(B.10)

Step 4: Writing the singular value decomposition of A
We now have all the information to write the singular value decomposition of
A. The three matrices U, V and ⌃ are formed thus: The ith column of U and
V are given by the ui and vi respectively, and the diagonal elements of ⌃ areordered such that ⌃11 = 1 , ⌃22 = 2 etc... Hence,
U=

"

p2
5
p1
5

p1
5
p2
5

#

,

V=

"

p1
5
p2
5

p2
5
p1
5

#

,

⌃=

"

#

3 0
0 2

.

(B.11)The singular value decomposition of A is thus
A = U⌃V

T

=

"

p2
5
p1
5

p1
5
p2
5

#"

3 0
0 2

#"

p1
5
p2
5

p2
5
p1
5

#

,

(B.12)and it can be easily verified that this representation is equivalent to Eq. (B.4).

107

B. SINGULAR VALUE DECOMPOSITION OF M ⇥ N MATRICES

B.1.2

Geometric picture of the SVD

Figure B.1: Action of the matrix A, Eq. (B.4), on the unit circle. The transformation of the vectors ~v1 and ~v2 , Eq. (B.9), is shown explicitly.

The SVD has a deep connection with geometry (see e.g.[100; 101]) which
we will touch upon here. The canonical example of how the SVD relates togeometry is to consider the e↵ect of a matrix A on the unit circle, by looking at
the contributions of each matrix of the SVD of A for A given by Eq. (B.4). The
idea is to think of V and U as rotations and reflections, and ⌃ as a stretching.
Fig. B.1 illustrates the sequence of transformations induced by acting the matrix
A on the unit vectors ~v1 and ~v2 , and all other vectors on the unit circle when
A. We can understand how ~v1 and ~v2 are transformed by solving the following
system of linear equations
A ~vi = U⌃VT ~vi = ~xi ,(B.13)

for i = (1, 2). We already know the solution to this equation, which is given in
Eq. (B.2), however it is useful to provide a geometric interpretation of the solution.
When ~v1 and ~v2 are multiplied by VT , the two vectors undergo a rotation and108become unit vectors along the x- and y-axes respectively:
VT ~v1 =VT ~v2 =

"
"

p1
5
p2
5

p2
5
p1
5

p1
5
p2
5

p2
5
p1
5

#"

p1
5
p2
5

#"

p2
5
p1
5

#
#

"

= ~e1 ="

= ~e2 =

1
0
0
1

#
#

,

(B.14).

(B.15)The matrix ⌃ then stretches ~e1 and ~e2 :⌃~e1 =⌃~e2 =
Finally,e1 and
1~

"
"

3 0
0 2
3 0
0 2

#
#

~e1 = 3 ~e1 =

e1 ,
1~

(B.16)~e2 = 2 ~e2 =

e2 ,
2~

(B.17)

e2 are projected onto2~

U 1~e1 =U 1~e1 =

"
"

p2
5
p1
5

p1
5
p2
5

p2
5
p1
5

p1
5
p2
5

#
#

e1
1~
e2
2~

u1 and
1~
=
=

1

2

"
"

u2 :
2~
p2
5
p1
5

#

p1
5
p2
5

=
#

u1 ,
1~
=

u2 .
2~

(B.18)(B.19)

Hence, A projects the vectors ~v1 and ~v2 onto 1~u1 and 2~u2 , respectively. It iseasily verified that 1~u1 and 2~u2 define the major and minor axes of an ellipse.
One can see how the matrix A transforms all the other vectors which describe
points on the unit circle, by realising that any other vector on the unit circle will
simply correspond to a rotation of ~v1 or ~v2 by a fixed angle. Hence any pointon the circle is projected onto a point on the ellipse with major and minor axes
defined by 1~u1 and 2~u2 .

109

B. SINGULAR VALUE DECOMPOSITION OF M ⇥ N MATRICES110

REFERENCES[1] D. Sigg and the LIGO Scientific Collaboration, “Status of the LIGO detectors,” Class. Quant. Grav., vol. 25, p. 114041, 2008. 1, 8, 63, 64, 65, 76, 82,
90, 92[2] G. M. Harry and the LIGO Scientific Collaboration, “Advanced LIGO: thenext generation of gravitational wave detectors,” Classical and Quantum
Gravity, vol. 27, p. 084006, 2010. 1, 7
[3] K. Somiya, “Detector configuration of KAGRA the Japanese cryogenic
gravitational-wave detector,” Classical and Quantum Gravity, vol. 29,
p. 124007, 2012. 1[4] B. Iyer, T. Souradeep, C. Unnikrishnan, S. Dhurandhar, S. Raja, A. Kumar,and A. S. Sengupta, “LIGO-India Tech. rep.,” 2011. 1[5] J. A. et. al The LIGO Scientific Collaboration, “Predictions for the rates
of compact binary coalescences observable by ground-based gravitationalwave detectors,” Class. Quant. Grav., vol. 27, p. 173001, 2010. 1, 2, 10, 21,
49[6] J. A. et al. The LIGO Scientific Collaboration and The Virgo Collaboration,“Search for gravitational waves from compact binary coalescence in LIGO
and Virgo data from S5 and VSR1,” Phys. Rev. D, vol. 82, p. 102001, 2010.
1, 10, 57[7] J. A. et al. LIGO Scientific Collaboration, Virgo Collaboration, “Search forgravitational waves from low mass compact binary coalescence in LIGOs

111

REFERENCES

sixth science run and Virgos science runs 2 and 3,” Phys. Rev. D, vol. 85,
p. 082002, 2012. 1, 10, 11, 57, 71[8] J. A. et al. The LIGO Scientific Collaboration, “Search for Gravitational
Waves from Binary Black Hole Inspiral, Merger and Ringdown in LIGOVirgo Data from 2009-2010,” arXiv:1209.6533, 2012. 1, 10, 11, 28, 57,
71[9] J. A. et al. The LIGO Scientific Collaboration and The Virgo Collaboration, “An upper limit on the stochastic gravitational-wave background of
cosmological origin,” Nature, vol. 460, pp. 990–994, 2009. 1[10] J. A. et al. The LIGO Scientific Collaboration and The Virgo Collaboration,“Searching for a stochastic background of gravitational waves with the laserinterferometer gravitational-wave observatory,” The Astrophysical Journal,
vol. 659, p. 918, 2007. 1[11] J. Taylor and J. Weisberg, “A new test of general relativity - gravitationalradiation and the binary pulsar psr 1913+16,” Ap. J., vol. 253, pp. 908–920,
Feb. 1982. 1[12] J. Veitch and A. Vecchio, “Bayesian coherent analysis of in-spiral gravitational wave signals with a detector network,” Phys Rev. D, vol. 81,
p. 062003, 2010. 2, 12, 15, 53, 64, 79, 86, 88, 90[13] M. van der Sluys, I. Mandel, V. Raymond, V. Kalogera, C. Roever, and
N. Christensen, “Parameter estimation for signals from compact binary
inspirals injected into LIGO data,” Class. Quant. Grav., vol. 26, p. 204010,
2009. 2, 15, 53, 55, 86[14] D. A. Brown, J. Brink, H. Fang, J. R. Gair, C. Li, G. Lovelace, I. Mandel,
and K. S. Thorne, “Prospects for detection of gravitational waves fromintermediate-mass-ratio inspirals,” Phys. Rev. Lett., vol. 99, p. 201102,
2007. 2, 21, 22[15] C. L. Rodriguez, I. Mandel, and J. R. Gair, “Verifying the no-hair property of massive compact objects with intermediate-mass-ratio inspirals in

112

REFERENCES

advanced gravitational-wave detectors,” Phys. Rev. D, vol. 85, p. 062002,
2012. 2, 22[16] R. OShaughnessy, C. Kim, V. Kalogera, and K. Belczynski, “Constraining
population synthesis models via empirical binary compact object merger
and supernova rates,” The Astrophysical Journal, vol. 672, p. 479, 2008. 2[17] B. Kiziltan, A. Kottas, and S. Thorsett, “The neutron star mass distribution,” arXiv:1011.4291, 2010. 2[18] J. A. et al. The LIGO Scientific Collaboration and The Virgo Collaboration,“Search for gravitational waves from binary black hole inspiral, merger, andringdown,” Phys. Rev. D, vol. 83, p. 122005, Jun 2011. 2[19] R. Smith, I. Mandel, and A. Vecchio, “Studies of waveform requirements
for intermediate mass-ratio coalescence searches with advanced detectors,”
arXiv:1302.6049, 2012. 2, 17[20] S. Farrell, N. Webb, D. Barret, O. Godet, and J. Rodrigues, “Anintermediate-mass black hole of over 500 solar masses in the galaxy ESO
243-49.,” Nature, vol. 460, p. 7251, 2009. 2, 21[21] T. J. Maccarone, A. Kundu, S. E. Zepf, and K. L. Rhode, “A black hole ina globular cluster,” Nature, vol. 445, p. 183185, 2007. 2[22] A. Patruno, S. Portegies Zwart, J. Dewi, and C. Hopman, “The ultraluminous X-ray source in M82: an intermediate-mass black hole with a giant
companion,” MNRAS Lett., vol. 370, pp. L6–L9, 2006. 2[23] C. D. Bailyn, R. K. Jain, P. Coppi, and J. A. Orosz, “The mass distributionof stellar black holes,” The Astrophysical Journal, vol. 499, no. 1, p. 367,
1998. 3[24] W. Farr, N. Sravan, A. Cantrell, L. Kreidberg, C. D. Bailyn, I. Mandel,
and V. Kalogera, “The mass distribution of stellar-mass black holes,” The
Astrophysical Journal, vol. 741, p. 103, 2010. 3

113

REFERENCES[25] V. Raymond, M. V. van der Sluys, I. Mandel, V. Kalogera, C. Rover,and N. Christensen, “Degeneracies in Sky Localisation Determination froma Spinning Coalescing Binary through Gravitational Wave Observations:a Markov-Chain Monte-Carlo Analysis for two Detectors,” Class. Quant.
Grav., vol. 26, p. 114007, 2009. 3, 14, 50, 55, 86, 93[26] A. Buonanno, B. R. Iyer, E. Ochsner, Y. Pan, and B. S. Sathyaprakash,“Comparison of post-newtonian templates for compact binary inspiral signals in gravitational-wave detectors,” Phys. Rev. D, vol. 80, p. 084043, 2009.
4, 6, 10, 12, 13, 17, 23, 24, 28, 29, 32, 69, 74, 81, 89, 98, 99, 103, 104[27] P. C. Peters and J. Mathews, “Gravitational Radiation from Point Masses
in a Keplerian Orbit,” Phys. Rev., vol. 131, pp. 435–440, 1963. 4[28] B. Schutz, A First Course in General Relativity. Cambridge University
Press, 2009. 4[29] M. Maggiore, Gravitational Waves Volume 1: Theory and Experiments.
Oxford University Press, 2008. 4, 5, 6, 8, 10, 11, 26, 34, 54[30] J. A. et al. The LIGO Scientific Collaboration and The Virgo Collaboration,“Search for gravitational wave ringdowns from perturbed black holes in
LIGO S4 data,” Phys. Rev. D, vol. 80, p. 062001, 2009. 7, 11[31] L. S. Finn, “Gravitational wave interferometer noise calculator,” 2013. 8[32] S. Hild, “Beyond the second generation of laser-interferometric gravitational
wave observatories,” Classical and Quantum Gravity, vol. 29, p. 124006,
2012. 9[33] B. Allen, J. D. E. Creighton, E. E. Flanagan, and J. D. Romano, “Robust statistics for deterministic and stochastic gravitational waves in nonGaussian noise: Frequentist analyses,” Phys. Rev. D, vol. 65, p. 122002,
Jun 2002. 8[34] R. Biswas, P. Brady, J. Burguet-Castell, K. Cannon, J. Clayton, A. Dietz,
N. Fotopoulos, L. Goggin, D. Keppel, C. Pankow, L. Price, and R. Vaulin,

114

REFERENCES“Detecting transient gravitational waves in non-Gaussian noise with partially redundant analysis methods,” arXiv:1201.2964, 2012. 8[35] R. Biswas, P. Brady, J. Burguet-Castell, K. Cannon, J. Clayton, A. Dietz, N. Fotopoulos, L. Goggin, D. Keppel, C. Pankow, L. Price, andR. Vaulin, “Likelihood-ratio ranking of gravitational-wave candidates ina non-Gaussian background,” arXiv:1201.2959, 2012. 8[36] P. Jaranowski and A. Krolak, “Gravitational-Wave Data Analysis. Formalism and Sample Applications: The Gaussian Case,” 2012. 12[37] I. Mandel, “Parameter estimation on gravitational waves from multiple coalescing binaries,” Phys. Rev. D, vol. 81, p. 084029, 2010. 12[38] E. Poisson, “Gravitational radiation from a particle in circular orbit arounda black hole. vi. accuracy of the post-newtonian expansion,” Phys. Rev. D,
vol. 52, pp. 5719–5723, Nov 1995. 13, 47[39] N. Yunes and E. Berti, “Accuracy of the post-newtonian approximation:
Optimal asymptotic expansion for quasicircular, extreme-mass ratio inspirals,” Phys. Rev. D, vol. 77, p. 124006, Jun 2008. 13, 47[40] T. W. Baumgarte and S. L. Shapiro, Numerical Relativity: Solving Einstein’s Equations on the Computer. Cambridge University Press, 2010. 13[41] Y. Pan, A. Buonanno, M. Boyle, L. T. Buchman, L. E. Kidder, H. P.
Pfei↵er, and M. A. Scheel, “Inspiral-merger-ringdown multipolar waveforms
of nonspinning black-hole binaries using the e↵ective-one-body formalism,”
Phys. Rev. D, vol. 84, p. 124052, Dec 2011. 13, 23, 28, 30, 82, 89[42] J. R. Gair and K. Glampedakis, “Improved approximate inspirals of testbodies into kerr black holes,” Phys. Rev. D, vol. 73, p. 064037, 2006. 13,
23[43] S. Babak, H. Fang, J. R. Gair, K. Glampedakis, and S. A. Hughes,“”kludge” gravitational waveforms for a test-body orbiting a kerr blackhole,” Phys. Rev. D, vol. 75, p. 024005, 2007. 13, 23

115

REFERENCES[44] E. Huerta and J. R. Gair, “Intermediate-mass-ratio-inspirals in the einstein telescope: I. signal-to-noise ratio calculations,” Phys Rev. D, vol. 83,
p. 044020, 2011. 13, 16, 23, 28, 44, 98[45] H. Nakano, Y. Zlochower, C. O. Lousto, and M.Campanelli, “Intermediatemass-ratio black-hole binaries: numerical relativity meets perturbation theory,” Phys. Rev. D, vol. 84, p. 124006, 2011. 13, 23, 48[46] B. J. Owen, “Search templates for gravitational waves from inspiraling binaries: Choice of template spacing,” Phys. Rev. D, vol. 53, pp. 6749–6761,
1996. 14, 88[47] I. W. Harry, B. Allen, and B. S. Sathyaprakash, “Stochastic template placement algorithm for gravitational wave data analysis,” Phys. Rev. D, vol. 80,
p. 104014, Nov 2009. 14[48] E. T. Jaynes, Probability Theory: The Logic of Science. Cambridge University Press, 2003. 14[49] K. Cannon, A. Chapman, C. Hanna, D. Keppel, A. C. Searle, and A. J.
Weinstein, “Singular value decomposition applied to compact binary coalescence gravitational-wave signals,” Phys. Rev. D, vol. 82, p. 044025, 2010.
15, 17, 50, 56, 57, 58, 59, 65, 66, 89, 98[50] K. Cannon, C. Hanna, and D. Keppel, “Interpolating compact binary
waveforms using the singular value decomposition,” Phys. Rev. D, vol. 85,
p. 081504, 2012. 15, 17, 50, 56, 57, 61, 62, 65, 66, 88, 89, 90[51] R. Smith, K. Cannon, C. Hanna, D. Keppel, and I. Mandel, “Towards rapid
parameter estimation on gravitational waves from compact binaries using
interpolated waveforms,” arXiv:1211.1254, 2012. 18, 86, 87, 88, 90, 93, 95[52] S. W. Davis, R. Narayan, Y. Zhu, D. Barret, S. A. Farrell, O. Godet,
M. Servillat, and N. A. Webb, “The Cool Accretion Disk in ESO 24349 HLX-1: Further Evidence of an Intermediate-mass Black Hole,” ApJ,
vol. 734, p. 111, June 2011. 21

116

REFERENCES[53] A. K. H. Kong, R. Di Stefano, and F. Yuan, “Evidence of an IntermediateMass Black Hole: Chandra and XMM-Newton Observations of the Ultraluminous Supersoft X-Ray Source in M101 during Its 2004 Outburst,” ApJ.
Lett., vol. 617, pp. L49–L52, Dec. 2004. 21[54] H. Feng and R. Soria, “Ultraluminous X-ray sources in the Chandra and
XMM-Newton era,” New Astronomy Reviews, vol. 55, pp. 166–183, Nov.
2011. 21[55] G. Fabbiano, “Populations of X-Ray Sources in Galaxies,” Annual Reviewof Astronomy & Astrophysics, vol. 44, pp. 323–366, Sept. 2006. 21[56] M. C. Miller, “Intermediate-mass black holes as LISA sources,” Class.
Quant. Grav., vol. 26, p. 094031, 2009. 21[57] M. C. Miller and E. J. M. Colbert, “Intermediate-mass black holes,” Int.
J. Mod. Phys. D, vol. 13, pp. 1–64, 2004. 21[58] G. Harry and The LIGO Scientific Collaboration, “Advanced LIGO: thenext generation of gravitational wave detectors,” Class. Quant, Grav.,
vol. 27, p. 084006, 2010. 21, 49, 72[59] The Virgo Collaboration, “Advanced Virgo Baseline Design note
VIR027A0,” tech. rep., 2009. 21, 49[60] K. Somiya, “Detector configuration of KAGRA the Japanese cryogenic
gravitational-wave detector,” Classical and Quantum Gravity, vol. 29,
no. 12, p. 124007, 2012. 21[61] I. Mandel, D. A. Brown, J. R.Gair, and M. C. Miller, “Rates and characteristics of intermediate-mass-ratio inspirals detectable by advanced ligo,”
Astrophys. J., vol. 681, p. 1431, 2008. 21, 22[62] E. E. Flanagan and S. A. Hughes, “Measuring gravitational waves frombinary black hole coalescences. i. signal to noise for inspiral, merger, andringdown,” Phys. Rev. D, vol. 57, pp. 4535–4565, 1998. 22

117

REFERENCES[63] E. A. Huerta and J. R. Gair, “Influence of conservative corrections on parameter estimation for extreme-mass-ratio inspirals,” Phys. Rev. D, vol. 79,
p. 084021, 2009. 23, 28, 31, 98[64] T. Futamase and Y. Itoh, “The Post-Newtonian Approximation for Relativistic Compact Binaries,” Living Rev. Relativity, vol. 10, 2007. 24[65] L. Blanchet, “Gravitational Radiation from Post-Newtonian Sources andInspiralling Compact Binaries,” Living Rev. Relativity, vol. 9, 2006. 24[66] I. Mandel and J. R. Gair, “Can we detect intermediate mass ratio inspirals?,” Class. Quant. Grav., vol. 26, p. 094036, 2009. 24[67] C. Cutler and E. Flanagan, “Gravitational waves from merging compact
binaries: How accurately can one extract the binarys parameters from theinspiral waveform?,” Phys. Rev. D, vol. 49, p. 2658, 1994. 27, 32[68] T. Damour, “Coalescence of two spinning black holes: An e↵ective onebody approach,” Phys. Rev. D, vol. 64, p. 124013, Nov 2001. 28[69] T. Damour, P. Jaranowski, and G. Schäfer, “E↵ective one body approachto the dynamics of two spinning black holes with next-to-leading orderspin-orbit coupling,” Phys. Rev. D, vol. 78, p. 024009, Jul 2008. 28[70] T. Damour and A. Gopakumar, “Gravitational recoil during binary blackhole coalescence using the e↵ective one body approach,” Phys. Rev. D,
vol. 73, p. 124006, 2006. 30[71] A. Buonanno and T. Damour, “Transition from inspiral to plunge in binary
black hole coalescences,” Phys. Rev. D, vol. 62, p. 064015, 2000. 30[72] A. Buonanno, Y. Pan, J. G. Baker, J. Centrella, B. J. Kelly, S. T.
McWilliams, and J. R. van Meter, “Toward faithful templates for nonspinning binary black holes using the e↵ective-one-body approach,” Phys.
Rev. D, vol. 76, p. 104049, 2007. 30[73] L. Blanchet, G. Faye, B. R. Iyer, and S. Sinha, “The third post-newtonian
gravitational wave polarizations and associated spherical harmonic modes

118

REFERENCES

for inspiralling compact binaries in quasi-circular orbits,” Classical andQuantum Gravity, vol. 25, p. 165003, 2008. 32[74] The LIGO Scientific Collaboration, “Tech. Rep. LIGO-T0900288-v3 LIGO
Project,” 2009. 32, 33[75] B. S. Sathyaprakash, B. F. Schutz, and C. V. D. Broeck, “Cosmography
with the einstein telescope,” Classical and Quantum Gravity, vol. 27, no. 21,
p. 215006, 2010. 34[76] T. A. Apostolatos, “Search templates for gravitational waves from precessing, inspiraling binaries,” Phys. Rev. D, vol. 52, pp. 605–620, Jul 1995.
38[77] B. Allen, W. G. Anderson, P. R. Brady, D. A. Brown, and J. D. E.
Creighton, “Findchirp: An algorithm for detection of gravitational waves
from inspiraling compact binaries,” Phys. Rev. D, vol. 85, p. 122006, 2012.
39[78] T. G. F. Li, W. D. Pozzo, S. Vitale, C. V. D. Broeck, M. Agathos, J. Veitch,
K. Grover, T. Sidery, R. Sturani, and A. Vecchio, “Towards a generic testof the strong field dynamics of general relativity using compact binary coalescence,” Phys. Rev. D, vol. 85, p. 082003, 2012. 49[79] I. Mandel and R. O’Shaughnessy, “Compact binary coalescences in theband of ground-based gravitational-wave detectors,” Class. Quant. Grav.,
vol. 27, p. 114007, 2010. 49[80] T. Bulik and K. Belczyski, “Constraints on the binary evolution from chirp
mass measurements,” ApJ., vol. 589, p. L37, 2003. 49[81] T. A. Apostolatos, C. Cutler, G. J. Sussman, and K. S. Thorne, “Spininduced orbital precession and its modulation of the gravitational waveforms from merging binaries,” Phys. Rev. D, vol. 49, pp. 6274–6297, Jun
1994. 52[82] B. A. Berg, Markov Chain Monte Carlo Simulations and Their Statistical
Analysis. Singapore, World Scientific, 2004. 54, 55

119

REFERENCES[83] N. Metropolis, A. Rosenbluth, M. Rosenbluth, A. Teller, and E. Teller,“Equation of state calculations by fast computing machines,” Journal ofChemical Physics, vol. 21, no. 6, p. 10871092, 1953. 56[84] W. Hastings, “Monte carlo sampling methods using markov chains and theirapplications,” Biometrika, vol. 51, no. 1, p. 97109, 1970. 56[85] The LIGO Scientific Collaboration, “LSC Algorithm Library (LAL),” 2013.
56, 89[86] A. Gelman and D. B. Rubin, “Inference from iterative simulation using multiple sequences,” Statistical Science 7, 4, 457-472 (1992), vol. 7, pp. 457–
472, 1992. 64, 79[87] P. Ajith, “Addressing the spin question in gravitational-wave searches:
Waveform templates for inspiralling compact binaries with nonprecessing
spins,” Phys. Rev. D, vol. 84, p. 084037, Oct 2011. 65, 72[88] L. Lindblom, B. J. Owen, and D. A. Brown, “Model waveform accuracy
standards for gravitational wave data analysis,” Phys. Rev. D, vol. 78,
p. 124020, 2008. 66[89] F. Ohme, “Analytical meets numerical relativity - status of complete gravitational waveform models for binary black holes,” Class. Quant. Grav.,
vol. 29, p. 124002, 2012. 66, 67, 69, 74, 78, 81, 99[90] W. J. Conover, Practical Nonparametric Statistics. John Wiley and Sons,
1971. 67[91] A. Taracchini, Y. Pan, A. Buonanno, E. Barausse, M. Boyle, T. Chu,
G. Lovelace, H. P. Pfei↵er, and M. A. Scheel, “Prototype e↵ective-one-bodymodel for nonprecessing spinning inspiral-merger-ringdown waveforms,”
Phys. Rev. D, vol. 86, p. 024011, Jul 2012. 72[92] B. A. et al., “Testing gravitational-wave searches with numerical relativity
waveforms: results from the first Numerical INJection Analysis (NINJA)
project,” Classical and Quantum Gravity, vol. 26, p. 165008, 2009. 74

120

REFERENCES[93] “NINJA BBH Project 2 Home Page,” 2013. 74[94] J. Skilling, “Nested sampling,” AIP Conf. Ser. Vol. 735, p. 395, 2004. 86[95] The LIGO Scientific Collaboration and The Virgo Collaboration, “Parameter estimation for compact binary coalescence signals with the first generation gravitational-wave detector network,” arXiv:1304.1775, 2013. 90[96] S. Field, C. Galley, F. Herrmann, J. Hesthaven, E. Ochsner, and M. Tiglio,“Reduced basis catalogs for gravitational wave templates,” Physical Review
Letters, vol. 106, no. 22, p. 221102, 2011. 93[97] P. Canizares, S. E. Field, J. R. Gair, and M. Tiglio, “Gravitational wave parameter estimation with compressed likelihood evaluations,”
arXiv:1304.0462, 2013. 93, 94[98] S. Mitra, S. V. Dhurandhar, and L. S. Finn, “Improving the efficiency of thedetection of gravitational wave signals from inspiraling compact binaries:
Chebyshev interpolation,” Phys. Rev. D, vol. 72, p. 102001, Nov 2005. 94[99] M. Galassi, J. Davies, J. Theiler, B. Gough, G. Jung-man, P. Alken,
M. Booth, and F. Rossi, GNU Scientific Library Reference Manual, 3rded., for version 1.12. Network Theory Ltd, United Kingdom, 2009. 105[100] G. Golub and C. Van Loan, Matrix Computations. The Johns Hopkins
University Press (3rd edition), 1996. 105, 108[101] J. Demmel, Applied Numerical Linear Algebra. Society for Industrial andApplied Mathematics, 1997. 108

121