MULTI-BASELINE SEARCHES FOR STOCHASTIC SOURCES AND
BLACK HOLE RINGDOWN SIGNALS IN LIGO-VIRGO DATABy
DIPONGKAR TALUKDER

A thesis submitted in partial fulfillment of
the requirements for the degree of
DOCTOR OF PHILOSOPHYWASHINGTON STATE UNIVERSITY
Department of Physics and Astronomy
AUGUST 2012

To the Faculty of Washington State University:
The members of the Committee appointed to examine the dissertationof DIPONGKAR TALUKDER find it satisfactory and recommend that it
be accepted.
—————————————Sukanta Bose, Ph.D., Chair—————————————Matthew Duez, Ph.D.
—————————————Gregory Mendell, Ph.D.
—————————————Guy Worthey, Ph.D.

ii

ACKNOWLEDGEMENTFirstly, I would like to express my gratitude and warm-hearted thanks to
Sukanta Bose for giving me the opportunity to contribute to gravitational-wave
science. I am indebted to him for his constant support and guidance throughout
my Ph.D. He certainly has been, without any doubt, a wonderful critic of mywork.
I enjoyed and benefited from the interactions I had over the years with themembers of the LIGO Scientific Collaboration and the Virgo Collaboration. I
would like to sincerely thank my research collaborators. First, I have had thegreat pleasure to work with Paul Baker, Sarah Caudill and Sanjit Mitra. My special thanks to Warren Anderson, Stefan Ballmer, Matthew Benacquista, Nelson
Christensen, Neil Cornish, Sanjeev Dhurandhar, Gabriela González, Kari Hodge,
Shivaraj Kandasamy, Albert Lazzarini, Vuk Mandic, Sean Morriss, Tania Regimbau, Emma Robinson, Joseph Romano, Eric Thrane, Alan Weinstein and many
more for helpful discussions.
I am deeply obliged to Bruce Allen for his support and hospitality during
my stay at Albert Einstein Institute, Hannover, Germany. Special thanks to myfriends and colleagues Benjamin Knispel, Chris Messenger, Ajith Parameswaran,
Holger Pletsch and Reinhard Prix for their help, friendship and entertainment.
I have benefited from enlightening discussions with Gregory Mendell andFred Raab in our weekly group meetings. I also want to express sincere thanks
to my teachers in Washington State University for their inspiration. I take thisopportunity to thank my friends in Washington State University. These peopleiii

have made my time at WSU productive. I especially enjoyed the company of
JiaJia Chang, Thilina Dayanga, Brett Deaton, Shaon Ghosh, Fatemeh H. Nouri,
and Szymon Steplewski. I thank Carl Brannen for reading the thesis carefully
and providing helpful comments.
Additionally, I would like to thank Matthew Duez, Gregory Mendell and
Guy Worthey for all of their time and being willing to serve on my thesis committee. I would like to acknowledge the National Science Foundation for funding my research in gravitational-wave data-analysis in support of LIGO. This
work has been supported by the National Science Foundation Grant Nos. PHY0758172 and PHY-0855679. I would also like to thank College of Sciences (COS)at Washington State University for awarding me a COS Research Assistantshipfor Spring 2012.
Finally, I cannot end without thanking my family for their sacrifice and encouragement. I owe my deepest thanks to Shampa Biswas for her unending loveand enormous patience. In February 2009, I lost my father ......, thanks to Steven
Tomsovic for letting me visit home in the middle of the semester, and Gordon
Johnson, Ranga Dias and Ashish Mishra for taking care of teaching my labs.

iv

MULTI-BASELINE SEARCHES FOR STOCHASTIC SOURCES AND BLACK
HOLE RINGDOWN SIGNALS IN LIGO-VIRGO DATA
Abstract

by Dipongkar Talukder, Ph.D.
Washington State University
August 2012
Chair: Sukanta Bose

We present a framework for the detection of stochastic gravitational-wave(GW) backgrounds, from cosmological and astrophysical sources, using radiometry with a network of gravitational-wave interferometers. The search statistic
itself is derived from the cross-correlations of the data across all possible baselines in a detector network, and reveals how much more sensitive a network isthan any of its component baselines. We model the neutron star distribution in
the Virgo cluster and apply the above framework to search for their stochastic
GW signature in LIGO-VIRGO data.
We also present a template-based multi-detector coherent search for perturbed black hole ringdown signals. Like the past “coincidence” ringdown
searches in LIGO data, our method incorporates knowledge of the ringdown
waveform in constructing the search templates. Additionally, it checks for consistency of signal amplitudes and phases in the different detectors with theirdifferent orientations and signal arrival times. We demonstrate the advantagesof implementing a coherent search in the ringdown search pipeline.
v

DEDICATION

There are many people in the world in poverty due to circumstances theycannot control. Even though many people fear the future, many have overcome
poverty, neglect, and difficult situations to become successful. To them I dedicate
my thesis.

vi

PREFACE
The thesis focuses on the efforts to detect gravitational waves from neutronstars and black holes. The work presented here is based on three peer-reviewed
journal publications [12, 13, 14], and my direct involvement in two LIGO-VIRGO
searches. My participation in all of these projects has been extensive. In thecontext of signal processing, I studied coherent contribution of gravitationalwave detectors to extracted signal power for different types of sources. In places
where I have presented results derived from LIGO data I have done so withpermission from the LIGO Scientific Collaboration (LSC). A part of the methods
and results presented here is under review and is potentially subject to change.
The opinions expressed here are my own and not necessarily those of the LSC.

vii

Contents
1

2

Introduction: Gravitational waves and their sources and detection1

1.1

Gravitational waves (GWs) in general relativity . . . . . . . . . . .

3

1.2

Sources of gravitational waves . . . . . . . . . . . . . . . . . . . . .

6

1.3

Principle of a GW detector . . . . . . . . . . . . . . . . . . . . . . . .

8Detecting an anisotropic stochastic background with multi-baseline GW
radiometry

15

2.1

The stochastic GW backgrounds (SGWBs) . . . . . . . . . . . . . . .

16

2.1.1

The detector strain due to a SGWB . . . . . . . . . . . . . .

19

2.1.2

Statistical properties of the signal and detector noise . . . .

22

2.2

Principle of a GW radiometer . . . . . . . . . . . . . . . . . . . . . .

23

2.3

Optimal search statistic . . . . . . . . . . . . . . . . . . . . . . . . . .

27

2.3.1

Cross correlation statistic . . . . . . . . . . . . . . . . . . . .

27

2.3.2

Detection statistic . . . . . . . . . . . . . . . . . . . . . . . . .

28

2.3.3

Performance of optimal detection statistic . . . . . . . . . .

35

Performance of multi-baseline radiometers . . . . . . . . . . . . . .

41

2.4.1

42

2.4

Sensitivity . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

viii

2.5
3

Sky coverage . . . . . . . . . . . . . . . . . . . . . . . . . . .

45

2.4.3

Parameter accuracy. . . . . . . . . . . . . . . . . . . . . . .

49

2.4.4

Map making . . . . . . . . . . . . . . . . . . . . . . . . . . . .

53

2.4.4.1

Dirty maps and clean maps . . . . . . . . . . . . .

55

2.4.4.2

Maximized-likelihood-ratio statistic . . . . . . . . .

57

Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

60

Searching for a SGWB from the population of neutron stars in the Virgocluster69

3.1

Motivation for the directed search . . . . . . . . . . . . . . . . . . .

70

3.2

Modeling the neutron stars in the Virgo cluster . . . . . . . . . . .

75

3.2.1

The gravitar model . . . . . . . . . . . . . . . . . . . . . . . .

76

3.2.2

A population synthesis model of millisecond pulsars . . . .

80

3.2.3

A simple power-law . . . . . . . . . . . . . . . . . . . . . . .

82

3.3

The directed SGWB search and Olbers’ paradox . . . . . . . . . . .

82

3.4

Setting up the analysis . . . . . . . . . . . . . . . . . . . . . . . . . .

84

3.4.1

Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . .

84

3.4.2

Angular extension of the Virgo cluster . . . . . . . . . . . .

86

3.4.3

The analysis pipeline . . . . . . . . . . . . . . . . . . . . . . .

86

3.4.4

Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . .

88Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

90

3.5
4

2.4.2

Searching for perturbed black hole ringdown signals

91

4.1

Motivation for the search . . . . . . . . . . . . . . . . . . . . . . . .

92

4.2

Quasi-normal modes of black holes . . . . . . . . . . . . . . . . . .

94

ix

4.3

4.4

4.2.1

Schwarzschild and Kerr black holes . . . . . . . . . . . . . .

95

4.2.2

Black hole perturbation . . . . . . . . . . . . . . . . . . . . .

95

4.2.3

The ringdown waveform . . . . . . . . . . . . . . . . . . . .

97

Ringdown searches with matched filtering . . . . . . . . . . . . . . 100
4.3.1

Matched filter . . . . . . . . . . . . . . . . . . . . . . . . . . . 100

4.3.2

Template banks . . . . . . . . . . . . . . . . . . . . . . . . . . 101

A coincident multi-detector search . . . . . . . . . . . . . . . . . . . 102
4.4.1

4.4.2
4.5

The coincidence search pipeline . . . . . . . . . . . . . . . . 104
4.4.1.1

Data set . . . . . . . . . . . . . . . . . . . . . . . . . 104

4.4.1.2

Template placement . . . . . . . . . . . . . . . . . . 106

4.4.1.3

Trigger generation . . . . . . . . . . . . . . . . . . . 109

4.4.1.4

Coincidence analysis . . . . . . . . . . . . . . . . . 113

Tuning the search . . . . . . . . . . . . . . . . . . . . . . . . . 115

The coherent multi-detector search . . . . . . . . . . . . . . . . . . . 125
4.5.1

Elliptically-polarized ringdown template . . . . . . . . . . . 125
4.5.1.1

Fourier transformation of a generic decaying sinusoid . . . . . . . . . . . . . . . . . . . . . . . . . . 129

4.5.1.2

Fourier transformation of an elliptically-polarized
ringdown template . . . . . . . . . . . . . . . . . . 131

4.5.1.3

The template normalization . . . . . . . . . . . . . 132

4.5.1.4

Numerical analysis . . . . . . . . . . . . . . . . . . 133

4.5.1.5

SNR statistic . . . . . . . . . . . . . . . . . . . . . . 136

4.5.2

Coherent statistics . . . . . . . . . . . . . . . . . . . . . . . . 137

4.5.3

Performance study . . . . . . . . . . . . . . . . . . . . . . . . 140

x

4.6
5

Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141

Aspects of a blind hierarchical coherent search for coalescing compact
binaries

145

5.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147

5.2

Multi-detector statistics . . . . . . . . . . . . . . . . . . . . . . . . . 149
5.2.1

Signal and noise . . . . . . . . . . . . . . . . . . . . . . . . . 150

5.2.2

The network detection statistic . . . . . . . . . . . . . . . . . 152

5.2.3

Alternative statistics . . . . . . . . . . . . . . . . . . . . . . . 159

5.3

Coherent hierarchical inspiral analysis pipeline . . . . . . . . . . . 162

5.4

Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 165

5.5Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 168

A Technical details

173

A.1 Antenna response . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 173
A.2 Parameter accuracy . . . . . . . . . . . . . . . . . . . . . . . . . . . . 176
B Advanced detectors studies

178

B.1 Performance of multi-baseline radiometers . . . . . . . . . . . . . . 178
C A directed SGWB search and Olbers’ paradox

183

C.1 Virgo cluster versus the isotropic background . . . . . . . . . . . . 183
D Notes on SNR threshold186

D.1 SNR thresholds for the S4 ringdown search . . . . . . . . . . . . . . 186
E List of abbreviations

188

xi

List of Figures
1.1

Effect of the two polarizations of GWs on test particles . . . . . . .

6

1.2

Schematic of an earth-based interferometric detector . . . . . . . .

10

1.3

Different noise sources that contribute to the sensitivity of the first
generation LIGO detectors . . . . . . . . . . . . . . . . . . . . . . . .

1.4

Strain amplitude spectral densities of LIGO detectors during their
S5 science runs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1.5

11

13

The noise in LIGO, Narrow-band (NB) LIGO-II, and Wide-band(WB) LIGO-II, and the estimated signal strengths plotted as afunction of gravitational wave frequency . . . . . . . . . . . . . . .

2.1

Gravitational wave spectrum as a function of frequency for SGWB
sources from inflationary models . . . . . . . . . . . . . . . . . . . .

2.2

17

The 95% Bayesian upper limit of the point estimate of a constant
SGWB spectrum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2.3

14

18

Gravitational wave spectrum as a function of observed frequency
for potential AGWB sources . . . . . . . . . . . . . . . . . . . . . . .

20

2.4

Geometry of an elementary GW radiometer . . . . . . . . . . . . .

24

2.5

Illustration of a GW radiometer beam pattern . . . . . . . . . . . .

30

xii

2.6

Designed noise amplitude spectrum of the initial and advanced
LIGO and VIRGO, and the Einstein Telescope . . . . . . . . . . . .

2.7

The injected, dirty, and clean sky maps of a weak polar-cap GWB
source . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2.8

38

Histograms of the dirty maps for the noise-only data set and weakpolar-cap signal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2.9

36

39

The distribution of the MLR statistic obtained for dirty maps generated from noise-only and weak polar-cap injection . . . . . . . .

40

2.10 The sensitivities and their area-weighted counterparts of three different baselines and their network are plotted as functions of thedeclination of a single-pixel SGWB source . . . . . . . . . . . . . . .

44

2.11 The sensitivities of three single baselines and their multibaseline
network plotted as functions of the central frequency of the sourceband . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

45

2.12 The standard deviation of dirty maps measured by the three LIGOVirgo baselines and the full network for a constant source power
spectrum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

47

2.13 The beam functions for the three LIGO-Virgo baselines for a constant source spectrum . . . . . . . . . . . . . . . . . . . . . . . . . .

49

2.14 Singular values of the Fisher matrices for individual baselines andthe whole network . . . . . . . . . . . . . . . . . . . . . . . . . . . .

50

2.15 The 1σ error and area-weighted 1σ error in the solid angle forlocating a source in the sky with three single baselines and their
network . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

xiii

54

2.16 Toy model of an extended source and the dirty maps made fromsimulated data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

56

2.17 Clean maps obtained by the deconvolution of the dirty maps ofFig. 2.16 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

57

2.18 Difference between the clean maps of Fig. 2.17 and the injected
map of Fig. 2.16(a) . . . . . . . . . . . . . . . . . . . . . . . . . . . .

58

2.19 The toy model of a localized source and the dirty maps made fromsimulated data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

59

2.20 Clean maps obtained by the deconvolution of the dirty maps ofFig. 2.19 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

60

2.21 Difference between the clean maps of Fig. 2.20 and the injected
map of Fig. 2.19(a) . . . . . . . . . . . . . . . . . . . . . . . . . . . .

61

2.22 Dirty maps made for simulated noise, without any injected signals

63

2.23 Clean maps obtained by the deconvolution of the dirty maps ofFig. 2.22 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

64

2.24 The toy model of a very weak extended source and the dirty maps
made from simulated data . . . . . . . . . . . . . . . . . . . . . . . .

66

2.25 Clean maps obtained by the deconvolution of the dirty maps ofFig. 2.24 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

67

3.1

The distribution of observed radio pulsars . . . . . . . . . . . . . .

72

3.2

The spectral density of gravitational radiation from NSs in the
Virgo cluster plotted as a function of GW frequency . . . . . . . . .

xiv

73

3.3

The noise amplitude spectral densities of various advanced detectors are compared with the effective source amplitude spectral
density . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

74

3.4

The distribution of observed millisecond pulsars . . . . . . . . . . .

80

3.5

Several model-dependent spectral profiles of SGWB from NSs inthe Virgo cluster . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

83

3.6

The map of a modeled signal-strength vector of the Virgo cluster .

86

3.7

The spherical harmonic coefficients of the modeled signal-strength
vector . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

87

3.8

A schematic diagram of the stochastic search pipeline . . . . . . . .

88

3.9

The MLR statistic of a simulated SGWB in the Virgo cluster isplotted for different harmonic cutoff . . . . . . . . . . . . . . . . . .

89

4.1

Ringdown waveform for a black hole of mass 100 M and spin 0.7

98

4.2

Central frequency of the ringdown waveform as a function ofmass for three spin values and quality factor as a function of spin
are plotted . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4.3

99

A schematic diagram of the coincidence and coherent stages of
the ringdown search pipeline. . . . . . . . . . . . . . . . . . . . . . . 105

4.4

The ringdown horizon distance as a function of mass for a spinning black hole . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107

4.5

The template bank for the LIGO-VIRGO S5/S6 ringdown search . 108

4.6

Contours of constant mismatch around templates of a constantfrequency and for a portion of the template bank . . . . . . . . . . 109

xv

4.7

Ringdown waveforms in the extremes of the parameter space,characterized by central frequency and quality factor . . . . . . . . 110

4.8

Decisive distance of missed and found injections is plotted against
chirp mass and frequency . . . . . . . . . . . . . . . . . . . . . . . . 119

4.9

Scatter plots of simulated and estimated effective distance of found
PHENOM injections . . . . . . . . . . . . . . . . . . . . . . . . . . . 120

4.10 Scatter plots of simulated mass 2 versus mass 1 of made and found
PHENOM injections and histogram of number of made and found
PHENOM injections per mass-ratio bin . . . . . . . . . . . . . . . . 121
4.11 Scatter plots of simulated and estimated effective distance of found
RINGDOWN injections . . . . . . . . . . . . . . . . . . . . . . . . . 122
4.12 Scatter plots of simulated and estimated frequency of found RINGDOWN injections, also scatter plots of simulated and estimatedmass . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123
4.13 Scatter plots of simulated quality factor and frequency, and spin
and mass of made and found RINGDOWN injections . . . . . . . . 124
4.14 Contour plots of H2 SNR versus H1 SNR, L1 SNR versus H1 SNR,
and L1 SNR versus H2 SNR for triple coincident time-slide events 126
4.15 Contour plots of H2 SNR versus H1 SNR, L1 SNR versus H1 SNR,
and L1 SNR versus H2 SNR for triple coincident injection events . 127
4.16 The H1-H2, H1-L1, and H2-L1 SNR distributions for time-slides
and injections . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128
4.17 Scatter plots of the old and new triple-coincidence ranking statistics129

xvi

4.18 The efficiency of finding injections and the ROC curve are compared between the old and new triple-coincidence ranking statistics 130
4.19 The integrands of the sine and cosine template variances are plotted as functions of frequency . . . . . . . . . . . . . . . . . . . . . . 135
4.20 Contours of the fractional difference of sine and cosine phases ofan elliptically-polarized ringdown template . . . . . . . . . . . . . . 136
4.21 The scatter plots for the null-stream and coherent SNR, and thebitten-L and effective-coherent SNR . . . . . . . . . . . . . . . . . . 142
4.22 Efficiency of finding signals using the coincidence search is compared with that of using the coherent search . . . . . . . . . . . . . 143
4.23 The ROC curve of the coincidence search is compared with the
ROC curve of the coherent search . . . . . . . . . . . . . . . . . . . 144
5.1

A schematic diagram of the coherent stage in the compact binary
coalescence search pipeline . . . . . . . . . . . . . . . . . . . . . . . 166

5.2

Scatter plots of the combined and coherent SNRs of the triggers
from simulated data . . . . . . . . . . . . . . . . . . . . . . . . . . . 171

5.3

The receiver operating characteristic curves of three CBC searches
are compared . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172

A.1 Definition of source location and orientation angles relative to theinertial frame of a GW detector . . . . . . . . . . . . . . . . . . . . . 174
A.2 The antenna response of interferometric GW detectors to GWs
with plus and cross polarizations . . . . . . . . . . . . . . . . . . . . 175

xvii

B.1 The sensitivities and area-weighted sensitivities of different LIGOVIRGO-INDIGO-KAGRA networks in their second and third generation configurations are plotted as functions of the declinationof a single-pixel SGWB source . . . . . . . . . . . . . . . . . . . . . 179
B.2 The beam functions of different second-generation LIGO-VIRGOINDIGO-KAGRA baselines for a constant source PSD . . . . . . . . 180
B.3 Singular values of the Fisher matrices for different second-generation
LIGO-VIRGO-INDIGO-KAGRA networks . . . . . . . . . . . . . . 181
B.4 The 1σ error and area-weighted 1σ error in the solid angle of locating a source in the sky with different LIGO-VIRGO-INDIGOKAGRA networks in their second and third generation configurations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182

xviii

List of Tables
1.1

International network of gravitational-wave interferometers . . . .

2.1

MLR statistic of dirty maps versus clean maps for the simulated
maps in Fig.s 2.16 and 2.17 . . . . . . . . . . . . . . . . . . . . . . .

2.2

65

MLR statistic of dirty maps versus clean maps for the simulated
maps in Fig.s 2.24 and 2.25 . . . . . . . . . . . . . . . . . . . . . . .

4.1

62

MLR statistic of dirty maps versus clean maps for simulated noisein Fig.s 2.22 and 2.23 . . . . . . . . . . . . . . . . . . . . . . . . . . .

2.4

62

MLR statistic of dirty maps versus clean maps for the simulated
maps in Fig.s 2.19 and 2.20 . . . . . . . . . . . . . . . . . . . . . . .

2.3

13

65

Extreme values for frequency and quality factor of the ringdown
template bank . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111

xixChapter 1
Introduction: Gravitational waves
and their sources and detection
In 1916, Albert Einstein published his famous paper entitled “Die Grundlage
der allgemeinen Relativitätstheorie (The Foundation of the General Theory ofRelativity)” [1]. General Relativity (GR) describes how matter and energy affects
the curvature of space-time. The existence of gravitational waves is predicted in
GR. A gravitational wave is a fluctuation in the curvature of space-time that
propagates with the speed of light. It transports energy and momentum awayfrom the source. When two massive objects, like neutron stars, orbit each other,
space-time is stirred by their motion, and gravitational radiation ripples outward
into the universe. The key possible sources of gravitational waves are colliding
compact binaries, rotating non-axisymmetric neutron stars, stellar core collapse,
and stochastic background.
In General Relativity, the lowest mass multipole that produces gravitational

1

waves is the quadrupole. The conservation of energy and momentum forbids
gravitational wave emission associated with the monopole and dipole moments.
Gravitational waves (GWs) decrease in amplitude as the inverse of the distancefrom the source. Owing to their extraordinarily weak interaction with matter,
GWs travel cosmological distances with negligible absorption or distortion. Thishas made direct detection of gravitational waves extremely challenging. Nevertheless, the weakness of interaction enables us to obtain precise information on
their sources.
Indirect evidence for the existence of gravitational waves was found by Hulse
and Taylor [2]. They discovered a binary system, PSR 1913+16, that is composedof a pulsar and a companion neutron star (NS). From the long term precise observations of the pulsar frequency, the decay rate of the orbital period was foundto be consistent with the predicted energy loss due to the emission of gravitational waves. For this discovery, they were awarded the Nobel Prize in physics
in 1993. However, no one has yet succeeded in directly detecting gravitational
waves. Several kilometer-scale gravitational-wave detectors such as LIGO andVIRGO have been constructed and a few more are in the process of being built
world-wide. Detection of gravitational waves is not only important to test GR
but also because it promises to open a new window for GW astronomy. It willtell us about the dynamics of large-scale events in the universe like the death ofstars, and the birth of black holes (BHs).

2

1.1

Gravitational waves in general relativity

In this work, we adopt the following notation: The space-time coordinates aredenoted by x µ := ( x0 , xi ) = ( x0 , x1 , x2 , x3 ) = (ct, x, y, z). While the Greek letters(e.g. µ, ν) are space-time indices (0, 1, 2, 3), the Latin indices (e.g. i, j) denotethe three spatial coordinate labels, 1, 2, 3. Partial derivatives with respect to thecoordinates are abbreviated using the symbols ∂µ ≡ ∂/∂x µ and ∂µ ≡ ∂/∂xµ .
We also define a comma-notation to indicate derivatives, i.e., ∂ β f ≡ f ,β and∂ β f ≡ f ,β . Repeated indices are summed over unless otherwise indicated.
The General Theory of Relativity asserts that the curvature of space-time isgravity. Curved space-time is best described in Riemannian manifolds [3] by ametric gµν , which is a function of the space-time coordinates x µ . The infinitesimal proper distance dτ between two space-time events separated by coordinate
distance dx µ is given bydτ 2 = gµν dx µ dx ν ,(1.1)

where the Einstein summation convention is understood on the right hand side.
In the limit of flat space-time, the metric gµν becomes the Minkowski metricηµν = diag(−1, 1, 1, 1).
The curvature and dynamics of space-time governed by the distribution andkinematics of energy are expressed geometrically by the Einstein Field Equations(EFE):Gµν ( g) =

8πG
Tµν ,
c4

(1.2)

where Gµν is the Einstein tensor, G is Newton’s constant and Tµν is known as
the stress-energy tensor. Due to its inherent nonlinearity, a complete description3

of the solutions of EFE is not yet possible. One class of solutions to the EFE is
obtained in the linearized theory close to flat space-time, i.e., where fields are
everywhere weak.
In weak gravity, the metric is approximately Minkowskian and can be decomposed into the flat space-time metric plus a small perturbation,gµν = ηµν + hµν ,|hµν |  1.

(1.3)

The most general coordinate transformation of the form
xµ → xµ + ξ µ (x) ,µ|ξ ,ν |  1,

(1.4)

preserves the form of the metric given by (1.3) if the perturbations are transformed using the formulahµν → hµν − ξ µ,ν − ξ ν,µ .

(1.5)

To see the physical effect of a gravitational wave it is useful to choose a gauge.
The Lorentz gauge freedom allows us to choose a coordinate system where thetrace-reversed metric perturbations,1h̄µν := hµν − ηµν hαα ,
2

(1.6)

h̄µν,ν = 0.

(1.7)

are divergence free,This reduces the Einstein field equation to a simple wave equation:h̄µν = −

16πG
Tµν ,
c4

4

(1.8)

where  denotes the d’Alembertian operator:
1 ∂2
 : = ∂ µ ∂ µ = − 2 2 + ∇2 .
c ∂t

(1.9)

We now focus on a region of space outside the source where all the componentsof the stress-energy tensor are zero, i.e.,h̄µν = 0.

(1.10)

This allows us to impose the transverse-traceless (TT) gauge, defined by theconditions:hµν,ν = 0(transverse) ,hµ(traceless) ,µ= 0(1.11)

which assures that the trace-reversed perturbation h̄µν is identical to the physicalperturbation hµν .
The effect of a gravitational wave on a test mass in space can be found bylooking at a plane wave solution of Eq. (1.10) under TT gauge with wave vectorkσ :hµν ( x ) = Aµν eikα x ,α(1.12)

with k α kα = 0 and Aµν kν = Aµ = 0. Choosing the wave vector k i along theµ

z-axis, one can completely parameterize the amplitude Aµν by two independent
polarizations, h+ and h× :




0
 0 0


 0 h+ h×
Aµν :=  0 h× − h+


0 0
0
5

0 


0 
 .

0 


0

(1.13)

h+
hx
0











Figure 1.1: The effect of the two polarizations h+ (top panel) and h× (bottom
panel) of a complete cycle of sinusoidal gravitational wave propagating throughthe page on a ring of test particles.

A gravitational wave with h+ polarization stretches distance along x-axis andshortens distance along y-axis during the first half cycle and does the opposite
during the second half cycle. The h× polarization does the same thing, but
in a coordinate system rotated by 45 degree. Figure 1.1 shows the effect of agravitational wave on a ring of particles.

1.2

Sources of gravitational waves

There are different kinds of “known” and “unknown” sources which generate
GW with different frequency spectra. Interferometric gravitational wave detectors have achieved unprecedented sensitivity to gravitational waves, expectedfrom sources of different kinds. Depending on lifetimes and other characteristics, the sources can be categorized as follows:6Transient: Transient sources are those which produce a very short duration,a few milliseconds to a few minutes, gravitational wave signal. Brief transientsfrom violent events like a coalescing compact binary, composed of massive neutron stars and/or black holes and asymmetrical core collapse in supernovae arein this category. Core collapse in massive stars has long been regarded as likely
to be an important source of gravitational waves.
According to GR a binary system composed of either neutron stars and/or
black holes (NS-NS, BH-NS, BH-BH) emits gravitational radiation at twice theorbital frequency, causing the orbit to shrink. As a result, the components ofthe system spiral in towards one another eventually combining to form a single
object, most likely a black hole, which too radiates gravitational waves as it approaches a stable state. The time when the two objects orbit around one anotheris called “inspiral phase” and the gravitational wave emitted is reasonably wellmodeled. However, the final few orbits and plunge in towards one another isstill “poorly” understood. This phase is called the “merger phase”. Numerical
relativity plays an important role in modeling the waveform associated with themerger phase. The final phase is known as “ringdown phase”, and its gravitational radiation waveform is well modeled in black hole perturbation theory.
However, a large class of transients including the most violent and catastrophic
ones are unmodeled.
Continuous: Periodic sources like pulsars and inspiraling binaries whichemit GW continuously over a long duration without significant change in theircharacteristics are called continuous sources. A rapidly spinning neutron staremits gravitational wave if its shape deviates from axisymmetry or its rotation7

is wobbling. Vibrational modes of a neutron star, if driven by accretion or an
instability, also would emit continuous GW for a long period of time.
Stochastic: Stochastic backgrounds of gravitational waves are produced by asuperposition of a large number of independent, uncorrelated astrophysical andcosmological origin that are not individually resolvable.
Depending on source characterizations, there exists several search methods.
The GW signals from compact coalescing binaries have been “precisely” modeled using post-Newtonian approximations, hence matched filtering technique
is used for extracting the true GW signal buried in strong detector noise. For
the unmodeled transients, most approaches involve looking for excess power incertain frequency bands at certain times. A targeted search, by correlating signals from two detectors with a time dependent phase factor that accounts forthe light travel time delay between two detectors, is the best strategy for thedetection of GW from known radio pulsars. Since the stochastic signals are unmodeled, the best strategy to detect stochastic sources is by correlating outputsof two detectors with appropriate time delays.

1.3

Principle of a GW detector

In a nutshell a laser interferometric gravitational-wave observatory is a sensitive interferometer that measures the change in relative spatial distance of suspended test masses induced by an impinging gravitational wave. Figure 1.2 isa schematic of an earth-based laser interferometric detector. The basic principle
on which a detector works is that of interference of light waves. Each of the two8orthogonal arms of the interferometer contains a pair of mirrored test masses
which form a Fabry-Perot optical cavity [4]. Laser light is divided at the beam
splitter, and reflected from the cavities in each arm. This reflected light acquiresa phase shift which is very sensitive to the separations of the test masses. These
phase shifts are then compared when the light recombines at the beam splitter.
If there is an impinging gravitational wave it will act to stretch or squeeze thesurrounding space-time, effectively causing distances to shrink in one arm while
stretching in the other. This then leads to an overall phase difference in the lightof the two arms.
Suppose the interferometer in Fig. 1.2 is arranged such that its arms lie along
the x and y axes. Suppose further that test masses are separated by a distance L
along the arms. When a gravitational wave with polarization h+ passes throughthe interferometer, the separation between the two masses in the x arm is givenby
L x (t) =p

L gxx L ≈

h+ (t)
1+2

L,(1.14)

whereas, in the y arm,
Ly (t) ≈h+ (t)
1−2

L.

(1.15)

The wave is thus detectable by measuring the length difference between two
arms. The fractional difference in length of the two arms is given by
L x (t) − Ly (t)δL(t):== h+ (t) .
L
L(1.16)

Generally, both polarizations of the wave influence the test masses and hence,δL(t)= F + h + ( t ) + F × h × ( t ) ≡ h ( t ),L
9

(1.17)

test mass
light storagearm

test mass
light storagearm
test mass

test mass

beam
splitter
power recycling
mirrorphotodetector

laser

Figure 1.2: Simplified schematic of an earth-based interferometric detector,
showing its primary components.

where F + and F × are antenna response functions, which depend on the source
location and the detector orientation, which changes with time, as the Earthspins (see Appendix A.1). Since a gravitational wave induces a strain in the detector, h(t) is often referred to as the “gravitational wave strain”. As an example,
the typical strain induced by a GW from a NS-NS binary in the Virgo cluster ish ∼ 10−21 . To be sensitive to such a strain, we need to measure δL = hL ∼ 10−18
m with 4 km long arms.
The principal noise sources affecting interferometer sensitivity are background10

motions of the test masses, mainly due to seismic noise, thermal noise, and noisein sensing the test mass motion due to laser noise, quantum fluctuations of thelight, and fluctuations in the number of residual gas molecules traversing theoptical beams. Among them, the seismic noise and the photon shot noise limit
the sensitivity of the interferometers at low and high frequencies, respectively.
The thermal noise associated with thermal vibration of mirrors and their suspension limits the sensitivity of the detectors in the intermediate band. Figure 1.3
illustrates how the different noise sources contribute to the sensitivity of the first
generation LIGO detectors.

Figure 1.3: An illustration of the sensitivity of the first generation LIGO detectors
along with the different noise sources that impact the sensitivity [9].

Table 1.1 specifies the locations of both existing and planned detectors all11over the world. Currently, there are three LIGO detectors: Livingston, Louisiana
hosts one detector of arm length 4 km and Hanford, Washington hosts two detectors of arm lengths 2 km and 4 km. The French-Italian 3 km long VIRGO
detector is located at Cascina, Italy. In the past few years, these first generation
detectors have successfully operated at or near design sensitivity [5, 6]. Figure 1.4 shows the strain amplitude spectral densities of LIGO detectors during
their S5 science runs. The LIGO/VIRGO detectors are currently in the process of
being upgraded to advanced configurations that are due to come on line in 2015.
The new detectors will improve the strain sensitivity of current instruments bya factor of ten, while also lowering the start of the sensitive frequency range
from 40 Hz down to 10 Hz, resulting in a thousandfold increase in the observable volume of space. Construction of a Japanese second generation 3 km longcryogenic gravitational-wave detector, KAGRA, has been started [7]. KAGRA
was previously known as LCGT for Large-scale Cryogenic Gravitational-wave
Telescope. There is a proposal for an Indian detector (INDIGO), the site forwhich has not been selected yet. The British-German 600 m long GEO is locatednear Hannover, Germany and is mainly used for testing optics for LIGO/aLIGO
detectors [8].
Figure 1.5 shows, along with the noise ñ( f ) curves, the estimated signal
strengths h̃( f ) for various sources. The signal strength h̃( f ) is defined in such a
way that, wherever a signal point or curve lies above the interferometer’s noise
curve, the signal, coming from a random direction in the sky and from a sourcewith a random orientation, is detectable with a false alarm probability of lessthan one percent.

12

Table 1.1: The geographic location and orientation of an international network ofgravitational-wave interferometers. The orientation is the geographic compass
angle, measured clockwise from north, of the line bisecting the arms of thedetector [10].
Detector

LabelLongitude

Latitude

Orientation

LIGO Livingston, LA

L (L1)

90◦ 460 27.300 W

30◦ 330 46.400 N

208.0◦LIGO Hanford, WA

H (H1 & H2)

119◦ 240 27.600 W

46◦ 270 18.500 N

279.0◦VIRGO, Italy

V (V1)

10◦ 300 1600 E

43◦ 370 5300 N

333.5◦KAGRA, JapanK (K1)

137◦ 180 3600 E

36◦ 240 3600 N

16.69◦

INDIGO, India

I74◦ 020 5900 E

19◦ 050 4700 N

270.0◦

Final S5 Strain Sensitivity of the LIGO Interferometers

−16

10

LHO 4km−(2007.03.18)
LLO 4km−(2007.08.30)
LHO 2km−(2007.05.14)
LIGO I SRD Goal, 4km−18√h(f) [1/ Hz]

10

−20

10

−22

10

−24

10

1

10

2

3

10

10

4

10

Frequency [Hz]Figure 1.4: Strain amplitude spectral densities of LIGO detectors during their
S5 science runs. For comparison, the solid magenta curve shows the designreference noise budget of first generation LIGO interferometers.

13

Vela Spindown
Upper Limit

Crab Spindown
Upper Limit

M
o /10

10-23

-I

ins

GO

pir

LI

al 1

00M

pc

Sco X-1

=1
Ω

Mer
ger
50M
B
/
H/B
BH/ o 50MHI
BH
NB LIGO-II
Mero
nsp
ger,
NS
ira
/NS
l 40
z=2
BH
0M M
/BH Inspi
pc erger
Ins ral 3
00M
pir
al,
p
z=0 c; N
WB
.4
S/B
LIG
HI
nsp
O-I
ira
I
l 65
0M
pc

ε=
r= 10 -7
10
kp
c

-9
0

=1
Ω

M
o BH
/BH

Kn
ε= own10 -6 P
, r= ulsa
10 r
kp
c

-1/2
1/2
~
h(f) = Sh , Hz

-7
0

=1
Ω

10

10-22

1
-1
0

10

-24

LMXBs

20Mo/20Mo
BH/BH Merger, z=1

10

20

50

100
200
frequency, Hz

500

1000

Figure 1.5: The noise ñ( f ) in LIGO, Narrow-band (NB) LIGO-II, and Wideband (WB) LIGO-II plotted as a function of gravitational wave frequency f , andcompared with the estimated signal strengths h̃( f ) from various sources [11].

14Chapter 2
Detecting an anisotropic stochastic
background with multi-baseline GW
radiometry
The Cosmic Microwave Background (CMB) was first detected, serendipitously, in
1964 by Penzias and Wilson [15]. This discovery supports that the universe had abeginning permeated by the cosmic microwave background as predicted by earlier theories [16]. The universe is also expected to have Stochastic Gravitational
Wave Backgrounds (SGWBs) of cosmological and astrophysical origin. Incoherent superposition of unmodeled and/or unresolved sources of GWs constitutea SGWB [17]. The GWs from these sources could be too weak to be detected
separately. But if the number of sources that contribute to each frequency bin ismuch larger than one, then searching their GW profiles in LIGO and Virgo data
is an interesting prospect.

15

The LIGO Scientific Collaboration and the Virgo Collaboration are searching for transients, continuous signals and stochastic backgrounds in the datafrom their detectors. The search for an isotropic stochastic GW background hascaught significant attention due to its cosmological significance. This primordial GW background is a direct probe of cosmological inflation [18]. It is wellestablished that the cosmic microwave background is isotropic, which indicatesthat the early universe was highly isotropic. But the assumption of isotropy
for SGWB may not be true. For example, if the dominant source of stochastic
gravitational waves in the frequency band of the earth-based detectors is an ensemble of nearby astrophysical sources, such as low-mass x-ray binaries or, even,coalescing compact objects, then the stochastic background will have a distinct
anisotropy. The astrophysical GW background (AGWB) is of major interest in
this chapter.

2.1

The stochastic GW backgrounds

In the early universe, inflation produced a stochastic background of gravitational waves through the parametric amplification of primordial quantum fluctuations [19]. Unfortunately this background itself is too weak to be directly
detected with existing ground-based detectors. In the period following inflation
there are a number of mechanisms that may have operated to produce additional SGWB from inhomogeneities in the fields that populate the early universe.
For example, first-order phase transitions in the early universe can generate astochastic background of gravitational waves that may be detectable today [20].

16Stochastic Background from Inflation
-6

-8

Ω gw
0.9 K
Blackbody

-10

-12

Inflation

Slope -2

-14
LISA

LIGO

-16

-18
-15

-10

-5

0

5

10

15

f/Hz

Figure 2.1: The spectrum of stochastic gravitational waves in inflationary models
is flat over a wide range of frequencies, and is shown as the solid curve [17]. The
horizontal axis is log10 of frequency, in Hz. The vertical axis is log10 ΩGW definedin Eq. (2.8).

Fig. 2.1 shows the spectrum of a SGWB predicted by inflationary models thatexist today over a wide range of frequencies. The inflationary spectrum rises
rapidly at low frequencies and falls off above the frequency scale f max associated with the fastest characteristic time of the phase transition at the end of
inflation. Because the spectrum falls off as the inverse of the square of frequency
at low frequencies, this is certainly too small to be detectable with either initialor advanced LIGO. Although LIGO has not yet detected a gravitational wave
signal, it has ruled out some early-universe models that predict a relatively large
stochastic GW background as shown in Fig. 2.2. The data rule out models of17early universe evolution with a relatively large equation-of-state parameter, as
well as cosmic (super)string models with relatively small string tension that are
favored in some string theory models. This search for the stochastic background
improves on the indirect limits from Big Bang nucleosynthesis and cosmic microwave background at 100 Hz [18].

Figure 2.2: The 95% Bayesian upper limit of the point estimate of a constant
SGWB spectrum, Ω0 < 6.9 × 10−6 is shown by a solid bar denoted as “LIGO
S5”. The upper limit is calculated from LIGO S5 data, using LIGO S4 result
as a prior for Ω0 and is valid for the frequency band 41.5-169.25 Hz [18]. The
projected advanced LIGO sensitivity (aLIGO) is also shown here.

In addition to the cosmological background, an astrophysical contribution
may have resulted from the superposition of a large number of unresolved astrophysical sources, whose frequency is expected to evolve very slowly compared to the observation time [21]. These can be either short-lived burst sources,18

such as stellar core collapses to neutron stars or black holes, phase transitionor quickly damped oscillation modes in young neutron stars, the final stage ofcompact binary mergers, or long-lived sources, typically non-axisymmetric rotating neutron stars, the early inspiral phase of compact binaries or captures ofstellar mass compact objects by supermassive black holes.
Fig. 2.3 shows the gravitational wave spectra plotted as a function of observed frequency for potentially significant AGWBs and their detectability using
the two upcoming advanced LIGO detectors. The detectability of an AGWB is
determined by the detector sensitivity, which sets the detection horizon. As thesensitivity of LIGO improves and the next generation of GW observatories comeon line, it is plausible that AGWBs can be observed for the first time, providing
new insights into both the sources generating them and the evolutionary historyof the universe.

2.1.1

The detector strain due to a SGWBAs we saw in Chapter 1, a gravitational wave is the propagation of a metricperturbation. In the transverse traceless gauge, the spatial part of the metricperturbations due to a SGWB can be written as a superposition of plane waves
h ab (t, r) =

Z ∞
−∞

df

Z
S2

A
dΩ̂ eab(Ω̂) h̃ A ( f , Ω̂)ei2π f (t+Ω̂·r/c) ,(2.1)

A ( Ω̂ ) are the components of the gravitational
where a and b are spatial indices, eab

wave polarization tensors, A = {+, ×} is the polarization index, and Ω̂ is aunit vector pointing in the direction of wave propagation. Since the GW strain
h ab (t, r) is real, the complex Fourier amplitudes h̃ A ( f , Ω̂) satisfy the reality con19

Figure 2.3: Gravitational wave spectrum as a function of observed frequency
for potential AGWB sources [21]. Theoretical limits for several primordial GW
background models are plotted as horizontal lines and the detection sensitivity from correlating two advanced LIGO type interferometers is shown as thedashed curve (labeled by “2 LIGO ad”).

dition, h̃ A (− f , Ω̂) = h̃∗A ( f , Ω̂).
The polarization tensors can be defined in terms of the spherical polar coordinates, θ ∈ [0, π ], φ ∈ [0, 2π ], on the sky. Let us defineΩ̂ = sin θ cos φ x̂ + sin θ sin φ ŷ + cos θ ẑ ,m̂ = cos θ cos φ x̂ + cos θ sin φ ŷ − sin θ ẑ ,n̂ = − sin φ x̂ + cos φ ŷ ,
20

(2.2)

so that {m̂, n̂, Ω̂} form a right-handed system of unit vectors. The axes are defined as follows: For a fixed but arbitrarily chosen origin of time t = 0, x̂ is directed towards the intersection of the equator and the longitude φ = 0, ẑ points
at the Celestial North Pole, and ŷ is chosen orthogonal to the previous two axes,A ( Ω̂ )forming a right-handed triad. Then, the polarization tensor components eab

are defined ase+ab ( Ω̂ ) = m̂ a m̂b − n̂ a n̂b ,e×ab ( Ω̂ ) = m̂ a n̂b + n̂ a m̂b .

(2.3)

in the aforementioned right-handed orthogonal basis.
Understanding the signal excited in an interferometric detector by a SGWB
is helped by specifying the detector’s location and orientation in the above orthogonal basis. Let the I th GW detector be located at r I (t), and let X̂ I (t) andŶI (t) be the unit vectors pointing along its arms. These three detector locationand orientation vectors are all time-dependent due to the Earth’s rotation. Then,the components of the I th detector tensor are given by1d ab
I (t) =
2hX̂ Ia (t) X̂ bI (t) − ŶIa (t) ŶIb (t)

i

,

(2.4)

and
h I (t) = h ab (t, r I (t)) d ab
I (t) ,

(2.5)

is the strain in it due to the SGWB.
The response of a detector to the polarization component A of a wave incident from direction Ω̂ is given by the antenna response functionA
FIA (Ω̂, t) = d ab
I ( t ) e ab ( Ω̂ ) ,

21

(2.6)

where we assumed the Einstein summation convention over the repeated indicesa and b. Contracting Eq. (2.1) with the detector tensors d I , the GW strain signal
in the I th detector can be expressed ash(t) ≡ h I (t) =

Z ∞
−∞

df

Z
S2

dΩ̂ h̃ A ( f , Ω̂) FIA (Ω̂, t)ei2π f (t+Ω̂·r I (t)/c) ,(2.7)

in terms of the antenna response functions.

2.1.2

Statistical properties of the signal and detector noise

The intensity of a stochastic GW background is usually characterized by a dimensionless quantity, GW spectrumΩGW ( f ) :=1 dρGW,ρc d ln f(2.8)

where dρGW is the energy density of the gravitational radiation in the frequency
range f to f + d f and ρc is the critical energy density follows from the Friedman
equation,c2hḣ (t, r)ḣ ab (t, r)i ,32πG ab3c2 H02
,
=
8πG

ρGW =
ρc

(2.9)

where H0 is the Hubble expansion rate for today. Here, ΩGW is not to be confused with the unit vector Ω̂ defined in Eq. (2.2).
The Fourier components of the strain h̃ A ( f , Ω̂) describing a stochastic GW
background are random variables whose expectation values define the statisticalproperties of the background. Without loss of generality we assume that thesecomponents have zero mean:hh̃ A ( f , Ω̂)i = 0,22(2.10)

where the angular brackets denote statistical average. In the presence of a signal,
the time series of the I th detector’s output x I (t) is a sum of the GW signal h I (t)and the detector noise n I (t):x I ( t ) = h I ( t ) + n I ( t ).

(2.11)

Statistically, the gravitational wave strain h I (t) is uncorrelated with the detector
noise, implying thathh I (t) n J (t0 )i = 0 ,(2.12)

for all I and J. We also assume that the noise is Gaussian with zero mean, i.e.,hn I (t)i = 0, and is uncorrelated in different detectors, namely,1
hn I ( f ) n J ( f 0 )i = δ( f − f 0 )δI J ξ ( I ) ( f ) ,2(2.13)

where ξ ( I ) is the one-sided noise power spectral density (PSD) of the I th detector.
The last assumption is not unreasonable when the detectors are widely separated
across the globe.

2.2

Principle of a GW radiometer

Radiometry or aperture synthesis is a well-known technique in radio astronomy
and CMB experiments [22, 23]. The principle of earth rotation aperture synthesis
can be used well for extracting the anisotropies of GW background using pairsof GW detectors [24]. Since the noise streams in different detectors are uncorrelated, the cross correlation between the outputs from a pair of detectors with anappropriate direction dependent filter is the best statistic for the estimation of
the strength of stochastic signals that may be present in those outputs.
23

A gravitational wave radiometer can be thought of as a pair of GW detectors separated by a distance, termed as the baseline (see Fig. 2.4). Owing to theearth’s rotation, the baseline vector changes direction while keeping its magnitude fixed and hence allows a time-dependent phase delay for a particular
direction of the sky. This delay corresponds to the difference between the timesof arrival between two detector sites of GW signals from that direction. The
cross correlation of the data from the two detectors with appropriate time delays would cause potential GW signals interfere constructively from the same
direction. Whereas signals from other directions will tend to cancel out becauseof destructive interference.

^
−Ω
^
−Ω
^ . ∆ r (t)−Ω∆ r (t)Detector 1

Detector 2

Figure 2.4: Geometry of an elementary GW radiometer [24].

In Fig. 2.4 two detectors are at locations r I (where, I = 1, 2), and the baseline
vector joining the two sites is 4r := r2 − r1 . The unit vector Ω̂ is the direction to the source, which is fixed in the barycentric frame. The phase difference−Ω̂ · 4r(t) between signals arriving at two detector sites from the same direction is shown in the figure. The baseline vector 4r rotates with the rotationof the earth, while keeping its magnitude fixed. The modulation of the crosscorrelation of the SGWB signals in the detector pair can be modelled for each24Ω̂. By searching for this modulation pattern in the detectors’ data one can infer
the nature and sky position of the source. A map of the SGWB can be, thus,
constructed by performing this synthesis for each location in the sky, patch by
patch. The sky resolution can be naively estimated from the diffraction formula4θ ∼ λGW /(|4r| sin θ ), where λGW is the GW wavelength and θ is the angle ofincidence. So, e.g., for λGW = c/1000Hz, |4r| = 3000km, and θ = 900 (overhead source), 4θ ≈ 0.1 radians. A better estimate of the resolution follows fromconsiderations involving the pixel-to-pixel Fisher information matrix [26, 13], inwhich case the solid-angle resolution scales inversely proportional to the squareof the signal-to-noise ratio.
A variety of data analysis strategies to search for an anisotropic GW background have been proposed and implemented in the past [27, 28, 29, 30, 31, 32].
These searches are usually performed in one of two types of bases in the sky,namely, the pixel or spherical harmonic basis. Use of the radiometer technique
for searching a GW background was proposed in Ref. [33] and was implemented
in the pixel basis on data from LIGO’s fourth science run [34]. An elaborate
study of this method, including the maximum-likelihood (ML) estimation of thetrue anisotropy of GW background by deconvolving the observed sky map, waspresented in [24]. Even though the pixel-based search is promising and simpler
to understand, it is not the best basis for probing sources with angular spreads
greater than the angular resolution of the GW radiometer. The spherical harmonic basis is expected to be better suited for detecting such sources [14]. Past
attempts at probing the GW anisotropy in the spherical harmonic basis were essentially studies of the periodic modulation of the observed background in the25detector baselines. Recently, a general ML formalism was developed to search
for the GW anisotropies in any basis, including the spherical harmonic basis, using a network of ground-based GW interferometers [14]. The pixel-based search
is a specific application of this formalism. One of the main goals of this chapter
is to perform a thorough comparison of the expected performances of individual
baselines and the whole network in detecting an astrophysical stochastic gravitational wave background (SGWB) and in estimating its parameters. The pixel
basis is used for this study.
Even though a pixel-based search is optimal for a localized source, the resolution of the source is limited by the length of the radiometer baselines, theorientation of the detectors, and their individual sensitivities. Probing a stochastic GW background with energy distributed across the pixelated sky demands astatistically meaningful integration of the energies received in every pixel. In order to accomplish this, we extend the maximized-likelihood ratio (MLR) statistic
for a single baseline to incorporate a network of detectors or, equivalently, multiple baselines. The rest of the chapter is based on Ref. [13] and devoted to
studying the performance of individual GW detector baselines and the wholenetwork by comparing different figures of merit for their performance, e.g., sensitivity, accuracy in localizing sources, sky coverage, and faithful extraction fromthe data of the sky distribution of a stochastic background.

26

2.3

Optimal search statistic

2.3.1

Cross correlation statistic

Since the targeted source is stochastic, we search for its GW signal by looking
for correlated patterns in the data of two or more detectors after accounting fortime delays and detector responses consistent with a given sky location. Thisis done by cross-correlating the data x I (t) from the detectors, taken in pairs,with a sky-position-dependent time-frequency filter Q̃k (t; f ), labeled by the skyposition index k. The cross correlation statistic combined for the observation
period T for the data x1,2 (t) from two detectors or, equivalently, for a baseline isdefined asT Z ∞S = 4 ∆t ∑k

t =0 − ∞d f x̃1∗ (t; f ) x̃2 (t; f ) Q̃k (t; f ) ,(2.14)

where x̃ I (t; f ) is the short-term Fourier transform of x I (t), over time interval ∆t,and is defined as in Ref. [24] as
xeI (t; f ) :=Z t+∆t/2t−∆t/20

dt0 x I (t0 ) e−2πi f t .

(2.15)

The filter that maximizes the signal-to-noise ratio (SNR) associated with thisstatistic is a scalar, square-integrable function on the sky [24] and, hence, can be
resolved linearly in an appropriate basis, such as a pixel basis or the sphericalharmonic basis. In the former case, k is the pixel index. The SNR can be expressed in terms of the filter function and a multiplicative factor [35, 36]. Oneof the multiplicative factors is the square root of observation time. This multiplicative factor means that we can dig out a weak signal from noise by allowing
a long enough observation time. The observation period could be from one
27

sidereal day to several months.
Let an astrophysical GW background be modeled such that the Fourier components of its GW strain h̃kA ( f ) of polarization A from the kth sky-position obey00hh̃kA∗ ( f )h̃kA0 ( f 0 )i = δAA0 δ( f − f 0 )δkk P(kA) H ( f ) ,(2.16)

where P(kA) is a dimensionless measure of the signal strength, and H ( f ) is its
two-sided power spectral density, with units of Hz−1 [24]. Here, we assume thesignal to be stochastic and uncorrelated in the two polarizations, different frequencies, and different sky locations. In the presence of a signal in the detector
data, the cross correlation statistic is
0

0

k
k
k
k
k
Sk = B+
k0 P(+) + B×k0 P(×) + n ,(2.17)

where the beam function B kA k0 is analogous to the point-spread function that maps
the power in the object (or sky) plane to that in the image plane. The exact formof the beam function is given in Eq. (2.22). The noise in the kth sky position is nk ,and Sk is termed as the dirty map [24]. We define P ( A) as a vector, with P(kA) as
its kth component, and B A as a matrix, with B kA k0 as its (k, k0 )th element.

2.3.2

Detection statistic

To get a single detection statistic, one must combine the measurements of Sk for
all k. When the detector noises are Gaussian and uncorrelated, an assumption
borne out to sufficient approximation for our purposes, the nk are Gaussian witha nontrivial covariance matrix, N, determined by the beam functions. The exact
form of N is discussed below.
28

If an astrophysical GW background signal, characterized by the pixel-strength
vector P , is present in the data, then the probability density function of the radiometer output S is given by
1
p(S|P ) = (2π )− Npix /2 exp[− ((S − B · P ) T · N−1 · (S − B · P )2+Tr[ln N])] ,(2.18)

whereas in the absence of a signal it is1p(n) = (2π )− Npix /2 exp[− (n T · N−1 · n + Tr[ln N])] .
2(2.19)

By the Neyman-Pearson criterion, the optimal detection statistic is the likelihood
ratio p(S|P )/p(n) [26].
For an unpolarized background from a source distributed across multiple
pixels and quantified by the signal-strength vector P = P (+) = P (×) , the loglikelihood ratio maximized over P ≡ kP k is
λ = q

= q

Sk (N−1 )kk0 (B · P̂ )k0(B · P̂ )q (N−1 )qr (B · P̂ )r
Sk P̂ kP̂ q B

qr

P̂ r

,

,

(2.20)

where P̂ is the unit vector along P . The beam matrix for an unpolarized source(see Fig. 2.5) is given by

B pq = B+ pq + B× pq ,(2.21)

T= 8 ∆ f ∆t ∑ Γ(Ω̂q , t)Γ(Ω̂ p , t)
"

×<

t =0
fu
2πi f (Ω̂q −Ω̂ p )·∆r(t)/c∑ef = fl29#
G (t, f ) ,(2.22)

where Γ(Ω̂q , t) is the time-varying baseline antenna-pattern, and G (t, f ) is ameasure of the spectral strength of the source relative to the baseline’s noisePSDs:
Γ(Ω̂, t) := F1+ (Ω̂, t) F2+ (Ω̂, t) + F1× (Ω̂, t) F2× (Ω̂, t) ,G (t, f ) :=H2 ( f ).
ξ (1) (t, f ) ξ (2) (t, f )(2.23)(2.24)

Therefore the fractional energy density in gravitational waves ΩGW ( f ) is related

Figure 2.5: Illustration of a GW radiometer beam pattern at declination +120 for
the LIGO detectors at Hanford and Livingston (with white noise, upper cutoff
frequency of 1024 Hz, H ( f ) = constant and observation time of one siderealday). For low declination, the beam is shaped like the figure 8, while in thehigher declination, the 8 smoothly turns into a tear drop [24].

to the strain power spectrum through 1
8π 2 3ΩGW (| f |) =| f | H( f )
3H02

Z
S2

dΩ̂ P (Ω̂) .

(2.25)

1 For an isotropic background the directionality is set to unity, i.e., P ( Ω̂ ) = 1, which simplifies

Eq. (2.25) byR

S2 dΩ̂ P ( Ω̂ ) = 4π.

30

In the weak-signal limit the noise-covariance matrix is approximately equal
to the beam matrix,
Npq ≈ B pq .

(2.26)

Its diagonal elements inform us about the sensitivity of the network to the different pixels in the sky for an SGWB with PSD H ( f ).
The statistic λ is the maximized (log-)likelihood ratio for a single-baseline
SGWB search and is the same statistic introduced in Appendix C of Ref. [14].
Here, it is expressed specifically in terms of quantities defined in the pixel basis.
It has zero mean and unit variance in the absence of a signal. When a signal ispresent in the data and its parameters are matched exactly by the template’s, themean of the statistic is

hλi = Pq

(B · P̂ )k (N−1 )kk0 (B · P̂ )k0 .

(2.27)The variance of the statistic remains unchanged. One can extend this singlebaseline statistic to the case of a multibaseline network. That statistic arises
directly from maximizing the log-likelihood ratio for a network and is given by
N

λN

= q

−1
k
k
b
∑I=
1 SI ( NI )kk0 (B I · P̂ )
N

0

−1
b
q
r
∑I=
1 (B I · P̂ ) ( NI )qr (B I · P̂ )

,

N

= q

I k
b
∑I=
1 Sk P̂

Nb
q I r
∑I=
1 P̂ Bqr P̂

,(2.28)

where I is the baseline index and the subscript N highlights that this MLR
statistic is for a network of baselines.
The MLR statistic is a detection statistic for SGWBs in the same manner as
the standard matched-filter statistic is for deterministic GW sources. The latter31

is also obtained by maximizing the likelihood ratio with respect to the strengthof the deterministic source. Searching for a signal from a deterministic source
involves maximizing the matched-filter statistic over a bank of templates definedon the signal’s parameter space. For SGWBs, as well, the detection statistic
can be the MLR, maximized further with respect to different SGWB models
given by P̂ , perhaps parametrized by a smaller number of parameters than
the number of components of P̂ . This is in contrast to the existing searches
for anisotropic GW backgrounds. Past dirty-map-based searches precluded thepresence of a signal by demonstrating that the map is consistent with a Gaussian
distribution, up to statistical fluctuations allowed by the number of independent“samples” on the sky [34]. However, they did not provide a confidence level
for the presence or absence of a broadband or spatially extended signal. A
better approach is to solve the inverse problem in an orthogonal basis, namely,
the pixel [24] or spherical harmonic [14] basis. This yields an estimate of the
background, i.e., a “clean” (deconvolved) map and the corresponding noisecovariance matrix. However, this approach depends heavily on how well-posed
the inverse problem is and how accurately it can be implemented numerically.
Consequently, a detection statistic constructed on the deconvolved data can beaffected by similar maladies. To work well, the inverse problem requires that the
network of interferometers is sufficiently nondegenerate, which is not always thecase. Indeed, the deconvolution procedure can enhance spatial noise correlations
and, sometimes, even introduce artifacts, thereby adversely affecting parameter
estimation and signal detection by such a procedure.
As we prove here, the detection problem does not require a well-posed in-32

verse problem and exists even for a degenerate network. A detection statistic
is best defined on the dirty map, as opposed to the clean map. As an added
advantage, a dirty-map-based statistic is faster to compute, since it obviates thecomputational overhead required for obtaining the clean map. While it is possible to use an arbitrary sky model, such as the one-dimensional basis P̂ , andestimate the strength of the SGWB for that particular model, the MLR statistic
in Eq. (2.20) provides a well-understood construct that can be maximized over a
set of parameters for selecting the model that best fits the data.
To elaborate further on the way the new statistic works, let us consider theexample of a directed pixel-space search, which is performed for only one source
and assumes that the angular extent of the source is, at most, one pixel. In
the standard radiometer search [25, 24], the dirty map Sk is computed for eachpixel k in the sky. The (signal part of the) dirty map is generally peaked at thesource pixel and has broad structures, including large negative patches, around
it. One way of inferring the presence or absence of a source in this image requires
deconvolving it. However, as we show later, deconvolution of a relatively weaksource can result in a clean map with significant errors, especially when the sky
is divided into around 3000 pixels or more. (A network resolution of several
square degrees requires a few thousand pixels across the sky.) Also, computing
the noise-covariance matrix can be numerically challenging. The MLR is a good
choice in this situation, since it combines all the pixel values to provide a single
number for the detection statistic that is simple to use in drawing inferences on
the presence or absence of a signal in the network data.
If a parametrized model of the background is available, one can construct the33likelihood-ratio statistic from the dirty map and maximize it over the parameter
space. The maximized likelihood-ratio statistic can also be used to perform amore advanced blind search, where no prior information is available about angular distribution of the power in the SGWB. For each basis component, one canassert that only that basis component is present in the signal and compute thestatistic with the corresponding sky model. Thus, using the dirty maps of “point
estimate” or SNR obtained by the existing radiometer search, our prescription
takes one step forward and can provide a map of likelihood ratios, which is statistically a more robust and meaningful quantity, given a set of highly correlated
observations.
The construction of the MLR statistic on a dirty map is simple. Eq. (2.20)
shows that it is the scalar product of the observed map Sk and a sky-modeldependent normalized “template.” The template is proportional to B · P , whichis the expected signal in the dirty map for a sky model P̂ . The inverse of thenoise-covariance matrix is the metric in the pixel space. The sky model canbe defined in a straightforward way. For instance, to search for a point source
localized to a single pixel, one would use a P̂ with all but one component,
namely, the component corresponding to that pixel, set to zero. Indeed, B · P issimply the point-spread function of the pixel with the nonzero component of P̂ .
Also notice that the inverse of the noise-covariance matrix, being proportional
to the beam matrix, cancels out algebraically in the expression for the MLR.
Therefore, unlike for deconvolution, for MLR construction the computation ofthis matrix is not needed. Otherwise, the latter procedure would have been
computationally similar to solving the inverse problem, avoiding which is one

34

of the main motivations for this work.
To complete the discussion, we note that the construction of a MLR statistic
is not limited to dirty maps and can be implemented for clean maps. A clean
map can be expressed as
e = P + nc ,P

(2.29)

where P is the true sky map, and nc is Gaussian noise with covariance Σ, whichis related to the dirty-map noise-covariance matrix through the relation Σ =(B T N−1 B)−1 . Therefore, following the same procedure as that for the dirty
map, one can write the MLR statistic for a clean map ase · Σ −1 · PP
λc = √,P · Σ −1 · P(2.30)

and, thereby, obtain model-based or blind likelihood-ratio maps.

2.3.3

Performance of optimal detection statistic

We numerically study the performance of the optimal statistic and compare withthe existing method. We use the LIGO 4km detectors located in Hanford (H1)and Livingston (L1). Unless otherwise stated, the noise PSDs of all detectors
are taken to be their (smoothed) first-generation design sensitivities as shownin the solid curves of Fig. 2.6. The frequency band considered here spans 40 −1024Hz, with a bin size ∆ f = 1Hz. The source PSD is taken to be a constant,H ( f ) = 1.516 × 10−48 /Hz. Note that the spectral index of the source PSD has asignificant effect on the resolution of the network. Predictions from astrophysicaland cosmological models suggest the nominal range of the spectral index to be
between −3 and 1. The higher the spectral index, the higher the resolution, and35

the more computationally expensive the directed search. We take the spectral
index to be zero here by setting H ( f ) as a constant.
LIGO
VIRGO
aLIGO
aVIRGOET−21√h(f) [1/ Hz]

10

−22

10

−23

10

−24

10

−25

10

1

10

2

3

10

10

4

10

Frequency [Hz]Figure 2.6: Designed noise amplitude spectrum of the initial LIGO and VIRGO(solid curves), advanced LIGO and VIRGO (dashed curves), and the Einstein
Telescope (dotted-dashed curve) [39, 40, 41].

The directed search is performed by dividing the (simulated) strain data fromall detectors into segments with a duration of 192 sec. The noise is taken to
be stationary. The sky is tessellated into 3072 pixels by using the HierarchicalEqual Area isoLatitude Pixelization [37, 38]. The choice of the signal integration
duration is taken to be a sidereal day, which leads to the azimuthal symmetry of
the baseline sensitivities and sky resolutions. The justification for choosing theabove parameter values can be found in [24].
We first construct simulated data sets of two kinds, one with only noise andthe other with a weak signal from a “polar-cap” source added to that noise. The
sky map of the latter case is shown in the first plot in Fig. 2.7. We make dirty
36maps for these two cases using single-pixel source templates for each of the 3072
pixels. These maps are essentially the maps of SNR for the directed search, as
can be seen by substitutingP̂ k =
 1 for target pixel ,(2.31)

 0 for remaining pixels ,into Eq. (2.20) [24]. The dirty maps for both cases look very similar, and only
one of them, namely, the one for the polar-cap source, is shown in the second
plot of Fig. 2.7. Not surprisingly, it is also similar to the dirty map presented in[34] for real data. Indeed, the MLR values over the dirty-map pixels appear tofollow a normal distribution, as shown in Fig. 2.8. This is consistent with the
distribution presented in [34]. Following that reference, we also plot the 1σ error
envelope around the Gaussian fit for 400 degrees of freedom and observe thatthe tops of every bar in the histogram lie within that envelope. The important
point to note here is that, for both noise-only and weak-signal (polar-cap) data
sets the distributions in Fig. 2.8 are very similar and consistent with a normal
distribution.
If we now pretend that we know the broad shape of the GWB sky and use
the P̂ k of the polar-cap signal as our template for computing the MLR statistic,
we find that the above two cases can be distinguished better: In the noise-only
case, the MLR statistic equals −0.023, while with the weak polar-cap signal it is1.400, which is significantly larger than the former, as explained below. To corroborate this claim, we computed the MLRs on an ensemble of 4000 realizationsof noise, with and without the weak polar-cap signal. In Fig. 2.9, we show the
distribution of the MLR statistic (λ) for noise-only (top) and weak-signal (bot37

(a) Injected map(b) Dirty map(c) Clean map

Figure 2.7: The P k sky map of a weak polar-cap GWB source is shown in (a).
The dirty (b) and clean (c) maps for this source were constructed using theradiometer algorithm for the LIGO H1L1 baseline. The last two maps for thisweak source are visually very similar to those for the noise-only case (which is
not shown here).

tom) cases. Clearly, the noise-only λ values are normally distributed with a zeromean, and the weak-signal λ values are normally distributed with a mean of∼ 1.4σ, where σ ≈ 1. This experiment confirms that, given our assumptions
on the detector noise and the signal, the MLR statistic can considerably enhancethe detectability of a weak diffuse stochastic background, if a reasonable model38

Histogram of λ (40 bins)
3500

sky area (deg2)

3000
2500
2000
1500
1000
500
0
−4

−3

−2

−1
0
1
2
λ
Data
Ideal Gaussian (mean=0,sigma=1)
1−sigma error for 400 indep. points

3

4

3

4

(a) Noise
Histogram of λ (40 bins)
3500

sky area (deg2)

3000
2500
2000
1500
1000
500
0
−4

−3

−2

−1
0
1
2
λ
Data
Ideal Gaussian (mean=0,sigma=1)
1−sigma error for 400 indep. points(b) CAP

Figure 2.8: Histograms of the dirty maps for the noise-only data set (top) andweak polar-cap signal (bottom), which is depicted in Fig. 2.7, are shown here.
These two histograms are consistent with that of zero-mean Gaussian data (solid
curved line), up to 1σ errors: Following [34], the 1σ error boundaries for 400
degrees of freedom have also been overlaid for a consistency check.

39

Histogram of λ (40 bins)
4000

sky area (deg2)

3500
3000
2500
2000
1500
1000
500
0
−4

−3

−2

−1
0
1
2
3
λ
Data
Ideal Gaussian (mean=0,sigma=1)
3−sigma error for 4000 indep. points

4

5

4

5

(a) Noise
Histogram of λ (40 bins)
4000

sky area (deg2)

3500
3000
2500
2000
1500
1000
500
0
−4

−3

−2

−1

0

λ

1

2

3

Data
Gaussian fit (mean=1.4045,sigma=1)
3−sigma error for 4000 indep. points(b) CAP

Figure 2.9: Monte Carlo simulations with 4000 noise realizations were performed
to study the performance of the MLR statistic. The distribution of the MLR
statistic obtained for dirty maps generated from noise-only (top) and weak polarcap injection (bottom) are shown. Clearly, the MLR statistic detects the signal at∼ 1.4σ level. Notice that the “CAP” histogram has shifted to the right.

40

of the background is available. This is true even when the distribution of thedirty-map pixel values for that source is close to zero-mean Gaussian.

2.4

Performance of multi-baseline radiometers

In this section we define a set of figures of merit to compare the performance ofa network of baselines with that of its individual baselines for a directed searchof a SGWB.
A single figure of merit may not suffice in capturing all the attributes of abaseline or a network of detectors benefiting a SGWB search. A certain baseline
or network configuration can have good sensitivity if the detectors are optimallyoriented, but have poor resolution if they are proximally located. Optimallyoriented detectors may be very sensitive to certain anisotropy modes, but insensitive to others, making the estimation problem highly degenerate. On the other
hand, a network of detectors that are oriented differently may have moderate,
yet uniform, sensitivity to all spherical harmonic modes of a SGWB, thereby mitigating the ill-posedness of the estimation problem. Such a network, however,will perform worse than one where all the detectors are aligned similarly in a
low-frequency, all-sky isotropic search. Therefore, the relevance of a figure of
merit is determined by the kind of search one is undertaking. Here, we propose a set of figures of merit that are relevant to current searches of anisotropic
stochastic background and that are special cases of the general ML frameworkpresented in [14].
For the numerical simulation studies below, we use the same detector char-

41

acteristics as mentioned in Scc. 2.3.3, but we now include the Virgo detector (V1)in Cascina, Italy to construct a three-baseline network. The baselines and their
network are named by concatenating the symbols for the participating detectors;e.g., the Hanford-Livingston baseline is termed as H1L1, and the network of theabove three detectors is termed as H1L1V1. A similar study for the second andthird generation detectors is presented in Appendix B.1.

2.4.1

Sensitivity

The first figure of merit is the “sensitivity” of a network and is motivated bya similar quantity defined in [42] for the all-sky isotropic search. In practice, agreater sensitivity implies a better confidence level, at which detection can bemade or upper limits can be inferred.
We define the single-baseline sensitivity for a directed search as the expectation value of the MLR in Eq. (2.27) for a SGWB source with P set to unity,
Sensitivity ==qq(B · P̂ )k (N−1 )kk0 (B · P̂ )k0 ,P̂ k Bkk0 P̂ k0 .

(2.32)The sensitivity can be expressed in the spherical harmonic basis as follows:
Sensitivity =q

Plm B lm l 0 m0 Pl 0 m0 ,(2.33)

where

Plm =
Blm l 0 m0 =

Z Z

Z

∗
dΩ̂ P̂ (Ω̂) Ylm(Ω̂) ,∗dΩ̂ dΩ̂0 Ylm(Ω̂) B(Ω̂, Ω̂0 ) Yl 0 m0 (Ω̂0 ) .

42(2.34)(2.35)

Owing to the statistical independence of the baselines, the multibaseline sensitivity squared is the sum of squares of the individual baseline sensitivities, as
was also noted for the isotropic-background baseline sensitivities in [42]:Sensitivity2N = ∑ Sensitivity2I .

(2.36)I

For an unpolarized background from a single pixel, say, labeled k, and withP̂ r = δr(k) , the sensitivity expression simplifies to
Sensitivity(k) =q

B q ( k ) ( N −1 )qr

B r (k) =p

B (k)(k) .

(2.37)

Unless otherwise mentioned, there is no sum over the repeated parenthetic indices in this chapter.
In the top panel of Fig. 2.10, we compare the sensitivities of the baselines andthe whole network as a function of declination. (As noted above, the sensitivities
are azimuthally symmetric.) For a fair comparison, we also replot them after
weighting them with the cosine of the latitude, in effect, to assign equal weightto every pixel on the sky. It is clear that the H1L1 baseline has much better
sensitivity due to the similar orientations of the two detectors. Still, inclusionof Virgo, which is oriented quite differently relative to H1 and L1, improves
the sensitivities of the network by ∼ 10% (which corresponds to an increase in
the observational volume by ∼ 30%), especially in the regions where the H1L1
baseline does not perform well. However, this network improvement is highly
superseded by all other performance improvements indicated by corresponding
figures of merit introduced in this section.
In Fig. 2.11 we plot the narrowband (5Hz) sensitivities at two locations,
namely, the Celestial North Pole (top) and the equator (bottom). Performance43

improvement of a network for a narrow band search at high frequencies is better
than the (frequency integrated) broadband search.

12

Sensitivity

10

8

H1L1V1
H1L1
L1V1

6

H1V1
4

2
−80

−60

−40

−20

0

20

40

60

80

60

80

Declination (degrees)
10

Area weighted sensitivity

98
H1L1V1

7

H1L1

6

L1V1
H1V1

5
4
3
2
1
0

−80

−60

−40

−20

0

20

40

Declination (degrees)

Figure 2.10: The sensitivities (top) and their area-weighted counterparts (bottom)of three different baselines and their network are plotted as functions of thedeclination of a single-pixel SGWB source. The source PSD (H ( f ) = 1.516 ×
10−48 strain/Hz) is chosen such that it has maximum SNR= 10 in the H1L1
baseline. (Note that the source parameter P is set to unity for these plots.) The
signal band considered here is 40-1024 Hz.

44

H1L1V1
H1L1
0

10

L1V1

Sensitivity

H1V1

−1

10

200

400

600

800

1000

Frequency (Hz)H1L1V1
H1L1
0

10

L1V1

Sensitivity

H1V1

−1

10

200

400

600

800

1000

Frequency (Hz)Figure 2.11: The sensitivities of three single baselines and their multibaseline
network plotted as functions of the central frequency of the source band. The
source is chosen to have a constant H ( f ) = 1.516 × 10−48 strain/Hz and a band
width of 5Hz. The top panel represents the sensitivities at the celestial poles,
and the bottom panel represents those at the celestial equator.

2.4.2

Sky coverage

In a directional search, the main advantage of a network lies in the fact that it
vastly improves sky coverage, which, in turn, leads to better parameter estima45

tion, including localization and sky-map reconstruction considered later in thissection. In this subsection, we illustrate the advantage of using a network ofdetectors, as compared to using its individual baselines, to this end.
In general, one radiometer baseline cannot sample the whole sky uniformly;
the measurement errors in some parts of the sky are much worse than thosein the other parts. Introduction of new baselines with different orientations improves filling in these “holes” by scanning the sky with different antenna-patternfunctions. In the first three plots in Fig. 2.12, we show the standard deviation inmeasuring the dirty map by the three individual baselines, namely, H1L1, H1V1,
L1V1, and their network H1L1V1 2 . The azimuthal symmetry mentioned before
is explicitly observed here. The H1L1 baseline has the least deviation at mostdeclinations, due to the optimal orientations of the H1 and L1 detectors. Again,
H1V1 and L1V1 baselines have low deviation in the regions where H1L1 does
not perform well. Since the dirty maps from different baselines are operationallycombined with an inverse-noise-variance weight, the harmonic mean of the variances provides the effective variance of the combined dirty map. The last plot
in Fig. 2.12 shows the effective deviation for a network of detectors. Clearly, thedeviation now has smaller spread and also, by construction (harmonic mean),the deviations are smaller than those for the individual baselines.
Most importantly, a network also complements the single-baseline observations in terms of angular resolution. The beam functions for each radiometerbaseline are highly asymmetric, which means that a given position on the sky is
2 The absolute scales of the plots are not important since the emphasis here is on relative

performance of the different baselines and their network.

46

(a) H1L1(b) L1V1(c) H1V1(d) H1L1V1Figure 2.12: The standard deviation of dirty maps measured by the three LIGOVirgo baselines and the full network for a constant H ( f ) are plotted here. The
aim of this figure is to show how individual baselines complement each other,
thereby making the “effective deviation” (i.e., the square root of the harmonicmean of variances of the individual baselines) of the combined map observedby the network nearly uniform. Note that the color scale in the network plot haslesser spread than the individual baselines. The absolute scale of the maps depends on the normalization of the filter, and only the relative scale is important
here. The azimuthal symmetry is present because we are considering a wholesidereal day’s observation, with stationary noise.

47

probed with quite different angular resolutions in the tangential directions. To
illustrate this aspect, the typical beam functions for the three LIGO-Virgo baselines in the direction of the Virgo cluster are shown in Fig. 2.13. If we consider
the beam for the H1L1 baseline, the sensitive part of the beam is similar to ahighly eccentric ellipse, suggesting that the angular resolution along the minor
axis is much finer than that along the major axis. The beams for the baselines in a
network involving the Virgo detector provide better resolution due to the longer
baselines: The beams are finer along the major axis of the H1L1 beam, thus
complementing the H1L1 observation, which is a major motivation for using a
network. This, in turn, improves the condition number of the Fisher information
matrix, thereby, reducing the numerical errors in the anisotropy estimation problem at “high” resolution, i.e., near or beyond the diffraction limited resolution,and significantly improves source localization accuracy.
Singular value decomposition of the Fisher information matrices provides amore quantitative verification of the above claim. Figure 2.14 shows the singularvalues of the Fisher matrices for the individual baselines and the whole network.
The LIGO baseline has very small singular values at higher resolutions (dashedcurved line), which implies that estimation of anisotropy at those resolutions is
an ill-defined problem. The network reduces the difference between high andlow singular values and regularizes the inverse problem at high resolution (solid
curved line).

48

Figure 2.13: The beam functions for the three LIGO-Virgo baselines for H ( f ) =constant are shown here. This figure illustrates that different baselines also complement each other in terms of angular resolution along different tangential sky
directions.

2.4.3

Parameter accuracy

An important figure of merit for a directed search is how well a point source
can be localized or its other parameters be constrained. In a noise-dominatedmapping experiment, it is not easy to identify sources in the observed images.
If there are candidate sources that have been modeled by other astronomical
observations, one can utilize that information to detect or constrain parameters

49Singular Values of the Fisher Matrices of SGWB Directed Search
6

10

4

singular value

10

2

10

0

10

−2

10

H1L1
H1V1
L1V1
H1L1V1

−4

10

−6

10

0

10

1

2

10

10

3

10

index

Figure 2.14: Singular values of the Fisher matrices for individual baselines andthe whole network are plotted in this figure. The LIGO H1L1 baseline (dashedcurved line) has very small singular values at high resolutions; therefore, estimation of anisotropy at those resolutions is an ill-defined problem. The network(solid curved line) makes the singular values much more uniform, thereby regularizing the inverse problem.

of such sources. If the source was very accurately modeled, the optimal strategy
would be to design a specific search focused on that source. But, in practice,with very limited knowledge of sources, the optimal strategy would be to vary
the parameters within a reasonable range and maximize the log-likelihood ratio.
One of the main advantages of the MLR statistic is that it allows estimationof parameters of the source, given a model. For example, if there is an extended
source, such as a cluster of galaxies, with an angular scale comparable to theresolution of the radiometer, and there exists a reasonable model for its mass

50

distribution, one can maximize the log-likelihood ratio to find the center of thecluster 3 . Even for a blind search, this method may prove to be advantageous
to perform a finer search around the poorly estimated parameters of a potential
candidate source.
In this section, we assess the accuracy with which a pointlike (single-pixel)
source can be located using a network of GW detectors as compared to its individual baselines. The parameter estimation accuracy is deduced from the elements of the Fisher information matrix [26].
For an unpolarized background from a single pixel, labeled k, and with P̂ r =δr(k) , the single-baseline detection statistic follows from Eq. (2.20) to beq

S p (N−1 ) pq B(k),λ(k) = q
B(rk) (N−1 )rs B(sk)(2.38)

which can be interpreted as the inner product of the data, S, and a unit-norm
template B̂ k . Hence, the match [43] between the unit-norm templates for the kth
and the k0 th pixels ispq

B(k) (N−1 ) pq B(k0 )
q 0
,
M = q
0
B(rk) (N−1 )rs B(sk) B(rk0 ) (N−1 )r0 s0 B(sk0 )= q

B(k)(k0 )q,
B(k)(k) B(k0 )(k0 )(2.39)

where the inner products are all defined in terms of N−1 . Define Θ(k) ≡ {µk , φk }as the pixel coordinates, where µk ≡ cos θk , with θk and φk being the declination
and right ascension of the kth pixel, respectively. Since the match has a maximum
3 Note that the image of the source may be too faint to be visually prominent in the dirty or

clean maps.

51

value of unity at k0 = k, one can expand M in a Taylor series about ∆µ(k) = 0
and ∆φ(k) = 0 as
∂2 Mµ∂Θ(k0 ) ∂Θν(k0 )1
M ≈ 1+2≈ 1−

Γ(k)µν

(SNR)2(k)!
∆Θ(k) ∆Θν(k) ,µΘ(k0 ) = Θ(k)∆Θ(k) ∆Θν(k) ,µ(2.40)

where Γ(k)µν ≡ Γµν (Θ(k) ) are the components of the Fisher information matrix
Γ(k) , and (SNR)(k) is the signal-to-noise ratio in the kth pixel.
For large SNR, the error variance-covariance matrix obeys


1
Γ−
(k)

µν

≈

D

µ
∆Θ(k) ∆Θν(k)E.

(2.41)The estimation error in the measurement of the sky-position solid angle (in steradians) is given by [44] 4
r
D
ED
E D
E2
∆Ω(k) = 2π(∆ cos θ(k) )2 (∆φ(k) )2 − ∆ cos θ(k) ∆φ(k) .

(2.42)The Fisher information matrix for multiple baselines is just the sum of the Fisher
matrices for the individual baselines,hΓ(k)µν

iN= ∑ ΓI(k)µν ,(2.43)

I

where I is the baseline index, and ΓI(k)µν is the Fisher information matrix of the

I th baseline, as given in Eq. (2.40). Hence, the error variance-covariance matrix
for the network is
h

Γ(k)

i−1 µνNDEµ≈ ∆Θ(k) ∆Θν(k)N,(2.44)

4 See Appendix A.2 for the detailed expression in terms of beam function and its derivatives.

52

for large SNR. Therefore, the 1σ estimation error in solid angle for locating apixel source with the multibaseline network is expressed as
rD
E
h
i
E D
D
∆Ω(k)

N

= 2π

(∆ cos θ(k) )2N(∆φ(k) )2N− ∆ cos θ(k) ∆φ(k)E2N. (2.45)

Note that this error diminishes with SNR as 1/SNR2 , i.e., localization is moreaccurate at higher SNR.
We present the source-localization errors for the individual LIGO-Virgo baselines and the network in the top panel of Fig. 2.15. We also show the corresponding area-weighted plots obtained by multiplying these errors with the cosine of
the latitudinal angle in the bottom panel of Fig. 2.15. The network clearly outperforms individual baselines by about 1 order of magnitude or more for almost
all declination angles.
The primary focus of this analysis was to obtain statistics that are based on
dirty-map constructs. However, it is straightforward to extend it to be applicable to clean maps. Similarly, although we considered broadband signals, it ispossible to easily extend our study to narrow band signals.

2.4.4

Map making

Finally, we compare the quality of sky maps made by the individual baselines
and their network since they are among the primary products of anisotropic
searches. Here, we consider two figures of merit, namely, the MLR statistic andthe normalized mean square error (NMSE) for comparing maps. To compare thedirty maps, we use the MLR statistic introduced earlier, and, to compare cleanmaps, we use both the MLR statistic and the NMSE (defined below).
53

−1

10

H1L1V1
H1L1
L1V1

Estimation error (sr)

−2

10

H1V1

−3

10

−4

10

−5

10

−80

−60

−40

−20

0

20

40

60

80

Declination (degrees)−1

10Area weighted estimation error (sr)H1L1V1
H1L1
−2

L1V1

10

H1V1
−3

10

−4

10

−5

10

−6

10

−80

−60

−40

−20

0

20

40

60

80

Declination (degrees)

Figure 2.15: The 1σ error (top) and area-weighted 1σ error (bottom) in the solid
angle for locating a source in the sky with three single baselines and the wholenetwork. The network accuracy is better by about 1 order of magnitude or moreat most of the declinations. Note that an error of 1sr ' 3282.81 sq-degrees, andthat the error here decreases as 1/SNR2 .

For an unpolarized and anisotropic gravitational wave background, the maximum likelihood estimators of the signal-strength vector are given by0
Pe k = (B −1 )kk0 Sk ,54

(2.46)

e k are components of thewhere Sk are components of the dirty map (2.17) and P
e , are the values of P that
deconvolved (clean) map. Note that the clean map, P
maximize the statistic. We extend this single-baseline analysis to a multibaseline
N

k =
k
b
one by simply adding the dirty maps and beam matrices as SN∑I=
1 SI andNbkkBNk0 = ∑I=1 BI k0 . So, the maximum-likelihood estimators for a multibaseline

network are given byk0
Pe k = (B N −1 )kk0 SN.

(2.47)

We simulate the data with signal as [24]x̃1∗ (t, f ) x̃2 (t, f ) = hh̃1∗ (t, f )h̃2 (t, f )i + ñ1∗ (t, f )ñ2 (t, f ) ,hh̃1∗ (t, f )h̃2 (t, f 0 )i = δ f f 0 H (| f |) ∑ P i γ(Ω̂i , t, | f |) ,(2.48)iwhere P i is the injected source strength at the ith pixel, and γ is the directiondependent overlap reduction function. We use a conjugate gradient (CG) method
to solve the set of linear equations (2.46) and (2.47).
The mismatch between two maps, injected and estimated, is measured using
NMSE :=

e − P |2
|P
,
|P |2

(2.49)

which is called normalized mean square error.

2.4.4.1

Dirty maps and clean maps

We performed numerical comparison of map-making performance for two typesof toy sky patterns – (i) extended, multideclination sky (Fig.s 2.16 and 2.17),
mimicking the (partially masked) image of the sky constructed by the WMAP
satellite [45], where essentially a modified galactic structure stands out; and (ii) a
relatively localized source peaked at the north pole (Fig.s 2.19 and 2.20).
55

(a) Injected map(b) H1L1(c) L1V1(d) H1V1

Figure 2.16: Toy model of an extended source is shown in (a). Dirty maps madefrom simulated data containing signal from that source are shown in the last
three panels for the three LIGO-Virgo baselines.

In both cases, the dirty maps from different baselines are quite distorted compared to the injected maps. (Compare the last three plots with the first one in
Fig. 2.16 and in Fig. 2.19). However, the deconvolution procedure yields reasonably resolved maps for all the baselines (Fig.s 2.17 and 2.20), signifying thatnone of the beam matrices are completely degenerate. The clean maps from the
network are, however, of better quality, as can be seen from the corresponding
low NMSE. To demonstrate this visually, we also show the difference between
56

(a) H1L1 (NMSE=0.4311)(b) L1V1 (NMSE=0.4171)(c) H1V1 (NMSE=0.4735)(d) H1L1V1 (NMSE=0.2206)

Figure 2.17: Clean maps obtained by the deconvolution of the dirty maps ofFig. 2.16, using 20 CG iterations, are shown here.

the clean and injected maps for the respective cases in Fig.s 2.18 and 2.21. As expected, the difference maps for the network look less noisy and uniform acrossthe sky than the individual baselines.

2.4.4.2

Maximized-likelihood-ratio statistic

We finally compute the MLR statistic introduced in Scc. 2.3.2 as a figure of merit
and, also, to demonstrate how this statistic can be powerful in identifying signal
in noisy maps. The MLR statistic for both dirty and clean maps for the two types
57

(a) H1L1(b) L1V1(c) H1V1(d) H1L1V1Figure 2.18: Difference between the clean maps of Fig. 2.17 and the injected mapof Fig. 2.16(a).

of sources considered here have been listed in Tables 2.1 and 2.2.
It is intriguing to note that, for both dirty and clean maps, one obtains similar values. This suggests that a deconvolution effected with only a few tens ofconjugate-gradient basis vectors does not cause a significant amount of information loss.
To understand the significance of the MLR statistic in the present context,
we perform two more exercises. First, we study the no-injection case; that is,
we make dirty maps of simulated noise (Fig. 2.22), deconvolve it (Fig. 2.23), and58

(a) Injected map(b) H1L1(c) L1V1(d) H1V1

Figure 2.19: The toy model of a localized source is shown in (a). Dirty maps
made from simulated data from three LIGO-Virgo baselines are shown in thelast three panels.

obtain its MLR (Table 2.3). The similarity of the values of this statistic for the
dirty and clean maps proves the unitarity of our deconvolution method.

One

can see that the MLR statistic (Table 2.3) is small (≈ 1) in all these cases. We
then introduce a small signal – the same as the extended (galaxylike) source
considered before, but at a much reduced strength.
Visually, the dirty maps are now weaker (Fig. 2.24), and clean maps almost
do not show the obvious presence of any source (Fig. 2.25), but the MLR statistic
59

(a) H1L1 (NMSE=0.6642)(b) L1V1 (NMSE=1.2939)(c) H1V1 (NMSE=1.4750)(d) H1L1V1 (NMSE=0.4961)

Figure 2.20: Clean maps obtained by the deconvolution of the dirty maps ofFig. 2.19, using 20 CG iterations, are shown here.

(Table 2.4) provides a clear and reliable indication of the presence of a signal,thus proving its usefulness in the search for signal in a noisy map.

2.5

Discussion

The search for an anisotropic stochastic gravitational wave background plays animportant role in present GW research. In addition to setting interesting upper
limits on astrophysical and cosmological backgrounds, the simplicity of the con60

(a) H1L1(b) L1V1(c) H1V1(d) H1L1V1Figure 2.21: Difference between the clean maps of Fig. 2.20 and the injected mapof Fig. 2.19(a).

comitant analysis reveals invaluable knowledge about the coherent performanceof the GW detector network.
So far, detailed analysis strategies have been developed to search for anisotropic background in pixel and spherical harmonic spaces, and a general maximumlikelihood-based framework has been established to search in any convenientbasis. The spherical harmonic search has been demonstrated using a networkof detectors [14]. In this paper, for the first time, we numerically implement
the directed radiometer search, including deconvolution, for a network of detec61

Table 2.1: MLR statistic of dirty maps (λ) versus clean maps (λc ) for the simulated maps in Fig.s 2.16 and 2.17.
Baseline

λ

λc

H1L1

785.555 783.271

L1V1

359.004 358.940

H1V1

315.717 315.662

H1L1V1 919.594 917.600

Table 2.2: MLR statistic of dirty maps (λ) versus clean maps (λc ) for the simulated maps in Fig.s 2.19 and 2.20.
Baseline

λ

λc

H1L1

284.652 284.173

L1V1

39.30839.377

H1V1

64.12964.113

H1L1V1 294.419 293.961

tors. These methods, in the past, were focused primarily on showing that theobserved map is consistent with Gaussian noise or in estimating sky maps. The
latter required the inversion of the convolution equation, which itself assumed
the network of detectors to be nondegenerate. Neither of these methods may
work in the presence of excessive noise and weak signal. Most importantly, astatistically meaningful, all-sky combined statistic, in the form of an optimal“detection statistic,” was needed in order to make precise statements about thepresence or absence of a given background model in a map. Here, we proposed62

(a) H1L1(b) L1V1(c) H1V1

Figure 2.22: Dirty maps made for simulated noise, without any injected signals.

a MLR statistic, which yields a single number when computed on the dirty or theclean map and can be used as a detection statistic. By computing the MLR statistic for a couple of toy models of the background, we observe that the detection
statistic is much larger than the noise-only case, even in the presence of weaksignals that are barely visible in dirty or clean maps. We corroborated thesestatements with results obtained from extensive Monte Carlo simulations of adiffuse background of known shape in an ensemble of noise realizations. However, a more detailed study using signals from a variety of background models
is surely worth pursuing in order to determine how accurate the templates need63

(a) H1L1(b) L1V1(c) H1V1(d) H1L1V1Figure 2.23: Clean maps obtained by the deconvolution of the dirty maps ofFig. 2.22, using 20 CG iterations.

to be in order to extract meaningful information from weak backgrounds.
We also compared the performance of individual baselines and the wholenetwork for the directed radiometer search using different figures of merit. Evaluating the performance of a network of GW detectors in SGWB searches is relatively straightforward compared to other GW signal searches [46, 47, 48, 49, 50,51]. This exercise was useful in drawing insights about the characteristics of a
network that are particularly helpful in boosting its performance. Our overall observation, not surprisingly, is that the network improves performance in mainly
64

Table 2.3: MLR statistic of dirty maps (λ) versus clean maps (λc ) for simulated
noise (which actually has an extremely weak signal added) in Fig.s 2.22 and 2.23.
Baseline

λ

λc

H1L1

0.512

0.433

L1V1

-1.549 -1.542

H1V1

1.105

1.120

H1L1V1

0.208

0.149

Table 2.4: MLR statistic of dirty maps (λ) versus clean maps (λc ) for the simulated maps in Fig.s 2.24 and 2.25.
Baseline

λ

λc

H1L1

98.64398.670

L1V1

43.52043.536

H1V1

40.43240.475

H1L1V1 115.132 115.176

three ways, namely, (1) by increasing the sensitivity by observing each direction
a greater number of times, (2) by observing the sky more uniformly, and (3) byprobing each direction on the sky with additional detectors on the globe. The
latter two enhancements lead to better localization of pointlike sources. This canbe understood via the behavior of the Fisher information matrix: More detectors
reduce its degeneracy and improve the well-posedness of the inverse problem.
This, in turn, leads to a more accurate production of clean maps.
Another question worth addressing in the future is about how closely spaced
65

(a) Injected map(b) H1L1(c) L1V1(d) H1V1

Figure 2.24: The toy model of a very weak extended source is shown in (a). Dirty
maps made from simulated data from three LIGO-Virgo baselines are shown in
the last three panels.

must the templates be on the parameter space to maximize the chances of detection with available computational resources. Indeed, the proposal for templated
searches for SGWB signals is not new to this paper. For example, it has beenaddressed earlier in the context of isotropic searches (see Ref. [53] and the references therein). Ref. [53] also introduced a metric on the parameter space of
those signals so as to enable an experimenter to infer what the principle axes areon that space and how fine a template bank one can afford based on the com66

(a) H1L1 (NMSE=1.0694)(b) L1V1 (NMSE=3.3557)(c) H1V1 (NMSE=3.8406)(d) H1L1V1 (NMSE=1.3503)

Figure 2.25: Clean maps obtained by the deconvolution of the dirty maps ofFig. 2.24, using 20 CG iterations.

putational resources available. A similar study can be carried out for finding amore optimal spacing of templates for directed searches than the one used here.
Whereas results presented here were derived for Gaussian noise, the codesused can be applied to real data as well. Indeed, the performance of the proposed statistic in real data sets from the LIGO and Virgo detectors can be determined through hardware injections that were done in the recent science runs,
such as the ones described in Ref. [54], and supplementing them with multiple software injections to improve the statistics. The expected improvement67

of network sensitivity over individual baselines, as demonstrated here, merits
the investment required for extending the current single-baseline analysis efforts [24, 25] to a multibaseline one. This conclusion is strengthened by the factthat adding a detector to a baseline can potentially mitigate the contribution ofcross correlated environmental noise that affects only one of the three resulting baselines. Including V1, which is on a different continental plate than theH1L1 baseline, can serve this purpose. Employing a null-stream statistic [52, 51]
to complement the detection statistic might also help in discriminating against
such noise.

68

Chapter 3
Searching for a SGWB from the
population of neutron stars in the
Virgo cluster
One of the promising applications of the methods presented in the previous
chapter is to search for a stochastic background from the population of rotating
non-axisymmetric neutron stars (NSs) in the Virgo cluster. The Virgo cluster isthe nearest and best-studied rich cluster of galaxies whose center is at a distanceof about 16.5 Mpc away in the constellation Virgo. At present, more than 1300
member galaxies of the Virgo cluster are known. Estimates suggest that there
are approximately 1011 neutron stars in the Virgo cluster [55], of which a small
fraction would emit GWs with frequencies within the sensitivity band of earthbased detectors. As such, it offers the best hope to detect a SGWB produced by
the population of the neutron stars in the Virgo cluster, if such a background69

exists and is reachable by the current and advanced detectors.
A recent study [55] suggests that the stochastic nature of gravitational waves
from the population of rotating non-axisymmetric neutron stars/pulsars in the
Virgo cluster can reveal itself as a GW hot spot in our searches. We begin thischapter by reviewing a series of assumptions made in this reference. We thendetail the steps that are taken towards a meaningful model-based search for
SGWB from NSs in the Virgo cluster.

3.1

Motivation for the directed search

The GW strain due to a rotating neutron star can be expressed as [56]h = 4π 2 α

GeI 2
f ,
c4 r

(3.1)

where α(≤ 1) is the orientation factor, e is the (equatorial) ellipticity of theneutron star, I its principal moment of inertia, r is the distance to the source,
G the Newton’s gravitational constant, c the speed of light, and f is the GW
frequency. If the rotation axis is optimally oriented with respect to the detector,
then α becomes unity.
The total number of Galactic neutron stars is estimated to be about 108 −
109 [57]. This is obtained by assuming the mean neutron star birth rate of 10−2
yr−1 and the Galactic age of about 10 Gyr. What is more uncertain for the study
reported in Ref. [55] is the number of neutron stars whose rotation period is of
the order of milliseconds. From the survey of radio pulsars in our Galactic disk,
the population of millisecond pulsars is estimated to be at least 40000 [58]. Thisimplies a mean birth rate of 4 × 10−6 yr−1 . To date, roughly 2000 NSs have been70

observed, not all of which are millisecond pulsars and the majority of which are
identified as isolated radio pulsars (PSRs) with ages . 100 Myr.
In Fig. 3.1, we show the distribution of observed radio pulsars. This distribution indicates two different populations: Each population is Gaussian distributed and has mean and standard deviation different from the other’s. Hence,the distribution is approximated as a bimodal Gaussian with the components asfollows:
P1 (log f r )d(log f r ) = √1e−((log fr −log µ1 ) /2σ1 ) d(log f r ) ,
2

2πσ1

2

(3.2)

(for f r > 50Hz) ,P2 (log f r )d(log f r ) = √1e−((log fr −log µ2 ) /2σ2 ) d(log f r ) ,
2

2πσ2

2

(3.3)

(for f r < 50Hz) ,where µ1 = 219 Hz, σ1 = 0.238, µ2 = 1.71 Hz and σ2 = 0.42, and f r = f /2. Here,P1 and P2 are normalized to unity when integrated from f r = 0 to infinity. We
assume a similar bimodal form of distribution of NSs in the Virgo cluster. We
assume that the total number of neutron stars in our Galaxy is 108 for f r < 50
Hz, and 40000 for f r > 50 Hz. Since there are of the order of 103 galaxies in the
Virgo cluster, the total number of NSs in the Virgo cluster is Nlow ∼ 1011 forf r < 50 Hz, Nhigh ∼ 4 × 107 for f r > 50 Hz. So, the distribution of NSs in the
Virgo cluster can be expressed asN ( f )d f = Nhigh P1 (log f r ) + Nlow P2 (log f r )



d fr
.
f r ln 10(3.4)

For a bandwidth of ∼ 103 Hz, the number of NSs in each mHz frequency bin
is ∼ 10 for f > 100 Hz. In the low frequency regime this number is much larger.
71

180
160

Number of pulsars

140
120
100
80
60
40
20
0

0.1

0.32

1
3.16
10
31.62
Pulsar spin frequency [Hz]100316.23 1000

Figure 3.1: The distribution of observed radio pulsars [59]. The horizontal axis
is log( f r ) where f r is the rotational frequency of pulsars. The red solid line is
the two component Gaussian-fit to the distribution [55].

Thus, it is not possible to resolve the signal from each NS and, therefore, the GW
signal from the NS population in the Virgo cluster is expected to be stochastic innature.
The spectral density of gravitational radiation from NSs in the Virgo cluster
is given by
H ( f ) = h h2 i N ( f ) ,
(3.5)



2
 e 
I
hα2 i f 4 N ( f ) , (3.6)= 7.05 × 10−34−5
10
1.1 × 1045 g cm2

where the unit of H ( f ) is Hz−1 , hα2 i represents the average of the orientation

factor with respect to the inclination and polarization angles. Here, we used
r = 16.5 Mpc, the distance to the Virgo cluster. Assuming a uniform distributionof the sources over those angles, we have hα2 i = 0.4. Figure 3.2 shows the
72

−48

10

ε =10−5
ε =10−6

S5 radiometer upper−limit

ε =10−7−50

10H(f) [strain 2 /Hz]

ε =10−8
−52

10

−54

10

−56

10

−58

10

1

2

10

10
Frequency [Hz]3

10

Figure 3.2: The spectral density of gravitational radiation from NSs in the Virgocluster is plotted as a function of GW frequency for different choices of average
ellipticity, e ∈ {10−5 , 10−6 , 10−7 , 10−8 }. The shaded area labeled “S5 radiometer
upper-limit” illustrates the radiometer upper-limit of GW strain power using S5
LIGO data [60].

spectral density as a function of GW frequency for different choices of average
ellipticity.
In order to obtain a rough idea of how large the H ( f ) is compared with
the noise PSD, a figure of merit “effective source spectral density”, Heff ( f ), is
definedHeff ( f ) =q

8 Tobs hΓ2 i1day f H ( f ) ,such that the signal-to-noise ratio can be expressed asv
uZu fu d f
H2 ( f )teffρ=,f ξ (1) ( f ) ξ (2) ( f )
fl
73

(3.7)

(3.8)

where Tobs is the observation time and hΓ2 i1day is the average of the squareof Γ (see Eq. (2.23)) over a sidereal day. The noise amplitude spectral densities(ASDs) of various detectors as well as the effective source ASD are plotted inqFig. 3.3. Here, we assume Tobs = 1 yr and hΓ2 i1day = 0.2.
LIGO
VIRGO
aLIGO
aVIRGOET
KAGRA
H1/2 (f)(ε =10−5)

−21

10

−22

√
h(f) [1/ Hz]

10

eff
1/2

−6

Heff (f)(ε =10 )1/2−7Heff (f)(ε =10 )1/2
Heff (f)(ε =10−8)

−23

10

−24

10

−25

10

−26

10

0

10

1

10

2

10
Frequency [Hz]

3

10

4

10

Figure 3.3: The noise amplitude spectral densities of various first, second andpthird generation detectors, as well as the effective source ASD, Heff , are plotted here. In this plot, Tobs = 1 yr and hΓ2 i1day = 0.04. This figure is reproduced
from Ref. [55].

The analysis results shown here strongly depend on the value of average ellipticity of the NSs. If e ≤ 10−6 , it will be difficult to observe stochastic GW
profiles of the NSs in the Virgo cluster with aLIGO, aVIRGO, and KAGRA detectors. Only when detectors have ET like sensitivity will it be possible to detectsuch a signal for average e = 10−6 . The detectability improves if we employ
several baselines of a network of detectors and a longer observation time [13].
74

3.2

Modeling the neutron stars in the Virgo cluster

Radio pulsars are highly magnetized, rapidly spinning neutron stars which emit
beams of radio waves. The wave is observed only when the beam crosses theline of sight. It is important to note that the observed distribution of the periodof radio pulsars is not necessarily the true representation of the entire population [61]. This is because various selection effects cause the number of detectedpulsars to be much smaller. For example, we can not observe the pulsars whoseradio beams do not cross the earth or whose radio emission is too weak to be
detected. The aperture of the emission cone decreases with the period as well as
the luminosity.
The optimal way of searching for stochastic GW profiles from the neutronstars in the Virgo cluster is to cross-correlate the data from the detectors, taken
in pairs, with a sky-position dependent time-frequency filter. To do that we needto assume what the spectral shape is of the targeted source. If the true spectrum
does not match with the modeled spectrum, the statistic becomes sub-optimal.
However, the source spectrum derived in Eq. (3.6) is severely affected by theradio selection effects. To model such a spectrum a set of corrections needs tobe applied to the known pulsar sample [62]. Here we propose to simulate theexpected spin frequency distribution assuming a power law, and then analyze
the data with a few different spectral index choices. If no detection is made, wewill present the upper limits on the average ellipticity of neutron stars in the
Virgo cluster as a function of the spectral index. In the following sections, wejustify three spectral index choices for the search.

75

There are several mechanisms by which neutron stars generate gravitational
waves. These include free precession, spinning non-axisymmetric deformation,and unstable oscillation. Free precession occurs if a neutron star’s rotation axis
is misaligned with respect to its symmetry axis. There are a number of mechanisms that may lead to a neutron star being deformed away from axisymmetry.
The non-axisymmetric rotating neutron stars are the most promising sources forcontinuous GWs. Unstable oscillation modes in the fluid part of a neutron starcan occur during the birth of the star, close encounter with another star or an
accretion phase.

3.2.1

The gravitar model

We consider a population of non-axisymmetric rotating neutron stars having
negligible magnetic field so that the dominant energy loss, due to a deviationfrom axisymmetry, goes into the production of gravitational waves. These objects
are often called “gravitars” [63, 64, 65]. There might exist such a population,
which is electromagnetically invisible. Detecting its GW profile with groundbased GW detectors is an interesting prospect. Even if gravitars do not exist,
they provide a relevant upper bound on the GW emission by objects that doexist, such as rapidly-rotating NSs. The GW spin-down due to a gravitar’s nonaxisymmetric shape is given by
32π GI 2 5
f˙ = −
e f ,
5 c5
4

(3.9)

where e = I1 −I I2 is the ellipticity of the gravitar, with I1 , I2 and I being the moment of inertia with respect to the three principal axes. Note that the gravitational76

wave frequency is twice the NS’s spin frequency. Integrating Eq. (3.9) gives theGW frequency at time t as
f (t) =



f 0−4 + β−1 e2 t

−1/4

,

(3.10)

where f 0 is the initial (birth) GW frequency and
β=

5 c5
.
128π 4 GI(3.11)

The constant β approximates to
β

1/3



= 5.3

1038 kg m2
I

1/3
kHz .

(3.12)

The frequency evolution of gravitars is determined by the distribution of spin
frequencies at birth, ellipticity, and the age since birth.
We considered a model motivated from the observed pulsars for the distribution of the initial spin frequencies. It is the one found as the maximum likelihood
model by population synthesis [66, 63]. It is described by a standard log-normal
distribution2
1
P f 0 ( f 0 )d f 0 = √d f0 ,
exp − 2 ln( f 0 ) − ln( f¯0 )2σ
2πσ f 0

1

(3.13)

where the unit of f 0 is Hz and the values σ = 0.69 and f¯0 = 400 Hz are takenfrom the Ref. [66]. Gravitars with f 0 > 4 kHz are excluded from the simulation.
The value of the ellipticity depends on the neutron star properties. The typical ellipticity of neutron stars is unknown, and the possible maximum value isvery uncertain. At present, the values 10−9 < e < 10−6 are considered plausible.
Following [63], we use a probability distribution which satisfies the principle ofmaximum entropy and is given by
Pe (e) = a(τ )e−e/τ ,77(3.14)

wherea(τ ) =1
τ (1 − e−emax /τ ),(3.15)

with emax being the maximum possible value for e. The relation between theparameter τ and the mean value of the distribution, ē, is given byτ = ē +emax.
exp(emax /τ ) − 1(3.16)

It is straightforward to derive an evolved distribution in frequency, given an
initial distribution of frequency, birth rate and ellipticity [65]. We define theconditional frequency distribution p f ( f |e, t) at time t such thatdPf = p f ( f |e, t) d f(3.17)

is the probability of finding a gravitar in the frequency range f to f + d f , assuming an ellipticity e and an age t. Let dPf0 be the probability of the birthfrequency in the frequency range f 0 to f 0 + d f 0 . The corresponding probability
density p f0 ( f 0 ) is defined by
dPf0 = p f0 ( f 0 ) d f 0 .

(3.18)

A frequency change by redshift is neglected since all the gravitars consideredhere are within the local universe (a few Mpc). To compute p f ( f |e, t) fromthe initial p f0 ( f 0 ), consider a single gravitar with ellipticity e whose current
frequency is f , and let f 0 (e, f , t) denote the gravitar’s frequency at time t in thepast. Solving Eq. (3.10) for the birth frequency yields
f 0 ( f , e, t) =

f −4 − β −1 e 2 t

78

−1/4

.

(3.19)

The fraction of gravitars at birth in the frequency range f 0 to f 0 + d f 0 is the sameas the fraction at present in the frequency range f to f + d f , so the intensity
p f d f = p f0 d f 0 yields∂ f 0 ( f , e, t)
df ,∂f
f 5 ( f , e, t)= p f0 ( f 0 ( f , e, t)) 0 5
df .
f

p f ( f |e, t) d f = p f0 ( f 0 ( f , e, t))(3.20)(3.21)The probability dP of finding a gravitar in a volume characterized by d f de dt is
dP = p f ( f |e, t) pt (t) pe (e) d f de dt ,(3.22)

where pt (t) and pe (e) are the probability distributions for the time and ellipticity, respectively. Note that the distributions of the initial frequency, time andR∞
ellipticity are assumed to be normalized, such that 0 p x ( x ) dx = 1. Then the
evolved frequency distribution is obtained by marginalizing dP over t and e aspev ( f ) d f =

Z T

f

0

dt

Z ∞
0

de p f ( f |e, t) pt (t) pe (e) d f .

(3.23)

Some of these expressions, between Eqs. 3.14 and 3.23, are re-derived from Benjamin Knispel’s Ph.D. thesis [65]. For a fixed age of all sources, which is obtainedby marginalizing Eq. (3.22) over ellipticity only, the evolved frequency distribution can be expressed aspev
f ( f |t) d f =Z ∞0

de p f ( f |e, t) pe (e) d f ,(3.24)
2

Z ∞ −1.05(− 14 ln( f −4 − β−1 e2 t)−5.992)
e
a(τ )
=
e−e/τ de (3.25)
51.73 f 0
f −4 − β −1 e 2 ta(τ ) τ −1.05(ln( f )−5.992)2=e,(3.26)
1.73 fwhere the last expression is obtained by setting t = 0 Myr.
79

3.2.2

A population synthesis model of millisecond pulsars

The population of the millisecond pulsars (MSPs) was the primary focus of thestudy reported in Sec. 3.1. Since MSPs are observed to spin down relatively
slowly, it is unrealistic to consider high ellipticity (e ≈ 10−5 ) for such a population. Figure 3.4 shows the distribution of observed MSPs and it certainly does
not represent the entire population since only a small fraction of the estimatedGalactic MSPs are observed. Also, the Virgo cluster is an elliptical rich cluster. These facts should be taken into consideration when deriving the estimatednumber of MSPs in the Virgo cluster. To take these into account, we begin witha population synthesis of radio pulsars in the Galactic disk [67, 68].
Bin width = 1 Hz.
4
3.5

Number of pulsars

3
2.5
2
1.5
1
0.5
0

100

200

300
400
500
Pulsar spin frequency [Hz]

600

700

Figure 3.4: The distribution of observed millisecond pulsars [59]. The frequency
bin width considered here is 1 Hz.

In a recent study, S. A. Story et al. used a population synthesis to investigate
the initial properties of MSP [67]. This study does not account for the populationof MSP observed in globular clusters. The spatial distribution of millisecond
80pulsars is modeled by assuming their birth in the Galactic disk with a random
kick velocity and by evolving them to the present within the Galactic potential.
Following [67], we derive the distribution of the actual period, by evolving a
set of simulated MSPs from the end of the spin-up phase to present 1 . For eachpulsar, we proceed as follows:• We assume a uniform birth rate distributed between 0 and 12 Gyr.
• The initial magnetic field B0 is drawn from a uniform distribution in therange from 108 to 1012 Gauss. In this model, there is no dissipation of themagnetic field.
• The initial period is derived from the formulaP0 = 0.18 × 103δ/7 B06/7 ,(3.27)

were P0 is in ms and B0 is in units of 108 G. The dithering parameter δvaries from 0 to 2.8.
• Finally, having the initial period, magnetic field and age, we calculate theactual period from the usual dipole formula [71]:P(t) =q

P02 + 0.154 B02 t ,(3.28)

where P and P0 are in ms, B0 in units of 108 G, and t in Gyr.
The Monte Carlo simulations for 10000 MSPs with P < 50 ms are done by
Tania Regimbau [72]. One could improve the model considered here by allowing
1 Similar studies for normal pulsars can be found in Refs. [69, 70].

81

magnetic field decay and simulating them in an elliptical rich galaxy. However,
if we account that one third of the total luminosity of the Virgo cluster comesfrom spiral type galaxies and two third from ellipticals, then the total numberof MSPs in the Virgo cluster is about 275 times larger than that in the Milky
Way [73].

3.2.3

A simple power-law

Here we consider a crude approximation that there are equal numbers of NSs inevery frequency bin, i.e, N ( f ) = constant.
In Fig. 3.5, we show different SGWB spectra calculated using the models
described above. For comparison, the SGWB spectrum obtained from observed
radio pulsars (without correcting for radio selection effects) is also plotted here.
Note that the quantity H̄ ( f ) is defined such that H̄ ( f ) = H ( f )/H ( f = 100). For
the gravitar model we assumed the average ellipticity to be ē = 10−7 , maximum
ellipticity emax = 2.5 × 10−6 , and an age in the range [0, 10] Myr. Spectra withdifferent ages in the range [0, 10] Myr are not much different. This is because we
assume a small average ellipticity for the gravitars. It is interesting to note thatthese model-based spectra are not much different from each other. But each ofthese three spectra is different from the spectrum of the observed population.

3.3

The directed SGWB search and Olbers’ paradox

Following Olbers’ paradox, for a putative SGWB from NSs in the Virgo cluster,
one may not need to invoke a directional search. An all-sky search should be
82

6

10

4

10

Observed radio pulsars
Simple power law
Gravitar modelPop. synthesis of MSP

2

H̄(f)

10

0

10

−2

10

−4

10

2

3

10

10
f [Hz]

Figure 3.5: Several model-dependent spectral profiles of SGWB from NSs in the
Virgo cluster. Note that H̄ ( f ) = H ( f )/H ( f = 100 Hz).

enough to add up coherent cross-power from all such sources thought the universe, out to a redshift beyond which stars were not yet forming. To address thisconcern, let us begin by stating the Olbers’ paradox: Light from the stars shouldmake the sky as bright as the sun in all directions. The darkness of the night
sky conflicts with the assumption of an infinite and static universe. In otherwords, the number of stars per unit solid angle increases with distance-squared.
Whereas, their intensity at the earth decreases as inverse distance-squared. Thismeans that the intensity per unit area per unit solid angle should be the same inevery direction.
There are many solutions to this paradox depending on the models about theuniverse’s structure and/or its evolution. For our case, the solution to Olbers’
paradox becomes obvious when we include the expansion of the universe. The

83

gravitational waves themselves are redshifted, which introduces the additional
decrease in flux and makes the distant universe darker. One can compare therelative strengths of an anisotropic SGWB and the isotropic SGWB created respectively by the local and distant universe. By counting the number of Galaxies
per unit solid angle per unit distance, one finds that the population in the direction to the Virgo cluster is roughly 8 times denser than that from all-sky up to 4
Gpc [74]. However, a better study considering the expansion of the universe isdone in Ref. [75]. For a universe which is currently dust dominated, it is foundthat the Virgo cluster is about 50 times brighter than the background created by
the distant sources (see Appendix C.1 for details). Hence, it makes sense to do a
directed search for GW from the Virgo cluster.

3.4

Setting up the analysis

The choice of basis for an analysis, in principle, should not affect the physicalsearch results. However, in practice, such a choice can bear on computationalcosts of a search and also on the systematic errors affecting observational results.
For these reasons, we expect that while searching for GW from a point source,
pixel basis is the natural choice. Whereas searching for a diffuse background, aspherical harmonic basis is a better choice.

3.4.1

Methodology

Since the angular extension of the Virgo cluster is about 8 degrees, we carryout the data analysis using a spherical harmonic (SpH) basis [14]. An extended84source with an arbitrary angular distribution can be characterized by sphericalharmonic coefficients Plm such that

P (Ω̂) = ∑ Plm Ylm (Ω̂) ,(3.29)

lm

which is consistent with Eq. (2.34). The series is truncated at l = lmax , allowing
for angular scale ∼ 2π/lmax . In SpH basis, the MLR statistic introduced in
Eq. (2.20) can be written as
λ= pP µ† XµP ν† Γντ P τ,(3.30)

where µ ≡ {l, m}, Xν is the dirty map, and Γµν is the beam matrix or Fisher
matrix and are given by [14]
Xν =
Γµν =

∑ γν∗ ( f , t) ξ

(3.31)ft

H( f )C ( f , t) ,(1) ( f , t ) ξ (2) ( f , t )∑ γµ∗ ( f , t) ξ

H2 ( f )γν ( f , t ) ,(1) ( f , t ) ξ (2) ( f , t )(3.32)

ft

where γµ ( f , t) is the angular decomposition of the overlap reduction functionγ(Ω̂, f , t), which characterizes the orientation and frequency response of thedetectors:
γlm ( f , t) ≡Zγ(Ω̂, f , t) =1 A
F (Ω̂, t) F2,A (Ω̂, t)e2πi f Ω̂·∆r(t)/c ,2 1

S2dΩ̂ γ(Ω̂, f , t) Ylm (Ω̂) ,(3.33)(3.34)

2 ∗and C ( f , t) = ∆tx̃1 ( f , t) x̃2 ( f , t) is the cross spectral density generated for eachinterferometer pair.

85

3.4.2

Angular extension of the Virgo cluster

We model the signal-strength vector of the Virgo cluster such that its non-zerocomponents follow a Gaussian distribution, centered at 12h 26m 32s RA and00+12◦ 430 24 Dec, extended over 12h−13h RA and 5◦ − 20◦ Dec. Figure 3.6 shows
the map of the signal-strength vector P ≡ Pmodel . The spherical-harmoniccoefficients Plm are set to zero for l > (lmax = 30) as shown in Fig. 3.7. This
corresponds to an angular scale about 12 degrees.

Figure 3.6: The map of a modeled signal-strength vector of the Virgo cluster. In
this map, the spherical harmonic coefficients are set to zero for l > (lmax = 30).

3.4.3

The analysis pipeline

We analyze data from LIGO’s 4 and 2 km detectors (H1, H2 & L1, excluding H1H2 pair) and VIRGO’s 3 km detector during their S5 and S6 and VSR1, VSR2,
and VSR3 science runs. The HLV network lifetime is about 3.2 years. The search
bandwidth considered here is 40-1500 Hz. We parse the time series into 60
second intervals, Hann-windowed, 50%-overlapping segments, coarse-grained
86

−48

3

x 10

−48

2

x 10

1.5

2

1

1
ℑ(Plm )ℜ(Plm )

0.5

0

0

−0.5

−1

−1

−2
−1.5

−3
0

5

10

15
l

20

25

−2
0

30

5

10

15
l

20

25

30

−20

−10

0
m

10

20

30

−48

3

x 10

−48

2

x 10

1.5

2

1

1
ℑ(Plm )ℜ(Plm )

0.5

0

0

−0.5

−1

−1

−2
−1.5

−3
−30

−20

−10

0
m

10

20

30

−2
−30

Figure 3.7: The spherical harmonic coefficients of the modeled signal-strength
vector in Fig. 3.6. In the top panel, the real (left) and imaginary (right) parts of

Plm are plotted for different m as a function of l. Whereas in the bottom panel,
the real (left) and imaginary (right) parts of Plm are plotted for different l as afunction of m.

to achieve 0.25 Hz resolution. We also mask frequency bins associated with
instrumental lines. The cross-correlation analysis is performed on each segment.
The results from the segments are then combined into the final result. Figure 3.8
shows a schematic of the stochastic search pipeline used for this search. For eachbaseline, given a model P , the MLR statistic is computed based on the dirty
map and beam matrix. The MLR statistics from the different baselines are then87

combined to find a network MLR statistic. The whole process is repeated butwith time-shifted data to estimate the background. The results are collected at
the end and post-processed if needed. Note that the upper stream of the pipeline
is common to both all-sky and directed searches. This pipeline has been tested
with hardware and software injections, and reviewed on many occasions.

Detector 1
Data­findDetector 2
Data­find

Condition data,
FFT

Baseline 1
Beam,
Beam
Dirty­map

Detector 3
Data­find

Condition data, Condition data,
FFT
FFT

Baseline 2
Beam,
Dirty­map

Baseline 3
Beam,
Dirty­map

Template­Bank

Detection
statistic

Detection
statistic

Detection
statistic

Network
Detection statistic,
Post­processing

Figure 3.8: A schematic diagram of the stochastic search pipeline.

3.4.4

Simulations

It is simple, using the spherical harmonics basis, to specify a resolution cutoff byonly allowing l ≤ lmax . This avoids oversampling and reduces the number of88required basis vectors. Also, since extending this cutoff to a larger lmax does notaffect the original basis vectors, it is straightforward to run the analysis with ahigher resolution, and later follow up at a lower resolution. The computationallydominant part of the analysis is the calculation of the Fisher matrix. The Fisher
matrix has (l + 1)4 elements, with (l + 1)2 the number of basis vectors. The symmetries in the spherical harmonic basis help to reduce the computational load.
However, reading, writing and saving the Fisher matrices and the dirty maps
are memory and disk space limited. The larger the l, the higher the memoryand disk space requirements.
35
30
25

λ

20
15
10
5
0
−5
0

5

10

15
lmax

20

25

30

Figure 3.9: The MLR statistic of a simulated SGWB in the Virgo cluster is computed and plotted for different harmonic cutoff.

Here we do a Monte Carlo simulation to choose the harmonic cutoff. Following the diffraction formula, the larger the baseline, the higher the resolution.
However, we consider 15 days of data from S5 science run with coincident between H1 and L1. We inject a signal, following the above Pmodel and simple
89

power-law spectrum, into the data stream. The data were first processed withlmax = 30, and later with different lower lmax . The MLR statistic for different
choices of lmax cutoff is presented in Fig. 3.9. It is obvious form the plot that
λ saturates at lmax = 25. So, we adopt lmax = 25 to search for the SGWB
from NSs in the Virgo cluster. It is also found in a similar study that the known
sky-distribution of the source contributes about 10% more to the signal-to-noise
ratio than the unknown counterpart.

3.5

Discussion

Since the search is ongoing and it requires internal review, we do not report thedata analysis results here. The actual population of NSs in the Virgo cluster is
unknown. Searching for their SGWB profile with current and advanced detectors has a merit if the actual number exceeds our expectation. The analysis itself
is nevertheless interesting, because even if we do not detect the signal, we canput an upper limit on the average ellipticity of the NSs in the Virgo cluster. In
addition to searching for SGWB targeting the Virgo cluster, we perform an allsky search using the same spectra. This will allow us to compare the results fromtargeting the Virgo cluster with those from an all-sky search, and hence address
Olbers’ paradox. Extending this analysis to other targets (e.g., NSs in the Comacluster) and searching for them in advanced detectors’ data are straightforward.

90

Chapter 4
Searching for perturbed black hole
ringdown signals
In this chapter, we study the methods of detecting the ringdown phase of a GW
signal that arises from a perturbed black hole. Such a black hole can result fromthe coalescence of a compact binary, a massive object falling into a BH or theformation of a BH through the asymmetric core collapse of a massive star. Thissignal is initially in the form of a superposition of quasi-normal modes [122].
However, at late times the waveform, which is known as a ringdown, is expected
to be dominated by a single mode. The existence of intermediate-mass black
holes (IMBHs), which have masses from hundreds to tens of thousands of solar
masses, has not yet been corroborated observationally. However, these objects
are of great interest for astrophysics and their discovery by LIGO and VIRGO in
the GW spectrum is an interesting possibility.
The optimal method for searching such a signal buried in detector noise is91

to match-filter the detectors’ output with theoretically modeled waveforms. The
coherent network statistic is optimal for detecting these signals in stationary,Gaussian noise [46, 47]. But in real noise, which is non-Gaussian and nonstationary, that statistic is sub-optimal unless additional discriminators of noise
artifacts are employed [12]. Here, we describe a hierarchical method for coherently searching ringdown signals in a network of detectors that is aided by suchdiscriminators.

4.1

Motivation for the search

A black hole itself does not emit any signals from within its horizon classically,
thus making it invisible and difficult to detect directly. However, astronomers
can infer the presence of a black hole by electromagnetically observing the influence of its strong gravity on nearby matter. In recent years, evidence has
mounted for both stellar-mass black holes and supermassive black holes. Stellarmass BHs are born in the collapse of massive stars and have masses in therange 3 M

≤ M ≤ 20 M . From the discovery of a number of bright X-

ray sources, it is evident that there are many binaries in our Galaxy and othergalaxies that consist of a stellar-mass BH, which accretes matter from a stellar companion [76, 77, 78, 79]. Supermassive BHs, which are believed to be
the engines behind quasars, have masses of 106 M

to 1010 M . Observa-

tions of the stellar environment near galactic centers, especially the high luminosity and small size of active galactic nuclei, provide the strongest empiricalevidence for the existence of massive black holes in the center of most galax-

92

ies [80, 81, 82, 83, 84, 85, 86, 87, 88, 89]. The formation of supermassive BHs is
not as well established as the formation of stellar-mass BHs, but it is theorizedthat they are produced during collisions between galaxies, or in a series of collisions between BHs that have masses greater than those of stellar-mass BHs, but
significantly less than those of supermassive BHs.
Observations of ultraluminous X-ray sources in galaxies and presence of excess dark mass in globular-cluster cores suggest that there may exist a population of intermediate-mass black holes with masses M ∼ 102 − 104 M [90]. But
the existence of IMBHs has never been confirmed [91]. Several IMBH formation
mechanisms have been proposed: (a) Formation of IMBHs in young clusters
via runaway collapse of stars [92, 93, 94]. (b) Formation of IMBHs in globular clusters through repeated mergers of stellar-mass BHs in binaries [95, 96].
(c) Formation of lower mass single IMBHs by the stalled supernova of early
Population III stars [90, 97]. (d) Formation of IMBHs from progressive accumulation of mass into a large seed black hole via coalescence of a population ofsmaller-mass BHs [95]. (e) IMBHs could be the relics of very massive metal-freestars [98]. The existence of binaries with IMBH components remains uncertain
since stellar winds may stall the growth of the IMBH progenitors in the runaway
collision scenario [99], or the merger recoil may also eject a newly formed blackhole from the globular cluster [100, 101].
Some recent observations indicate that IMBHs could exist in the core of globular clusters [102, 103, 104, 105]. These objects are thought to grow from accretion of smaller compact objects [106]. It has been suggested that IMBHs are theengines powering ultraluminous X-ray sources, such as M82 X-1 [107, 108, 109]

93

and NGC 1313 X-2 [110, 111]. Ultraluminous X-ray sources are found off-centerin young star clusters in starburst galaxies: X41.4 + 60 in M82 has approximate
mass range of 500 − 105 M[112, 113] and HLX-1 in ESO 243-49 has a lower

mass limit of approximately 500 M [114].
The detection of GWs from IMBHs could have important consequences fortheories about the formation of supermassive black holes and the dynamics andevolution of globular clusters [115, 116] and could ascertain information about
their masses and spin parameters. The merger and ringdown phases of theGW signal are important for detection of IMBH sources because for massive
systems the characteristic frequencies of the inspiral phase are usually outside a
detector’s frequency band.

4.2

Quasi-normal modes of black holes

According to GR, a black hole (BH) is a region of space-time where gravity is sostrong that nothing, not even light, can escape it. At the center of a black hole,space-time has a curvature singularity. That singularity is enclosed by an “event
horizon”, a mathematically defined surface that marks the point of no returnfor an object falling in. Black holes can have masses varying from a few solar
masses to millions of solar masses. In general, a black hole can be completely
characterized by three parameters, namely, its mass, spin and charge. And only
the mass and spin are likely to be significant for any astrophysical black hole.

94

4.2.1

Schwarzschild and Kerr black holes

The Kerr-Newman metric is a stationary and axisymmetric solution of the Einstein-Maxwell equations, describing the space-time geometry in the region surrounding a charged, rotating black hole [117]. The other stationary black holesolutions can be obtained as particular cases of the Kerr-Newman solutions. The
space-time geometry outside an uncharged non-spinning spherical black hole ofmass M is given by the Schwarzschild metric2GM
ds = − 1 − 2
c r
2





2GM
c dt + 1 − 2
c r
2

2

 −1

dr2 + r2 dθ 2 + r2 sin2 θdφ2 . (4.1)

And the space-time geometry outside an uncharged spinning black hole of mass
M and angular momentum J is described, in terms of the Boyer-Lindquist coordinates (r, θ, φ), by the Kerr metricds2Σ
4GM r a sin2 θdφ dt + dr2 + Σ dθ 2cΣ∆!
2GM r a2 sin2 θ+ r 2 + a2 +sin2 θdφ2 ,
(4.2)
c2 Σ


2GMr
= − 1− 2c Σc2 dt2 −where a ≡ J/cM is the spin parameter, Σ ≡ r2 + a2 cos2 θ and ∆ ≡ r2 −2GMr/c2 + a2 . From this point on, unless otherwise mentioned, we refer tothe spin by a dimensionless spin parameter â = Jc/GM2 . This parameter cantake values ranging from 0 (Schwarzschild BH) to 1 (extreme Kerr BH).

4.2.2

Black hole perturbation

The theory of black hole perturbations and the associated radiation has a long
history. The study of BH perturbations was initiated by the pioneering work of95

Regge and Wheeler [118] and was continued by Zerilli [119] and Chandrasekhar
and Detweiler [120, 121]. Their study found that black holes remain stable if subjected to a small perturbation. Their study also showed that the solution to thelinearized Einstein Equation can be expressed in terms of spherical harmonics
Ylm . The existence of quasi-normal modes (QNMs) of a BH was first pointed out
by Vishveshwara [122] in calculations of the scattering of gravitational waves by
a Schwarzschild black hole. In GR, quasi-normal modes arise as perturbationsof stellar or black hole space-time. QNM oscillations have been found in perturbation calculations of particles falling into Schwarzschild [123] and Kerr black
holes [124, 125] and in the collapse of a star to form a black hole [126, 127, 128].
A nice review of quasi-normal modes can be found in Ref. [129]. A BH canbe perturbed in a variety of ways, e.g., by the incidence of GWs, by an object
falling into it, by the interaction with a companion, by the accretion of mattersurrounding it, or by the formation process in a gravitational collapse. There
are no normal mode oscillations associated with non-radial perturbations. This
is in contrast with the normal modes of Newtonian gravity, because in GR they
are damped by the emission of GWs and, hence, they are called quasi-normal
modes. Each mode has a characteristic complex frequency, whose real part represents the actual frequency of the oscillation and whose imaginary part represents the inverse of the damping time. We express the angular frequency ωlm interms of the oscillation frequency f lm and the quality factor Qlm as−1ωlm = 2π f lm − iτlm
= 2π f lm − iwhere τ is the damping time.

96

π f lm
,
Qlm

(4.3)

4.2.3

The ringdown waveform

Numerical simulations (see for example, Refs. [130, 131]) have demonstrated
that the fundamental mode, l = m = 2, dominates the GW emission. Far fromthe source, the plus and cross polarizations of a ringdown waveform can be expressed in terms of the central frequency f 0 ≡ f 22 and the quality factor Q ≡ Q22
as [132, 133, 134, 135]
π f0 th+ (t) = A (1 + cos2 ι) e− Q cos(2π f 0 t) ,(4.4)

h× (t) = A (2 cos ι) e− Q sin(2π f 0 t) ,(4.5)

π f0 twhere A is the amplitude of the l = m = 2 mode, and ι is the inclination angleof the source. The strain produced in the detector is thenh(t) = h+ (t) F+ (θ, φ, ψ) + h× (t) F× (θ, φ, ψ) ,(4.6)

where F+,× are the detector antenna-pattern functions (see Appendix A.1 formore details), with ψ being the wave-polarization angle and (θ, φ) being thesky-position of the source. For a given source at a distance r, one defines theeffective distance as
Deff = qr
2 (1 + cos2 ι )2 /4 + F2 cos2 ι
F+
×

.

(4.7)

The effective distance gives the distance to a BH as if it were optimally locatedand oriented.
Following Refs. [136, 133, 135], the ringdown waveform is approximated by
π f0 th0 (t − t0 ) = A e− Q cos(2π f 0 t − ϕ0 ) ,
97

(4.8)

were ϕ0 is the signal’s initial phase and t0 is the time of arrival at the detector.
If e is the fraction of the black hole’s mass radiated as gravitational waves, thenthe amplitude is given by [134]
r


GM
5
A=
e
Q−1/2 F ( Q)−1/2 g( â)−1/2 ,
2
c2 r

(4.9)



7
where g( â) = 1.5251 − 1.1568(1 − â)0.1292 and F ( Q) = 1 + 24Q
2 , M is the blackhole mass, c is the speed of light, G is the gravitational constant, and r is the
distance to the source. An example of a ringdown waveform is shown in Fig. 4.1.
The frequency and quality factor for each quasi-normal mode can be related
−21

x 10

r = 100 Mpc [f0 = 172.8109 Hz, Q = 3.2871]
M = 100 M⊙ , â = 0.7

543h0 (t)

2
1
0
−1
−2
−3
0

0.01

0.02

0.03
t [s]

0.04

0.05

0.06

Figure 4.1: Ringdown waveform for a black hole of mass M = 100 M and spinâ = 0.7 corresponding to frequency f 0 = 172.8109 Hz and quality Q = 3.2871, at
a distance of r = 100 Mpc.
to the black hole mass M and dimensionless spin parameter â through fitting98formula [141, 142, 143, 144, 145]. For the l = m = 2 mode, we haveQ = 0.7000 + 1.4187(1 − â)−0.4990 ,i
1 c3 h0.1292,f0 =1.5251 − 1.1568(1 − â)
2π GM(4.10)(4.11)and the inverse relationship is given byâ = 1 −

Q − 0.7000
1.4187"

−2.0040

1 c3
1.5251 − 1.1568
M =2π G f 0,(4.12)
Q − 0.7000
1.4187

−0.2589 #.

(4.13)

Note that Q is determined by â and f 0 is determined by M and â. Figure 4.2
shows how frequency and quality factors are related to mass and spin.
4

10

16

â = 0â = 0.495â = 0.99

14
12

3

10

Q

f0 [Hz]

10
8
2

10

6
4
1

10

0

200

400
M [M⊙ ]

600

800

2
0

0.2

0.4

0.6

0.8

1

â

Figure 4.2: Central frequency of the ringdown waveform is plotted as a functionof mass for three spin values (left). The mass range (above a few solar masses)
reflects the sources that earth-based detectors are most sensitive to. Quality
factor is plotted as a function of the dimensionless spin parameter (right).

99

4.3

Ringdown searches with matched filtering

When the signal is known, the optimal method of extracting the signal fromGaussian noise is matched filtering [26]. We filter the strain data from a detector
with each template characterized by either mass and spin or frequency and quality factor. A trigger is generated when the signal-to-noise ratio (SNR) yieldedby matched filter crosses a pre-defined threshold for a given filter. Of course,
the threshold is chosen carefully so as to minimize the false dismissal rate for agiven false-alarm rate. The higher the threshold, the more significant the trigger.
However, as the noise in the data stream is non-stationary and non-Gaussian,matched filtering alone is not enough to establish the significance of a trigger asa GW signal. Noise can often mimic the signal we are searching for, and so alarge effort goes into characterizing the noise to best separate it from a potential
gravitational wave signal.

4.3.1

Matched filter

In GW data analysis, the data from multiple detectors is match-filtered with
templates derived from theoretical waveforms to test the presence or absence of
signals in the data. Filtering the data x (t) with an elliptically-polarized template
hep (t; µi ) characterized by the source parameters µi yields the SNR statistic givenby|h x, hep i||h x, hep i|
ρep (hep ) = q,≡
σep
hhep , hep i(4.14)

where h x, hep i denotes the noise-weighted scalar product of the data and thetemplate defined in Eq. (5.10). The data x (t) may or may not contain a signal
100

and is expressed as
x (t) = n(t) ,if signal is absent.

 n(t) + h(t) ,if signal is present.

(4.15)

Far from the source, the ringdown template can be expressed as
π f0 thep (t) = e− Q {cos(2π f 0 t) − i sin(2π f 0 t)} .

(4.16)

π f0 tThis choice is equivalent to modeling the template as e− Q cos(2π f 0 t − ϕ0 ) inmatched filter searches. For each template, triggers that have SNRs greater thana pre-defined threshold are retained. These triggers are used for determining
coincidence across different detectors.

4.3.2

Template banks

Although the form of the signal is known, the exact values of the intrinsic parameters are unknown. For the ringdown waveform, these parameters are thecentral frequency f 0 and the quality factor Q. We create a bank of templates tosearch over the range of intrinsic parameters of interest. We then filter the datawith each of the templates in the bank. The templates are placed discretely inthe parameter space. If the signal parameters match those of a template and the
SNR crosses the pre-defined threshold, we record the template parameter, thetime of the trigger and other relevant information. In general, more than one
template in the bank can generate a trigger for a given signal. We then obtain an
estimate of the intrinsic parameters for the source of the signal from the templatethat provided the largest SNR.

101

In order to understand how to place templates in the parameter space, we define two templates ũ( f ; µ, λ) and ũ( f ; µ + ∆µ, λ + ∆λ) characterized by intrinsic
parameters λ and extrinsic parameters µ, such as time of arrival and phase. The
match, M, between two templates is defined as the inner product maximized
over the extrinsic parameters:
M (λ, ∆λ) ≡ maxµ,∆µ hu(µ, λ), u(µ + ∆µ, λ + ∆λ)i .

(4.17)

This is the fraction of the SNR achieved by filtering a signal with parameters λagainst a template with parameters λ + ∆λ [43]. Expanding M in a power series
about ∆λ gives1
M (λ, ∆λ) ≈ 1 +
2



∂2 M
∂∆λi ∂∆λ j
∆λk =0∆λi ∆λ j .

(4.18)The metric in the template bank is defined by the mismatch, (1 − M), which isjust the square of the proper distance between templates
ds2 ≡ 1 − M = gij ∆λi ∆λ j ,(4.19)

where the metric is given by1
gij (λ) = −
2



∂2 M
∂∆λi ∂∆λ j
∆λk =0.

(4.20)

Having defined this metric on the intrinsic parameter space, we can now use it
to calculate the spacing of the discrete template family required to retain a given
fraction of the ideal event rate.

4.4

A coincident multi-detector search

The first LSC ringdown search was carried out by Creighton [136] with a single
filter on data from the LIGO 40 m prototype in 1994 using the GRASP soft102

ware [137]. In 2004, Adhikari [138] performed a coincident ringdown searchon 300-hour-long data from the LIGO S2 science run using LAL software [176].
The TAMA collaboration also carried out a search for ringdown signals using
data from their 300 m interferometer [139]. In 2009, a 90%-confidence upper
limit was placed on the rate of ringdowns from black holes with mass between
85 M

and 390 M

in the local universe, assuming a uniform distribution ofsources, of 3.2 × 10−5 yr−1 Mpc−3 [133, 134]. This search was carried out on datafrom LIGO S4 science run, which took place between February 22 and March24, 2005. Let us name it the “LIGO S4 ringdown search”. This includes a total of 567.4 hours of analyzable data from H1, 571.3 hours from H2, and 514.7
hours from L1. This yielded 360 hours of data in coincidence among H1, H2,and L1. To check whether a signal is in coincidence between detectors, a test isapplied separately on time, frequency and quality factor. The upper limit waspresented based on the study of simulated ringdown only signals added in thedata stream by software. A weakly modeled burst search for GWs from mergers
of non-spinning intermediate-mass binary black holes was performed on datafrom LIGO S5 and VIRGO VSR1 science runs. A 90%-confidence upper limitof 0.13 Myr−1 Mpc−3 is placed on the rate of non-spinning sources with component masses m1 = m2 = 88 M [140]. In the following sections, we present themethod for searching for GWs from perturbed BH ringdown signals in the datafrom LIGO S5/S6 and VIRGO VSR2/VSR3 science runs. We refer to it as the“LIGO-VIRGO S5/S6 ringdown search”.

103

4.4.1

The coincidence search pipeline

The “ringdown search pipeline”, illustrated in Fig. 4.3, is a multi-detector data
analysis pipeline designed for searching the l = m = 2 quasi-normal mode ofgravitational waves from perturbed black holes [134, 133]. Here we summarize
the main steps. The first stage of the pipeline involves reading in and conditioning the data from each of the detectors. We read in uncalibrated data at asampling rate of 16384 Hz, low-pass filter it to remove any power above 4096 Hz
as a part of downsampling it to 8192 Hz to reduce the computational cost. The
data is then high-pass filtered to remove power below 40 Hz. This is converted tostrain by applying the detector response function. The one-sided power spectral
density is calculated for each 2176s long segment of the calibrated data. The datais then broken further into sets of 16 overlapping blocks, each 256s in length, andfiltered using a bank of ringdown templates. This segmentation is discussed in
more detail in Ref. [185]. Once the triggers are generated, we require coincidenceof a trigger from at least two detectors in waveform parameters to increase theconfidence level for the presence of a signal. At this stage we also veto triggers
occurring during times when data quality flags are on.

4.4.1.1

Data set

The search uses data from LIGO S5 and S6 science runs and from VIRGO VSR2
and VSR3 science runs. LIGO S5 science run took place between 2005 November
4 and 2007 October 1, during which time the detectors reached the initial LIGO(iLIGO) design sensitivity. The detectors had duty factors of 78% for H1, 79%
for H2, and 66% for L1. We exclude data from VIRGO VSR1 science run. That104

Data­find
H1

H2

L1

Generate bankof templates

Generate bankof templates

Generate bankof templates

Match­filter
data

Match­filter
data

Match­filter
data

Coincidence test: time,
waveform parameters

DiscardFail

Pass
Candidateevent

Coherent stage: phase andamplitude consistency check

Figure 4.3: A schematic diagram of the coincidence and coherent (see Section 4.5)
stages of the ringdown search pipeline. The LIGO detectors (H1, H2 and L1) areconsidered as an example of this diagram.

means we only analyze data from H1, H2, and L1 detectors during S5/VSR1.
We sometimes refer to this as the ”LIGO S5 ringdown search”. LIGO S6 sciencerun was held from 2009 July 07 to 2010 October 20. During this run, the 4 km
H1 and L1 detectors were operated at sensitivities that surpassed that of theprevious S5 run, with duty factors of 52% and 47%. The 2 km H2 detector
was not operated during S6. VIRGO VSR2 science run was held from 2009 July
07 to 2010 January 08 with an improvement in sensitivity roughly a factor of105

2 over VIRGO VSR1 science run. VIRGO VSR3 was held from 2010 August
11 to 2010 October 20. The overall Virgo duty cycle over VSR2 and VSR3 was
78%. Figure 1.4 shows the best sensitivities, in terms of noise spectral density, ofthe LIGO interferometers during S5 science run. Using the power spectrum we
can estimate a detector’s horizon distance, DH , which is the distance at whicha specified source with optimal location and orientation produces a SNR of 8
in that detector. We consider a spinning black hole with â = 0.7 and e = 0.01
(where e is the fraction of BH mass emitted as GWs). Figure 4.4 shows thehorizon distance as a function of mass for LIGO/VIRGO detectors at the best
sensitivities during their joint S5/VSR1 science run.

4.4.1.2

Template placement

The ringdown analysis pipeline constructs a two-dimensional lattice templatebank in the ( f 0 , Q) space where the mismatch between two templates differing in ringdown frequency by d f 0 and in quality factor by dQ is given by the
metric [133].
"
#
4
2
2
1
3
+
4Q
3
+
8Q
3
+
16Q
ds2 =
dQ2 − 2dQd f 0 +d f 02 ,2
8 Q2 (1 + 4Q2 )2
f 0 Q(1 + 4Q2 )f0(4.21)

which assumes white noise. The mismatch is defined in terms of φ ≡ log( f 0 ) sothat the metric coefficients in the ( Q, φ) variables are independent of f 0 :1
3 + 16Q4
3 + 4Q22222dQ − 2dQdφ + (3 + 8Q )dφ , (4.22)
ds =8 Q2 (1 + 4Q2 )2
Q(1 + 4Q2 )= gQQ dQ2 + gQφ dQdφ + gφφ dφ2 .

(4.23)

Currently, the population algorithm uses only the diagonal terms containingd f 02 and dQ2 . Starting with the smallest central frequency and quality factor, the106Figure 4.4: The ringdown horizon distance as a function of mass for a black holewith spin â = 0.7 and e = 0.01 plotted for LIGO and VIRGO detectors. For
LIGO 4 km detectors during S5 science run, the ringdown horizon distance for a
black hole of mass 210 M reached about 600 Mpc. The quantity in parentheses
represents date in YYMMDD format when the horizon distance shown here wasmeasured.
algorithm places templates along Q in steps of [145, 146]dQ =dseff Q(1 + 4Q2 )p,
3 + 16Q4(4.24)

and along dφ = d log( f ) in steps ofdφ = pwhere dseff = 4

p

dseff
3 + 8Q2,(4.25)(1 − MM) and MM is the specified minimal match between

the template and signal. The templates are placed in log( f 0 ) and Q space suchthat the mismatch between any signal and its closest template is 3% or less. Figure 4.5 shows the ringdown template bank for f 0 ∈ [50, 2000] Hz, Q ∈ [2.1187, 20]
107

used in this study. For this bank, ds2 = 0.03 and it has a total of 615 templates.
The frequency spacing of the templates is finer for lower frequencies and forhigher quality factors. The template placement can be well understood by looking at the contours of constant mismatch around each template in the bank.
These contours form overlapping ellipses in the template parameter space. Eachcontour covers a region in the template parameter space so that any given signal
in that region will match a template by no less than 97%. Figure 4.6 illustrates
the contours around templates for a single frequency and for a snippet of thetemplate bank.
20

Q

15

10

5

0

2

3

10

10
f [Hz]

Figure 4.5: The template bank for the LIGO-VIRGO S5/S6 ringdown search.
Following Refs. [133, 136, 134], the template used in this search is given by
π f0 thc (t) = e− Q cos(2π f 0 t) ,0 ≤ t ≤ tmax

(4.26)

with a length of 10 e-folding times, tmax = 10τ, where τ = Q/π f 0 . In order
to illustrate the parameter space covered in this search, we plot expected signals
108

30
25

30

ds2=0.03
template

25
20
Q

Q

20
15

15

10

10

5

5

0
39.8

ds2=0.03
template

44.7

50.1
f [Hz]

56.2

63.1

0
39.8

44.7

50.1
f [Hz]

56.2

63.1

Figure 4.6: Contours of ds2 = 0.03 around templates of a constant frequency(left) and for a portion of the template bank (right).

in time domain (see Fig. 4.7), in the extremes of the parameter space. Figure 4.7
shows the signals in the extremes of the parameter space characterized by frequency and quality factor. Using the relations given in Eqs. (4.12) and (4.13),
we tabulate the masses and spins, corresponding to the extreme frequencies andquality factors, in Table 4.1.

4.4.1.3

Trigger generation

For the ringdown template given in Eq. (4.26), the SNR is|h x, hc i|| Zc |ρc ( hc ) = p≡.
σch hc , hc i

(4.27)

In practice, the detector output is not continuous but is a discretely sampled
quantity with sampling interval ∆t seconds,
x j ≡ x (t j ) ,109(4.28)

−20

1.5

x 10

r = 100 Mpc [M = 605.1196 M⊙ , â = 0.9947]−22

3

f0 = 50 Hz, Q = 20

x 10

r = 100 Mpc [M = 15.1280 M⊙ , â = 0.9947]
f0 = 2000 Hz, Q = 2021

1h0 (t)

h0 (t)

0.5
0

0
−1
−0.5

−2

0

0.2
−20

x 10

0.4

0.6
t [s]

0.8

1

1.2

0

r = 100 Mpc [M = 237.9536 M⊙ , â = 0]

−22

x 10

f0 = 50 Hz, Q = 2.1187

0.01

0.015
t [s]

0.02

0.025

0.03

r = 100 Mpc [M = 5.9488 M⊙ , â = 0]f0 = 2000 Hz, Q = 2.1187

5

1.5

4

1

3

0.5

2
h0 (t)

h0 (t)

0.005

1

0
0
−0.5
−1
−1
−1.5
0

−2
0.02

0.04

0.06
0.08
t [s]

0.1

0

0.12

0.5

1

1.5
t [s]

2

2.5

3
−3
x 10

Figure 4.7: Ringdown waveforms in the extremes of the parameter space, characterized by central frequency and quality factor. To illustrate the short natureof these waveforms, the e-folding time is also marked with a dashed line.

where t j = j∆t. Thus, in order to filter the data we need to recast the expressions
described above. The digital matched filter operates on a single data segment
consisting of N consecutive samples of x (t j ). The length of this data segment
is T = N∆t s. The frequency domain quantity x̃ ( f k ) denotes the value of thecontinuous function x̃ ( f ) at a particular frequency, labeled
fk =

k
.
N∆t

110

(4.29)

Table 4.1: Extreme (pair) values for frequency and quality factor of the ringdown
template bank, together with their corresponding masses and spins.
f 0 (Hz)Q

M (M )â

50

2.1187237.9536

02000

2.11875.9488050

20605.1196 0.9947

2000

20

15.1280

0.9947

If the units of x j are counts, then x̃ ( f k ) has units of counts/Hz. We define x̃k byx̃k =x̃ ( f k ),∆t(4.30)

which has units of counts. Negative k defines negative frequencies. Since x (t j ) is
sampled at intervals of ∆t, the sampling theorem implies that x (t j ) is bandwidth
limited to the frequency range − f Ny ≤ f ≤ f Ny , where f Ny = 1/(2∆t) is the
Nyquist frequency. Any power in x (t) at frequencies above f Ny will be aliased
into the range − f Ny ≤ f ≤ f Ny , thus, corrupting the data. To prevent this, dataof frequency higher than f Ny are removed using analog low-pass filters before
the data is digitized. Therefore, x (t j ) completely determines the data x (t) in theband of interest. Thus the Fourier transformed data is expressed as
N −1x̃k = ∑ x j e−i2πjk/N ,(4.31)

j =0

and its inverse Fourier transform is given byxj =1 N −1 i2πjk/Nx̃ j e
.
N k∑
=0
111

(4.32)

Similarly, the template used in this search is given in its discrete form by
hc j ≡ hc ( t j ) = e −π f 0 j∆tQcos(2π f 0 j∆t) ,(4.33)The template time series is Fourier transformed to get h̃c k ≡ h̃c ( f k ).
The filtering is done in Fourier domain template by template for each of thesixteen 256s overlapping block of data in a 2176s analysis segment. Then thefilter output is inverse Fourier transformed to get the time series
)
(
N/2−1

2
Sv−1
<
| Zc j | =∑ dyn R (∆t x̃k∗ ) ∆t h̃c k e−i2πjk/N ,
N∆t
k =1

(4.34)

where the inverse calibrated power spectrum
1
1
=
Sv−1
2
2
ξ( f )
dyn R(4.35)

is expressed in terms of uncalibrated spectral density Sv ( f ), response function
R( f ), and the “dynamical range factor” dyn = 1020 [134]. The numerical valueof the response function is very small, so to save the computational cost of extra
precision we scale this quantity by the dynamical range factor. Note that in
Eq. (4.34) we set the DC (k = 0) term and the Nyquist frequency (k = N/2)
term to zero. The variance of the template is evaluated as
σc2 =

2

2 i
4 N/2−1 Sv−1 h 
<
∆t
h̃
+
=
∆t
h̃
dyn2 .
ck
ck
2 | R |2N∆t k∑dyn=1(4.36)

Using the pre-defined threshold on the SNR, ρc∗ , we calculate the equivalentthreshold on the filter output, | Zc |∗ . For a given template we compute| Zc |∗ = ρc∗

σc
2 dyn

(4.37)and compare it with | Zc |. We cluster those triggers above threshold in time of 1s
sliding window [134]. The trigger with the loudest SNR in a cluster of triggers
112

is marked as the trigger to be recorded. The SNR of the clustered triggers is
calculated asρc (t j ) = | Zc (t j )|

2 dyn
.
σc

(4.38)

The amplitude A is calculated (according to Eq. (4.9)) from the template parameters f 0 and Q, and with e = 0.01 and the distance to the source r = 0.5
Mpc. The effective distance is a derived quantity based on the amplitude andthe sensitivity of the detector to a source at 1 Mpc. This sensitivity is calculatedas
σ2 = (1 Mpc) σc2 A2 ,(4.39)

from which the effective distance of the trigger is estimatedσ
Deff = (1 Mpc) .
ρc(4.40)All of the parameters mentioned here are recorded for each trigger.

4.4.1.4

Coincidence analysis

Once triggers are found in one detector they are checked for parameter consistency and time delay with triggers from other detectors, operating concurrently.
This is commonly known as a coincidence test. Triggers failing this test are discarded. Whereas, triggers fulfilling the test are recorded as coincidences and arefollowed up further because they can be signals.
Random noise is considered to be uncorrelated between geographically separated detectors. Thus, requiring coincidence between triggers of separated detectors lowers the false alarm rate significantly [150]. In LIGO S4 ringdown search,
the coincidence test was composed of two separate stages. The time requirement113

is enforced first, only allowing triggers that appear within 4 ms of each other for
co-located detectors (H1-H2 pair) and 14 ms for the H1-L1 pair. Note that H1and L1 observatories are separated by 10 ms of light travel time. The second
coincidence test is based on the metric used for template placement to measure
the distance between the triggers’ template parameters [133]. However, it was
found that the time of arrival difference of recovered coincidences is actually frequency dependent, which led to having a frequency dependent time-difference
cut, instead of a flat dt cut. In this S5 search, we apply a coincidence test thatuses the 3D metric in ( f 0 , Q, t) space to compute the distance between coincident
triggers [146]:
ds2 =

1 + 6Q2 + 16Q4 2 1 + 28Q4 + 128Q6 + 64Q8
dQ2
d f0 +
4Q2 (1 + 6Q2 + 8Q4 )2
4 f 02 (1 + 2Q2 )π 2 f 02 (1 + 4Q2 ) 2
1 + 2Q2 + 8Q4
dt−
d f 0 dQ
Q2
2 f 0 Q(1 + 6Q2 + 8Q4 )2π f 0 (1 − 2Q2 )2πQ(1 + 4Q2 )dtdf+
dt dQ .
−
0
1 + 2Q2(1 + 2Q2 )2+(4.41)

Using this metric allows us to specify a single threshold that simultaneously
covers variations in f 0 , Q, and t. By using the 3D coincidence test, we recover
the same number of signal triggers, compared to S4 coincident test, while reducing the level of false triggers. Simulation experiments were used to set thecoincidence threshold at ds2 = 0.4.
At this stage of the coincident search pipeline we also veto triggers occurring
during times when data quality flags were on. Data quality flags are intended
to indicate periods of data taking which suffer from environmental and instrumental effects, inducing noise into the data [147, 148, 149]. We differentiate
among three categories of environmental and instrumental vetoes: a) CAT1 114

Data marked as CAT1 indicate severe malfunctions of the detector. Data segments surviving the CAT1 veto are used in the analysis. b) CAT2 - Data marked
as CAT2 indicate well-understood short duration noise transients. Data segments flagged by this category should not be used for upper limit calculations
but are used for searching detection candidates. An example of a CAT2 flag iswhen simulated waveforms are present in the data in the form of “hardware
injections”, which are made for calibration purposes. c) CAT3 - Data marked as
CAT3 denote periods when the GW strain channel is only weakly correlated toenvironmental and instrumental artifacts. Detection candidates found at thesetimes should be followed up with care.

4.4.2

Tuning the search

In any experiment, there will always be false triggers that mimic the behaviorof signals. Typically, a detector’s false alarm rate due to the background canbe determined simply by turning off the source or changing the orientation of
the detector. However, these are not options for GW detectors, which are moreor less omni-directional and cannot be shielded from gravitational waves. So,
we estimate the background due to accidental coincidences of noise by shifting
triggers in one detector relative to those in another by durations greater thanthe light-travel-time between those detectors. The shifts minimize the chanceof putting an actual GW trigger in the background sample. We refer to theseas background triggers or slide triggers, as opposed to the in-time coincident
triggers (i.e., zero-lag triggers or foreground triggers) obtained without suchtime shifts. Background triggers corresponding to times which are flagged by115

data quality studies are discarded. In LIGO S5 ringdown search, H2 data is not
slid with respect to H1 data because the correlation of noise in these co-locateddetectors under-estimates the rate of false alarms [146]. However, we slide L1
data on a ring fifty times in steps of 5s and fifty times in steps of -5s, forming
100 experimental trials with no true signals [185, 134].
In order to test the sensitivity of our search to GWs from black holes, we add
various simulated signals to the data stream. These are known as ”software injections”. We then run the search pipeline to find them in the noise. To measure
the efficiency of the search we simulate three different populations of waveforms [146]. These include (a) PHENOM - phenomenologically motivated postNewtonian inspiral-merger-ringdown (IMR) gravitational waveforms for binary
black holes [152, 153, 154, 155, 156], (b) EOBNR - effective-one-body (EOB) postNewtonian model fitted to numerical-relativity simulations [157, 158], and (c)RINGDOWN - quasi-normal modes of the perturbed black hole as discussed inSec. 4.2.3. The injection parameters are chosen in such a way that they covera wide range of signal parameter space. The injections are placed at random
times within an interval of 300s distributed both uniformly in distance and in
logarithmic distance between 1 Kpc and 1 Gpc and uniformly in initial phase,sky position, inclination and polarization angles.
We utilize the injection and the slide triggers to find a balance between recovering as many simulated signals in coincidence between multiple detectors
as possible while keeping the rate of false coincidences to a minimum. This
process is known as tuning a search pipeline. We finished tuning the LIGO S5
ringdown search recently. We find the SNR thresholds to be 5.5, 4.0, and 5.5 for116H1, H2, and L1 detectors, respectively, and the 3D coincidence threshold to be
ds2 = 0.4 for each detector. We found that the SNR thresholds for the LIGO
S4 ringdown search were set at higher values (see Appendix D.1). This was
one of the main reasons why no triple-coincident background trigger was foundin H1, H2 and L1 detectors. It is also found that the number of slide triggers
grows exponentially as the SNR threshold is lowered. So, to strike a balance
between the computational cost and the detection efficiency, we only lowered
the SNR threshold for the H2 detector. This is also motivated by the fact that thesensitivity of H2 detector is roughly half of the sensitivity of H1 or L1 detector.
We have not opened any boxes on the analysis yet. So, here we present somepreliminary results from the tuning exercise. These results are produced using
a portion of the LIGO S5 data. The results presented here contain only software injections and time-slides after CAT1, CAT2 and CAT3 vetoes have been
applied. In Fig. 4.8, we show the sensitivity of the search pipeline to GWs fromblack holes by comparing made and found software injections in about 30 daysof LIGO S5 data when H1, H2 and L1 detectors were operating in science mode.
These plots include RINGDOWN as well as IMR (i.e., EOBNR and PHENOM)
injections. Note that in the following plots, "Sim." means simulated parameters,and "Sngl." means estimated signal parameters. In Figs. 4.9 and 4.10, we separately study the PHENOM injections. Estimated effective distance is higher thansimulated effective distance for a large number of IMR injections. We understand that the parameter estimation is bad for unequal mass binaries and theestimated ringdown mass is consistently higher than the simulated ringdown
mass. And in Figs. 4.11, 4.12, and 4.14, we study the RINGDOWN injections.

117

We find bias in frequency recovery for higher mass black holes. Better mass
recovery is found for signals with higher quality factor. This is because signals
with higher quality factor have more cycles than signals with lower quality factor. The ringdown templates match earlier with the ringdown phase of the IMR
signal, i.e, they match with the merger. Thus, effectively, they pick up a longer”ringdown” signal with a larger quality factor and lower frequency. We also findthat inclusion of IMR waveforms improves detection efficiency. This is because
the ringdown only signals at high frequency are more difficult to recover thanthose at low frequency because detector sensitivity decreases at higher frequencies. However, the inspiral and merger phases preceding the ringdown phasesweep up through the detector band, finding a significant match with ringdown
templates.
The coincidence statistic is a ranking mechanism using the SNRs of coincident triggers. The exact form of the statistic for a given population depends on
the properties of the SNR distributions of simulated and background signals (seeFigs. 4.14, 4.15, and 4.16). In the LIGO S4 ringdown search, no triple-coincident
event was found. However, a ranking statistic was proposed, which was motivated by Gaussian distribution of triggers, and is given by [133]ρD,S4 =

q

ρ2H1 + ρ2H2 + ρ2L1 .

(4.42)

Note that the original equation in Ref. [133] does not have the square root andwas named the detection statistic there. This was done in order to enhance thelikelihood of ranking the triple-coincident (injection) triggers higher than the118

Figure 4.8: Decisive distance of missed and found injections is plotted against
chirp mass (top) and frequency (bottom). The stars are injections found withzero combined false-alarm rate (FAR), and the circles are injections found with
non-zero combined FAR. The red crosses are missed injections (vetoed injections
are excluded). Note that the decisive distance is the second smallest effective
distance. Chirp mass is defined in Sec. 5.2.1.
119Figure 4.9: Scatter plots of simulated and estimated effective distance of found
PHENOM injections. The color-bar represents estimated frequency (top left),
quality factor (top right), mass (middle left), and spin (middle right). The colorbar of the bottom plot represents simulated mass-ratio. These parameters are
estimated in H1 detector.
120

Figure 4.10: Scatter plots of simulated mass 2 versus mass 1 of made and found
PHENOM injections (left). Histogram of number of made PHENOM injections
per mass-ratio bin is compared with histogram of number of found PHENOM
injections per mass-ratio bin (right).

double-coincident ones; the double coincident statistic was introduced asρDS = min ρifo1 + ρifo2 , 2ρifo1 + 2.2 , 2ρifo2 + 2.2 .

(4.43)The optimal statistic is the sum of SNR squares if the distribution is Gaussian.
However, the distribution in real data was found to have long "tails", i.e, coincidences with a very loud SNR in one detector and a much lower SNR in theother, which motivated the form in Eq. (4.43). For the LIGO S5 ringdown search,the distribution of double-coincident triggers follow the same distribution as in
the LIGO S4 ringdown search. So, we use this statistic to rank double-coincident
triggers. However, by lowering the SNR threshold and with better understanding of systematics, sources and the pipeline, we managed to generate hundredsof triple-coincident background triggers. Their SNR distributions are shown in
121

Figure 4.11: Scatter plots of simulated and estimated effective distance of found
RINGDOWN injections. The color-bar represents difference between simulated
and estimated parameters: Frequency (top left), quality factor (top right), mass(bottom left), and spin (bottom right). These parameters are estimated in H1
detector.

122

Figure 4.12: Scatter plots of simulated and estimated frequency of found RINGDOWN injections. The color-bar represents simulated mass (top left) and difference between simulated and estimated mass (top right). The scatter plots ofsimulated and estimated mass of found RINGDOWN injections (bottom). Here,the color-bar represents simulated quality factor. These parameters are estimatedin H1 detector.

123

Figure 4.13: Scatter plots of simulated quality factor and frequency (left), andspin and mass (right) of made and found RINGDOWN injections.

Figs. 4.14 and 4.16, along with the injection distributions in Figs. 4.15, and 4.16.
We find the background triggers tend to tail along L1 SNR. Since we do nothave any signal-based veto for this pipeline, we down-weight triggers in the tailof one IFO significantly. This was done by defining a new triple-coincidence
ranking statistic (see Fig. 4.16):
ρTR =

h

nq

ρ2H1 + ρ2H2 + ρ2L1 , ρH1 + ρH2 + 0.75 ,oiρH1 + ρL1 + 0.75 , ρH2 + ρL1 + 0.75

min.

(4.44)

ρH1 /ρH2 >1This new triple-coincidence ranking statistic is compared with the old triplecoincidence ranking statistic in Fig. 4.17. There are hundreds of injection triggers
which become louder than the loudest background trigger in this new ranking
statistic, that otherwise would not have been achieved. This shows significant
improvement in detection efficiency. More quantitatively, the detection probability of this new statistic is higher than that of the old statistic at low false-alarm124

probability (see Fig. 4.18).

4.5

The coherent multi-detector search

A study performed on data from LIGO S4 science run shows that the Coherent
WaveBurst (cWB) search is more sensitive to ringdown signals than the ringdown search. Note that the cWB pipeline is designed to search for unmodeled
burst signals [151]. We understand that this gap in performance is due to thefact that the templated ringdown search described above did not have the infrastructure for checking the coherence of signals from a common astrophysicalsource in multiple detectors. To improve this scenario we implement several
consistency checks. For example, we check for consistency of signal amplitudes
and phases in the different detectors with their different orientations and withthe signal arrival times in them. In the following sections, we discuss the utilities
for the coherent stage, namely the coherent SNR and null-stream statistic, of the
ringdown pipeline.

4.5.1

Elliptically-polarized ringdown template

Traditionally, the coincidence search pipeline has not searched for the initial
phase of the ringdown signal. In other words, it uses only single-phase ringdown templates (see Eq. 4.26) to filter the data. However, this makes the pipeline
limited from having enough information to check for signal amplitude andphase consistency across the detectors. A convenient way to search for the unknown phase of a signal is to filter the data with both phases of the template125

H1H2L1 in H1H2L1−time: CAT3

H1H2L1 in H1H2L1−time: CAT3

2

1.3

1.3
2

1.8

1.2

1.6

1.2

1.4

1.5

1
1

0.9

log10(ρL1)

log10(ρH2)

1.1

0.8

1.1

1.2

1

1
0.8

0.6

0.9
0.5

0.7

0.4

0.8
0.8

0.9

1
1.1
log10(ρH1)

1.2

1.3

0.2

0

0.8

0.9

1
1.1
log10(ρH1)

1.2

1.3

0

H1H2L1 in H1H2L1−time: CAT3
1.3
2

1.2

log10(ρL1)

1.5

1.1

1

1

0.9
0.5

0.8
0.7

0.8

0.9
1
log10(ρH2)

1.1

1.2

1.3

0

Figure 4.14: Contour plots of H2 SNR versus H1 SNR (top let), L1 SNR versus
H1 SNR (top right), and L1 SNR versus H2 SNR (bottom) for triple coincident
time-slide events. The color-bar represents log10 ( N ), where N is the number oftriggers.

126

H1H2L1 in H1H2L1−time: CAT3

H1H2L1 in H1H2L1−time: CAT3

2

1.6
2

1.8

1.4

1.8

1.2

1.6

1.6

1.5

1

1.2

log10(ρL1)

log10(ρH2)

1

1.4

1.4

0.8

0.6

1.2
1

0.4

0.5

1

0.8

0.2

0.8
0.8

1

1.2

1.4
log10(ρH1)

1.6

0

1.8

0.8

1

1.2

1.4
log10(ρH1)

1.6

1.8

0

H1H2L1 in H1H2L1−time: CAT32

1.8

log10(ρL1)

1.6

1.5

1.4
1

1.2
0.5

1

0.8
0.8

1

1.2
1.4
log10(ρH2)

1.6

1.8

2

0

Figure 4.15: Contour plots of H2 SNR versus H1 SNR (top let), L1 SNR versus
H1 SNR (top right), and L1 SNR versus H2 SNR (bottom) for triple coincident
injection events. The color-bar represents log10 ( N ), where N is the number oftriggers.

127

H1H2L1 in H1H2L1−time: CAT3

24

24

22

22

20

20

18

18

16

16
ρL1

ρH2

H1H2L1 in H1H2L1−time: CAT3

14

14

12

12

10

10

8

8

6

6

4

5

10

15
ρH1

20

4

25

5

10

15
ρH1

20

25

H1H2L1 in H1H2L1−time: CAT3
24
22
20
18

ρL1

16
14
12
10
8
6
4

5

10

15
ρH2

20

25

Figure 4.16: The H1-H2 (top left), H1-L1 (top right), and H2-L1 (bottom)
SNR distributions for time-slide (black crosses) and injection triggers (redpluses). The curves represent the contours of constant values of the new triplecoincidence ranking statistic defined in Eq. 4.44.

128

H1H2L1 in H1H2L1−time: CAT3
30
28
26
24

ρTR

22
20
18
16
14
12
10
8

10

15

20
ρD,S4

25

30

Figure 4.17: Scatter plots of the old and new triple-coincidence ranking statistics.
Time-slide triggers are represented by the black crosses and injection triggers are
represented by the red pluses.

separately [47]. This straightforwardly applies to the frequency domain inspiralsignal in the stationary phase approximation [185]. This is because the sine andthe cosine phases of the inspiral template are exactly orthogonal. But this is nottrue for a generic damped sinusoid. In the following sections, we justify thisstatement.

4.5.1.1

Fourier transformation of a generic decaying sinusoid

In general, a decaying sinusoid can be expressed as a real function
x (t) = A e−αt cos(ωt − ϕ) u(t)129

(0 ≤ t ≤ ∞ ) ,(4.45)

Fraction of injections found louder than loudest background

H1H2L1 in H1H2L1−time: CAT3
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1

ρD,S4
ρTR
InjFindEff

0 1
10

2

3

10
Deff
(Mpc)
H1

10

H1H2L1 in H1H2L1−time: CAT3
0.57
0.565

Detection probability

0.56
0.555
0.55
0.545
0.54
ρD,S4

0.535

ρTR
0.53
0

0.05

0.1
False−alarm probability

0.15

0.2

Figure 4.18: The efficiency of finding injections using the old triple-coincidence
ranking statistic is compared with that of using the new triple-coincidence ranking statistic (top). The maximum achievable efficiency, for the injections considered here, is shown in "InjFindEff”. In this case, we assume that all foundinjections are louder than the loudest background. The ROC curve of the old
triple-coincidence ranking statistic is compared with the ROC curve of the new
triple-coincidence ranking statistic (bottom).

130where u(t) is the Heaviside step function
 1for t ≥ 0 .
u(t) = 0otherwise .

(4.46)Let us define a complex function
z(t) = A e−αt ei(ωt− ϕ) u(t)(4.47)

such that x (t) is the real part of z(t). The Fourier transformation of z(t) is givenbyz̃( f ) ==Z ∞z(t) e−2πi f t dt

(4.48)

A e−iϕ,α + i (2π f − ω )(4.49)A eiϕ.
α + i (2π f + ω )(4.50)−∞and hence,z̃∗ (− f ) =Now, the Fourier transform of the real part of z(t) is given by
1 ∗[z̃ (− f ) + z̃( f )]2
A eiϕ [α + i (2π f − ω ) + e−iϕ [α + i (2π f + ω )]=.
2
α2 + 4iπ f α − 4π 2 f 2 + ω 2

x̃ ( f ) =(4.51)(4.52)In the next section, we show how the above relations are used to get the Fouriertransform of an elliptically-polarized ringdown template.

4.5.1.2

Fourier transformation of an elliptically-polarized ringdown template

The elliptically-polarized ringdown template is expressed ashep (t) = hc (t) − ihs (t) ,131(4.53)

where hc (t) and hs (t) are two damped sinusoid functions given byhc (t) = e−π f0 t/Q cos(2π f 0 t) ,(4.54)

hs (t) = e−π f0 t/Q sin(2π f 0 t) .

(4.55)The real part of hep (t) can be found from Eq. (4.45) by setting A = 1, α =π f 0 /Q, ω = 2π f 0 and ϕ = 0. Hence, the Fourier transformation of hc (t) isfound Eq. (4.52), and is given by
π f0
Q + i2π f.
h̃c ( f ) =π2 f 24π 2 f f
4π 2 f 02 − 4π 2 f 2 + Q20 + i Q 0(4.56)

Similarly, the imaginary part of hep (t) is found from Eq. (4.45) by setting A =1, α = π f 0 /Q, ω = 2π f 0 and ϕ = π/2. Hence, the Fourier transform of hs (t) isfound from Eq. (4.52), and is given byh̃s ( f ) =2π f 0π2 f 24π 2 f 02 − 4π 2 f 2 + Q20 + i

4π 2 f f 0
Q.

(4.57)

By combining Eqs. (4.56) and (4.57), one gets the Fourier transformation of hep (t).
4.5.1.3

The template normalization

Filtering the data x (t) with an elliptically-polarized template hep (t; µi ) characterized by the source parameters µi yields the SNR statistic|h x, hep i|,ρ(hep ) = q
hhep , hep i|h x, {hc − ihs }i|q,hhep , hep iph x, hc i2 + h x, hs i2q=.
hhep , hep i=

132(4.58)(4.59)(4.60)

Note that the noise-weighted scalar product introduced above is defined for thetwo vectors X (t) and Y (t) by

h X, Y i ≡ 2

Z ∞
−∞

df

X̃ ( f )Ỹ ∗ ( f ).
ξ (| f |)(4.61)The detailed calculation for the template normalization is as follows:2
σep= hhep , hep i ,(4.62)= h{hc − ihs }, {hc − ihs }i ,



Z ∞
h̃c ( f ) − i h̃s ( f ) h̃c∗ ( f ) + i h̃s∗ ( f )= 2df,ξ (| f |)
−∞
Z ∞
Z ∞
h̃c ( f )h̃c∗ ( f )h̃s ( f )h̃s∗ ( f )
= 2
df
+2
df
ξ (| f |)ξ (| f |)
−∞
−∞
Z ∞
Z ∞
∗
h̃s ( f )h̃c ( f )h̃c ( f )h̃s∗ ( f )
−2i
df
+ 2i
df
,
ξ (| f |)ξ (| f |)
−∞
−∞
Z ∞
Z ∞
h̃c ( f )h̃c∗ ( f )h̃s ( f )h̃s∗ ( f )
= 2
df
+2
df
,
ξ (| f |)ξ (| f |)−∞−∞= h hc , hc i + h hs , hs i .

(4.63)(4.64)(4.65)(4.66)(4.67)The justification for canceling the last two terms in Eq. (4.65) is numericallyshown in Fig. 4.19. It is important to note thath hc , hc i 6 = h hs , hs i ,for a range of µi .

(4.68)

This is because hc and hs are not orthogonal due to the damping factor. In thefollowing section, we numerically justify this statement.

4.5.1.4

Numerical analysis

To discretize the elliptically-polarized template we consider N data points (assume N to be even) sampled over a given time. When one Fourier transforms a133

function one seeks a discrete frequency array
fk =

k
;
N∆t

k=−

N
N
, ...., ,22(4.69)

where k is an integer and ∆t is the inverse of the sampling frequency f s , i.e,∆t = 1/ f s . We choose the sampling frequency 4096 Hz, with 256 data points. For
simplicity, we consider noise to be white Gaussian, and set its power spectrum,
which is constant, to unity, i.e, ξ (| f |) = 1. In Fig. 4.19 we consider an ellipticallypolarized template characterized by the central frequency f 0 = 1200 Hz and thequality factor Q = 4. We plot the integrands of the last two terms of Eq. (4.65) asfunctions of frequency: The top left panel represents their real parts and the top
right panel represents their imaginary parts. Their sum (integration) over thefrequencies is also noted. So, the last two terms of Eq. (4.65) together contributezero (except for negligible numerical errors) in the template norm expression. In
the bottom panel, we plot and compare the integrands of the first two terms of
Eq. (4.65). It proves that the first two terms are not equal.
We now consider a bank of templates that includes both sine and cosine
phases and characterized by a range of central frequencies f 0 ∈ [50 − 2000] Hz
and quality factors Q ∈ [2.1187 − 20]. For each template, we define a quantity
called “fractional difference” asFractional difference ≡h hs , hs i − h hc , hc i.
h hc , hc i

(4.70)

Figure 4.20 shows the contours of the fractional difference in the template parameter space. For higher central frequencies and lower quality factors, thefractional difference is significant. Therefore, in order to coherently search for134f0 = 1200 Hz, Q=4

−8

4

x 10

3.5

sum:
2.0052x10−6
2.0052x10−6

1.5Im(h̃s (f)h̃∗c (f))
Im(h̃c (f)h̃∗s (f))

1
Integrand

IntegrandRe(h̃c (f)h̃∗s (f))
2.5
2

x 10

2

Re(h̃s (f)h̃∗c (f))

3

f0 = 1200 Hz, Q=4

−7

3

0
sum:
−4.7976x10−23
4.7976x10−23

−1

1
−2

0.5
0
−2000

−1000

0
Frequency (Hz)
3

−3
−2000

1000

2000

−7

f0 = 1200 Hz, Q=4

x 10

−1000

0
Frequency (Hz)

1000

2000

h̃c (f)h̃∗c (f)2.5h̃s (f)h̃∗s (f)Integrand2sum:
1.4799x10−5
1.6042x10−5

1.5
1
0.5
0
−2000

−1000

0
Frequency (Hz)

1000

2000

Figure 4.19: The integrands of the sine and cosine template variances (Eq. (4.65))
are plotted as functions of frequency for a given point in the parameter space.
Here, the frequency plotted in the x-axes is meant to be the frequency f k definedin Eq. (4.69). The top left and the top right panels depict the frequency behaviorof the last two terms in expression (4.65), whereas the bottom panel depicts thesame aspect of the first two terms in that expression.

ringdown signals in a network of detectors, it is necessary to model the templateby incorporating both phases.

135Fractional difference
0.6
18.4972
0.5
16.6672
14.8372

0.4

13.0072Q

0.3
11.1772
9.34720.2

7.51720.1

5.687203.8572540

1040
f0 (Hz)1540

Figure 4.20: Contours of the fractional difference of sine and cosine phases of theelliptically-polarized ringdown template. The color bar represents the fractional
difference as defined in Eq. (4.70).

4.5.1.5

SNR statistic

The coherent stage (see Fig. 4.3) follows up on triggers that survive the coincident stage by computing their coherent and null-stream statistics. We constructsmall snippets of C-data around the peak-time of every trigger that is found to
be coincident in multiple detectors in a network. The detailed prescription is
given in the next chapter. However, this construction is done by filtering thedata around the trigger peak-time with elliptically-polarized templates from the
ringdown template bank, and the SNR is given by|h x, hep i|| Zep |≡.
ρep (hep ) = q
σep
hhep , hep i

136(4.71)

Following the same procedure as in Sec. 4.4.1.3, we find)!2
 −i2πjk/N
Sv−1∗<(4.72)∑ dyn R (∆t x̃k ) ∆t h̃ck e
k =1)!2 #1/2
(
N/2−1

Sv−1
, (4.73)<∑ dyn R (∆t x̃k∗ ) ∆t h̃sk e−i2πjk/N
k =1

2| Zep j | =
N∆t

+

"

(

N/2−1

and the template variance is given by
2
σep
=

N/2−1

2

2 i
Sv−1 h 
<
∆t
h̃
+
=
∆t
h̃
(4.74)ck
ck∑ dyn2 | R|2
k =1
!
N/2−1
2

2 i
Sv−1 h < ∆t h̃sk+ = ∆t h̃sk
dyn2 . (4.75)+ ∑2 | R |2
dyn
k =1
4
N∆t

The SNR-threshold procedure introduced in Sec. 4.4.1.3 is also followed here,
but with the quantities defined in terms of Zep , ρep , and σep .

4.5.2

Coherent statistics

In this section we briefly outline the basic expressions for the coherent SNR andthe null-stream statistics for use in the ringdown search pipeline. We postpone
the more detailed discussion of those statistics (see, e.g., Eqs. (5.31) and (5.39)) to
the next chapter. Specifically, we focus on three cases of detector networks most
relevant to the S5 coherent ringdown searches:
Case-1: Two co-aligned detectors with different noise PSDs.
The coherent SNR for two co-aligned detectors with different noise power
spectral density is given by|C σ + C2 σ2 |$12 = p 1 1,(σ1 )2 + (σ2 )2137(4.76)

where σI is the template-norm
σI ≡ σep I =qhhep I , hep I i ,(4.77)

h x I , hep I i,
σep I(4.78)andC I ≡ ρ I eiΦ I =is the matched-filter output against an elliptically-polarized template in the Ithdetector. The null-stream statistic is given by
C1
C2
σ1 − σ2η12 = r 
2
1
σ1

+

 2 .

(4.79)1
σ2

Case-2: Three detectors with two of them coaligned and colocated at one site
and the third one located at a second site, and all with different noise PSDs.
Let detectors I = 1, 2 be at the same site. The coherent SNR for this networkis given by$123 =q($12 )2 + (ρ3 )2 .

(4.80)

And the null-stream η123 is the same as the null-stream defined for Case-1.
Case-3: Three detectors at different sites with different orientations and noise
PSDs.
The antenna response functions are expressed as

 
 
 F+   cos 2ψ sin 2ψ   u =  ,
F×− sin 2ψ cos 2ψ
v(4.81)

with u(θ, φ) and v(θ, φ) being detector orientation (and source-location) dependent functions. For a network of three detectors, we introduce the following138

shorthand notations for the quantities involving these functions:A12 ≡ (u1 v2 − u2 v1 ) ,(4.82)

A23 ≡ (u2 v3 − u3 v2 ) ,(4.83)

A31 ≡ (u3 v1 − u1 v3 ) ,(4.84)

where, u I and v I are the Ith detector orientation-dependent functions. We also
define Aij = −A ji . Then the coherent and null-stream statistics for this network
are given byvuu B C +B C 2+ B C +B C 2+ B C +B C 2
32 3
23 2
31 3
21 2
13 1
t 12 1,$123 =
2
2
2(B12 ) + (B23 ) + (B31 )(4.85)

and
v
u
2
2
2
u C1
C
C
C
C
C
u B12 − B323 + B313 − B212 + B232 − B131
η123 = u, 2  2  2
t1
1
1
+ B23 + B31
B12(4.86)
where

B12 ≡ A12 σ1 σ2 ,(4.87)

B23 ≡ A23 σ2 σ3 ,(4.88)

B31 ≡ A31 σ3 σ1 ,(4.89)

are the shorthand notations for the quantities involving detector orientation,
source-location, and template normalization. Like Aij , we also have Bij = −B ji .

139

4.5.3

Performance study

The coherent stage for the LIGO-VIRGO S5/S6 ringdown search is in the processof tuning, and the results discussed here are preliminary and may change. So,
we do not present them here. However, there is preliminary evidence that thecoherent statistic performs better than the coincident statistic. To study the utility of the coherent statistics, we ran the ringdown search pipeline (see Fig. 4.3)
on a month of data from LIGO S4 science run. For this example, we consider
only H1 and H2 detectors. A total of 225 simulated ringdown signals were
injected into the data, of which 86 were found by the pipeline. The injection
finding efficiency is low because most of the injections were too weak. A total of 143 background triggers, obtained through time-slide experiments, were
found. Figure 4.21 scatter plots the coincident statistic (bitten-L) against thecoherent statistic (effective-coherent SNR) for the found injection and the slide
triggers. The effective-coherent SNR defined here is the quadrature sum of thecoherent SNR and null-stream. The curved lines in the top panel of Fig. 4.21
represent constant effective-coherent SNR values. In Fig. 4.22, we compare theefficiency of finding injection triggers using the coincidence statistic with the coherent statistic. Note that both the amplitude and phase consistency are appliedin this analysis. In Fig. 4.23 we plot and compare the ROC (receiver operative
characteristic, see Ref. [26]) curve for these two statistics. We find that the coherent statistic performs better than the coincidence statistic in the low false-alarm
regime.

1404.6Discussion

We have demonstrated using simulated signals that the ringdown pipeline is anefficient search pipeline for detecting gravitational waves from perturbed black
holes. As discussed above, we also show that the added-coherent-stage ringdown pipeline performs better than the coincidence-only ringdown pipeline.
We expect its performance to be boosted for triple-site searches (LIGO-VIRGO S6
ringdown search). It is decided to use multivariate statistical classifier (MVSC) asa detection statistic for the LIGO-VIRGO S5/S6 ringdown searches [146]. It is an
implementation of the random forest of bagged decision trees algorithm. MVSC
combines information from many statistics. So, the amalgamated information
should be better than any individual statistic alone. MVSC’s performance is thebest for a best set of parameters that best separate background from signals.
The effective-coherent SNR and the new triple-coincidence ranking statistic arefound to be the members of the best set of MVSC’s parameters.

141

1

Null−stream

10

0

10

−1

10

SlidesInjections
−1

0

10

1

10

2

10

10

Coherent SNR
2

10

Bitten−L

SlidesInjections

1

10

1

2

10

10
Effective−coherent SNR

Figure 4.21: The scatter plots for the null-stream and coherent SNR (top panel),and the chopped-L and effective-coherent SNR (bottom panel). The injection
triggers are denoted by red pluses, and the background triggers (slide triggers)
are denoted by black crosses. All background triggers have been retained. The
curved lines in the top figure are the contours of constant effective-coherent
SNR, which is a combination of the null-stream and coherent SNR. There are
10 found injection triggers, marked by green circles, which are quieter than theloudest background in the coincidence search whereas louder than the loudest
background in the coherent search.

142

Fraction of injections found louder than loudest background
Fraction of injections found louder than loudest background

Efficiency Plots
0.4
Coherent search using effective−coherent SNR
Coincident search using bitten−L

0.35
0.3
0.25
0.2
0.15
0.1
0.05
0
10

20

30

40
50
60
70
Injected distance (Mpc)

80

90

100

Efficiency Plots
0.12
Coherent search using effective−coherent SNR
Coincident search using bitten−L
0.1

0.08

0.06

0.04

0.02

0
50

60

70
80
Injected distance (Mpc)

90

100

Figure 4.22: Here we compare the efficiency of finding signals using the coincidence search with that of using the coherent search. In the top panel, we focus
on a set of stronger injections, with distance from 10 Mpc to 100 Mpc. In the bottom panel, a similar study was performed but with weaker signals, i.e, injections
were made allowing distance from 50 Mpc to 100 Mpc. The overall performanceof the coherent search is quite a bit better than the coincidence search. In the last
distance bin, from 85 Mpc to 100 Mpc, the coherent search above does not show
any improvement over the coincidence search simply because no coincidence
trigger was found in that bin.
143

ROC curves
0.31
0.3

Detection probability

0.29
0.28
0.27
0.26
0.25
0.24
0.23
0.22
0

Coherent search using effective−coherent SNR
Coincident search using bitten−L
0.05

0.1
0.15
False−alarm probability

0.2

Figure 4.23: The ROC curve of the coincidence search is compared with the ROC
curve of the coherent search. The overall performance of the coherent searchis better than that of the coincident search. This is because there are 10 foundinjections that are weaker than the loudest background trigger in the coincidence
search as shown in Fig. 4.21. In the low false-alarm region, detection probabilityof coherent search is higher than the coincident search.

144

Chapter 5
Aspects of a blind hierarchicalcoherent search for coalescing
compact binaries
We describe a hierarchical data analysis pipeline for coherently searching for
GW signals from non-spinning compact binary coalescences (CBCs) in the dataof multiple earth-based detectors. This search assumes no prior information on
the sky position of the source or the time of occurrence of its transient signals
and, hence, is termed “blind”. The pipeline computes the coherent network
search statistic that is optimal in stationary, Gaussian noise and more advanced
than the current state-of-the-art approach, called coincident statistics. More importantly, it allows for the computation of a suite of alternative multi-detector
coherent search statistics and signal-based discriminators that can improve theperformance of CBC searches in real data, which can be both non-stationary

145

and non-Gaussian. Also, unlike the coincident multi-detector search statistics
that have been employed so far, the coherent statistics are different in the sensethat they check for the consistency of the signal amplitudes and phases in the
different detectors with their different orientations and with the signal arrival
times in them. Since the computation of coherent statistics entails searching inthe sky, it is more expensive than that of the coincident statistics that do not require it. To reduce computational costs, the first stage of the hierarchical pipeline
constructs coincidences of triggers from the multiple interferometers, by requiring their proximity in time and component masses. The second stage follows up
on these coincident triggers by computing the coherent statistics. Here, we compare the performances of this hierarchical pipeline with and without the second(or coherent) stage in Gaussian noise. Whereas introducing hierarchy can be
expected to cause some degradation in the detection efficiency compared to thatof a single-stage coherent pipeline, nevertheless it improves the computational
speed of the search considerably. The two main results of this work are: (1) The
performance of the hierarchical coherent pipeline on Gaussian data is shown tobe better than the pipeline with just the coincident stage. (2) The three-site network of LIGO detectors, in Hanford and Livingston (USA), and VIRGO detector
in Cascina (Italy) cannot resolve the polarization of waves arriving from certain
parts of the sky. This can cause the three-site coherent statistic at those sky positions to become singular. Regularized versions of the statistic can avoid thatproblem, but can be expected to be sub-optimal. The aforementioned improvement in the pipeline’s performance due to the coherent stage is true in spite ofthis handicap. This work (see Ref. [12]) was the result of a collaboration that in-146

cluded two other graduate students, Thilina Dayanga and Shaon Ghosh, at the
Washington State University.

5.1

Introduction

Signals from binaries of neutron stars and black holes (BHs) enjoy the prospectof being the first signals to be detected by GW detectors [159]. They are among
the best understood of all GW sources and a large enough number of them are
expected to appear in the data of second generation detectors [160]. The last
several science runs at LIGO [161], GEO600 [162], and VIRGO [163] revealed
that searches for signals from these compact binary coalescences (CBCs) benefit
from the networking of multiple detectors because of the reduction in the rateof accidentals or false alarms, especially, from non-stationary and non-Gaussian
noise artifacts. Further, studies with injection of simulated signals show that theestimation of source parameters, such as sky position and wave polarization,is also helped by networks involving detectors at three or more sites aroundthe globe [41, 164, 165]. This is important since CBCs may not always emit
electromagnetic signals that are strong enough to be observable.
This chapter describes blind CBC search strategies, which must be contrasted
with a targeted search method [166]. An example of the latter case is a search
for GW signal triggered by a short-duration gamma-ray burst (GRB). Short GRBs
have been conjectured to be associated with NS-NS or NS-BH coalescences [167,168]. Owing to an electromagnetic counterpart, the sky-position of the short
GRB and the time of arrival of its gamma-ray signal are known in advance for147

offline searches. This implies that searches for GW signals from these sources
require three less parameters to scan, and are, therefore, computationally lessexpensive. (In reality, one searches over a several-second window around thearrival time of the gamma-ray signal because it is not clear yet how separated theemission of the gamma-ray burst and the binary-object merger are in time [169,170].) Perhaps more significantly, it reduces the probability of false-alarms and,
therefore, increases our detection confidence.
In this chapter, we address how one tackles both these issues, namely, ofincreased computational costs and false-alarm rates, affecting a blind search forsignals from CBCs with non-spinning components. To reduce the excess computational cost arising from scanning the arrival time, one introduces hierarchicalstages in the search pipeline, whereby, first, the triggers of interest are identified
in the detectors individually. This is followed by recognizing triggers that arecoincident in multiple detectors and then computing network-based statistics forthem that reveal their significance as GW candidates. (These hierarchical steps
were introduced in Ref. [171] and have been used in multiple CBC searches eversince.) The final stage is used to compute the coherent network statistics for thesecoincident triggers. To address the second problem of increased false-alarms, especially, from non-stationary noise transients, we introduce signal-based multidetector discriminators that check for consistency of the signals appearing inindividual detectors with a CBC source, after accounting for the different orientations of the detectors and the delays in their times of arrival in them.
Past experiments with multi-detector searches for gravitational-wave signals
from CBCs have shown that the statistics that are optimal in Gaussian and sta-148

tionary noise (OGSN) cease to be so in real data, in general [171, 172, 173]. Instead a function of the chi-square-weighted [174] matched-filter [175] outputs
has been found to deliver a better performance [171, 172]. This function isarrived at empirically by comparing the distribution of the matched-filter andchi-square statistics for simulated CBC signal injections with that of the background. These statistics did not, however, use the phase of the matched-filter
output to discriminate signals from noise, which a coherent statistic [46, 47] isequipped to do. We will call the former coincident statistics. Their construction
has nevertheless helped inspire techniques for obtaining empirically an effective
coherent statistic that performs better in real data than the coherent statistic ofRefs. [46, 47]. It is this statistic and its variants, which can be useful in searching for non-spinning CBC signals in real data, that we discuss in detail in thischapter.

5.2

Multi-detector statistics

We begin by describing the statistic that is optimal for coherently searching fornon-spinning CBC signals in data from multiple detectors when their noise is
Gaussian and stationary. The first part of this section gives an alternative derivation of this statistic, as compared to that available in the literature [46]. In theprocess, it introduces notation and convention followed here. It also introduces
signal parameters and variables used in the coherent search code available in
the LIGO (Scientific Collaboration) Algorithm Library LAL [176]. We then compare that statistic with the aforementioned empirically-motivated multi-detector

149

coincident statistics, which have been applied in real data.

5.2.1

Signal and noise

Consider a non-spinning coalescing compact binary with component masses
m1,2 , such that its total mass is M = m1 + m2 and its reduced mass is µ =m1 m2 /M. In the restricted post-Newtonian approximation, the two polarizations
determining the GW strain are:GMh+ (t; r, M, µ, ι, ϕc , tc ) =c2 r

tc − t
5G M/c3

−1/4

1 + cos2 ι2× cos[ ϕ(t; tc , M, µ) + ϕc ] ,

−1/4
GM
tc − th× (t; r, M, µ, ι, ϕc , tc ) =
cos ι
c2 r
5G M/c3
× sin[ ϕ(t; tc , M, µ) + ϕc ] ,

(5.1)

(5.2)

which depend on M, µ, the luminosity distance to the source r, the inclination
angle of the source’s orbital-momentum vector to the line of sight ι, the time of
coalescence of the signal tc , and the coalescence phase of the signal ϕc . Above,ϕ(t; tc , M, µ) is the orbital phase of the binary [177, 178], M = µ3/5 M2/5 is the
chirp mass, G is the gravitational constant and c is the speed of light in vacuum.
The GW strain in a detector can then be modeled as,h(t) = F+ h+ (t) + F× h× (t) ,(5.3)

where F+,× are antenna response functions that quantify the sensitivity of thedetector to the sky-position and polarization of the source,

 
 
 F+   cos 2ψ sin 2ψ   u =  ,
F×− sin 2ψ cos 2ψ
v
150

(5.4)

with ψ being the wave-polarization angle and u(α, δ) and v(α, δ) being detectororientation dependent functions of the source sky-position angles (α, δ): [179, 47]u(α, δ) ≡ −
1
1 + cos2 α cos 2δ ,2v(α, δ) ≡ − cos α sin 2δ .

(5.5)
(5.6)

Following Ref. [41], let us map the CBC signal parameters (r, ψ, ι, ϕc ), into
new parameters, a(k) , with k =1,...,4, such that the strain in any given detector
has a linear dependence on them:4h ( t ) = ∑ a ( k ) hk ( t ) ,

(5.7)

k =1

where the hk (t)’s are completely independent of those four parameters. By
comparing the above expression for the GW strain with that defined throughEqs. (5.1), (5.2), and (5.3), we findh1 (t) ∝ u(α, δ) cos[ ϕ(t; M, µ, α, δ, tc )] ,h2 (t) ∝ v(α, δ) cos[ ϕ(t; M, µ, α, δ, tc )] ,h3 (t) ∝ u(α, δ) sin[ ϕ(t; M, µ, α, δ, tc )] ,h4 (t) ∝ v(α, δ) sin[ ϕ(t; M, µ, α, δ, tc )] ,(5.8)

where the proportionality factor is [ G M/c2 ][(tc − t)/(5G M/c3 )]−1/4 . This method of resolving the GW strain signal in a basis of four time-varying functions
was first found in Ref. [56] for pulsar signals.
The new parameters, a(k) , with the index k taking four values, are defined in151terms of (r, ψ, ι, ϕc ) as,
1
1 + cos2 ι(1)a=− sin 2ψ sin ϕc cos ι ,cos 2ψ cos ϕc
r
2


1
1 + cos2 ι(2)+ cos 2ψ sin ϕc cos ι ,a=sin 2ψ cos ϕc
r
2


1
1 + cos2 ι(3)+ sin 2ψ cos ϕc cos ι ,a= −cos 2ψ sin ϕcr2
1 + cos2 ι1(4)
sin 2ψ sin ϕc− cos 2ψ cos ϕc cos ι .
a
= −
r
2

(5.9)

These constitute an alternative set of parameters that define the likelihood ratio. We used parenthetic indices above to avoid confusing them with numerical
exponents.

5.2.2

The network detection statistic

Let the inner-product of two temporal functions p(t) and q(t) be defined as

h p, qi( I ) = 4<

Z ∞
0

df

p̃∗ ( f ) q̃( f ),ξ(I)( f )(5.10)

where p̃( f ) and q̃( f ) are the Fourier transforms of p(t) and q(t), respectively,and ξ ( I ) ( f ) is the one-sided noise PSD of the Ith detector [26], with I = 1, ..., M
for a network of M detectors. The angular brackets denoting the inner-productare subscripted with the detector index since that product depends on the noise
PSD of the detector. Assuming that detector noise n I (t) is additive, the strain ina detector in the presence of a CBC signal isx I (t) = n I (t) + h I (t) ,(5.11)

where h I (t) is given by Eq. (5.3), but now with the antenna response functions
superscripted with the detector index. (The polarization components h+,× (t)152also depend on I through the coalescence time, as explained below.) Moreover,
if the noise is zero-mean Gaussian and stationary, the log-likelihood ratio (LLR)
is [26]1log Λ I = h x I , h I i( I ) − hh I , h I i( I ) ,2(5.12)

which can serve as a statistic for detecting signals in a single detector.
To explore the properties of the LLR, it will be useful to define the (complex)
unit-norm template S I (t) associated with the circular-polarization component of
a GW, namely, h+ (t) + ih× (t). It can be shown [47] thathi−1/4
S I (t) = g(−I1) τ I (tc − t)eiϕ(t) ,where g( I ) (with units of√(5.13)Hz) is a normalization factor, such that hS I , S I i = 1,

and
5
τ =
256 f sII

G M f sIc3

−5/3

(5.14)

is the time spent by the signal in the detector band, in the Newtonian approximation. Above, f sI is the seismic cut-off frequency of the Ith detector below
which it has little sensitivity for GW signals. The single detector matched-filter
output against S I (t) can then be defined as


I
I
I
C I = hS I , x I i ≡ c+
+ ic−
= ρ I eiφ ,(5.15)I , ρ I and φ I are all real; ρ I = |C I | is often termed as the signal-to-noise
where c±ratio (SNR) in the Ith detector. Since the detector strain due to a GW signal is
expected to be tiny, one has g( I )  1. Therefore, for computational efficiency,
we define a new factor that is closer to unity,
G M/c2
5G Mτ 1/4σ( I ) ≡g( I ) ,
1 Mpc
c3
153

(5.16)

with τ computed for a reference detector selected from one of those in the network. This is convenient since, as explained below, the detection statistics andthe parameters {ψ, ι, ϕc } are all independent of the above parenthetic scale factors; only the source distance depends on them, and is computed after accounting for them.
Using the strain expression in Eq. (5.7), the LLR for a network of multiple
detectors can be recast in terms of a(k) , provided one knows how the strain fromthe same CBC signal varies from one detector to the other. This was explainedin Refs. [46, 47]. Here, it suffices to note that this dependence arises owing to:(a) The spatial separation of the detectors, which can cause relative delays in thearrival of the signal. These delays are determined by the source’s sky-positionand can be accounted for in Eqs. (5.1) and (5.2) by adding those delays to tc .
(b) The different orientations of the detectors, which change u and v. Assumingthat the noise in the different detectors are statistically independent, the joint
log-likelihood ratio for a network of M detectors is
log



( M)

Λ



M

=

∑ log Λ I

I =1

1
= Nk a(k) − Mij a(i) a( j) ,2

(5.17)

where, in the last expression, the sum over detectors has been absorbed in Nk
and Mij , as defined below:






M
I
 N1 
 ∑ I =1 σ( I ) u I c+  uσ · c+ 








 M


I 
 N2 
 ∑ I =1 σ( I ) v I c+ vσ · c+ 
 = χ
 = χ
.


 M



I
 N3  ∑ uσ · c− 
 I =1 σ( I ) u I c− 








M
I
N4
vσ · c−∑ I =1 σ( I ) v I c−
154

(5.18)

Above, uσ and vσ are network vectors with components σ( I ) u I and σ( I ) v I , reI , andspectively, c± are network vectors with components c±

χ≡π

2/3



GM /c2
1Mpc

3/4
Mpc

is a normalization factor with dimensions of length. Also, A B 0 0  B C 0 0 


M=

 0 0 A B 
0 0 B C
with







(5.20)





 A 




2
 B  = χ  u ·v  , σ σ 






2
kvσ k
C

k u σ k2

(5.19)(5.21)

which define the network template-norm, namely, twice the second term on
the right-hand side of Eq. (5.17); the first term there can be interpreted as thematched-filter output of the network data-vector, x ≡ { x1 , x2 , ..., x M } [47].
Maximizing 2 log ( M) Λ with respect to a = { a(1) , a(2) , a(3) , a(4) } yields
2 log ( M) Λā

= N T · M −1 · N ,(5.22)

which is still a function of { M, µ, α, δ, tc }. (Note that the above statistic is independent of χ.) The concomitant maximum likelihood estimates (MLEs) of thecomplementary set of four parameters are denoted with an overline:ā = M−1 · N .
155(5.23)

These estimates are also functions of { M, µ, α, δ, tc }, and are determined by the

I as follows:data through c±




(
1
)
2
 ā

 kvσ k (uσ · c+ ) − (uσ · vσ ) (vσ · c+ )  (2) 2 (v · c ) 
−
u
·
v
u
·
c
+
k
u
k
(
)
(
)
 ā


σ
σ
σ
+
σ
σ
+ 

= χ
,
 (3) 


∆
2
 ā

 kvσ k (uσ · c− ) − (uσ · vσ ) (vσ · c− ) 








ā(4)
− ( u σ · v σ ) ( u σ · c − ) + k u σ k2 ( v σ · c − )(5.24)

where ∆ ≡ AC − B2 . The MLE of a parameter will be denoted by placing an
overline on its symbol.
It is important to note that the maximization in Eq. (5.22) assumes that the
network matrix M is invertible. This is not true, in general. Indeed, M is singular when uσ is aligned with vσ . These two vectors are determined by how
the interferometers in the network are oriented with respect to the wave propagation vector, but are not affected by the polarization angle ψ. In addition tothis singularity, M can be rank deficient, thus, making the problem of inverting it ill-posed [180]. Physically, this implies that the network does not haveenough linearly independent basis detectors to be able to resolve the source parameters a. Note that these maladies of M are dependent on the sky-position
angles. This means that a network that is able to resolve the signal parameters for certain source sky-positions may not be able to do so for others. These
problems can be tackled by regularizing M in a variety of ways that have been
explored in the context of searches of transient signals from unmodeled sources,
also called “burst” searches [180, 181, 182]. These methods obviate the rankdeficiency problem at the cost of making the search statistic sub-optimal. Thus,
any deficiencies arising from potential singularities in M or its regularization156

method adopted by a search pipeline will affect its performance. Since M isindependent of the detector strain data, such effects will arise in searches insimulated Gaussian data sets as well, such as the ones studied here. Since ourresults below are devoid of these maladies, we are confident that they will notarise in real data searches as well.
The maximum-likelihood estimates for the four physical parameters (r, ψ, ι, ϕc )can now be expressed in terms of the above estimates by inverting Eq. (5.9) andreplacing a with ā. Specifically, for the luminosity distance we get:√1 + 6 cos2 ῑ + cos4 ῑr̄ =,(5.25)2kākq
2where kāk ≡ ∑4i=1 ā(i) is the norm of the four-parameter vector MLE, andῑ is defined below along with the other MLEs. Since those angular parameter
estimates should not depend on an overall scaling of ā, it helps to define thedimensionless unit-norm components āˆ (k) ≡ ā(k) /kāk. In terms of the āˆ (k) , themaximum-likelihood estimates for the three angular parameters are, (1) āˆ (2) + āˆ (3) āˆ (4)
ˆ
2
ā
1
,
p
ψ̄ =
sin−1 
2
4
1−ζ

 
(1) āˆ (3) + āˆ (2) āˆ (4)
ˆ
2
ā
1
,
p
φ̄c = − sin−1 
2
1 − ζ2
!
√
2
1
−
1
−
κ
,
ῑ = cos−1
κ(5.26)




(
1
)
(
4
)
(
2
)
(
3
)
ˆ
ˆ
ˆ
ˆ
where ζ ≡ 2 ā ā − ā ā
and
κ=

ζ
p
.
1 + 1 − ζ2
157

(5.27)

Note that the expression for ψ̄ goes over to that of φ̄c under the transformationψ̄ −→ (−φ̄c )/2 and āˆ (2) ↔ āˆ (3) . This relation arises from a similar symmetry

exhibited by the a(k) defined in Eq. (5.9). Expressions for the CBC MLEs and thecoherent statistic were first obtained in Refs. [46, 47]. Above, we reexpress themin terms of the four parameters a(k) since the search code in LAL uses them [176].
Substituting for M and N, the MLR can be expanded as,
2 log Λā

= ( w + · c + )2 + ( w − · c + )2 + ( w + · c − )2 + ( w − · c − )2 ,I ,where w± are network vectors with components w±

 


 w I +   O11 O12   σ( I ) u I 

=

,
wI−
O21 O22σ( I ) v I(5.28)(5.29)

and






√

1  C + A + D/G1
 O11 O12 
 √

= √
2∆
O21 O22
C + A − D/G2

with D ≡p( A − C )2 + 4B2 and G1,2 ≡

q

√
√



C + A + D (C − A − D )/(2BG1 ) ,
C + A − D (C − A + D )/(2BG2 )(5.30)(C − A ∓ D )2 + 4B2 /(2B). The above

matrix diagonalizes M and, in so doing, identifies the dominant polarizationbasis, first identified in [46] and named as such in [181].
The coherent search statistic is just 2 log Λmaximized over { M, µ, α, δ, tc },ānamely,ρ2coh = 2 log Λϑ̄,(5.31)

where ϑ = { a(1) , a(2) , a(3) , a(4) , M, µ, α, δ, tc } is a set of nine parameters for the
non-spinning CBC signal. The last five parameters are searched for numerically,by using a grid for the masses and the sky-position and by using the fast Fourier
158

transform [183] to search for the coalescence time. ϑ̄ denotes the MLE values
Iof these parameters. Searching over (α, δ) requires the flexibility to delay c±Jby an interval that can be anywhere between zero and the lightrelative to c±

travel-time between the locations of the Ith and Jth detectors or the negative of
it. This is why we construct small snippets of C I (t) called C-data around theend-time of every trigger that is found to be coincident in multiple detectors in
a network. The statistic defined above will be termed as the coherent network
SNR and is the detection statistic optimal in stationary, Gaussian noise [47].
On the other hand, the combined signal-to-noise ratio, which was used asa detection statistic in the past and is used here in Fig. 5.2 for comparison, isdefined asM  2
ρ2comb = ∑ ρ I = kρk2 ,(5.32)

I =1

which is devoid of two significant pieces of information present in the coherent
search statistic in Eq. (5.31). The first piece of information is in the form of
the w I ± factors, which assign more weight to the matched-filter output of thedetector that is more sensitive to a given sky-position and has a lower noise
PSD (or bigger σ( I ) ). The second piece of information is in the form of thecross-detector terms that check for the consistency of the phases φ I with those
expected of a real signal.

5.2.3

Alternative statistics

The last several science runs at LIGO, GEO600, TAMA, and VIRGO have shown
time and again that real detector data is both non-stationary and non-Gaussian.
Consequently, neither the single-detector matched-filter-based SNR nor the co159

herent network SNR are optimal in that data. It is also known that empirically constructed search statistics perform better there. These alternative searchstatistics are based on signal discriminators such as the chi-square [174] and rhosquare tests [184], and their performances are compared against the statistics thatare optimal in Gaussian and stationary noise. These performances are evaluatedin terms of their receiver-operating characteristics, which in turn are constructedfrom detection efficiencies of simulated signals injected into network data andfrom the background rates obtained through multiple time-slide experiments.
The statistic that performs better in single-interferometric detector searches
is the matched-filtered output weighted by a function of the χ2 (or chi-square)statistic [171],ρeff ≡ ρ  χ2



2pχ2 − 2



1+


ρ2
ρ2o

−1/4


,

(5.33)

where ρ0 and pχ2 are empirical parameters that are deduced by examining theperformance of ρeff in real data. In the latest low-mass LIGO search, with 2M <M < 35M , they were chosen to be 250 and 16, respectively [173]. For the highmass (25M < M < 100M ) search studied below, these choices are 50 and 10,
respectively. Here, pχ2 is the number of degrees of freedom of the chi-squarestatistic, and ρo is chosen so that for small ρ and average chi-square values,ρeff ≈ ρ. A large chi-square value indicates that the disagreement between the
PSDs of the search template and the putative signal (or noise artifact) in the datais large, and imparts a greater penalty on ρeff by reducing its value relative to ρ.
The network equivalent of the effective SNR is
v
uM
u
2
( M)
ρI
ρ =t
eff

∑

I =1

160

eff

(5.34)

and is defined this way simply because it works in real data in discriminating
signal injections from background. A coherent statistic that can perform betterIin real data than its OGSN kin is constructed straightforwardly by replacing c±

with

I
I 
c±
eff ≡ c±χ2I
2pχ2 − 2



I

1+

2
ρI

(ρoI )2!−1/4,(5.35)

in Eq. (5.28). Since the ρ I and χ2I statistics are computed in the CBC search
pipeline when the data from the individual detectors are filtered, their values
are available to the coherent stage for computing the chi-square-weighted coherent
statistic defined above at little additional computational cost.
Scrutinizing expression (5.28) of ρcoh , one finds that it can be decomposed
into two parts. The first part isρ2auto−coh =

M 

∑

I =1

w2I + + w2I −



CI

2

(5.36)and is a sum of auto-correlation terms in each detector. This part of the coherent
statistic is less discriminatory between signal and noise triggers. The second
part,

M

ρ2cross−coh = ∑

M

h
i
I J
I J
w
w
+
w
w
c
c
+
c
c
(
)
∑ I+ J+ I− J− + + − − ,(5.37)

I =1 J =1( J 6= I )by contrast, is a sum of cross-correlation terms across pairs of detectors, or baselines, and is critical in checking for phase consistency among signals appearing
in the detectors from a GW source. Once again, both of the above statistics can
I with c I
be made more robust against noise glitches by replacing c±± eff to obtain

their chi-square-weighted counterparts.
Another statistic that is helpful in discriminating signals from noise glitches161in multi-detector data is the null-stream [51]. If C̃ I ( f ) is the Fourier transformof C I (t), then one can show that for GW signals in the data, the mean ofM(I)Y ≡ ∑ K I σinv Sh( I ) ( f )C̃ I ( f )(5.38)

I =1

is zero. Above, K I = e I JK F+J F×K , with e I JK being the Levi-Civita symbol, and(I)

σinv ≡ (σ( I ) )−1 . For non-stationary artifacts, however, this need not be true,
thereby, motivating the following discriminator:
η= p

h|Y |i,Var (|Y |)(5.39)

where h x i and Var( x ) denote the statistical average and variance of x, respectively. The above construct is called the null-stream statistic. Just like thecoherent SNR, it can be decomposed into two parts as well, comprising autocorrelation and cross-correlation terms, respectively. The former is akin to theincoherent energy defined in Ref. [52] for burst searches and will be denotedas ηauto . For GW signals one expects ηauto to be large while η itself is small.
On the other hand, for noise artifacts, η is expected to be large, on the average,
even when ηauto itself is not very strong. This analysis argues for a new statistic,
namely,R = ηauto /η ,(5.40)

which we call the ratio-statistic.

5.3

Coherent hierarchical inspiral analysis pipeline

The coherent hierarchical inspiral analysis (CHIA) pipeline mainly comprises
two stages, namely, the coincident and coherent stages, respectively. Both in162

volve multiple steps. The coincident stage has been discussed in the past in
Refs. [171, 172] and is described here briefly for completeness. It includes thefollowing steps: (a) Compute noise PSDs and generate template-banks of thetwo component masses for each detector in the network. The noise PSDs varyfrom one detector to another, and in time. A template bank is constructed forevery 2048s chunk of data from every detector [185]. (b) Use the template bankfor each detector to filter the data from that detector and output the parametersof triggers crossing the chosen SNR threshold. For the injection studies, simulated software-injections are added in software to the data in this step, before thedata are match-filtered. (c) Parameters of the triggers from the participating detectors are then compared to identify coincidences [186]. Before these coincident
triggers are considered as detection candidates, in real data one usually applies
data-quality vetoes. For our study in simulated data, we forego this stage of thepipeline and, instead, apply the coherent stage directly to the triple-coincident
triggers. For the computation of the coherent and null-stream statistics the Cdata time-series, which include both the amplitude and the phase time-series ofthe matched-filter outputs, are required. These time-series are computed in thecoherent stage and not upstream in the pipeline since it is computationally lessexpensive to identify coincidences and construct the C-data only for them.
The coherent stage in the CBC search pipeline is constituted of 4 steps. In
the first step, a “coherent bank” of templates is constructed from the parametersof the coincident triggers. Triggers in different detectors that are coincident andarise from the same GW source can have different mass pairs owing to the possibility that the noise PSDs of the detectors they arise in are somewhat different

163and because of the random nature of noise. For every coincident trigger we construct a network template with a single mass-pair, namely, the one corresponding to the loudest SNR among all the detectors, to search coherently around theend-time of that putative signal. This mass-pair will be termed as the max-SNR
pair and the corresponding detector the max-SNR detector. For example, consider a triple-coincident trigger with {ρ, m1 /M , m2 /M } = {10.0, 1.43, 1.39},{10.9, 1.40, 1.36}, and {8.9, 1.51, 1.32} in the first, second, and third interferometric detector (or IFO), respectively. Then the max-SNR detector is IFO-2 and themax-SNR mass-pair is {m1 /M , m2 /M } = {1.40, 1.36}, which is the templateincluded in the coherent bank to represent this coincident trigger in the coherent
stage.
While this mass pair will not necessarily give the loudest SNR in the twoother detectors, it has been found to yield a better performance for the coherentstatistic and null-stream than when they are computed using the original and,often, non-identical mass pairs in the different detectors. (Note, however, thatsimulated software injections in real data must be used to determine empiricallyif the detection efficiency is helped by using the same mass pair across all detectors in any given science run.) Also, since error-covariances are known to exist
between the mass parameters and the trigger end-time, we search at and aroundthe end-times of the single-detector triggers that constitute a given network trigger.
The second step in the coherent stage is the construction of trigger-banks,
whereby the coherent-bank template for every coincident trigger is copied asa single-detector template. (See Fig. 5.1.) In the subsequent step, the single-

164detector templates are used to filter the data from the individual IFOs. This
step outputs the time-series of C-data around the trigger end-times in that detector. Additionally, this step computes the template normalization factor andchi-square for the max-SNR mass-pair across all detectors per coincident trigger.
Note that the values of these constructs are not available earlier in the pipeline
for the triggers in the detectors complementary to the max-SNR detector since, in
general, the mass-pairs would be somewhat different in the preceding coincident
stage of the search pipeline. In summary, this step outputs a C-data time-series
and the corresponding signal parameters, such as the template normalization
factor, for every trigger listed in the coherent-bank output file.
The final step of the coherent stage is the coherent-statistics step, which matches
the parameters of each triple-coincident trigger to the C-data time-series output
by the matched-filtering step and uses them and the corresponding templatenorms, chi-square values for the respective detectors to compute a variety ofmulti-detector statistics, such as the coherent SNR, null-stream, the chi-squareweighted coherent SNR, and other alternative statistics.

5.4

Results

To study the performance gain arising from using the coherent stage, we ran
the CBC search pipeline with and without that stage on simulated Gaussian
noise, with LIGO-I noise PSD [176] in the 4km LIGO detectors in Hanford (H1),
Livingston (L1), and in the VIRGO detector (V1), for the duration of approximately a month. Specifically, this search pipeline was run once with signal in-

165

Coincident stage

Coherent bankIFO­1
trigger bank

IFO­M
trigger bank

IFO­1matchedfiltering

IFO­M
matchedfiltering

Coherent statistics

Figure 5.1: A schematic diagram of the coherent stage in the compact binary
coalescence search pipeline.
jections and again (parallelly) without injections but with time-slid data so thatthe background could be estimated. The top plot in Fig. 5.2 compares the performance of the coherent statistics and the combined effective SNR. The bottom
plot compares the coherent SNR and null-stream statistics. For these simulations, 1051 signals were injected in software in all three detectors. The source
distances of all injections were between 100-500 Mpc. The total masses of thesesources were chosen to be in the range 25-100 M , and component masses between 1-99 M . A total of 55 of those injections were found, above the single-

166

interferometer detection thresholds of 5.0 and coherent SNR threshold of 3.751 .
The latter threshold was intentionally chosen to be lower since we anticipatedthat some coincident background triggers will have negative cross-terms owing
to incoherent phases, thereby, yielding lower coherent SNRs.
All injections recovered by the coincident stage were also found by the coherent stage, and are symbolized by red pluses. The black crosses depict the background triggers that are found by the coincident stage and survive the coherent
stage. The blue circles, on the other hand, denote background triggers in the coincident stage that got vetoed by the choice of the threshold on the coherent SNR
in the coherent stage. To include them in the top plot, we arbitrarily assign all ofthem ρcoh = 3.0. Comparing the sets of black crosses and blue circles reveals thatthe coherent stage not only reduces the number of background triggers but, in
this case, also vetoes some of the loudest ones (in combined-effective SNR). Furthermore, whereas all found injections have coherent SNR greater than that of
the loudest background trigger, 13 of them have combined-effective-SNR weaker
than that of the loudest background trigger (shown in blue circles). When compared to the loudest black cross, that number drops to 7. It drops further whensome of the background triggers with the loudest null-stream (as shown in thebottom plot) are vetoed. The resulting performance improvement is depicted in
the blue dash-dotted Receiver-Operating-Characteristic (ROC) curve in Fig. 5.3;
1 The detection probabilities are small because, first, all injections made were weak and, sec-

ond, here we focused only on triggers that are coincident in all three detectors. Owing tosensitivity disparities, it is more likely to find injection trigger coincidences in two of the threedetectors. Only weak injections were made since that is where the coherent code can help improve the performance of current searches.

167

its performance is better than that of the coincident stage (shown in red), without the null-stream vetoes. The former is asymptotic to the ROC curve of thecoherent stage (shown in black dashes) for higher false-alarm probabilities.
Finally, Fig. 5.2 reveals the existence of a gap between the loudest backgroundand the weakest injection ρcoh values. One might argue that this is owing to thelack of a sufficient number of weak signal injections made into the data. We haveverified that, indeed, one can get some injection triggers to show up in that gap
by making multiple weak injections (say, with source distances between 500-750Mpc) in the data. Those studies also reveal that the detection efficiency in that
region is very low (i.e., less than 1 in 250). We believe that this low efficiency
is partly caused by the coincident stage, in the way it has been designed andtuned, acting as a bottleneck for the coherent stage.

5.5

Discussion

The main advantage of implementing a blind coherent search in the hierarchical manner explained above is that it has a lower computational cost comparedto that of a fully coherent search pipeline. This is primarily because it reduces
the number of time-of-arrival values for the coherent code to search for, andbecause recognizing coincidences is relatively cheaper computationally. There
are additional reasons, such the inherent detector-bound nature of data-quality
cuts, which are best implemented in the matched-filtering stage. This in turncan reduce an otherwise triple-coincident trigger into a double-coincident one if
the third IFO data-points around the concurrent time get vetoed. Since the co-168incident and coherent statistics are the same for two-site CBC searches, it makes
sense to not follow them up with the coherent stage.
There are, however, some demerits of searching hierarchically. The first one
came to the fore in the results presented above, where the coincident stage ispotentially affecting the efficiency of the coherent stage in finding injections. Indeed, it may be possible to improve the injection finding efficiency by reducing
the SNR thresholds in the matched-filtering step of the coherent stage. Whilethat may happen, it is also likely that the overall performance of the pipeline
will be hurt since it will tend to increase the background rate as well. An alternative solution is to retain the original mass-pairs of the coincident triggers
in the coherent stage instead of replacing them with max-SNR mass-pairs. This
will ensure that injection-finding efficiency of the matched-filtering stage is unaffected, but may hurt the coherence of the triggers and, therefore, ultimately
affect the injection finding efficiency of the coherent stage. It may also cause thefalse-alarm rate to rise, owing to the less stringent requirements on the agreement of the mass-pair values across the network of detectors.
A more optimal solution that addresses the drawbacks of the last two solutions is to assign to every coincident trigger multiple mass-pair templates tosearch the data with in the coherent stage. This approach makes sense since statistical errors alone are known to cause substantially different mass-templates
to be triggered by signals in different detectors arising from the same (injected
simulated) source. However, as was shown by the work in Ref. [186] on identifying coincidences, the separation in the mass parameter-space between triggers in two detectors from the same source can be wide enough to allow for

169

multiple other mass templates to fit in between. Some of these intermediate
mass-templates can have a greater chance of not only passing the SNR thresholdin individual detectors but also appearing as coherent. The main problem toattack here is to find what the optimal density and size are of these relatively
small template banks localized around the coincident mass-pairs. Too small adensity or size can hurt signal-finding efficiency and too large a density or size
can increase the background rate.

170

Coherent SNR vs null−stream
14

Coherent null−stream

12

Background
Injections

10
8
6
4
2
0
0

10

20
30
Coherent SNR

40

50

Figure 5.2: These are scatter plots of the combined and coherent SNRs of injection triggers, represented by red plus symbols, and background (or “slide”)
triggers, represented by the black crosses. The coherent SNR was used to cluster
the triggers, from both injections and slides. The coherent SNR performs noticeably better than the combined effective SNR in discriminating signals frombackground: In the top plot, at a detection threshold of a little above 6 in thecoherent SNR all the injections found in the coincident stage are recovered witha vanishing false-alarm probability. For the same false-alarm probability, thecombined effective SNR detects a lesser number of injected signals.
171

ROC Curves
0.024

Coincident stage
Coincident stage with null−stream veto
Coincident followed by coherent stage

Detection probability

0.022

0.02

0.018

0.016

0.014

0.012
0

0.2

0.4
0.6
False−alarm probability

0.8

1Figure 5.3: The receiver operating characteristic (ROC) curves of three CBC
searches are compared above. The ROC of the search with the coincident stage
alone is plotted in solid red line, and has the weakest performance owing to the
13 found injections that are weaker than the loudest background trigger in thatsearch. On the other hand, the ROC curve for the hierarchical pipeline, withcoherent stage included, is shown in black dash-dotted line and has the best
performance. It has a constant detection probability because all found injections
are louder than the loudest background trigger for this pipeline. Finally, thethird ROC curve, shown as a blue dashed line is the coincident stage, with thenull-stream veto applied. This veto improves the performance of the coincident
pipeline, so much so that for low detection-thresholds (or high false-alarm probability) its ROC curve rises to match that of the pipeline with the coherent stage.
The average error in the detection probabilities plotted here is less than 3 × 10−4 .

172Appendix A
Technical details
A.1

Antenna response

The gravitational wave strain is the linear combination of the strains induced byplus and cross polarizations. The antenna response functions F+ and F× encode
the detector responses to the plus and cross polarizations defined in the radiationbasis. Their functional forms are given by [187]1
F+ (θ, φ, ψ) = − (1 + cos2 θ ) cos 2φ cos 2ψ − cos θ sin 2φ sin 2ψ ,2
1
(1 + cos2 θ ) cos 2φ sin 2ψ − cos θ sin 2φ cos 2ψ ,F× (θ, φ, ψ) =
2

(A.1)
(A.2)

where θ and φ are the spherical polar coordinates of the source in the detector
frame, and ψ is the polarization angle defined in Fig. A.1. Since the polarization
angle ψ is unknown a priori, it is often useful to use the quadrature sum ofantenna response functions, also called antenna-pattern:
Fsum (θ, φ) =q

2 ( θ, φ, ψ ) + F2 ( θ, φ, ψ ) ,
F+
×

173

(A.3)

which is independent of ψ. The antenna response functions and the antennapattern are shown in Fig. A.2.

y'source
frame

x'
ι
y''

ψ

Z
z''

x''

z'
radiation
frame

θ
Y
detector
frame

ϕ
X

Figure A.1: Definition of source location and orientation angles relative to theinertial frame of a GW detector with arms along X and Y axes. Wave propagatesalong z00 -axis from a source at an arbitrary sky position defined by (θ, φ) andsubtending an inclination angle ι. The polarization angle ψ, subtended by thewave, is also shown here.

174

Figure A.2: The antenna response of interferometric gravitational wave detectors
to gravitational waves with plus (top left), and cross (top right) polarizations.
For making these plots, we have chosen ψ = 0. The bottom plot depicts theantenna-pattern (for an unpolarized wave.)

175

A.2

Parameter accuracy

The match can be rewritten as
M = 1 − gαβ ∆Θα(k) ∆Θ(k) ,
β

(A.4)

where
Γ(k)αβgαβ= (SNR(k) )2 gαβ (Θ(k) ) ,


2
∂ M
1

= gαβ (Θ(k) ) ,= − 
2 ∂Θα 0 ∂Θ β 0
Θ(k0 ) = Θ(k)
(k )
(k )


 gµ(k) µ(k) gµ(k) φ(k) := .
gφ(k) µ(k) gφ(k) φ(k)

(A.5)

(A.6)

The components of the above gαβ matrix are obtained from the derivatives of thebeam matrix:gµ (k ) µ (k ) =1"(B(k)(k) )2(B(k)(k) )2
2 #

∂
B0
,
−
∂µk0 k k k0 =k∂2
B0∂µk0 ∂µk k k k0 =k

gµ(k) φ(k) = gφ(k) µ(k) =
"


1
∂2
(B(k)(k) )B0∂µk0 ∂φk k k k0 =k
2(B(k)(k) )2


#
∂
∂
−
B0
B0
,
∂µk0 k k k0 =k∂φk0 k k k0 =k

176

(A.7)

(A.8)

gφ(k) φ(k) =1"(B(k)(k) )2(B(k)(k) )2

2 #
∂
−
.
B0
∂φk0 k k k0 =k∂2
B0∂φk0 ∂φk k k k0 =k



(A.9)

The estimation error (2.42) is obtained from Eqs. (A.5) and (2.41) by utilizing thefact that the inverse of the determinant of a matrix is the same as the determinantof the inverse of that matrix. Above, we assumed that B is nonsingular.

177

Appendix B
Advanced detectors studies
B.1

Performance of multi-baseline radiometers

We assess the improvement in performance of a LIGO-VIRGO network gainedby adding detectors in Japan and India, namely, KAGRA and INDIGO, respectively. The source PSD is taken to be a constant, H ( f ) = 1.516 × 10−48 /Hz. The
signal band considered here is 40-1024 Hz. For the second-generation configuration, the noise PSDs of LIGO, INDIGO, and KAGRA are taken to be the aLIGO
noise PSD. The noise PSD of VIRGO is chosen to be the aVIRGO one. And forthe third-generation configuration, the noise PSDs of all detectors are taken to
be that of ET, as shown in Figure 2.6.

178

4000

4000
3500
Area−weighted sensitivity

Sensitivity

3500

3000

2500
HLV
HLVI
HLVK
HLVIK

2000

1500

−50

3000
2500
2000
1500
1000

HLV
HLVI
HLVK
HLVIK

500
0
Declination (degrees)

0

50

−50

0
Declination (degrees)

500Declination (degrees)

50

5

5

6.5

7

x 10

6
Area−weighted sensitivity

6

Sensitivity

5.5
5
4.5
4
3.5
3

x 10

HLV
HLVI
HLVK
HLVIK
−50

5
4
3
2
1

0
Declination (degrees)

0

50

HLV
HLVI
HLVK
HLVIK
−50

Figure B.1: The sensitivities (left) and their area-weighted counterparts (right)of different LIGO-VIRGO-INDIGO-KAGRA networks in their second (top) andthird (bottom) generation configurations are plotted as functions of the declination of a single-pixel SGWB source.

179

Figure B.2: The beam functions of different second-generation LIGO-VIRGOINDIGO-KAGRA baselines for a source at declination +0.25 radian and with aconstant PSD are shown here.

180

7

Singular value

10

6

10

5

HLV
HLVI
HLVK
HLVIK

10

0

10

1

2

10

10
Index

Figure B.3: Singular values of the Fisher matrices for different second-generation
LIGO-VIRGO-INDIGO-KAGRA networks are plotted in this figure.

181

−8

−8

10
Area−weighted estimation error (sr)Estimation error (sr)

10

−9

10

−10

10

−11

10

HLV
HLVI
HLVK
HLVIK
−50

10

50

−50

0
Declination (degrees)

500Declination (degrees)

50

10Area−weighted estimation error (sr)Estimation error (sr)

HLV
HLVI
HLVK
HLVIK

−12

10

−13

10

−14

10

−15

−10

10

−11

0
Declination (degrees)

−12

10

−9

10

HLV
HLVI
HLVK
HLVIK
−50

−13

10

−14

10

−15

10

−16

0
Declination (degrees)

10

50

HLV
HLVI
HLVK
HLVIK
−50

Figure B.4: The 1σ error (left) and area-weighted 1σ error (right) in the solid angle of locating a source in the sky with different LIGO-VIRGO-INDIGO-KAGRA
networks in their second (top) and third (bottom) generation configurations.

182

Appendix C
A directed SGWB search and Olbers’
paradoxC.1

Virgo cluster versus the isotropic backgroundWe explicitly compute the ratio between the effective number of sources in the
nearby Virgo cluster to that in the distant universe in the same solid angle subtended by the Virgo cluster. This is taken from the Ref. [75].
Consider the flat (k = 0) model of the universe which is currently dust dominated. Then the scale factor is given bya(t) = α t2/3 .

(C.1)

The Hubble parameter H := ȧ(t)/a(t) at the current epoch t = t0 from Eq. (C.1)
becomes
H0 := H (t0 ) =ȧ(t0 )2=,a ( t0 )
3t0

183

(C.2)

which gives t0 = (2/3) H0−1 . Setting the present scale factor to 1, i.e., a(t0 ) = 1,
one gets
a(t) = α t2/3

 2/3 t
3H0 t 2/3
= a ( t0 )
=
.
t0
2

(C.3)

The cosmological redshift z, defined as 1 + z := a(t0 )/a(t) = 1/a(t), is avery useful measure of distance. From the Friedmann-Robertson-Walker metricfor flat space-time, the proper radius r, from which light emitted at the epoch treaches us at the epoch t0 , can be written as
t0
dt
cdt
=
,
t ( t/t0 )2/3
t a(t)"
 1/3 #qt2c= 3c t0 1 −=1 − a(t) .
t0
H0

r =

Z t0

Z

(C.4)
(C.5)

Then, introducing the definition of redshift, one can write,


2c
1
r=
1− √
.
H0
1+z

(C.6)

c
(1 + z)−3/2 dz .
H0

(C.7)

Therefore,
dr =Since there are n(1 + z)3 galaxies per unit volume at a redshift z, the effective
number of sources in a solid angle ∆Ω can be expressed as
Z ∞n(1 + z)3 r2 ∆Ω dr,
r 2 (1 + z )5
0
Z ∞
dz
2c
c
= n∆Ω
= n∆Ω
,
7/2
H0 0 (1 + z)
5H0

∆neff =

(C.8)
(C.9)

where n is the average number density of sources in a homogeneous universe,the (1 + z)5 factor signifies the decrease in radiation energy density by a factorof (1 + z)4 and increase in time interval by a factor of 1 + z (which reduces the184effective luminosity of the source). Assuming cH0−1 = 3 × 103 h−1 Mpc, whereh ∼ 0.72, and n ∼ 0.004 Mpc−3 , the effective number of galaxies at 1 Mpc
distance per unit solid angle becomes ∼ 7.
The number of sources (galaxies) at the Virgo cluster is nVirgo ∼ 1500, spreadover a solid angle of ∆ΩVirgo ≈ 0.015 (that is 50 square degrees). The Virgocluster is rVirgo ∼ 16.5 Mpc away. The redshift to the Virgo cluster is negligible.
Thus, for the Virgo cluster the effective number of sources at 1 Mpc distance perunit solid angle is
Virgo

∆neff

=

nVirgo∼ 360 .
2∆ΩVirgo rVirgo(C.10)

This is of course a very rough estimate. However, the Virgo cluster is about 50
times brighter than the background created by the distant sources.

185

Appendix D
Notes on SNR threshold
D.1

SNR thresholds for the S4 ringdown search

In the low-mass CBC search, the two orthogonal phases of the binary inspiralwaveform are expressed as [185]2µhc ( t ) = 2[πGM f (t)]2/3 cos[2φ(t) − 2φ0 ] ,
M
c


2
µ
hs ( t ) = 2[πGM f (t)]2/3 sin[2φ(t) − 2φ0 ] .
M
c

(D.1)
(D.2)

The normalized filter-output is found by filtering the data x (t) against the elliptically-polarized template hep (t) = hc (t) − ihs (t) and is expressed as|h x, hep i||h x, hc i + i h x, hs i|=p,ρep (hep ) ≡ q
hhc − ihs , hc − ihs ihhep , hep iph x, hc i2 + h x, hs i2= p,
h hc , hc i + h hs , hs ip1
h x, hc i2 + h x, hs i2
p
= √
,
2
h hc , hc i1
≡ √ ρD (hep ) .
2
186

(D.3)

(D.4)
(D.5)
(D.6)

The quantity ρD defined above is termed as SNR in the low-mass CBC search.
Any triggers having ρD ≥ 5.5 are kept, which is equivalent to having ρep ≥ 3.9.
Contrastingly, for the ringdown search the normalized filter-output is given
by Eq. (4.27) as|h x, hc i|
ρc ( hc ) = p.
h hc , hc i

(D.7)

The quantity ρc is termed as SNR in the ringdown search. That means thequantity ρc in the ringdown search is equivalent to the quantity ρep in the lowmass CBC search.

187

Appendix E
List of abbreviations

AGWB: Astrophysical Gravitational-Wave BackgroundaLIGO: Advanced Laser Interferometer Gravitational-Wave Observatory, the second generation LIGO detectors
ASD: Amplitude Spectral Density
aVIRGO: Advanced VIRGO, a second generation European GW observatory
BH: Black Hole
CBC: Compact Binary Coalescence
CG: Conjugate Gradient
CHIA: Coherent Hierarchical Inspiral Analysis
CMB: Cosmic Microwave BackgroundcWB: Coherent WaveBurst
Dec: Declination
EFE: Einstein Field Equations

188

EOB: Effective-One-Body
EP: Equatorial Plane
ET: Einstein Telescope, a third generation gravitational-wave observatory
FAR: False Alarm Rate
GEO: British-German 600 m long gravitational-wave observatory
GR: General Relativity
GRB: Gamma-Ray Burst
GW: Gravitational Wave
HEALPix: Hierarchical Equal Area iso-Latitude Pixelization
H1: The LIGO detector with 4 km arm-length at Hanford
H1L1: The baseline constructed from H1 and L1
H1L1V1: The triple-baseline network constructed from H1, L1 and V1
H1V1: The baseline constructed from H1 and V1
IFO: Interferometric detector
IMBH: Intermediate-Mass Black Hole
IMR: Inspiral-Merger-Ringdown
INDIGO: Proposed Indian Gravitational-Wave Observatory
KAGRA: Japanese second generation 3 km long cryogenic gravitational-wave
detector (previously known as LCGT for Large-scale Cryogenic Gravitationalwave Telescope)
LIGO: Laser Interferometer Gravitational-wave Observatory
LLR: Log-Likelihood Ratio
LMXB: Low-Mass X-ray Binary
L1: The LIGO detector with 4 km arm-length at Livingston

189

LSC: LIGO Scientific Collaboration
ML: Maximum Likelihood
MLE: Maximum Likelihood Estimate
MLR: Maximized-Likelihood Ratio
MSP: Millisecond Pulsar
NMSE: Normalized Mean Square Error
NS: Neutron Star
OGSN: Optimal in Gaussian and Stationary Noise
PSD: Power Spectral Density
PSR: Pulsar
RA: Right Ascension
ROC: Receiver Operating Characteristic
S5: The fifth LIGO science run
S4: The fourth LIGO science run
SGWB: Stochastic Gravitational-Wave Background
SNR: Signal-to-Noise Ratio
S6: The sixth LIGO science run
TT: Transverse-Traceless
VIRGO: European gravitational-wave observatory located at Cascina, Italy
V1: The “European gravitational-wave observatory" with 3 km arm-length at
Cascina
VSR: VIRGO Science Run

190

Bibliography[1] A. Einstein, Annalen der Physik, 354, 769 (1996).
[2] R. A. Hulse, and J. H. Taylor, Astrophys. J. 195, L51 (1975).
[3] B. F. Schutz, A first course in general relativity (Cambridge University Press,
1985).
[4] P. R. Saulson, Fundamentals of Interferometric Gravitational Wave Detectors(World Scientific, 1994).
[5] http://www.ligo.caltech.edu/.
[6] https://wwwcascina.virgo.infn.it/.
[7] http://gwcenter.icrr.u-tokyo.ac.jp/en/.
[8] http://www.geo600.org/.
[9] D. Sigg, Internal working note, LIGO-P980007-00-D.
[10] B. F Schutz, Class. Quant. Grav. 28, 125023 (2011).
[11] C. Cutler, and K. S. Thorne, An Overview of Gravitational-Wave Sources[arXiv:gr-qc/0204090].
191

[12] S. Bose, T. Dayanga, S. Ghosh, and D. Talukder, Class. Quant. Grav. 28,
134009 (2011).
[13] D. Talukder, S. Mitra, and S. Bose, Phys. Rev. D 83, 063002 (2011).
[14] E. Thrane, S. Ballmer, J. Romano, S. Mitra, D. Talukder, S. Bose, and V.
Mandic, Phys. Rev. D 80, 122002 (2009).
[15] A. A. Penzias, and R. W. Wilson, Astrophys. J. 142, 419 (1965).
[16] R. H. Dicke, P. J. E. Peebles, P. G. Roll, and D. T. Wilkinson, Astrophys. J.
142, 414 (1965).
[17] B. Allen, The stochastic gravity-wave background:sources and detection[arXiv:gr-qc/9604033].
[18] B. P. Abbott et al. (LIGO Scientific), Nature 460, 990 (2009), 0910.5772.
[19] A. A. Starobinsky, JETP Lett. 30, 682 (1979).
[20] T. Kahniashvili, A. Kosowsky, G. Gogoberidze, and Y. Maravin, Detectabilityof Gravitational Waves from Phase Transitions [arXiv:astro-ph/0806.0293].
[21] D. Coward, and T. Regimbau, New Astron. Rev. 50, 461 (2006).
[22] B. F. Burke, and F. Graham-Smith, An Introduction to Radio Astronomy (Cambridge University Press, 1997).
[23] R. B. Partridge, 3K: The Cosmic Microwave Background Radiation (Cambridge
University Press, 1995).

192

[24] S. Mitra, S. Dhurandhar, T. Souradeep, A. Lazzarini, V. Mandic, S. Bose,and S. Ballmer, Phys. Rev. D 77, 042002 (2008), arXiv:0708.2728.
[25] S. W. Ballmer, Class. Quant. Grav. 23, S179 (2006), gr-qc/0510096.
[26] C. W. Helstrom, Statistical Theory of Signal Detection (Pergamon Press, London, England, 1968), 2nd ed.
[27] B. Allen and A. C. Ottewill, Phys. Rev. D 56, 545 (1997).
[28] N. J. Cornish, Class. Quant. Grav. 18, 4277 (2001).
[29] H. Kudoh and A. Taruya, Phys. Rev. D71, 024025 (2005), gr-qc/0411017.
[30] A. Taruya and H. Kudoh, Phys. Rev. D72, 104015 (2005), gr-qc/0507114.
[31] A. Taruya, Phys. Rev. D74, 104022 (2006), gr-qc/0607080.
[32] Kipp C. Cannon, Phys. Rev. D75, 123003 (2007).
[33] A. Lazzarini and R. Weiss, Internal working note LIGO-T040140-00-Z, Laser
Interferometer Gravitational Wave Observatory (LIGO) (2004).
[34] B. Abbott et al., Phys. Rev. D 76, 082003 (2007), astro-ph/0703234.
[35] E. Flanagan, Phys. Rev. D 48, 2389 (1993).
[36] B. Allen, and J. D. Romano, Phys. Rev. D 59, 102001 (1999).
[37] Hierarchical Equal Area isoLatitude Pixelization (HEALPix),
http://healpix.jpl.nasa.gov/.
[38] K. M. Górski et al., Astrophys. J. 622, 759 (2005), astro-ph/0409513.
193

[39] B. P. Abbott et al. (LIGO Scientific), Rep. Prog. Phys. 72, 076901 (2009),
0711.3041.
[40] F. Acernese et al., Optics and Lasers in Engineering 45, 478 (2007).
[41] P. Ajith, and S. Bose, Phys. Rev. D 79, 084032 (2009).
[42] G. Cella, C. N. Colacino, E. Cuoco, A. Di Virgilio, T. Regimbau, E. L. Robinson, and J. T. Whelan, Class. Quant. Grav. 24, S639 (2007), 0704.2983.
[43] B. J. Owen, Phys. Rev. D 53, 6749 (1996), arXiv:gr-qc/9511032.
[44] C. Cutler, Phys. Rev. D 57, 7089 (1998), arXiv:gr-qc/9703068.
[45] G. Hinshaw, M. R. Nolta, C. L. Bennett, R. Bean, O. Doré, M. R. Greason,
M. Halpern, R. S. Hill, N. Jarosik, A. Kogut, et al., Astrophys. J. Suppl 170,
288 (2007), arXiv:astro-ph/0603451.
[46] S. Bose, A. Pai, and S. Dhurandhar, International Journal of Modern Physics
D 9, 325 (2000), arXiv:gr-qc/0002010.
[47] A. Pai, S. Dhurandhar, and S. Bose, Phys. Rev. D 64, 042004 (2001),
arXiv:gr-qc/0009078.
[48] L. Wen and Y. Chen, Phys. Rev. D 81, 082001 (2010), 1003.2504.
[49] K. Hayama, S. D. Mohanty, M. Rakhmanov, and S. Desai, Class. Quant.
Grav. 24, S681 (2007), 0709.0940.
[50] S. Dhurandhar, B. Krishnan, H. Mukhopadhyay, and J. T. Whelan, Phys.
Rev. D 77, 082001 (2008), 0712.1578.
194

[51] Y. Gürsel and M. Tinto, Phys. Rev. D 40, 3884 (1989).
[52] S. Chatterji, A. Lazzarini, L. Stein, P. J. Sutton, A. Searle, and M. Tinto, Phys.
Rev. D 74, 082005 (2006), arXiv:gr-qc/0605002.
[53] S. Bose, Phys. Rev. D71, 082001 (2005), astro-ph/0504048.
[54] S. Bose et al., Class. Quant. Grav. 20, S677-S687 (2003).
[55] S. Dhurandhar, H. Tagoshi, Y. Okada, N. Kanda, and H. Takahashi, Phys.
Rev. D 84, 083007 (2011).
[56] P. Jaranowski, A. Królak, and B. F. Schutz, Phys. Rev. D 58, 063001 (1998).
[57] N. Sartore, E. Ripamonti, A. Treves, and R. Turolla, Astron. Astrophys. 510,
A23 (2010).
[58] D. R. Lorimer et al., Astrophys. J. 439, 933 (1995).
[59] http://www.atnf.csiro.au/people/pulsar/psrcat/.
[60] J. Abadie et al. (LIGO Scientific), Phys. Rev. Lett. 107, 271102 (2011).
[61] D. R. Lorimer et al., Mon. Not. R. Astron. Soc. 372, 777 (2006).
[62] D. R. Lorimer, Binary and Millisecond Pulsars, Living Rev. Relativity 11, 8(2008).
[63] C. Palomba, Mon. Not. R. Astron. Soc. 359, 1050 (2005).
[64] B. Knispel, and B. Allen, Phys. Rev. D 78, 044031 (2008).

195

[65] B. Knispel (2011), Ph.D. Thesis, Albert Einstein Institute, Hannover,
http://edok01.tib.uni-hannover.de/edoks/e01dh11/66410357X.pdf.
[66] Z. Arzoumanian, D. F. Chernoff, and J. M. Cordes, Astrophys. J. 568, 289(2002).
[67] S. A. Story, P. L. Gonthier, and A. K. Harding, Astrophys. J. 671, 713 (2007).
[68] P. L. Gonthier, C. Billman, A. K. Harding, I. A. Grenier, and M. Pierbattista, Population Synthesis of Radio and Gamma-ray Pulsars in the Fermi Era[arXiv:astro-ph/1103.1856].
[69] T. Regimbau, and J. A. de Freitas Pacheco, Astron. Astrophys. 359, 242(2000).
[70] T. Regimbau, and J. A. de Freitas Pacheco, Astron. Astrophys. 374, 182(2001).
[71] J. P. Ostrikder, and J. E. Gunn, Astrophys. J. 157, 1395 (1969).
[72] T. Regimbau (Private communication, 2012).
[73] T. Regimbau, Internal working note, LIGO-VIRGO Scientific Collaboration(2012).
[74] S. Ballmer, Internal working note, LIGO Scientific Collaboration (2012).
[75] N. Mazumder, S. Mitra, and S. Dhurandhar, Internal working note, LIGO
Scientific Collaboration (2012).

196

[76] R. E. E. Rothschild, E. A. Boldt, S. S. Holt, and P. J. Serlemitsos, Astrophys.
J. 189, L13 (1974).
[77] H. V. D. Bradt, and J. E. McClintock, Ann. Rev. Astron. Astrophys. 21, 13(1983).
[78] J. A. Orosz et al., Nature 449, 872 (2007).
[79] T. Bulik, Nature 449, 799 (2007).
[80] J. Kormendy, and D. Richstone, Annu. Rev. Astron. Astrophys. 33, 581(1995).
[81] R. Schöde et al., Nature 419, 694 (2002).
[82] J. H. Lacy, C. H. Townes, T. R. Geballe, and D. J. Hollenbach, Astrophys. J.
241, 132 (1980).
[83] E. Serabyn, and J. H. Lacy, Astrophys. J. 293, 445 (1985).
[84] R. Genzel, D. M. Watson, M. K. Crawford, and C. H. Townes, Astrophys. J.
297, 766 (1985).
[85] K. Sellgren, M. T. McGinn, E. E. Becklin, and D. N. B. Hall, Astrophys. J.
359, 112 (1990).
[86] A. Krabbe et al., Astrophys. J. 447, L95 (1995).
[87] J. W. Haller et al., Astrophys. J. 456, 194 (1996).
[88] A. Eckart, and R. Genzel, Nature 383, 415 (1996).

197

[89] K. Gebhardt, Nature 419, 675 (2002).
[90] M. C. Miller, and E. J. M. Colbert, Int. J. Mod. Phys. D 13, 1 (2004).
[91] G. Fabbiano, Science 307, 533 (2005).
[92] S. F. P. Zwart, and S. L. W. McMillan, Astrophys. J. 576, 899 (2002).
[93] D. Vanbeveren, H. Belkus, J. V. Bever, and N. Mennekens, Astrophys. Space
Sci. 324, 271 (2009).
[94] T. Pan, A. Loeb, and D. Kasen, Mon. Not. R. Astron. Soc. (2012),
doi:10.1111/j.1365-2966.2012.21030.x .
[95] M. C. Miller, and D. P. Hamilton, Mon. Not. R. Astron. Soc. 330, 232 (2002).
[96] S. Goswami, S. Umbreit, M. Bierbaum, and F. A. Rasio, (2011),[arXiv:1105.5884].
[97] T. Ohkubo, K. Nomoto, H. Umeda , N. Yoshida, and S. Tsuruta, Astrophys.
J. 706, 1184 (2009).
[98] A. Heger, and S. E. Woosley, Astrophys. J. 567, 532 (2002).
[99] E. Glebbeek, E. Gaburov, S. E. de Mink, O. R. Pols, and S. F. P. Zwart, Astron.
Astrophys. 497, 255 (2009).
[100] H. Baumgardt, J. Makino, and T. Ebisuzaki, Astrophys. J. 613, 1143 (2004).
[101] R. M. O’Leary, F. A. Rasio, J. M. Fregeau, N. Ivanova, and R.
O’Shaughnessy, Astrophys. J. 637, 937 (2006).

198

[102] K. Gebhardt, R. M. Rich, and L. C. Ho, Astrophys. J. 578, L41 (2002).
[103] J. Gerssen et al., Astron. J. 124, 3270 (2002).
[104] H. Baumgardt, J. Makino, and P. Hut, Astrophys. J. 620, 238 (2005).
[105] M. Safonova, and P. Shastri, Astrophys. Space Sci. 325, 47 (2010).
[106] K. Gültekin, M. C. Miller, and D. P. Hamilton, Astrophys. J. 616, 221 (2004).
[107] T. Okajima, K. Ebisawa, and T. Kawaguchi, Astrophys. J. 652, L105 (2006).
[108] P. Casella, G. Ponti, A. Patruno, T. Belloni, G. Miniutti, and L. Zampieri,
Mon. Not. R. Astron. Soc. 387, 1707 (2008).
[109] T. E. Strohmayer, and R. F. Mushotzky, Astrophys. J. 703, 1386 (2009).
[110] A. Patruno, and L. Zampieri, Mon. Not. R. Astron. Soc. 403, L69 (2010).
[111] F. R. Ferraro, A. Possenti, E. Sabbi, P. Lagani, R. T. Rood, N. D’Amico, and
L. Origlia, Astrophys. J. 595, 179 (2003).
[112] S. F. P. Zwart, and E. P. J. van den Heuvel, Nature 450, 388 (2007).
[113] P. Kaaret et al., Mon. Not. R. Astron. Soc. 321, L29 (2001).
[114] S. A. Farrell, N. A. Webb, D. Barret, O. Godet, and J. M. Rodrigues, Nature
460, 73 (2009).
[115] E. Vesperini, S. L. W. McMillan, A. D’Ercole, and F. D’Antona, Astrophys.
J. 713, L41 (2010).

199

[116] S. Umbreit, J. M. Fregeau, S. Chatterjee, and F. A. Rasio, Astrophys. J. 750,
31 (2012).
[117] E. T. Newman, E. Couch, K. Chinnapared, A. Exton, A. Prakash, and R.
Torrence, J. Math. Phys. 6, 918 (1965).
[118] T. Regge, and J. A. Wheeler, Phys. Rev. 108, 1063 (1957).
[119] F. J. Zerilli, Phys. Rev. D 2, 2141 (1970).
[120] S. Chandrasekhar, and S. Detweiler, Proc. R. Soc. Lond. A 344, 441 (1975).
[121] S. Chandrasekhar, Proc. R. Soc. Lond. A 358, 421 (1978).
[122] C. V. Vishveshwara, Nature 227, 936 (1970).
[123] M. Davis, R. Ruffini, W. H. Press, and R. H. Price, Phys. Rev. Lett. 27, 1466
(1971).
[124] S. L. Detweiler, Proc. R. Soc. Lond. A 352, 381 (1977).
[125] S. L. Detweiler, and E. Szedenits, Astrophys. J. 231, 211 (1979).
[126] C. T. Cunningham, R. H. Price, and V. Moncrief, Astrophys. J. 224, 643(1978).
[127] C. T. Cunningham, R. H. Price, and V. Moncrief, Astrophys. J. 230, 870(1979).
[128] C. T. Cunningham, R. H. Price, and V. Moncrief, Astrophys. J. 236, 674(1980).

200

[129] K. D. Kokkotas, and B. G. Schmidt, Quasi-Normal Modes of Stars and Black
Holes, Living Rev. Relativity 2, 2 (1999).
[130] A. Buonanno, G. B. Cook, and F. Pretorius, Phys. Rev. D 75, 124018 (2007).
[131] E. Berti, V. Cardoso, J. A. Gonzalez, U. Sperhake, M. Hannam, S. Husa,
and B. Brügmann, Phys. Rev. D 76, 064034 (2007).
[132] L. S. Finn, Phys. Rev. D 46, 12 (1992).
[133] B. P. Abbott et al. (LIGO Scientific), Phys. Rev. D 80, 062001 (2009).
[134] L. M. Goggin (2008), Ph.D. Thesis, California Institute of Technology,[arXiv:0908.2085].
[135] L. M. Goggin, Class. Quant. Grav. 23 S709 (2006).
[136] J. D. E. Creighton, Phys. Rev. D 60, 022001 (1999).
[137] B. Allen, “GRASP users manual”, (1997),
http://www.lsc-group.phys.uwm.edu .
[138] R. Adhikari (2004), Ph.D. Thesis, Massachusetts Institute of Technology,
http://dspace.mit.edu/bitstream/handle/1721.1/28646/1/58964895.pdf .
[139] Y. Tsunesada et al. (TAMA Collaboration), Class. Quant. Grav. 21, S703
(2004).
[140] J. Abadie et al. (LIGO Scientific), (2012), [arXiv:1201.5999].
[141] E. W. Leaver, Proc. R. Soc. Lond. A 402, 285 (1985).

201

[142] F. Echeverria, Phys. Rev. D 40, 3194 (1989).
[143] E. Berti, V. Cardoso, and C. M. Will, Phys. Rev. D 73, 064030 (2006).
[144] E. Berti, J. Cardoso, V. Cardoso, and Marco Cavaglià, Phys. Rev. D 76,
104044 (2007).
[145] S. Caudill, S. E. Field, C. R. Galley, F. Herrmann, and M. Tiglio, Class.
Quantum Grav. 29, 095016 (2012).
[146] S. Caudill et al., The Ringdown Search Group Internal Working Note
(LIGO) (2012).
[147] LSC Internal working note LIGO-T1000056-v5, Laser Interferometer Gravitational Wave Observatory (LIGO).
[148] N. Leroy (for the LIGO Scientific Collaboration and Virgo Collaboration),
Class. Quantum Grav. 26, 204007 (2009).
[149] J. Aasi et al., (2009), [arXiv:1203.5613].
[150] B. P. Abbott et al., Phys. Rev. D 79, 122001 (2009).
[151] S. Klimenko, I. Yakushin, A. Mercer, and G. Mitselmakher, Class. Quant.
Grav. 25, 114029 (2008).
[152] P. Ajith et al., Class. Quant. Grav. 24 S689 (2007).
[153] R. Sturani, S. Fischetti, L. Cadonati, G. M. Guidi, J. Healy, D. Shoemaker,
and A. Viceré, JPCS 243, 012007 (2010).
[154] P. Ajith et al., Phys. Rev. D 77, 104017 (2008).
202

[155] P. Ajith et al., Phys. Rev. Lett. 106, 241101 (2011).
[156] L. Santamaría et al., Phys. Rev. D 82, 064016 (2010).
[157] T. Damour, A. Nagar, M. Hannam, S. Husa, and B. Brugmann, Phys. Rev.
D 78, 044039 (2008).
[158] Y. Pan, A. Buonanno, M. Boyle, L. T. Buchman, L. E. Kidder, H. P. Pfeiffer,and M. A. Scheel, Phys. Rev. D 84, 124052 (2011).
[159] K. Thorne, in Three Hundred Years of Gravitation, edited by S. Hawking and
W. Israel (Cambridge University Press, Cambridge, U.K.; New York, U.S.A.,
1987), pp. 330–458.
[160] J. Abadie et al. (LIGO Scientific), Class. Quant. Grav. 27, 173001 (2010),
1003.2480.
[161] D. Sigg, Class. Quant. Grav. 25, 114041 (2008).
[162] H. Grote, Class. Quant. Grav. 25, 114043 (2008).
[163] F. Acernese et al., Class. Quant. Grav. 25, 184001 (2008).
[164] S. Nissanke, D. E. Holz, S. A. Hughes, N. Dalal, and J. L. Sievers, Astrophys. J. 725, 496 (2010), 0904.1017.
[165] S. Fairhurst, New J. Phys. 11, 123006 (2009), 0908.2356.
[166] I. W. Harry, and S. Fairhurst (2010), LIGO Document No. P1000104-v6.
[167] D. Eichler, M. Livio, T. Piran, and D. N. Schramm, Nature 340, 126 (1989).

203

[168] B. Paczynski, Acta Astron. 41, 257 (1991).
[169] B. Abbott et al. (LIGO Scientific Collaboration), Phys. Rev. D 77, 062004(2008), 0709.0766.
[170] J. Abadie et al. (LIGO Scientific Collaboration and Virgo Collaboration),
Astrophys. J. 715, 1453 (2010), 1001.0165.
[171] B. Abbott et al. (LIGO Scientific Collaboration), Phys. Rev. D 77, 062002(2008), 0704.3368.
[172] B. P. Abbott et al. (LIGO Scientific Collaboration), Phys. Rev. D 79, 122001(2009), 0901.0302.
[173] B. P. Abbott et al. (LIGO Scientific Collaboration), Phys. Rev. D 80, 047101(2009), 0905.3710.
[174] B. Allen, Phys. Rev. D 71, 062001 (2005), gr-qc/0405045.
[175] B. S. Sathyaprakash, and S. V. Dhurandhar, Phys. Rev. D 44, 3819 (1991).
[176] LSC Algorithms Library,https://www.lsc-group.phys.uwm.edu/daswg/projects/lalsuite.ht.
[177] L. Blanchet, T. Damour, B. R. Iyer, C. M. Will, and A. G. Wiseman, Phys.
Rev. Lett. 74, 3515 (1995), gr-qc/9501027.
[178] L. Blanchet, B. R. Iyer, C. M. Will, and A. G. Wiseman, Class. Quant. Grav.
13, 575 (1996), gr-qc/9602024.

204[179] P. Jaranowski, A. Krolak, K. D. Kokkotas, and G. Tsegas, Class. Quant.
Grav. 13, 1279 (1996).
[180] M. Rakhmanov, Class. Quant. Grav. 23, S673 (2006), gr-qc/0604005.
[181] S. Klimenko, S. Mohanty, M. Rakhmanov, and G. Mitselmakher, Phys. Rev.
D 72, 122002 (2005), gr-qc/0508068.
[182] S. D. Mohanty, M. Rakhmanov, S. Klimenko, and G. Mitselmakher, Class.
Quant. Grav. 23, 4799 (2006), gr-qc/0601076.
[183] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical
Recipes 3rd Edition: The Art of Scientific Computing (Cambridge University
Press, 2007).
[184] A. Rodrígues (2007), Masters Thesis, Lousiana State University, Baton
Rouge.
[185] D. A. Brown (2004), Ph.D. Thesis, University of Wisconsin, Milwaukee,[arXiv:0705.1514].
[186] C. A. K. Robinson, B. S. Sathyaprakash, and A. S. Sengupta, Phys. Rev. D
78, 062002 (2008), 0804.4816.
[187] W. G. Anderson, P. R. Brady, J. D. E. Creighton, and É. É. Flanagan, Phys.
Rev. D 63, 042003 (2001).

205