Signatures and Dynamics of Compact Binary Coalescences
and a Search in LIGO’s S5 Data

Thesis byDrew Garvin Keppel

In Partial Fulfillment of the Requirements
for the Degree ofDoctor of Philosophy

California Institute of Technology
Pasadena, California

2009(Defended April 20, 2009)

ii

c 2009�
Drew Garvin Keppel
All Rights Reserved

iii

To my everloving wife, Lisa

iv

Acknowledgements
Without the help of a multitude of people, all of the work presented here would not have beenpossible. I would like to thank my advisor, Alan Weinstein, for taking me under his wing as a
student and sharing his passion for LIGO with me. He has given me guidance as a developing
scientist, as well as the freedom to explore the paths I choose.
Additionally, I would like to thank my committee members, Yanbei Chen, Jay Marx, Tom Prince,
and Kip Thorne for all of their time and advice with regard to my thesis.
I would like to thank the review committee members for the S5 first year search, Marie-AnneBizouard, Laura Cadonati, Eric Gustafson, Vicky Kalogera, Greg Mendell, Peter Saulson, and John
Whelan, who have been an invaluable resource for validating the S5 first year search.
I would like to thank the many postdoctoral scholars and research scientists who have walked
the halls of the 6th floor of Millikan during my time at Caltech. Your examples have driven me tobecome a better scientist. I would especially like to thank Duncan Brown, Stephen Fairhurst, andEirini Messaritaki for answering my endless questions about CBC searches and helping to guide me
through the process; Kipp Cannon, Chad Hanna, and Antony Searle for sharing their visions of the
future of CBC searches with me; and Parameswaran Ajith for his insight in IMR waveforms andhelp with parameter estimation.
I would also like to thank my fellow graduate students at Caltech: Diego Fazi, Lisa Goggin, Kari
Hodge, Pinkesh Patel, and Becky Tucker; and throughout the LSC: Rahul Biswas, Collin Capano,
Nick Fotopolus, Ian Harry, Gareth Jones, Dave McKechan, Craig Robinson, Andy Rodriguez, Lucia
Santamaria, Jake Slutsky; as well as other members of the CBC group: Sukanta Bose, Patrick Brady,
Jessica Clayton, Thomas Cokelaer, Jolien Creighton, Alex Dietz, Gaby Gonzalez, Romain Gouaty,

v
Anand Sengupta, Ruslan Vaulin; and collaborators in the LSC-Virgo Collaboration (and any one
else I have failed to mention) for the friendships we have built and the many conversations we have
had throughout the years while working, traveling, and playing together.
Pertaining to my work in relativity, I would like to thank Kip Thorne for leading us on afascinating exploration that started from a homework problem; Yanbei Chen for his endless insights;
and David Nichols, Jaﬀ Kaplan, and the rest of the TAPIR group for the many discussions we have
had.
I would also like to thank all of the administrative assistants: including Irene Baldon, JoAnn
Boyd, and Marylin Wright; and LIGO computer staﬀ: Phil Ehrens, Dan Kozak, and Stuart Anderson, with whom I have interacted, for all of their aid and support, which allowed me to focus morestrongly on my research.
Last, but certainly not least, I would like to thank my parents and siblings for their love andsupport over the years and for encouraging me to pursue my dreams; my wife, Lisa, for her unending
love, patience with my work, and support on so many levels; and Maxwell’s Equation for his perpetualenergy and enthusiasm with which he tackles life. I also thank God for all of the opportunities andabilities He has given me to explore His creation.

vi

Abstract
In the chapters that follow, we probe several aspects of compact binary environments, focusing
on results of orbits and collisions of compact objects. First, we describe a search for low-masscompact-binary-coalescence gravitational-wave signals in data from the LIGO detectors’ most sensitive, longest-running science run to date (S5). Next, we go into detail on the interpretation of theresults including its development. We then investigate the bounds on the mass of the graviton that
could be achieved from the detection gravitational waves from a binary black hole merger. Last, westudy the flow of momentum in compact binaries using the Landau-Lifshitz formalism.

vii

Contents
Acknowledgementsiv

Abstractvi

1 Introduction

1

2 Gravitational Waves and Their Sources, Including Compact Binary Coalescences3

2.1

Pillars of Relativity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3

2.2

Einstein’s Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

2.3

Gravitational Wave Solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

2.3.1

Linearized Gravity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5

2.3.2

Geodesic Deviation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

2.3.3

Gravitational-Wave Polarizations . . . . . . . . . . . . . . . . . . . . . . . . .

7

Gravitational Wave Sources . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

9

2.4.1

Compact Binaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

10

2.4.1.1

Compact Binary Inspiral Waveforms . . . . . . . . . . . . . . . . . .

11

2.4.1.2

Inspiral Waveform Approximants . . . . . . . . . . . . . . . . . . . .

16

2.4.1.3

Inspiral Waveform Stationary Phase Approximation . . . . . . . . .

17

Other Types of Gravitational Waves . . . . . . . . . . . . . . . . . . . . . . .

17

2.4

2.4.2

3 Introduction to LIGO

19

3.1

History of Gravitational-Wave Detectors . . . . . . . . . . . . . . . . . . . . . . . . .

19

3.2

Design of LIGO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

22

viii
3.3

Dominant Noise Sources . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

26

3.4

Calibration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

28

3.5

Interferometer Antenna Response . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

30

3.6

Networks of Detectors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

35

3.7

LIGO’s Past and Future . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

38

4 Astrophysics of Compact Binaries

39

4.1

Evidence for Neutron Stars and Black Holes . . . . . . . . . . . . . . . . . . . . . . .

39

4.2

Predicted Compact Binary Coalescence Rates . . . . . . . . . . . . . . . . . . . . . .

40

4.2.1

Extrapolations from Merging Binary-Neutron-Star Observations . . . . . . .

40

4.2.2

Population Synthesis of Field Binaries . . . . . . . . . . . . . . . . . . . . . .

41

4.2.3

Dynamical Simulations of Star Clusters . . . . . . . . . . . . . . . . . . . . .

42

4.2.4

Extrapolations from Short-Gamma-Ray-Burst Observations . . . . . . . . . .

43

5 Gravitational Wave Signal Processing

46

5.1

The Optimal Filter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

46

5.2

The Waveform Overlap . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

49

5.2.1

Matched Filtering for Compact Binary Signals . . . . . . . . . . . . . . . . .

50

5.2.2

Template Bank Construction . . . . . . . . . . . . . . . . . . . . . . . . . . .

52

5.2.3

Signal-Based Vetoes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

54

5.2.3.1

The χ2 Veto Calculation . . . . . . . . . . . . . . . . . . . . . . . . .

55

5.2.3.2

The r2 Veto Calculation . . . . . . . . . . . . . . . . . . . . . . . . .

56

5.2.3.3

The Eﬀective Distance Cut . . . . . . . . . . . . . . . . . . . . . . .

56

5.2.3.4

The Amplitude Consistency Check . . . . . . . . . . . . . . . . . . .

58

Previous Detection Statistic: Eﬀective SNR . . . . . . . . . . . . . . . . . . . . . . .

59

5.3

6 Data Analysis Pipeline and Tuning

60

6.1

Blinding the Analysis Pipeline. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

64

6.2

Simulated Waveform Injections . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

65

ix
6.3

Template Bank Generation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

66

6.4

First-Stage Single-Detector Trigger Generation . . . . . . . . . . . . . . . . . . . . .

68

6.5

First-Stage Coincident Trigger Generation . . . . . . . . . . . . . . . . . . . . . . . .

71

6.5.1

Coincidence Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

72

6.5.2

Background Estimation: Time Shifts . . . . . . . . . . . . . . . . . . . . . . .

74

6.5.3

Coincidence Tuning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

75

6.6

Trigger to Template . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

77

6.7

Second-Stage Single-Detector Trigger Generation . . . . . . . . . . . . . . . . . . . .

77

6.7.1

Signal-Based Veto Tuning . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

78

Second-Stage Coincident Trigger Generation . . . . . . . . . . . . . . . . . . . . . . .

84

6.8.1

Application of Data Quality Vetoes . . . . . . . . . . . . . . . . . . . . . . . .

84

6.8.2

Multidetector Signal-Based Veto Tuning . . . . . . . . . . . . . . . . . . . . .

88

6.8

7 False Alarm Rate Detection Statistic
7.1

90

Definition of Trigger Categories . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

91

7.1.1

Coincidence Types . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

91

7.1.2

Veto Diﬀerences . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

92

7.1.3

Mass Diﬀerences . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

95

7.1.4

Trigger Categories . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

98

7.2

FAR: False Alarm Rate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

98

7.3

FARc: Combined False Alarm Rate . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100

7.4

Comparison of Detection Statistics . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103

7.5

Plotting Background Trials: Lightning Bolt Plots . . . . . . . . . . . . . . . . . . . . 105

7.6

False Alarm Probability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107

8 Ranking of Detection Candidates

108

8.1

Separating Triggers by Data Quality Flag Vetoes . . . . . . . . . . . . . . . . . . . . 110

8.2

Calculate Analyzed Times . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113

x
8.3

Trigger Clustering and Separation by Category . . . . . . . . . . . . . . . . . . . . . 115

8.4

FAR Calculation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116

9 Rate (Upper Limit) Calculation

117

9.1

Posterior and Upper Limit Calculation . . . . . . . . . . . . . . . . . . . . . . . . . . 120

9.2

Incorporating Systematic Errors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121

9.3

Cumulative Luminosity and Background Probability . . . . . . . . . . . . . . . . . . 122

9.4

9.3.1

Cumulative Luminosity and its Derivative . . . . . . . . . . . . . . . . . . . . 122

9.3.2

Background Probability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123

Systematic Error Calculation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124
9.4.1

Monte Carlo Errors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124

9.4.2

Calibration Errors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125

9.4.3

Waveform Errors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126

9.4.4

Galaxy Errors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126
9.4.4.1

Distance Errors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126

9.4.4.2

Magnitude Errors . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127

10 Results of the Search for Low Mass Compact Binary Coalescences in the FirstYear of Data from LIGO’s Fifth Science Run

129

10.1 Follow-up Procedure for Coincident Triggers . . . . . . . . . . . . . . . . . . . . . . . 130
10.2 Loudest Events . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132
10.3 Coalescence Rate Upper Limits . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135
10.3.1 Upper Limits Neglecting Spin . . . . . . . . . . . . . . . . . . . . . . . . . . . 137
10.3.2 Upper Limits Including Spin . . . . . . . . . . . . . . . . . . . . . . . . . . . 139
10.4 Summary of the Search . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140
11 Future Improvements for Compact Binary Coalescence Searches

141

11.1 Detection Statistic . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141
11.2 Background Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143

xi
11.3 Coherent Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144
11.4 Data Cleaning and Noise Regression . . . . . . . . . . . . . . . . . . . . . . . . . . . 145
11.5 Low Latency Searches . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145
12 Testing Eﬀects of a Massive Graviton

147

12.1 Eﬀects of Massive Gravitons . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148
12.2 Dispersion of Gravitational Waves . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148
12.3 Dispersion of Compact Binary Coalescence Gravitational Waveforms . . . . . . . . . 151
12.4 Recovery of Dispersed Compact Binary Coalescence Gravitational Waveforms . . . . 153
12.5 Using the Dispersion of Gravitational Waves to Bound the Mass of the Graviton . . 154
12.6 Beyond the Fischer Matrix Approach . . . . . . . . . . . . . . . . . . . . . . . . . . . 162
12.7 Other Bounds on the Graviton’s Mass . . . . . . . . . . . . . . . . . . . . . . . . . . 162
13 Toward Understanding Black Hole Merger Dynamics

164

13.1 The Landau-Lifshitz Formalism in Brief . . . . . . . . . . . . . . . . . . . . . . . . . 167
13.2 Momentum Conservation for a Fully Nonlinear Compact Binary. . . . . . . . . . . 169

13.3 Post-Newtonian Momentum Flow in Generic Compact Binaries . . . . . . . . . . . . 171
13.3.1 Field Momentum Outside the Bodies . . . . . . . . . . . . . . . . . . . . . . . 171
13.3.1.1 Field Momentum Specialized to the Extreme-Kick Configuration . . 175
13.3.2 Centers of Mass and Equation of Motion for the Binary’s Compact Bodies . . 175
13.3.2.1 Equations of Motion Specialized to the Extreme-Kick Configuration 179
13.3.3 The Momenta of the Binary’s Bodies . . . . . . . . . . . . . . . . . . . . . . . 180
13.3.3.1 The Momentum of a Star . . . . . . . . . . . . . . . . . . . . . . . . 180
13.3.3.2 The Momentum of a Black Hole . . . . . . . . . . . . . . . . . . . . 182
13.4 Momentum Conservation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 186
13.4.1 Momentum Conservation for Extreme-Kick Configuration . . . . . . . . . . . 186
13.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 188
Bibliography

189

xii

List of Figures
2.1

Eﬀects of Passing Gravitational Waves . . . . . . . . . . . . . . . . . . . . . . . . . . .

8

2.2

Inspiral Waveform . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

15

3.1

LIGO Schematic . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

25

3.2

LIGO Noise Budget . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

27

3.3

Model Control Loop . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

28

3.4

Euler Angles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

31

3.5

Detector Antenna Patterns . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

34

3.6

Individual Detector Angle-Average Sky Maps . . . . . . . . . . . . . . . . . . . . . . .

36

3.7

Combined Detector Angle-Average Sky Maps . . . . . . . . . . . . . . . . . . . . . . .

37

3.8

LIGO Science Runs Timeline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

38

6.1

The HIPE Pipeline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

63

6.2

Data Preprocessing. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

68

6.3

Sample Template Bank . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

69

6.4

Trigscan Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

72

6.5

E-Thinca Tuning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

76

6.6

χ2 Tuning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

79

6.7

χ2 Threshold . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

80

6.8

r2 Tuning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

81

6.9

Diﬀerent r2 Veto Populations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

82

6.10

r2 Sample Time Series . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

83

xiii
6.11

Eﬀects of DQ Vetoes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

87

6.12

H1H2 Eﬀective Distance Cut . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

89

7.1

Eﬀect of Coincidence on Backround Trigger Rate . . . . . . . . . . . . . . . . . . . . .

92

7.2

Eﬀect of Amplitude Veto on Backround Trigger Rate . . . . . . . . . . . . . . . . . .

93

7.3

Eﬀects of Diﬀerent Amplitude Vetoes for H1 and H2 . . . . . . . . . . . . . . . . . . .

94

7.4

Combined ρ2eﬀ Distributions for Background Triggers for Diﬀerent Chirp Mass Bands

96

7.5

Combined ρ2eﬀ Distributions for Background Triggers for Diﬀerent Chirp Mass Bins .

97

7.6

Background IFAR Plots . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

99

7.7

Combined Background IFAR Plot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101

7.8

Background IFARc Plot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102

7.9

Comparison of Detection Statistics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104

7.10

In-Time IFARc Plot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106

8.1

Postprocessing Pipeline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109

8.2

septime Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112

8.3

Background Analyzed Time . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114

9.1

Upper Limit Calculation Pipeline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119

10.1

Comparison of Triggers from Loudest Event to an Injected Signal . . . . . . . . . . . 131

10.2

IFAR Loudest Events . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134

10.3

Upper Limits versus Mass . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138

12.1

Waveform Inversion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152

12.2

Massive Graviton Bounds . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161

13.1

Extreme-kick Configuration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 164

13.2

Bobbing Motion in Extreme-kick Configuration . . . . . . . . . . . . . . . . . . . . . . 166

13.3

Pretorius’ Explanation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 166

13.4

Division of Space . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 169

xiv
13.5

Field Momentum Density . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 176

xv

List of Tables
3.1

Combined Horizon Volumes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

35

10.1

Detailed Results of the BNS Upper Limit Calculation . . . . . . . . . . . . . . . . . . 135

10.2

Detailed Results of the BBH Upper Limit Calculation . . . . . . . . . . . . . . . . . . 135

10.3

Detailed Results of the BHNS Upper Limit Calculation . . . . . . . . . . . . . . . . . 136

12.1

Phenomenological Coeﬃcients . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156

12.2

Massive Graviton Bounds . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 160

13.1

Pieces of Total Momentum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 187

1Chapter 1

Introduction
General relativity is a very rich theory that makes many exciting predictions about the universe’s
most extreme environment. One of those predictions is the emission of gravitational waves (GWs)from accelerating objects (e.g., orbiting compact binaries), which have yet to be detected directly.
To that end, the Laser Interferomenter Graviational-wave Observatory (LIGO) detectors have been
built and measure these waves as they pass through the earth. Recently, the LIGO detectors reachedtheir design sensitivity and recorded data from their fifth science run, which accumulated a year’s
worth of coincident data where multiple detectors were in operation. The first part of this thesis will
focus on a search for GW signals from compact binary coalescence (CBC) (LIGO’s flagship search)in the first year of data from this science run. In those chapters we describe the signals we searchfor, the design of the GW detectors, how the search was performed, and its results. For this part of
the thesis, my role was to tune and run the analysis mentioned above. This included contributing to
the code base used for the analysis, implementing new features in many stages of the pipeline, anddesigning and implementing the postprocessing pipeline to calculate the detection statistic. Finally,
this included managing the writing of and acting as the corresponding author on the resulting paper.
In the second part of this thesis, we focus on testing the eﬀects of a massive graviton. That
chapter describes a massive graviton’s eﬀects on the propagation of GWs and uses those eﬀects toconstrain the graviton’s mass for a given CBC GW signal observation. The main result of my work is
an application of previous methods to a new class of CBC waveforms that contain more information(i.e., complete CBC waveforms) and can more stringently constrain the graviton’s mass.

2
The final part of this thesis focuses on the dynamics of compact objects as they interact andsubsequently merge. Momentum flow between the bodies and the surrounding gravitational field
is studied and found to explain interesting features seen in Numerical Relativity simulations oforbiting black holes that have nonzero spin angular momentum. For this part of the thesis, my role
was to perform and verify the surface and volume integral calculations associated the momentum
conservation results.

3

Chapter 2

Gravitational Waves and Their
Sources, Including Compact
Binary Coalescences
In this chapter we give a brief introduction to General Relativity, focusing on GW emission. We
then focus our attention on GWs from CBCs. The following presentation uses the conventions foundin [1].

2.1

Pillars of Relativity

Einstein developed the pillars of Relativity Theory in the early 20th Century, which can be summarized by two thought experiments, or “gedankenexperiments.” These gedankenexperiments reveal
diﬀerent implications of a single assumption, that the Laws of Physics are the same for everyone
observing an experiment.
The first application of this assumption was to the speed of light. Einstein assumed that the
speed of light was a constant, independent of an observer’s frame of reference. The conclusions
Einstein drew from this led him to develop the Special Theory of Relativity, which is useful inLorentz reference frames (i.e., those without gravity).
The second application of this assumption was the universality of the motion of objects in gravitational freefall. Einstein’ gedankenexperiment went as follows: suppose there are two observers
observing the motion on a object in freefall in their local reference frame. Over small distances and4times, they will observe the same motion of the object whether the motion is due to gravitational
acceleration towards a gravitating body (e.g., for an observer whose frame is resting on the Earth),or due to a uniform acceleration of the reference frame (e.g., for an observer whose frame is in a
rocket ship that is constantly accelerating in a specific direction). From this gedankenexperiment,
Einstein concluded that the acceleration around a gravitating body is produced by a curvature ofspacetime in the object’s vicinity. Using this framework, Einstein developed the General Theory ofRelativity.

2.2

Einstein’s Equations

In the General Theory of Relativity, the equations that link the curvature of spacetime to thepresence of matter and energy are called the Einstein Field Equation (EFE) given as

Gµν =

8πG
Tµν ,
c4

(2.1)

where Gµν is the Einstein Tensor given by Gµν = Rµν − 12 Rgµν , R = g µν Rµν , Rµν = Rγ µνγ , Rγ µνγ
is the Riemann Curvature Tensor, gµν is the spacetime metric, and Tµν is the stress-energy tensor.
These equations tell matter how to move due to the curvature of spacetime and spacetime how tocurve due to the presence of matter.

2.3

Gravitational Wave Solutions

The EFEs not only describe how spacetime curvature is influenced in the presence of matter, it also
predicts how the curvature should propagate away from accelerating objects. Just as electromagnetic
waves are emitted and travel away from accelerating electric charges, GWs are emitted and travel
away from accelerating masses. If we are far away from the accelerating source, the EFEs simplifytoGµν = 0 .

(2.2)

5
In this region, the spacetime is essentially Minkowski except for the perturbations do to GWs:gµν = ηµν + hµν , where ηµν is the Minkowski metric. Assuming these perturbations are small
so that all components of hµν � 1, we can calculate how these perturbations propagate using alinearized version of the EFEs. In the following discussion, we follow closely the description foundin reference [1].

2.3.1Linearized Gravity

To first order in hµν the EFE become−h̄µν,α α − ηµν h̄αβ, αβ + h̄µα, α ν + h̄να, α µ = 0 ,(2.3)

where h̄µν ≡ hµν − 12 ηµν h and h ≡ η αβ hαβ . We can eliminate all of the terms after the first by
imposing the tensor analogue of the Lorentz gauge condition of electromagnetic theory, namelyh̄µα, α = 0. This leaves us with the wave equationh̄µν,α α = 0 .

(2.4)

This wave equation has plane wave solutions��αh̄µν = � Aµν eikα x ,(2.5)

where � [. . .] means the real part of [. . .]. From the wave equation we find kα k α = 0, which meansk is a null vector and the waves travel at the speed of light. From the gauge condition h̄µα, α = 0
we find Aµα k α = 0, which means the amplitude tensor A is orthogonal to k.
At this point A has six degrees of freedom, however four of these can be pinned down with anadditional gauge specification. Let us choose Aµν uν = 0 where uν is a 4-velocity of a particular
frame. This fixes an additional three degrees of freedom. For the final gauge choice, let us choose6Aµ µ = 0. All of these constraints in a Lorentz frame with u0 = 1 and uj = 0 givehµ0 = 0 ,(2.6)

hkj, j = 0 ,(2.7)

hkk = 0 .

(2.8)

The first of these equations implies only the spatial components of hµν are nonzero. The second
implies that the spatial components are divergence free. Finally, the third implies that hµν is trace
free. This gauge is called the transverse-traceless (TT) gauge because in this gauge hµν is transverseto the time direction and the direction of propagation, and is trace free.
These gauge choices leave us with two degrees of freedom for the propagating waves. This can beinterpreted as two polarizations of GWs, the plus (+) and cross (×) polarizations, named for theiraﬀect on test particles. In terms of the these polarizations and in the TT gauge, Aµν takes the form


 0


 0

Aµν = 

 0



0

2.3.2

0

0

h+

h×

h×

−h+

0

0



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

(2.9)

Geodesic Deviation

Objects undergoing free motion will follow geodesics described by the geodesic equation

d2 xλ
dxµ dxν+ Λλ µν
=0.
2
dt
dt dt

(2.10)

Consider two particles, A and B, separated by nj in a coordinate system that is a local Lorentz
frame along A’s world line. In this coordinate system the separation vector becomes nj = xB j −xA j = xB j , since A is at the origin. In a TT coordinate system that moves with particle A andwith its proper reference frame to first order in the metric perturbation hTTjk , the geodesic deviation

7
equation is written as
d 2 njTT= −Rj0k0
nk .
∂t2(2.11)In the TT gauge the Riemann curvature tensor takes the simple form

1
Rj0k0 = − hTT.
2 jk,00(2.12)

Using this along with the separation vector’s definition, the geodesic deviation equation becomesd2 xB j1 ∂ 2 hTTjk=
xB k .
∂t2
2 ∂t2(2.13)

This equation of motion can be integrated, assuming particles initially at rest relative to eachother, to obtain the location of B as a function of time in the proper reference frame of A�1
xB k (t) = xB (0) δjk + hTT
2 jk
j

�

,

(2.14)

where hTT
jk should be evaluated at the location of A. The second in the parentheses on the right-handside can be viewed as a strain and is proportional to the metric perturbation hjk .

2.3.3

Gravitational-Wave Polarizations

As mentioned above, there are two degrees of freedom for a GW, which can be seen as distinct
polarizations. Let us consider a plane, monochromatic wave propagating in the z direction. In the
TT gauge of this wave, the gauge constraints (equation (2.6)) leave the only nonzero components of
hTT
µν
�
�
TT
−iω(t−z)hTT,
xx = −hyy = � A+ e��TT−iω(t−z)hTT,xy = hyx = � A× e

(2.15)(2.16)

8

Figure 2.1: Eﬀects of Passing Gravitational Waves
The eﬀect of the + and × polarization GW passing perpendicular to a ring of test particles.
where A+ and A× are the amplitudes of the + and × polarizations, respectively, and ω is thefrequency of the GW.
Defining the polarization tensors e+ = (ex ⊕ ex ) − (ey ⊕ ey ) and e× = (ex ⊕ ey ) + (ey ⊕ ex ) we
find a GW can be decomposed into its independent polarizations as

hTT
jk = h+ e+jk + h× e×jk .

(2.17)The eﬀects of each polarization can be seen independently by looking at a circular ring of test
particles surrounding a central particle: let all of the particles lie in the plane perpendicular to thedirection of propagation of the incident GW. For an incident +-polarized wave, the polarization
tensor e+ shows that the test particles will oscillate about their circular equilibrium positions insuch a way that as the distance from the origin along the x direction increases, the distance fromthe origin along the y direction decreases. This results in the circle being stretched and squeezed
along the x and y axes at one moment, and vice versa half a wavelength later. For a ×-polarized
wave, a similar thing happens, except rotated by 45◦ . This eﬀect is visualized in figure 2.1.

9

2.4

Gravitational Wave Sources

In the case of electromagnetic radiation, a multipolar analysis reveals that electromagnetic waves are
sourced by, in order of significance, accelerating electric-charge monopoles, electric-charge dipoles
and electric-current dipoles, electric-charge quadrupoles and electric-current quadrupoles, etc. In thecase of gravitational radiation, the conservation of energy prevents radiation due to the accelerationof mass monopoles, the conservation of linear momentum prevents radiation due to the acceleration ofmass dipoles, and the conservation of angular momentum prevents radiation due to the accelerationof current dipoles. This results in the leading-order radiation coming from the acceleration of mass
quadrupoles.
In reference [2], Thorne gives the relationships between GW radiation and its source. Radiation
can be related to the source through the multipolar expansion. For Newtonian sources the relevant
multipolar moments are given by

IAl =
SAl =

��

��

�STF
ρXAl d3 x,(2.18a)�STF(�al pq xp ρvq ) XAl−1 d3 x,

(2.18b)where XAl = xa1 xa2 . . . xal , ρ is the Newtonian mass density, v is the velocity, �ijk is the Levi-Civita
tensor, and STF stands for the symmetric trace-free part.
For these sources, the radiation hTT
jk is given by

hTT
jk

=

�∞
�4
l=2

∞
�
+
l=2

l!

rad
r−1(l) IjkA(t − r)NAl−2
l−2

8
rad
r−1 �pq(j (l) Sk)pA(t − r)nq NAl−2l−2(l + 1)!

�TT,(2.19)

where (l) I and (l) S represent the lth time derivative of these moments, and NAl = n1 n2 . . . nl .
The superscript “rad” denotes this as the moment in the radiation field, which is the same as themoment in the near zone if we ignore nonlinear corrections proportional to ln(λ/2πR) where R isthe background curvature of the local universe.

10
The energy and momentum carried away from a radiating system by GWs is given simply in theform of the stress-energy tensor for GWs:

GW
Tµν
=

1 � TT TT �hh.
32π jk,µ jk,ν(2.20)

Expanding the energy term (i.e., µ = ν = 0) in terms of the multipole moments, we find theenergy flux carried by the waves to be

dE
dt

=

∞
�
(l + 1)(l + 2)
l=2
∞
�

+

l=2

(l − 1)l�l+1�
1
IAl l+1 IAll!(2l + 1)!!

�l+1�
4l(l + 2)
1
SAl l+1 SAl ,(l − 1) (l + 1)!(2l + 1)!!

(2.21)

where l! = l(l − 1) · · · 2 · 1 and l!! = l(l − 2)(l − 4) · · · (2 or 1).
GWs can come any source with an accelerating mass or current multipole, however, in mostsituations the dominant term in the multipole expansion is the l = 2 mass quadrupole. Diﬀerent
types of sources can be characterized into the following categories: known-transient waveforms (e.g.,compact binary signals), unknown-transient waveforms (i.e., burst signals), discernible and narrowband continuous wave sources (e.g., pulsar and white-dwarf binary signals), and indiscernible orbroad-band continuous wave sources (e.g., stochastic and confused white-dwarf binary sources). In
section 2.4.1 we derive the waveforms associated with compact binary signals. In section 2.4.2 webriefly discuss additional examples and how each type of waveform is searched for.

2.4.1

Compact Binaries

One of the families of known transient waveforms is associated with compact binary systems (i.e.,binary systems whose objects are either neutron stars or black holes). As these objects orbit eachother, they act as a source for GWs.
Compact binaries are interesting signals, in part, because they are so well understood. Thesesignals can act as standard candles for measuring astrophysical distances due to the fact that their11amplitude is uniquely determined by their phase evolution and the luminosity distance to the source.
In addition, GW signals from compact binaries are considered some of the most extreme probes ofGeneral Relativity as the signals originate from objects that massively curve space-time and travelat speeds approaching the speed of light.

2.4.1.1

Compact Binary Inspiral Waveforms

In this section we go over the leading-order eﬀects in calculating a compact binary inspiral waveform.
To aid in this calculation, we assume the objects are Newtonian point particles with masses mA andmB . In this case, the mass density term in equation (2.18a) become Dirac delta functions and thequadrupole moment Ijk is given as

Ijk

=
=

�

�

mA xAj xAk

A

�

mA

A

�

�STF

1
2
xAj xAk − δjk rA
3�.

(2.22)The leading-order term in the GW radiation (equation (2.19)) is due to the mass quadrupole andgiven as
hTT
jk =

2 TT
Ï (t − D) ,D jk

(2.23)

where D is the distance from the source.
Let us assume these objects are in quasicircular orbits around each other described by Kepler’s
Laws where their orbital speed is related to their separation by

Ω=

�

M
,
r3

(2.24)

where M = mA + mB is the total mass of the system, and r is the separation.
As these objects orbit each other, in a coordinate system whose origin is at the center of mass

12with the plane of the orbit in the x-y plane, their locations are given byµr cos(Ωt) ,mAµ
yA (t) =r sin(Ωt) ,mAµ
xB (t) = −r cos(Ωt) ,mBµ
yB (t) = −r sin(Ωt) ,mB
xA (t) =(2.25a)(2.25b)(2.25c)(2.25d)where µ = mA mB /M is the reduced mass of the system.
Using these trajectories and equation (2.22) for Ijk , we find quadrupole moments for this systemas��1Ixx = r2 µ cos2 (Ωt) −
,
3
�
�
1
Iyy = r2 µ sin2 (Ωt) −
,
3
�
�
1
Ijz = r2 µ −δjz,
3(2.26b)

Ixy = Ixy = r2 µ sin(Ωt) cos(Ωt) .

(2.26d)(2.26a)(2.26c)As GWs leave the system, they carry away energy, which has a back-reaction eﬀect on theorbiting binary. The energy carried away by quadrupole radiation is given from equation (2.21) asdE
1 ... ...
= � I jk I jk � ,dt mass quadrupole
5(2.27)

which, substituting these moments in, we find to bedE
dt mass quadrupole

=
=

�
�2 �
32r4 µ2 Ω6 �
4 cos2 (Ωt) sin2 (Ωt) + 2 sin2 (Ωt) − 1
54 2 6
32r µ Ω.
5

(2.28)

13
Using the Keplerian relation equation (2.24), we find this to bedE
32µ2 M 3=.
dt mass quadrupole
5r5(2.29)The energy contained in the system as a function of their separation is given by the Newtonian
formula
E=−

µM
.
2r

(2.30)As energy leaves, the orbiting objects sink further into their respective gravitational potentials,
causing their orbital separation to decrease. Assuming this energy loss occurs adiabatically, we canfind the orbital separation as a function of time. First, we find the derivative of the separation withrespect to time as
�

� �
�
dE
dE
=
/
dt
dr
64µM 2
= −
.
5r3

dr
dt

(2.31)

Integrating this, we find the separation as a function of time to ber(t) =

�

256µM 25

� 14

1

(tc − t) 4 ,(2.32)

where tc denotes the time when r = 0.
We can use this result along with equation (2.24) to find the orbital frequency Ω(t) as a functionof timeΩ(t) =

�

5

256M 3
5

� −3
8

−3

(tc − t) 8 ,where M = η 3/5 M is the chirp mass and η = µ/M is the symmetric mass ratio.

(2.33)

14
We can now use equation (2.23) to find the quadrupole radiation asTT
hc (t) ≡ −hTT
xx = hyy�
4r2 (t)µΩ2 (t) � 2
cos (Ω(t)t) − sin2 (Ω(t)t)
D
4r2 (t)µΩ2 (t)cos(2Ω(t)t) ,
D

=
=

TT
hs (t) ≡ −hTT
xy = −hyx

8r2 (t)µΩ2 (t)sin(Ω(t)t) cos(Ω(t)t)
D
4r2 (t)µΩ2 (t)sin(2Ω(t)t) .
D==(2.34a)(2.34b)

We can split the radiation into separate parts for the amplitude evolution A(t) and the phaseevolution Φ(t). The radiation then takes the formhc (t) = A(t) cos(2Φ(t)) ,(2.35a)hs (t) = A(t) sin(2Φ(t)) .

(2.35b)Using the formulae for the separation and orbital frequency evolutions, we find A(t) and Φ(t) tobeA(t) ==

4r2 (t)µΩ2 (t)
D
4µM
Dr(t)
1

=

Φ(t) =−1
554
M 4 (tc − t) 4 ,D�(2.36a)Ω(t)dt= φ0 −

�tc − t5M

� 58,(2.36b)where φ0 is the orbital phase at coalescence.
The signal described by equations (2.36a and 2.36b) is characterized by a sinusoidal oscillation

15

Figure 2.2: Inspiral Waveform
The frequency f and strain h evolution of an inspiral waveform given in arbitrary units.
that increase in both amplitude and frequency as a function of time (i.e., a chirp signal). This is
illustrated in figure 2.2, which shows both the frequency evolution and the strain evolution over a
few cycles near the end of an inspiral waveform.
This continues as long as there is a local minimum in the eﬀective potential of the system. For
the Schwarzschild potential, which is the unique potential outside a spherically symmetric mass
distribution, there exists a local minimum for a test particle in a circular orbit around the object
down to a separation of 6M , called the innermost stable circular orbit (ISCO). At this point theabove derivation breaks down and we terminate the calculation of the waveform. This occurs at a(tc − t) value of

405M
,
16η

(2.37)

1
fISCO = √.
6 6πM

(2.38)(tc − t) =and a frequency of

The bodies will then plunge together, merging, leaving a final perturbed black hole, which relaxes
through quasi-normal mode oscillations. These final stages of the coalescence we call the “merger”
and “ringdown.”

16
2.4.1.2

Inspiral Waveform Approximants

The above derivation was a rather simple one using only the leading-order energy-loss terms and thequadrupole approximation. All of these arguments have been expanded in Post-Newtonian (PN)
theory where matching between the near-zone gravitational field and the wave-zone gravitational
field is done. In various studies diﬀerent assumptions and expansions are used to calculate thesewaveforms to varying orders in v/c, where an nPN order calculation is a controlled Taylor expansion
up to (v/c)n . The following are a few examples: reference [3] computes waveforms to Newtonianorder in amplitude and 2PN order in phase, references [4, 5, 6, 7] use the Eﬀective One Body (EOB)
approximation for the Hamiltonian to compute waveforms to Newtonian order in amplitude and
3PN order in phase, references [8, 9] use the Padé resummation technique to compute waveforms to
Newtonian order in amplitude and 3.5PN order in phase. All of the previous waveforms ignored spin
eﬀects associated with angular momentum of the bodies. The SpinTaylor approximant is computed
to Newtonian order in amplitude and 3.5PN order in phase using formulae from references [10] andbased on references [11, 12, 3, 13, 14, 15, 8, 16].
An example waveform constructed at Newtonian order in amplitude and 2PN order in phase isgiven similarly as above, however, in order to simply the expressions, we introduce a dimensionlesstime variable
Θ=

� η

5M

�1/8
(tc − t).

(2.39)In terms of this new variable, the orbital phase is given by���
13715 55
3π 2
Φ(t) = φ0 −
Θ5 ++ η Θ3 −Θη8064 964�� �
9 275 495
284 875
1855 2++η+η Θ1 .
14 450 688 258 048
2048(2.40)In this notation, equation (2.36b) is given by the first two terms on the right hand side of equation(2.40).

17
2.4.1.3

Inspiral Waveform Stationary Phase Approximation

When searching for these signals in GW detector data, it is convenient to convert these waveforms
to the frequency domain using the Stationary Phase Approximation (SPA)h̃(f ) =A(t) iΨ(f )
�
e
f˙(t)

(2.41)= A(f )eiΨ(f ) .

(2.42)In this approximation, the Newtonian amplitude A(f ) is given as
A(f ) =

�

5
24

�

M5/6
2/3
π Df 7/6

�,(2.43)and the 2PN phase Ψ(f ) is given in terms of a dimensionless frequency variable,

1/3τ = (πM f ),(2.44)

as���
33715 55
Ψ(f ) = 2πf tc + φ0 −
τ −5 ++ η τ −3 − 16πτ −2
128η
756
9
�
�
�
15 293 365 27 145
3085 2++η+η τ −1 .
508 032
504
72

2.4.2

(2.45)

Other Types of Gravitational Waves

Additional examples of transient situations in which the gravitational waveform is well modelled
include GWs from cusps traveling along cosmic strings, hyperbolic encounters between compact
objects, and ringdowns from perturbed black holes. Since the waveform is well modeled, GW signals
from these transients can be searched for using matched filtering, which we discuss in chapter 5, aswas used in references [17, 18, 19, 20, 21, 22, 23].
We call transient GW events “bursts” when their waveform is not well modeled. Examples of

18
situations generating these types of signals include GWs from supernova and neutron star glitches.
In order to search for these types of signals, since the waveform is unknown, GW data is searched fortransient excess power. In order to ensure that these transients are not just noise in the detectors,
these types of searches require either coherent or coincident signals to be seen in multiple detectors.
Examples of searches that have been performed from the LIGO Scientific Collaboration can be foundin [24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 23, 34].
Stochastic waves come from sources that are continuously emitting in either an incoherent fashionor in such a large number that the signals at the detectors are confused and inseparable. Examplesof these types of sources include relic GWs from inflation, GWs from Galactic binary white-dwarf
systems, and GWs from slow-spinning Galactic pulsars. These types of signals are searched for bylooking for correlations between multiple GW detectors. Examples of these searches can be foundin [35, 36, 37, 38, 39, 40, 41, 42, 43].
Continuous wave sources are are defined as signals that are quasimonochromatic and are separablefrom other sources. Such signals are commonly emitted by asymmetric spinning neutron stars, suchas pulsars. These signals are emitted from their source at a constant frequency, however when theydetected, they have been modified by Doppler shifts from motion due to the Earth’s rotation, theEarth’s orbit around the sun, the orbit of the source (e.g., when the source is part of a binary system),
etc. Searches for these types of signals then correct for these shifts, coherently or incoherently sum
the data over larger periods of time, and also correlate data from multiple detectors in order toreduce noise. Examples can be found in [44, 45, 46, 47, 48].
The remainder of this thesis focuses on GWs from CBCs.

19

Chapter 3

Introduction to LIGO
In this chapter we describe the history of GW detectors and the design of the LIGO GW detectors,
which have been built for the detection of GWs. This description is broken up into several parts:
a brief history of GW detectors in given in section 3.1, the design of the detectors is described in
section 3.2, the dominant noise sources for LIGO are described in section 3.3, the calibration of
LIGO is described in section 3.4, the response of the detectors to GWs is given in section 3.5, adescription of using multiple detectors in a global network is given in section 3.6, and a descriptionof previous LIGO science runs are given in section 3.7.

3.1

History of Gravitational-Wave Detectors

There exists a rich, and somewhat controversial, history of searching for GWs [49]. In the late 1950s
and early 1960s, Weber designed and built the first GW detector in the form of a large, suspended,
metal bar [50]. This bar was fit with piezoelectric transducers around its center to measure theoscillations of this driven harmonic oscillator. This bar was isolated from the surrounding environment using a combination of seismic isolation suspensions and a vacuum chamber. As GWs passed
through this bar, the dynamical Riemann tensor would induce oscillations of the bar’s fundamental
mode, giving it the name of a resonant mass detector. These oscillations would then be detectedby the transducers and, looking at the fundamental mode oscillations using a Fourier Transform,
recorded. Weber’s group went on to build several additional detectors so that they could look for coincident triggers (i.e., times when the detectors exhibited excess noise at their resonance frequency).

20
Results from these searches can be found in Refs. [51, 52] where they find rates of coincident triggers
well above the background.
After Weber built this detector and reported data from it, several groups around the world built
their own detectors. These group (with detectors located in e.g., Argonne, Glasgow, Moscow, ISR,
IBM, Reading, Rutherford, Tokyo, Munich, Frascati, Stanford, Bell Labs, and Rochester [53]) thenwent on to publish additional papers about the rate of coincident triggers not being consistent with
Weber’s observations.
The next set of advanced detectors were the ultralow temperature bars Allegro [54], AURIGA[55], Explorer [56], NAUTILUS [57], and Niobe [58]. These detectors further suppressed the noisein the detectors by reducing the temperature of the detectors to a few Kelvin using refrigeration,
allowing detector noise temperatures of a few to submilliKelvin levels to be reached. Searches for
GWs with these detectors also could not measure a rate above their background.
Another design for resonant mass detectors, a spherical mass rather than a cylindrical bar, wasconceived by Forward in Ref. [59]. If we look at the antenna response of the diﬀerent resonant mass
detectors, we see that the bars are sensitive to strains parallel to their cylindrical axis with the
response tensor
Rab = uaX ubX ,(3.1)

where uaX is the unit vector along the cylindrical axis of the bar. This implies that these detectors
will be relatively insensitive to GWs traveling along this axis as GWs are transverse(-traceless) (TT)waves (see chapter 2).
However, this is not the case for spherical resonant detectors. For spherical detectors, the response
tensor would be dependent on the number and locations of transducers, with their locations not being
limited to the locations conducive to measuring the fundamental mode. Since GWs are TTwaves,a basis of five tensors is all that is needed in order to determine all directions, polarizations, and21
magnitudes of incident GWs [59]. One independent set given by Forward is1
Aab = √ (uaX ubX − uaY ubY ) ,2
1
B ab = √ (−uaX ubX − uaY ubY + uaZ ubZ ) ,6
1
C ab = √ (uaX ubY + uaY ubX ) ,
2
1 a b
ab
D = √ (uX uZ + uaZ ubX ) ,2
1
E ab = √ (uaY ubZ + uaZ ubY ) .
2(3.2a)(3.2b)(3.2c)(3.2d)(3.2e)

Several detectors of this sort are being pursued around the world. These include the MiniGRAIL
project [60] in Leiden University, the Netherlands and the “Mario Schenberg” GW detector [61] in
the University of São Paulo, Brazil.
A diﬀerent detector technology that has also been designed and built is the interferometer detector. This type of detector, in its most basic form, consists of a Michelson Interferometer thatmeasures diﬀerential length changes along its two arms. The advantage of this type of detector isthat, as one is measuring actual length changes instead of the amplitude of oscillations on resonanceof a harmonic oscillator, the bandwidth of the detector can be much larger than that of the resonantmass detectors. If the arms are orthogonal and oriented along the X and Y directions respectively,
this type of detector has the response tensor

Rab = uaX ubX − uaY ubY .

(3.3)

A meter size version of this detector was first used for this purpose by Forward in the 1970s [62], ofcomparable sensitivity to the bar detectors operating at that time. Around the same time, Weiss
analyzed the noise sources of such a detector and developed the first ideas of kilometer-scale devices[63]. This inspiration is what has led to the modern day construction of TAMA300 [64], GEO600[65], Virgo [66], and the three LIGO detectors [67], all kilometer-scale GW detectors.
Resonant mass detectors are typically sensitive to GWs in the frequency range of 103 –104 Hz,

22
and ground-based interferometer detectors are sensitive to GWs in the frequency range of 101 –104 Hz. Other techniques are also currently being pursued to detect GWs in diﬀerent frequency
regions. These include space-based interferometers, such as LISA [68] and DECIGO [69], designed
to be sensitive to GWs with frequencies between 10−3 –10−2 Hz and 10−1 –101 Hz respectively,
pulsar timing experiments (first described in Refs. [70, 71]) aimed at detecting very low frequency
GWs in the range of 10−9 –10−7 Hz, and Cosmic Microwave Background polarization measurements[72, 73, 74, 75, 76, 77, 78, 79, 80, 81] designed to detect GWs in the extremely-low frequency bandof 10−18 –10−15 Hz.

3.2

Design of LIGO

The LIGO detectors were designed to detect diﬀerential strains from incident GWs, based on asimple Michelson Interferometer (see figure 3.1). This design starts with an input laser incident on
a beam splitter located at the origin with the face at a 45◦ angle between the x- and y-axes. The
beam splitter directs 50% of the light in two perpendicular directions, x and y. The light travels
toward mirrors located at the end of the arms, (Lx , 0) and (0, Ly ), which redirect the light backtoward the beam splitter. The phase the light accumulates as it travels from the beam splitter tothe end mirror is given by
φ=

� L
0

2π
2πL
dx =
,
λ
λ

(3.4)

where L is the distance traveled and λ is the wavelength of the light. The total phase of one roundtrip will then be the sum of the integral down and the integral back.
When the light returns to the beam splitter, the two beams interfere with each other in such a
way that, depending on the diﬀerential arm length, the light is either directed back toward the laser(the symmetric port), or toward a photodetector (the antisymmetric port), or some combination of
the two. The diﬀerential arm length Lx − Ly can then be measured, modulo one wavelength, by23
looking at the power of the light coming from the antisymmetric port, PAS . This is given byPAS= P0 sin2 (φx − φy )��
4π
2
= P0 sin(Lx − Ly ) ,
λ

(3.5)

where P0 is the input laser power, φx is the phase accumulated for the beam in the x-arm, andφy is the phase accumulated for the beam in the y-arm. This shows that PAS is dependent on thediﬀerential arm length.
Let us assume that the arms are of the same length L, and an incident GW with amplitude his stretching the detector along the x-axis and squeezing it along the y-axis. We see from chapter 2that this eﬀect causes the distance the light needs to travel to be given by L� ≈ L(1 ± h/2), giving
φx and φy as

φx

=
≈

� L�

� 0
2π
2π
dx −
dx
λ
� λ
0
�
�L
4πL
h
1+
,
λ
2

4πL
φy ≈
λ

�

h
1−
2

�

(3.6a)(3.6b).

With these phase shifts, we find that PAS is then given by�hL
PAS = P0 sin 4π
λ

�2

,

(3.7)

such that for the same change in PAS , the strain h is inversely proportional the length L of thearms. However, although there is a signal at the antisymmetric port for this interferometer, this is apoor way to detect GWs. For one, since the above power at the antisymmetric port is proportional
to sin2 (∆φ), which can be Taylor expanded about zero as ∆φ2 , the signal we are measuring isquadratically proportional to the small GW signal we are trying to sense. In addition, there is no
information on which direction we need to move the mirrors to bring the signal back to the dark

24fringe.
In order to get around these problems, LIGO has adopted the use of the “heterodyne detection”
technique, which combines a phase-modulated input laser with a Schnupp Asymmetry. In this setup,
before the laser light enters the interferometer, the light is phase modulated with a frequency Ω suchthat instead of all the light entering the detector with one frequency, there are eﬀectively three
superposed laser beams entering with three diﬀerent frequencies and electric field amplitudes: theoriginal beam with frequency ω and amplitude E0 , also known as the carrier, and two sidebands
with frequency ω ± Ω and amplitude E1 .
When the arm lengths are perfectly equal, all of these beams will come out at the symmetric port,
which results in the same answer as above. However, there is a diﬀerent response when a Schnupp
Asymmetry is introduced between the arms. This occurs when the arms of the interferometer are
unequal by an amount that is an integer number of wavelengths of the carrier light but not an integerof wavelengths for the sidebands. In this situation, in the absence of a GW signal, the carrier lightwill come out the symmetric port while the sidebands will be transmitted to the antisymmetric port.
When there is a GW signal present, a phase shift will be introduced on all three beams such that a
mixture of all three beams will exit the antisymmetric port. Processing the signal coming from theantisymmetric port with demodulation results in an error signal that is linear in the GW amplitude,
just as is needed.
In order to increase the sensitivity of the detector, LIGO uses two modifications to the standard
Michelson interferometer. The first is the use of Fabry-Perot cavities in the two arms in order toincrease the storage time of light within the arms, eﬀectively increasing the length of the arms, andthus the phase shift on the light, by the same factor. The second is to boost the light circulating
within the interferometer for a given laser input power with the addition of a power-recycling mirror. A schematic for the design of LIGO, including all of the additions to the standard Michelsoninterferometer, can be seen in figure 3.1.

25















Figure 3.1: LIGO Schematic
A simple schematic of the LIGO detectors. The mirror labels are as follows: PRM (power recycling
mirror), BS (beam splitter), ITMX/Y (x-/y- arm input test mass), ETMX/Y (x-/y- arm end testmass), AS (antisymmetric port). The width of the laser beam denotes power in diﬀerent portionsof the detector. A simple Michelson interferometer includes the mirrors BS and ETMX/Y. The
addition of the ITMX/Y mirrors makes use of Fabry-Perot cavities for the arms. The use of the
PRM boosts the power circulating in the whole interferometer.

26

3.3

Dominant Noise Sources

The design sensitivity of the LIGO detectors is limited by several fundamental noise sources fordiﬀerent frequency bands [82]: Seismic Noise for f < 40 Hz, suspension thermal noise for 40 Hz <f < 150 Hz, and photon shot noise for f > 150 Hz.
Seismic Noise is the result of ground motion coupling to the motion of the mirrors. This typeof noise is suppressed in multiple ways, both actively and passively. Passive isolation is achieved bysuspending the mirrors as pendula, on seismic isolation stacks. This suspension acts as a low passfilter for motion coming from the ground, resulting in a pendulum response function for the mirror
motion that dies as 1/f 2 for noise above the resonant frequency. Active isolation is also used bysensing ground motion and feeding this back into the controls of the Hydraulic External Pre-Isolator
for cancellation.
Suspension Thermal Noise results from the random thermal motions of the room-temperaturesuspension feeding into motions of the mirrors. This motion has been calculated by Saulson in
Ref. [83].
Photon Shot Noise results from a finite number of photons arriving each second, informing usof the motion of the mirrors. Since we are limited in the number of photons that arrive at thephotodetector per second by the power of the laser light, we end up not being able to detect strains
better than [82]∆h[Hz]−1/2 =

1
L

�

hcλ
π2 P

�1/2

,

(3.8)

where h is Planck’s constant, c is the speed of light, λ is the wavelength of the laser light, and P
is the power of the laser light. The arm cavity introduces an enhancement in P at the expense of
reduced sensitivity at high frequencies due to the cavity pole. This results in a shot-noise equivalentstrain that rises with f above the cavity pole frequency of ∼90 Hz.
These are not the only noise sources in the LIGO detectors. figure 3.2 shows all of the known
noise sources during the S5 run [84]. There the noise in channel that carries the GW signal (DARM)
is compared to the sum of all the known noise sources (Total). Above 100 Hz, the dominant noise

27

Figure 3.2: LIGO Noise Budget
The noise budget for the 4 km instrument at Hanford from Ref. [84]. The line labeled “Total” isthe sum of all known noise sources. The black line (measured DARM noise) should be compared tothis.
source is shot noise and DARM agrees very well with Total. Between 50 and 100 Hz there are noise
sources that have not been accounted for in addition to the expected dominant source, which issuspension thermal noise. Below 50 Hz the noise is mainly due to man-made activities and is highly
nonstationary, thus it is hard to estimate at any given point in time.

28

Figure 3.3: Model Control Loop
A simple model for the control loop of the LIGO length sensing and control system [85].

3.4Calibration

In this setup, the light incident on the antisymmetric port’s photodetector carries the informationof the strain the detector is experiencing. However, in order to keep the detectors in a stable datataking mode, this signal is fed back into detector such that the light arriving from the antisymmetricport is kept roughly constant. This means that the strain of the detector must be reconstructed fromthe signal arriving from the antisymmetric port, known as the error signal q, and the transformations
made on q to control the detector, what we call the response function R. This reconstruction iseasily understood in the frequency domain as [85]h(f ) = R(f )q(f ) .

(3.9)

The feedback loop diagram for this control system is shown in figure 3.3. In the frequency

29
domain, the relationships depicted in the diagram take the simple form

xr (f ) = h(f ) + xc (f ) ,(3.10a)q(f ) = α(t)C0 (f )xr (f ) ,(3.10b)d(f ) = β(t)D0 (f )q(f ) ,(3.10c)xc (f ) = A(f )d(f ) ,(3.10d)where C0 (f ) is a reference sensing function measured at some reference time whose frequency dependence is determined primarily by the Fabry-Perot arm cavities and is roughly (1 − f /fp )−1 withfp ≈ 90 Hz and α(t) is a real overall gain that depends on the light power stored in the arms and ismeasured every ∼ 10 seconds, A is an actuation function that is largely determined by the pendulum
suspension of the mirrors, and D0 is the feedback filter and β(t) is a real overall gain.
From these relations, R(f ) is found to beR(f ) =

1 + α(t)β(t)G0 (f ),α(t)C0 (f )(3.11)

where G0 (f ) is the reference open loop gain given by

G0 (f ) = A(f )C0 (f )D0 (f ) ,(3.12)

which has a unity gain frequency of ∼ 150 Hz. The time-dependent gain factors α(t) and β(t) aremeasured by tracking the amplitude of calibration lines that are injected into the detectors at around
1 kHz, 150 Hz, and a few tens of Hz.
This procedure tells us how to calibrate the error signal from the detector in order to obtain thestrain signal in the frequency domain. However, for LIGO’s fifth science run (S5), an additional step
was added to the calibration [85], giving the strain signal in the time domain h(t). This is the main
method used for analyses of the S5 data.

30

3.5

Interferometer Antenna Response

In chapter 2, we described the eﬀect of incident GWs on test particles perpendicular to the directionof propagation of the waves. We noted that interferometric detectors have a response tensor, in thereference frame of the beam splitter, of the form


 0


 0

αβ
R =

 0



0

0

0

1

0

0

−1

0

0



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

(3.13)The resulting response of the detector is given by∆l1= hαβ Rαβ .
l2(3.14)

GWs from arbitrary directions will have diﬀerent projections onto this response tensor. Here we
describe those projections as the antenna response of the detector to GWs incident from directions
given by θ and φ, and the polarization angle ψ, which bring equation (3.14) into the form∆l
= F+ h + + F× h × ,l(3.15)

where F+ and F× are the antenna factors we are evaluating. The above angles are identified infigure 3.4.
In this derivation we start with the assumption that the wavelength of the GW is much longer
than the storage time of the light in the interferometer’s arms times the speed of light λ � cT .

31

Figure 3.4: Euler Angles
The Euler angles {φ, θ, ψ} used in converting from the GW propagation frame {X, Y, Z} to thedetector frame {x, y, z}.
In the radiation frame, denoted with “� ”s above the indices, the GW hα� β � is given by


 0


 0

hα� β � = 

 0



0

0

0

h+

h×

h×

h+

0

0



0 


0 

 ,

0 



0

(3.16)

which can be split intohα� β � = Sα+� β � h+ + Sα×� β � h× .

(3.17)To transform this into frame of the detector (i.e., the frame of the beam splitter) we need totransform both indices of hα� β � . This transformation can be seen as the combination of diﬀerent

32
Euler angle rotations with the combined rotation matrix given by


�

α
Rα

 0


 0

= 

 0



0

 0


 0

= 

 0



0

0

0

− cos θ cos φ cos ψ − sin φ sin ψ− cos θ cos φ sin ψ + sin φ cos ψ

cos θ sin φ cos ψ − cos φ sin ψ

cos θ sin φ sin ψ + cos φ cos ψ− sin θ cos ψ− sin θ sin ψ

0

0

A

B

D

E

G

H

0

0 


C 

 .

F 



I






sin θ cos φ − sin θ sin φ − cos θ(3.18)

Applying this to B+ and B× , we find


 0


 0

+
Sαβ
=

 0



0


 0


 0

×
Sαβ
=

 0



0

0

0

A2 − D 2

AB − DE202AB − DE

B −E

AC − DFBC − EF002AD

AE + BD

BD + AE

2EB

CD + AF

CE + BFAC − DF  ,
BC − EF 



2
2
C −F
0

(3.19a)AF + CD  .

BF + CE 



2CF

(3.19b)+×
We can now calculate the antenna factors where F+ = (1/2)SαβRαβ and F× = (1/2)SαβRαβresulting in
F+ =
F× =

�
1�
1 + cos2 θ cos 2φ cos 2ψ + cos θ sin 2φ sin 2ψ ,2

�1�
1 + cos2 θ cos 2φ sin 2ψ + cos θ sin 2φ cos 2ψ .
2

(3.20a)(3.20b)

We can see the structure of these antenna factors when plotted with ψ = 0 in figure 3.5. From

33
these we find that there are several regions of the sky where the detectors are sensitive to neither
polarization, which we call nodes.
There is one additional angle that will aﬀect the amplitude of the signal that a interferometer
detector will receive from a quasi-circular inspiralling compact binary, the angle of inclination ι
between the orbital plane and vector connecting the origin of the detector and the origin of thesource. This will have an eﬀect on resulting signal seen by the detector h(t) such that��h(t) = h+ (t)F+ 1 + cos2 ι + h× (t)F× (2 cos ι) .

(3.21)

34

Figure 3.5: Detector Antenna Patterns
The antenna patterns for interferometer detectors for plus (top left), cross (top right), andcombined (bottom) polarizations with polarization angle ψ = 0. For this figure, the detectors arms
are oriented such that the first node in the bottom plot is found at 45◦ between the two arms ofthe detector and subsequent nodes occur every 90◦ as you rotate around the vertical axis.

35
Table 3.1: Combined Horizon Volumes
Detector Combination Volume (Mpc3 )H1H2
2.13 × 103
H2L1
2.11 × 103
H2V1
1.81 × 103
H2L1V1
8.84 × 103
H1L1
1.23 × 104
H1H2L1
H1V1
8.68 × 103
H1H2V1H1L1V1
1.54 × 104
H1H2L1V1
The combined horizon volume to which each detector combination is sensitive. When calculating
this, we compute the distance at which each detector would see a binary neutron star signal withinclination angle ι at an SNR of 8. For each point in the sky and each polarization ψ, we take thedistance from the second most sensitive detector (i.e. only double coicident observations are
required), volume-averaged over ψ and ι, and then integrate over the sky. These calculations aredone using the design sensitivity of each detector.

3.6

Networks of Detectors

As was discussed above, coincidence can be used to verify a GW signal when it is seen by multiple
detectors. However, this can only be done when at least two detectors see a signal. One common
way to quantify the sensitivity of a detector is by looking at the distance at which a given detector
can see a optimally located and oriented binary neutron star signal (i.e., an overhead, face on orbit;
θ = φ = ψ = ι = 0) with an SNR of 8. This is referred to as the horizon distance. Non-optimallylocated and oriented sources can be detected with SNR = 8 at shorter distances. Table 3.1 shows thevolumes of space that are screened after averaging over all of these angles when diﬀerent combinationsof interferometric GW detectors are operating. The angle-averaged sky maps associated with eachdetector and combination of detectors can be seen in figure 3.6 and 3.7 respectively, which are
Mollweide projections of the sky. As additional detectors are built around the world, we will be able
to more completely cover the sky.
Just as combining the signals from diﬀerent transducers on a spherical resonant mass detector
allows us to determine the sky location and polarization of a given signal, coherently combining
the output from multiple, misaligned interferometric detectors accomplishes the same thing [86].
Another advantage of coherent combinations of signals is the amplification of the signal relative to36

Figure 3.6: Individual Detector Angle-Average Sky Maps
The Hanford (top left), Livingston (top right), and Virgo (bottom) angle-averaged sky maps fortheir design sensitivities. The color scale shows the average distance at which each detector can see
a binary neutron star signal with an SNR of 8, volume-averaged over the polarization andinclination angles.
the noise in the detectors. Because the noise in the diﬀerent detectors is independent random noise,when the signals are combined, the random fluctuations are suppressed though incoherent addition
while external GW signals add coherently.
A final advantage of coherent analysis is the overdetermination of the signal from multiple data
streams. When at least three detectors are used, since there are only two independent polarizationsof a GW signal, the data from the detectors can be combined in such a way that the signal should
be cancelled out. This type of combination, called a null-combination, can be used to identify anaccidental coincidence of noise triggers in the detector network since this will cause excess energy inall the detectors that is incoherent.

37

Figure 3.7: Combined Detector Angle-Average Sky Maps
The Hanford-Livingston (top left), Hanford-Virgo (top right), Livingston-Virgo (bottom left), andHanford-Livingston-Virgo (bottom right) combined angle-averaged sky maps for their design
sensitivities. The color scale shows the average distance at which each that detector combination
can see a binary neutron star signal with an SNR of 8, volume-averaged over the polarization andinclination angles.

38

Figure 3.8: LIGO Science Runs Timeline
A timeline showing the progress made in the LIGO detectors strain sensitivity as well as thediﬀerent periods of science data taking [84].

3.7

LIGO’s Past and Future

After the LIGO detectors were built, there were several times when science data was taken in between
commissioning periods on their way toward their design sensitivity (figure 3.8). These science runs
(S1–S4) were analyzed for many diﬀerent types of signals [17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28,
29, 30, 31, 32, 33, 23, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48]. One of the important
aspects of these runs was finding times during which the detectors were not operating properly.
During these times, many auxiliary channels of the detectors were recorded and analyzed, leading to
the discovery of diﬀerent ways the interferometers glitch and the suppression of such glitches during
the next commissioning break.
In the upcoming months, the LIGO detectors will finish being upgraded from their initial design(Initial LIGO) to their enhanced design (called enhanced LIGO) and will start taking data for theirsixth science run with an expected factor of two increase in sensitivity over S5. A few years later,
the detectors will be taken oﬄine again to undergo a more extensive reworking, returning with theiradvanced design (Advanced LIGO) in place, which should further improve the sensitivity by anorder of magnitude with respect to Initial LIGO.

39

Chapter 4

Astrophysics of Compact Binaries
This chapter is broken up into two sections. The first (section 4.1) presents evidence for the existence of compact objects (i.e., neutron stars and black holes). The second (section 4.2) presents
predicted rates for the coalescence of binary compact objects based on observations and theoreticalconsiderations.

4.1

Evidence for Neutron Stars and Black Holes

To date, there have been many observations of compact objects, which confirm the existence ofneutron stars and strongly suggest the existence black holes. The theoretical maximum mass foran electron-degenerate star is given by the Chandrasekhar mass limit of ∼1.4 M⊙ . This limit is
implicated in multiple types of stellar explosions, including Type Ia/b/c and II supernovae, in which
either the mass of an electron-degenerate, white dwarf approaches this limit (Type Ia), or the coreof a massive a star approaches this limit (Type 1b/c II). For the latter situation, the observationsof the remnant object has been linked to compact object observations, as in the case of Cassiopeia
A [87].
Compact objects with a tightly orbiting companion are good testing grounds for diﬀerentiating
between neutron stars and stellar mass black holes. If the companion object is close enough, it fills
its Roche Lobe, the boundary at which the gas of the companion object is no longer gravitationallybound to its parent star, and leaks onto the compact object. This gas then forms an accretion disk
around the compact object as it is transferred. By measuring properties of the binary system as well

40
as of the accretion disk, we can learn many things about the compact object. Studying the orbitalparameters of the binary will tell us the mass of both the compact object and companion object. If
the mass of the compact object is larger than the Chandrasekhar Limit, then the compact object
must be either a neutron star or a black hole. If thermonuclear explosions are seen as matter falls
onto the compact object, it is determined to be a neutron star. In the case of black holes, studying
the profile of Fe-lines in the accretion disk can allow the spin of the black hole to be determined [88].
In addition, observation of binary pulsars such the Hulse-Taylor pulsar [89] have also yielded
evidence for compact objects. The two objects, roughly 1.44 and 1.39 M⊙ , in this system areorbiting each other once every 8 hours, corresponding to a separation of ∼1 R⊙ . Additionally, thissystem is seen to lose orbital energy at the rate predicted by GW emission in general relativity.
At this separation, if the objects were ordinary stars, there should be some other much stronger
force causing them to depart from this prediction. Electromagnetic observations have also found nocounterpart, which at a distance of 8 kpc, excludes these objects from being white dwarfs.

4.2

Predicted Compact Binary Coalescence Rates

In this section we summarize a number of diﬀerent ways the rate of CBC are predicted. Thesemethods are broken up into three categories, namely extrapolations from merging binary neutron star
(BNS) observations (section 4.2.1), population synthesis of field binaries (section 4.2.2), dynamical
simulations of star clusters (section 4.2.3), and extrapolations from short Gamma Ray Burst (GRB)
observations (section 4.2.4) [90].

4.2.1

Extrapolations from Merging Binary-Neutron-Star Observations

Observations of binary systems involving pulsars that will coalesce within a Hubble time can beextrapolated in order to obtain a BNS merger rate. There are currently four known systems ofbinary pulsars that will merge within a Hubble time, and one more possible system, all in ourgalaxy. These systems are PSR B1913+16 (the Hulse-Taylor pulsar), PSR B1534+12, J0737-3039A,
J1756-2251, and possibly J1906+0746.

41
In [91, 92], the authors use these binary pulsar systems to calculate a combined probability on
the rate of coalescences in our galaxy. To do this they combine the number Ntot of observed Galactic
pulsars, the fraction of Galactic pulsars with pulse and orbital characteristics similar to those of aparticular type, the lifetime τlife of each of the observed systems, and the upward correction factor
fb−1 due to pulsar beaming. For each type of pulsar population, with the ratio α of the mean numberof observed pulsars �Nobs � to the total number of pulsars Ntot for that populationα=�Nobs �
,
Ntot

(4.1)

they obtain a probability density function P (R) for the rate given by

2

P (R) = (ατlife fb ) Re−(ατlife fb )R .

(4.2)

These probability distributions can be combined (Appendix A of [93]) to obtain the most probablerate of R ≈ 71 MWEG−1 Myr−1 and a 90% confidence interval of ∼15–240 MWEG−1 Myr−1 [94]without the inclusion of J1906+0746 and a factor of 2 larger with its inclusion, where an MWEG is
a Milky Way Equivalent Galaxy. This can be converted to 7.1 × 10−1 Mpc−3 Myr−1 using a galaxy
number density of 10−2 MWEG Mpc−3 .

4.2.2

Population Synthesis of Field Binaries

Population synthesis simulations have been used to estimate the rate of binary coalescences for
BNS, binary black-hole neutron-star (BHNS), and binary black hole (BBH) systems in the binary
evolution scenario (i.e., stars formed as binaries during their stellar formation). These simulations
start with stars distributed according to an initial mass function. These systems are then evolved
keeping track of their evolutionary details. At the end of the simulation, observational constraints
can be imposed on the results such that the results are consistent with the observed sample ofmerging Galactic BNS, wide Galactic BNS, white dwarf–neutron star binaries, the observed rate ofType Ib/c, Type II supernovae, etc.

42With these constraints imposed, the final rates can be can be obtained by looking at the number of BNS, BHNS, and BBH with orbits tight enough to merge within a Hubble time. Results
suggest that the range of merger rates are 1–103 MWEG−1 Myr−1 for BNS [91, 92], 5 × 10−2 –1 × 102 MWEG−1 Myr−1 for BHNS [95], and 4 × 10−2 –1 × 102 MWEG−1 Myr−1 for BBH [96]systems.

4.2.3

Dynamical Simulations of Star Clusters

In this section we summarize arguments of [97, 98, 99, 100] for predicting the CBC rate from diﬀerent
types of star clusters. Star clusters may provide the necessary breeding ground for BBH due to theirincreased star formation rate and their increased density. Studies have shown [101] that in suchan environment, mass segregation occurs through dynamical interactions within the cluster, driving
the higher mass objects toward the center, increasing the possibility of forming stellar mass BBH
systems.
For a cluster of mass Mcl , since stars form according to a power-law mass distribution [102], thenumber of stars that are massive enough to form black holes after supernova (i.e., those with masses
greater than 20 M⊙ ) is roughly 3 × 10−3 (Mcl /M⊙ ). These objects undergo rapid evolution forming
their black holes quickly (on a timescale of about tSN ∼ (M )−2.5 × 1010 yr). Then, since these are
the most massive objects, they tend to sink toward the center with a timescale oftrelax � tcross × 0.1N/ ln N ,(4.3)

where tcross is the typical crossing time of the cluster [103]. This happens through a mechanism
known as equipartition, which is a statistical tendency for objects undergoing two-body interactions
to equilibrate their kinetic energy. This tendency causes the lower mass object to leave a two-bodyinteraction with a larger average speed than the higher mass object, thus causing the higher mass
object to sink further in the potential well of the cluster.
Another eﬀect that occurs in the mass segregation process is the accumulation of high mass

43
binaries. When a binary interacts with a third object, the two most massive objects tend to leave
the interaction as a binary [104]. In these interactions, the binding energy of the binary tends toincrease. These interactions thus favor the formation of tight binaries from the highest mass objects
in the cluster.
The combination of mass segregation and binary exchanges continue until the recoil speed fromthe three-body interactions is large enough such that even the BBH systems receive a recoil larger
than the escape velocity of the cluster [105, 106]. At that point, the binary has hardened enoughthat even though it leaves the cluster, the time until merger due to gravitational radiation is lessthan a Hubble time.
These arguments have been applied to globular clusters in [98]. They find that the rate ofmergers from globular clusters is given by ∼gcl gevap Mpc−3 Myr−1 , where gcl is the fraction of total
star formation that occurs in clusters and gevap is the fraction of cluster-forming mass that possesses
the birth conditions necessary for this process to occur.
These arguments have also been applied to nuclear star clusters both in the presence of a supermassive black hole (SMBH) [99], and not [100]. In [99], they found a wide range of merger rates in
the presence of a SMBH with rates varying between 1.5 × 10−6 and 2 × 10−4 Myr−1 per galaxy (notMilky Way equivalent galaxy) depending on the model chosen. In [100], they argue for merger rates
in the absence of a SMBH of a few times 10−3 Mpc−3 Myr−1 .

4.2.4

Extrapolations from Short-Gamma-Ray-Burst Observations

Here we present arguments used in [107] for estimating the rate of mergers of dynamically formed
BNS and BHNS systems in globular clusters from short GRB observations. Short GRBs are theorized
to originate in the merger of two compact objects forming a black hole surrounded by an accretion
disk [108, 109, 110, 111]. In order for an accretion disk to form, there must be matter present in thesystem from either one or both objects, which excludes the possibility of short GRB coming from
BBH systems.
BNS and BHNS systems can form in one of two ways. The objects were originally binary stars

44
that both underwent core collapse forming a binary of compact objects, or the objects separately
underwent core collapse and then dynamically formed through companion capture and possible
exchange interactions. The former we will refer to as “primordial” systems, which can be found in
the field. Results from population synthesis have shown that primordial systems that will merge
within a Hubble time do so soon after their formation, thus the redshift associated with such systems
should closely follow the star formation rate of the universe [112, 113, 114]. On the other hand,
dynamically formed systems will be delayed by the relaxation time trelax of the cluster, which canbe on the order of a Hubble time.
Using this time delay, the short GRB rate can be calculated from the star formation rate history
using
RGRB (z) ∝� t(z)0

RSFR (t − τ )P (τ )dτ ,(4.4)

where RSFR (t) is the star formation rate at time t, P (τ ) is the distribution of time delays τ before
merger, which goes as P (τ ) ∼ 1/τ for primordial systems [115, 116] while it increases for increasingtime delays for dynamically formed systems [117].
Using equation (4.4) and the distribution of observed luminosities for short GRB, [107] calculate
the local rate of events to be 1.3 × 10−3 and 4.0 × 10−3 Mpc−3 Myr−1 assuming all observed events
come from primordial or dynamically formed systems respectively. Given the observed distributionof redshifts associated with short GRB, the best fit of the data comes from a 60% contribution fromdynamically formed systems. Reducing the contribution from dynamically formed systems suchthat the Kolmogorov-Smirnov probability that the observed distribution comes from the expecteddistribution reaches a value of 0.1, the contribution from dynamically formed systems is foundto be at least 10%. These two combinations yield a local event rate of 2.9 × 10−3 and 1.6 ×
10−3 Mpc−3 Myr−1 for 60% and 10% fraction of dynamically formed mergers respectively.
The above rate calculations have only taken into account the observed rate, which ignores theeﬀects of beaming at the source. In [118] and [119], the authors have estimated the beaming factor to
be on the order of fb−1 ∼ 100, where fb is the fraction of the total solid angle within which the GRB isemitted. Taking this into account, the beaming-corrected rates are found to be 2.9(fb−1 /100) × 10−1

45
and 1.6(fb−1 /100) × 10−1 Mpc−3 Myr−1 for 60% and 10% fraction of dynamically formed mergers
respectively.

46

Chapter 5

Gravitational Wave Signal
Processing
In this chapter, we describe the matched filtering that is done in order to search for a signal of
known form within data that also contains noise, as is the case for analyzing GW data.
This problem has been well studied within the field of signal-processing, and the optimal tool for
extracting the signal from the noise has been found to be the Weiner Filter [120], also known as thematched or optimal filter.

5.1

The Optimal Filter

Let us assume we have a data stream s(t) that is the sum of a stationary, Gaussian noise n(t) anda signal h(t)s(t) = n(t) + h(t) .

(5.1)

Since n(t) is stationary, the mean of n(t) can be taken to zero by definingn(t) ≡ nraw (t) − �nraw � ,(5.2)

where the brackets denote the expectation value. Since this is true for any stationary series, fromhere we will assume n(t) has a mean of zero. The two-sided power spectral density (PSD) Sn (f ) of

47
this Gaussian noise is defined by�ñ(f )ñ(f � )� ≡ δ(f − f � )Sn (f ) .

(5.3)

Let us define a real-valued filter F (t) asA=�∞

F (t)a(t)dt =−∞�∞F̃ ∗ (f )ã(f )df ,

(5.4)

−∞

where A is the filtered value of a(t). Above we have used the convolution theorem and transformed
from the time domain to the frequency domain using the Fourier transform of a function, defined byã(f ) =�∞a(t)e−2πif t dt .

(5.5)

−∞

Using the above properties of the noise, we find

�

N2

�

=
=
=

� ∞� ∞−∞ −∞� ∞� ∞
−∞ −∞
� ∞�
−∞

F̃ ∗ (f )F̃ ∗ (f � ) �ñ(f )ñ(f � )�df � df
F̃ ∗ (f )F̃ ∗ (f � )δ(f − f � )Sn (f )df � df

�2���F̃ (f )� Sn (f )df ,(5.6)

where we have also used the property of real-valued functions that ã∗ (f ) = ã(−f ) in the last line.
We now wish to find the form of the filter F (t), or in this case F̃ (f ), which will optimally extract
the signal from the noise. We measure this by maximizing the ratio of filtered values2

H�N 2 �

=

=

��
�2
� ∞ ∗
�
� −∞ F̃ (f )h̃(f )df ��2
� ∞ ��
�F̃(f)
�� Sn (f )df
−∞
��
�2
�
� ∞ ∗ ��� −∞ F̃ (f ) Sn (f ) h̃(f )/ Sn (f )df �.
�2
� ∞ ��
�
F̃
(f
)
S
(f
)df
�
�
n
−∞

(5.7)

48
We can use the Cauchy-Schwarz inequality,
� ∞
�2
��
�
�∞
�∞
�
�
2
2
�
�
A(f )B(f )df � ≤|A(f )| df|B(f )| df ,
�
�
�
−∞

−∞

(5.8)

−∞

to argue that in order for the ratio of equation (5.7) to be maximized, the functions corresponding�to A(f ) and B(f ) must be equal up to a constant C, where A(f ) = F̃ ∗ (f ) Sn (f ) and B(f ) =�h̃(f )/ Sn (f ). Equation (5.7) then becomesH2
�N 2 �

C
=

= C

�

��
�
�2
�2
� ∞ ��
� ∞ ��
�
�
F̃
(f
)
S
(f
)df
h̃(f
)
/S
(f
)df
�
�
�
�
n
n
−∞
−∞
�
�
2
�∞ �
�
F̃ (f )� Sn (f )df
−∞ �

� ∞

h̃∗ (f )h̃(f )df .
Sn (f )−∞

(5.9)

From this exercise we find thatF̃ ∗ (f ) = Ch̃(f ),Sn (f )

(5.10)

which is the definition of the optimal filter for h(t) embedded in stationary, Gaussian noise.
Let us use the optimal filter to define an inner product of the data s with the template h:(h|s) =� ∞h̃∗ (f )s̃(f )df .
Sn (f )−∞

(5.11)

An interesting interpretation of equation (5.11) is to split the PSD Sn (f ) into two amplitude
spectral density (ASD) terms
rately,�Sn (f ) that can be associated with the data and the template sepa-(h|s) =� ∞h̃∗ (f )s̃(f )��
df .
Sn (f ) Sn (f )−∞(5.12)

This ends up weighting both the template and the data by the inverse of the ASD, which can be seen
as “whitening” both the template and the data. The reason this is known as whitening is because
for stationary, white, Gaussian noise the PSD is frequency independent resulting in an optimal filter

49of the form(h|s) =� ∞h̃∗ (f )s̃(f )df .

(5.13)−∞∗∗Sn (f )Sn (f )

If we define h̄(f ) = √h̃ (f ) and s̄(f ) = √s̃ (f ) , equation (5.12) comes into the same form as equation(5.13)(h|s) =� ∞h̄∗ (f )s̄(f )df .

(5.14)−∞

5.2

The Waveform Overlap

Now that we have an optimal filter for extracting a signal for noisy data, let us use it to define theoverlap M of two vectors a and b
M = (a� |b� ) ,(5.15)

where the “� ” denotes a normalization such thata,σa(5.16a)σa2 = (a|a) ,(5.16b)a� =andand similarly for b.
This definition of the overlap will take its largest value of 1 when a� and b� are the same function,
which can be seen by(a� |a� ) =
=

�

a a
|
σa σa
1
(a|a)
σa2

�

= 1.

Since we are interested in the strength of the signal that matches the template h in the data s,

50
we compute the signal-to-noise ratio (SNR) by normalizing the overlap with respect to the templatebut not the dataSNR =1(s|h) .
σh(5.17)The value of the SNR will then be proportional to the amplitude of the signal buried in the noise.
So far we have defined the overlap between s and h where we have assumed the template to beof the same length as our data. However, what we are actually interested in is the matched filter ofs using h where the data s(t) is an extended time series whose length is longer than the templatetime series h(t). The matched filter output as a time series is given by(s(t)|h) =� ∞s̃∗ h̃e−2πif t df .

(5.18)−∞For real-valued time series s and h, this filter outputs a real-valued time series that, in the absenceof signal, is χ2 distributed with one degree of freedom.

5.2.1

Matched Filtering for Compact Binary Signals

As derived in chapter 2, the GWs we are searching for come in two polarizations: the plus (+) andcross (×) polarizations. The actual signal seen by a detector is a combination of the two polarizations
that can be calculated using the antenna pattern of the detector and the parameters of the source.
Specifically, the signal seen by the detector h(t) is given byh(t) = F+ h+ (t) + F× h× (t) ,(5.19)

where h+ and h× are the plus and cross polarizations of the signal respectively, and F+ and F× are
the antenna patterns of the detector giving the sensitivity to the plus and cross polarizations of GW
signals respectively. Descriptions of these are given in chapter 3.
Due to the nature of GW signals from the inspiral phase of CBCs (see chapter 2 for a discussion), we find that the phase evolution of the cross polarization is 90◦ out of phase with the plus

51
polarization, given byh+ (t) = A(t) cos (φ(t)) ,(5.20)h× (t) = A(t) sin (φ(t)) ,(5.21)

where A(t) is the amplitude evolution of the signal and φ(t) is the phase evolution of the signal.
From this, we find a simple relation between the Fourier transform of the two polarizations:h̃+ = ih̃× .

(5.22)

When we filter s = (Xh+ /σh ) + (Y h× /σh ) with the template h+ , we obtain the matched-filter
real-time series z+
z+ = (s|h+ ) =XY(h+ |h+ ) +(h× |h+ ) .
σhσh(5.23)

Using equation (5.22), we find(h× |h+ ) =
=
=

� ∞

h̃∗× h+ + h̃∗+ h×
df
Sn
−∞
�
�∗
� ∞ ih̃+ h+ + h̃∗ ih+
+
−∞
� ∞

Sn

df

−ih̃∗+ h+ + h̃∗+ ih+
df
Sn
−∞

= 0,

which then impliesz+ = Xσh ,(5.24a)where we have used equation (5.16b) for the definition of σh . A similar procedure shows that when
we filter s with the template h× , we obtain the matched filter time series z×

z× = Y σ h .

(5.24b)

52

From equation (5.24) we construct our SNR ρ using the combination

ρ=

1
σh

�
2
2
|z+ | + |z× | .

(5.25)

This is referred to as a two-phase filter, which has twice the degrees of freedom of a single-phasefilter. The bonus of extracting information from both polarizations of the GW comes with the costof increasing the expectation value of ρ when there is no signal present.
By combining the definition of the SNR (equation (5.25)) with the template normalization σh(equation (5.16b)) we can define an eﬀective distance for a given trigger

Deﬀ =

σh
.
ρ

(5.26)We choose a normalization distance of 1 Mpc to set the templates’ amplitudes. This results in theunits of the eﬀective distance being Mpc.

5.2.2

Template Bank Construction

The manifold of waveforms is a continuous space in the component masses, of which we are only
able to search discrete points. In order to make sure we do not miss a signal because its parameters
are slightly diﬀerent from what we are searching for, we construct a bank of templates in such a way
as to minimize the loss of a signal’s SNR. This is done by computing the overlap of waveforms with
diﬀerent parameters using equation (5.15).
As noted above, the overlap of a waveform with itself is unity. Any mismatch, 1 − M, between
two waveforms will show up as a reduction of the recovered SNR when searching for one waveformwith the other. When constructing a template bank, it is useful to view this mismatch as a measureof the distance between two templates. With this in mind, we can create a metric that will tell us
the distance between two templates of diﬀerent parameters. This is done by defining the parameter
space in which the metric exists, and then computing the mismatch between templates. For two

53
templates infinitesimally far apart in the parameter space θµ , the mismatch is [121, 122]

1 − M ≈ 1 − (h (θµ + dθµ ) |h (θν ))≈≈1[1 + 1 − 2 (h (θµ + dθµ ) |h (θν ))]21[(h (θµ ) |h (θν )) + (h (θµ + dθµ ) |h (θν + dθν ))2(5.27)− (h (θµ + dθµ ) |h (θν )) − (h (θµ ) |h (θν + dθν ))]
≈
≈
≈
≈

1
[(h (θµ + dθµ ) − h (θµ ) |h (θν + dθν ) − h (θν ))]
2
��
��
1
∂h µ ∂h νdθ | ν dθ
2
∂θµ
∂θ
1
[(hµ |hν )] dθµ dθν21gµν dθµ dθν ,2(5.28)∂hwhere hµ ≡ ∂θ
µ , and we have used the symmetry (a|b) = (b|a). The above equation defines themetric gµν to be of the quadratic formgµν ≡ (hµ |hν ) .

(5.29)The metric (equation (5.29)) can be used to predict to quadratic order how quickly the mismatch
grows as we move in diﬀerent directions in parameter space away from a particular point. Using thismetric we place templates such that the furthest distance any point is from a template is less than atolerance value �. This will ensure that no more than � of the SNR is lost due to our discretizationof parameter space.
The question we must now ask is which parameter space to use. Answering this question specifies
which directions we will use when taking partial derivatives of the templates. In the end, it is usefulto choose those directions such that the metric is as flat as possible across the space. This helps indeveloping an algorithm that most eﬃciently covers the parameter space. Such a space has beenfound (equations (23) and (24) in Ref. [123]) whose directions correspond to the chirp times of a
waveform. Using 2.0 PN order waveforms for nonspinning binary objects in quasi-circular orbital

54
evolution, there are chirp times τi for i = 0, 2, 3, 4, corresponding to the chirp times proportionalito (v/c)) in the PN expansion of the phase. These chirp times are given in terms of the mass
parameters as [124, 125, 122]

τ0

=

τ2

=

τ3

=

τ4

=

5
8/3

256 (πf0 )5M−5/3 ,
2

192η 2/5 (πf0 )π5/3
8η 3/5 (πf0 )��
743 11+ η M−1 ,336
4

M−2/3 ,

5
4/3

128η 4/5 (πf0 )��
3 058 673 5429
617 2
+
η+
η M−1/3 ,1 016 064 1008
144

(5.30a)(5.30b)(5.30c)(5.30d)where M is the total mass of the binary, η = m1 m2 /M 2 is the symmetric mass ratio, M = η 3/5 M
is the chirp mass, and f0 is some fiducial frequency. If f0 is chosen to be the starting frequency for
the waveform generation, then the chirp times correspond to the 0PN length of the chirp for τ0 , andthe (i/2)PN correction to the chirp time for i = 2, 3, 4. For the search described in later chapters,
we have chosen to use (τ0 , τ3 ) as the parameter space in which we lay our templates, however wecompute a 3D metric that additionally includes the diﬀerence in coalescence times of two waveforms
for use in other portions of the search. Although the metric is not perfectly flat for PN order greaterthan 1.0PN, we assume it is slowly varying, and thus flat in the local vicinity of a point. When
constructing the bank, we use the hexagonal placement algorithm [123], which is the most eﬃcientplacement algorithm for a flat, 2D space.

5.2.3

Signal-Based Vetoes

Signal-based vetoes are discriminators that can be used to separate triggers arising from either arandom instantiation of the Gaussian noise or transient excess power in the data from triggers arisingfrom actual signals. These vetoes work by comparing how a particular aspect of the data should
behave in the presence of a signal and not.
In the following sections we discuss several signal-based vetoes that look at diﬀerent aspects of thedata. The first two signal-based vetoes described below, the χ2 and r2 vetoes, use information from55
a single detector. Since we require triggers to be coincident between multiple detectors, which wewill discuss in more detail in chapter 6, there are addition signal-based vetoes we can perform using
information from more than one detector. The eﬀective distance cut and the amplitude consistency
check are examples of this, using information from two coaligned detectors.

5.2.3.1

The χ2 Veto Calculation

The χ2 veto tests the consistency of the data with what we expect from a signal by looking at how the
SNR, ρ, is accumulated from diﬀerent parts a template. This is done by breaking up the template
into p continuous bins such that each is expected to provide ρ/p to the SNR calculation in thepresence of a real signal. This calculation is described by

χ2 =

p �
�
i=1

ρi −

ρ
p

�2

,

(5.31)

where p is the number of χ2 bins, ρi is the SNR from the ith bin of the template, and ρ is the total
SNR.
With this formulation, in the presence of Gaussian noise, χ2 is a quantity that is χ2 distributedwith 2p − 2 degrees of freedom:� 2�
χ N = 2p − 2 .

(5.32)

Normally such a quantity would be χ2 distributed with p − 1 degrees of freedom, however there is an
extra factor of 2 because ρ2 is calculated from a two-phase filter, which has an expectation value of� 2�
ρ = 2 in Gaussian noise.

The χ2 is not actually the quantity we threshold on when applying the χ2 veto. This is because

the expectation value of χ2 in the presence of signal is not the same as in Gaussian noise. Since weare using a discrete template bank to search for signals in a continuous parameter space, we will notrecover a signal with a template of those exact parameters. If we calculate the expectation valueof χ2 in the presence of signal without noise, evaluated at the time when the SNR is maximized, we

56
find
� 2�
χ S = δρ2 ,(5.33)

where δ is the mismatch between the template and the signal. Since a real signal can have a large
value of χ2 for high SNRs, we end up thresholding on a normalized χ2 , ξ 2 , where we veto triggers
with
ξ2 ≡

χ2
> ξ∗2 .
p + δρ2(5.34)For the χ2 veto, p, δ, and ξ∗ are tunable parameters.
It should be noted that in order to do this calculation, p additional inverse Fourier Transformsof the data need to be computed for the p diﬀerent parts of the template.

5.2.3.2

The r2 Veto Calculation

Another signal-based veto that has been developed is the r2 veto. This veto is also a measure of
how much the data looks like an actual signal. The diﬀerence between the r2 veto and the χ2 veto isthat the r2 veto looks over a stretch of data rather than at a single point in time. Specifically, ther2 veto measures the amount of time the χ2 time series spends above a particular threshold for the2
T seconds prior to an trigger. The duration rduration
is given as

2rduration=dfsampling,(5.35)

where d is the number of data points in the time window (t0 − T, T ) that have χ2 /p > r∗2 , andfsampling is the sampling rate of the data. For the r2 veto, T and r∗2 are tunable parameters.
5.2.3.3

The Eﬀective Distance Cut

When looking for coincident triggers between detectors that are coaligned, we apply an eﬀective
distance cut. Coaligned detectors, such as the two collocated Hanford detectors H1 and H2, havethe interesting property of sharing the same antenna patterns. The eﬀective distance seen by a

57
particular detector for a given CBC GW signal isD
Deﬀ = �,2
2
2
F+ (1 + cos ι) /4 + F×2 cos2 ι(5.36)

which is only dependent on the distance to the source D, the inclination angle ι between the orbital
plane and the line of site between the source and the detector, and the detector antenna patterns
F+ and F× . In terms of measured quantities, the eﬀective distance is defined in equation (5.26).
Since coaligned detectors have the same antenna patterns, they should see the same Deﬀ for triggers
coming from a real CBC GW signal.
Since noise fluctuations in the diﬀerent detectors can separately change the eﬀective distanceseen by each detector, we do not require the eﬀective distances to agree perfectly. Instead we allow
for some variation and veto coincident triggers that have a fractional eﬀective distance diﬀerence κ
greater than a given threshold κ∗ (i.e., κ > κ∗ ). We compute κ as

κ=

2 |Deﬀ,A − Deﬀ,B |,Deﬀ,A + Deﬀ,B(5.37)

where Deﬀ,A is the eﬀective distance from detector A, Deﬀ,B is the eﬀective distance from detector
B.
This veto is only applied between coaligned detectors due to the fact that the same signal canbe seen by two nonaligned detectors with very diﬀerent amplitudes, and thus eﬀective distances,since the blind spots of the detectors’ antenna patterns do not coincide. The blind spots are a verysmall portion of the sky and thus it is not unreasonable to apply a similar veto between nonaligned
detectors. In that case, the fraction of true signals being vetoed can be calculated by looking at thevolume of the universe in which such a signal would originate, bounded by the distance at which
each detector would see a signal with an SNR above a particular threshold.

58
5.2.3.4

The Amplitude Consistency Check

In this section we describe the calculation of an amplitude consistency check. This check is very
similar to the eﬀective distance cut, however instead of checking the consistency between two triggers
that make up a coincident trigger, we are checking whether the lack of a trigger from one coaligned
detector is consistent with the measured parameters of the trigger in the other detector.
As always, we use the PSD from each interferometer to calculate a template normalization σ(equation (5.16b)) for a given time. With this, we compute a horizon distance for a canonical mass
BNS template as the distance to which we would see an optimally located and oriented signal withan SNR of 8. Similar to the eﬀective distance (equation (5.26)), this is given as

Dhorizon =σ.
8(5.38)

We can rearrange the eﬀective distance cut (equation (5.37)) using the eﬀective distance (equation(5.26)) for the trigger from detector A and the horizon distance (equation (5.38)) from detector B.
Solving for ρA , we find the veto to be

ρA >

�

2 + κ∗
2 − κ∗

��

Dhorizon,ADhorizon,B�ρthreshold,B .

(5.39)In order to use this equation to veto triggers from detector A when there is no coincident trigger in
detector B, we have substituted the SNR threshold of detector B ρthreshold,B in the place of ρB . This
veto tells us to veto triggers from detector A when the ratio of the horizon distances in detectors A
and B is small enough such that there should have also been a trigger above the SNR threshold ofdetector B. The factor (2 + κ∗ ) / (2 − κ∗ ) allows for noise fluctuations to change ρA , Dhorizon,A , andDhorizon,B by some amount.
This veto is also only applied between coaligned detectors due to the fact that the there exist
locations in the sky and polarizations of a signal such that for two diﬀerent, nonaligned detectors,
one of them should see the signal above the SNR threshold, while the other should not. As in the59case for the eﬀective distance cut, a similar veto can be constructed for nonaligned detectors if one
is willing to veto a fraction of true signals.

5.3

Previous Detection Statistic: Eﬀective SNR

For the previous searches, a detection statistic was developed to separate signals from background.
This statistic is the combined eﬀective SNR and is constructed as follows.
The single-detector SNR ρ is produced by matched filtering the data against our templates. From
ρ we define the eﬀective SNR, ρeﬀ , as

ρ2eﬀ = ��

ρ2
�
�
�,
χ2
ρ2
1 + 2502p−2(5.40)

where p is the number of bins used by the χ2 veto, 2p − 2 is the number of degrees of freedom ofthe χ2 veto in Gaussian noise, and the 250 is a tunable parameter that helps to further separate
signals from background. This definition of the eﬀective SNR reduces the apparent significance of
nonGaussian instrumental artifacts since it weights the SNR by the χ2 . This eﬀectively reduces thesignificance of outliers to the expected SNR distribution due to Gaussian noise (i.e., nonGaussian
instrument artifacts) while minimally aﬀecting the apparent significance triggers of real signals.
We then combine the eﬀective SNRs for the single-detector triggers that form a coincident trigger
into the combined eﬀective SNR, ρc , for that coincident trigger using

ρ2c =

N
�
i=1

ρ2eﬀ,i .

(5.41)

60

Chapter 6

Data Analysis Pipeline and Tuning
In this chapter we describe the standard pipeline used to analyze data from GW detectors, such as
the LIGO detectors, for GW signals from CBC. In addition, we will also describe the tuning thatwas done for the Search for Low Mass CBCs in the First Year of LIGO’s Fifth Science Run (S5)
Data, referred to as “this search.”
The name of the executable that generates the standard pipeline is the Hierarchical Inspiral
Pipeline Executable (HIPE); and we colloquially refer to this portion of the pipeline as the HIPE
pipeline. The goal of the HIPE pipeline is to distill GW data in the form of a time series of strain
data into a collection of candidate GW triggers. The HIPE pipeline is generic enough that it is used
for many diﬀerent matched filtering searches including searches for GW signals from low mass CBC,
triggered CBC searches (i.e., GRBs), high mass CBC, CBC whose component objects are spinning,and black hole ringdowns.
The HIPE executable, lalapps inspiral hipe, generates a directed acyclic graph (DAG) ofcomputational jobs that are run by the Condor High Throughput Computing environment, wherea DAG, as described by the name, is a collection of jobs with interdependencies such that there
are no loops. HIPE reads an input file containing all of the arguments to invoke at diﬀerent stagesof the pipeline as well as the locations of files containing lists of times associated with individual
detectors’ operating status. The time files for each detector included for this search were times thedetector was in “science mode” (i.e., when the detectors were operating nominally), which we callanalyzed time files, and four sets of time files indicating increasingly severe excursions from nominal

61
operation, which we call category veto files.
Using these inputs, HIPE first computes all of the times an individual detector was in sciencemode for 2064 continuous seconds. Because of the way our filtering is done, HIPE throws out the
beginning and end 72 seconds in order to calculate the “analyzable times”. It then takes intersectionsof these times across diﬀerent detectors to obtain diﬀerent kinds of “coincident times” during which
two or more detectors were in operation. Finally, HIPE sets up job dependencies across the diﬀerent
stages of the pipeline and produces the DAG to be run by Condor.
The diﬀerent stages of the HIPE pipeline employed in this search were:• template bank generation (detailed in section 6.3);• first-stage single-detector trigger generation (detailed in section 6.4);• first-stage coincident trigger generation (detailed in section 6.5);• template bank regeneration (detailed in section 6.6)
• second-stage single-detector trigger generation (detailed in section 6.7);• second-stage coincident trigger generation (detailed in section 6.8).
The template bank generation stage generates a bank of templates in some specified portion of theparameter space. Those templates are then used for the first-stage single-detector trigger generation
by performing matched filtering of the templates against the data, producing a list of single-detector
triggers. The first-stage single-detector triggers from diﬀerent detectors are then combined andused to search for coincident triggers in the first-stage coincident trigger generation by identifying
triggers that are coincident across detectors in parameter space, including time. The first-stage
coincident triggers are converted back into single-detector template banks for a second round ofmatched filtering in the template bank regeneration stage. This new template bank is used for thesecond-stage single-detector trigger generation by matched filtering the data against the templates
again, this time also calculating and vetoing triggers based on signal-based vetoes. Finally, thesecond-stage single-detector triggers are used for second-stage coincident trigger generation, where

62
triggers occurring during times in the category veto files are thrown out, coincident triggers are found,and more signal-based vetoes are calculated and used. The output of this stage of the pipeline is
the final output of the HIPE pipeline. The pipeline flow is represented by figure 6.1.
Multiple stages of matched filtering and coincidence occur because the calculation of the signalbased vetoes in the second-stage single-detector trigger generation is very compute intensive. In
order to speed the search up, we only compute these signal-based vetoes on triggers that havealready passed coincidence.
The follow sections will go through each of these stages in detail, but first let us describe theblinding and tuning procedure we use and how we test our pipelines using software injections.

63

Figure 6.1: The HIPE Pipeline
A flowchart showing how GW detector data are analyzed in the HIPE pipeline.

64

6.1

Blinding the Analysis Pipeline

One of the goals of searching for GWs is to give a statistical probability that the triggers seen are
real signals or due to background. This requires that the search pipeline’s tuning not be informed by
the in-time triggers that are present in the data. This can be achieved by blinding the data during
the tuning process.
To do this, we only look at coincident triggers associated with 1) time-shifted data in which thedetector data streams streams have been slid against each other such that no true GW signal could
produce coincident triggers in that time shift (this data estimates the background of the pipeline),
2) software injected data in which GWs have been artificially introduced into the data in software(we use these to measure the eﬃciency of our pipeline at detecting real signals), and 3) playgroundin-time data that we have defined to be 600 of every 6370 seconds of data, which we use to checkour pipeline to make sure nothing has gone awry in our analysis (we search for GW signals in thisdata but exclude it from statistical statements resulting from the search). Using these triggers, wetune the parameters of the pipeline and construct a detection statistic to rank the importance of thecoincident triggers such that the time-shift coincident triggers are suppressed while the coincident
triggers associated with software injections are emphasized.
Only at this point do we look at the in-time coincident triggers (“open the box”), unblinding theanalysis to the part of the data in which putative GW detections can be found. In this way we avoid
any bias in establishing the detection statistic that can come from tuning the search parameters aftercandidates are examined, which could lead to subjectively elevating or “killing” such candidates.
However, we choose to apply an unavoidably subjective follow-up procedure (described in chapter10) to develop suﬃcient confidence in a potential detection candidate to warrant the announcementof the first discovery of GWs.

65

6.2

Simulated Waveform Injections

In order to measure the eﬃciency of our pipeline to recovering GW signals from CBC, we inject
several diﬀerent PN families of waveforms into the data and check to see the fraction of signals thatare recovered. The diﬀerent waveform families used for injections in this search include GeneratePPN
computed to Newtonian order in amplitude and 2PN order in phase using formulae from reference[3], EOB computed to Newtonian order in amplitude and 3PN order in phase using formulae fromreferences [4, 5, 6, 7], PadéT1 computed to Newtonian order in amplitude and 3.5PN order in phaseusing formulae from references [8, 9], and SpinTaylor computed to Newtonian order in amplitude
and 3.5PN order in phase using formulae from references [10] and based upon references [11, 12,
3, 13, 14, 15, 8, 16]; using code from reference [126]. Each of these families except for SpinTaylor
ignores the eﬀects of spin on the orbital evolution.
The parameters of the signals injected for each of the families are chosen from either an astrophysical distribution or a nonastrophysical distribution. The parameters that are astrophysically
distributed include:• the sky location (right ascension, declination) is distributed uniformly on the surface of asphere, since we expect no preferred direction in the sky for extragalactic GW sources;• the inclination angle (ι) where cos ι is distributed uniformly between 0 and 1, since there is nopreferred orientation of sources with respect to the observer’s line of sight;• the polarization azimuthal angle (ψ) distributed uniformly between 0 and 2π, since there is nopreferred orientation of sources with respect to the observer’s orientation;• (for the SpinTaylor waveform family) the component objects’ spin orientations relative to theinitial orbital angular momentum are individually distributed uniformly on the surface of asphere. No correlation is expected amongst the three orientations of the component objects’
spin orientations and the orientation of the orbital angular momentum for captured binaries.
However, there may be correlations for primordial binaries.
The parameters that are nonastrophysically distributed include:

66• the total mass is distributed uniformly between 2 and 35 M⊙ , since there is no reliable theory
for predicting the distribution in total mass, and we wish to measure the CBC rate for a rangeof masses;• the mass ratio q = m1 /m2 is distributed uniformly between 1 and the maximum value suchthat m1 , m2 > 1M⊙ , since there is no reliable theory for predicting the distribution in mass
ratio;• the physical distance D is distributed uniformly in log10 D, since we want to test our pipeline
on a large range of signal amplitudes;• (for the SpinTaylor waveform family) the component objects’ unitless spin magnitudes â ≡��(cS) / Gm2 are individually distributed uniformly from 0 to 1, since there is no consensus
on the physical distribution of the component objects’ spin magnitudes nor the correlations
between the component objects’ spin magnitudes.

6.3

Template Bank Generation

The search we are performing for GW signals from CBC is a matched filter search where we compare
the data against a model signal waveform to see how much the data looks like the model. In thiscase the model, which we call a “template,” is a gravitational waveform calculated from PN theory
for a CBC of a certain mass (as discussed in chapter 2). In section 5.2.2 we discussed theoreticallythe generation of a template bank with the use of a metric on the parameter space. In this section
we will go over the specific choices that were made in the construction of template banks used inthis search.
Based on the goals of the search, we know the masses of the templates should cover the range of
total mass from 2 to 35 M⊙ with a minimum component mass of 1 M⊙ . A factor that governs thetotal number of template is that we are willing to lose at most 10% of the signals. A mismatch dueto the discreteness of the template bank reduces the SNR we can recover, and thus the distance towhich we can see, by a factor of the mismatch. Assuming an signal distribution uniform in volume,

67
this translates the 10% tolerance in signal recovery to a match of (90%)

1/3= 97% or a 3% mismatch.

The next choice we make is which waveform approximation is used to generate the templates.
In interests of speed and simplicity of code, we choose to use 2.0 restricted PN templates, which
are second PN order in the phase evolution and Newtonian order in the amplitude, generated in thefrequency domain using the SPA. These choices are validated by studies that show these templates
can recover time domain signals of diﬀerent orders and diﬀerent approximant families [127, 124, 128].
The final items we must choose in order to generate a template bank are the frequency range ofour calculations and the calculation of the PSD of the data, Sn . The LIGO noise varies over manyorders of magnitude for the frequency range we are interested in. The lowest mass signals we aresearching for extend from arbitrarily low frequencies up to a few kHz. Because of the steep slope ofthe seismic noise in the detectors starting at 40Hz and below, we choose to suppress the data below
30Hz using two rounds of eighth-order Butterworth high-pass filters, and then analyze data above
40 Hz. This is shown in the spectra of the data after the diﬀerent stages of filtering in figure 6.2.
The gap between where our high-pass filters end and where we start analyzing data is to ensure weare not suppressing the signals in the band we analyze.
To choose our high frequency cutoﬀ, we look at the frequencies our signals end on, a few kHz
and below, and we choose 2048 Hz as our high frequency cutoﬀ. In order to reduce the amount ofdata we work with, we down-sample the data from 16384 samples per second to 4096 samples per
second, applying low pass filters above the Nyquist frequency of 2048 Hz in order to prevent aliasingof power from frequencies above this to low frequencies.
We calculate Sn using a 2048-second-long segment of data. This is broken up into 15 overlapping
chunks of 256 seconds, which we Fourier transform to calculate 15 PSDs. Sn is then constructed bytaking the median of these in each frequency bin. This PSD is used in the metric calculation of thechosen parameter space and, later on, in our matched filtering.
With the above procedures, we end up with ∼5000 templates covering the parameter space weare interested in. An example template bank can be seen in figure 6.3. The essentially flat metricin (τ0 , τ3 ) (discussed in chapter 5) produces an essentially uniform hexagonal spacing of templates

68

Figure 6.2: Data Preprocessing
Data spectra after the diﬀerent levels of filtering we perform prior to matched filtering.
in that space, and a nonuniform spacing in (m1 , m2 ) space. Conversely, the straight bounds of thebank in (m1 , m2 ) space results in current boundaries in (τ0 , τ3 ) space.

6.4

First-Stage Single-Detector Trigger Generation

The first-stage single-detector trigger generation is the first round of matched filtering we perform
on the data. In this stage we perform the same data conditioning as described in the template banksection with an addition step of injecting signals to the data when performing software injections.
In that case, after the first round of high-pass filtering, the injection waveforms are added to thedata, as noted in figure 6.2. The data conditioning continues after this as before with the second
round of high-pass filtering and then resampling.
Once we get to this stage, we perform the matched filtering as described in chapter 5 in orderto obtain an SNR time series for each template. Since the templates we are filtering with have a
specific duration, there will be wrap-around eﬀects near the boundaries of each chunk where the69

Figure 6.3: Sample Template Bank
A sample template bank is shown in both (m1 , m2 ) space and (τ0 , τ3 ) space.
SNR time series will be corrupted. We set the 256-second chunks to overlap by 128 seconds on eachend, enabling us to throw away triggers within 64 seconds of each boundary containing the corrupted
portion of the time series.
The SNR time series is then used to generate single-detector triggers by finding times when the
SNR exceeds a threshold value of 5.5 for this search. This may happen for many adjacent points intime, so to reduce the amount of triggers we generate, we cluster the time series over the length of
the template ttmplt using a sliding window. This is implemented in the following way:
1. start at the beginning of the time series;
2. find the next time the SNR crosses the threshold;3. generate a trigger, which records that time and SNR ρ∗ ;4. search the time series over the next ttmplt seconds for the first time when ρ > ρ∗ ;5. if there exists such a time, update the trigger to contain the new time and SNR, then repeat
step 4;
6. if there does not exist such a time, store the trigger in a trigger list and, if there is additional

70time in the time series, repeat step 2.
After this procedure, we have a list of single-detector triggers from a single template. This process
is repeated for each template in the bank. Once we are done filtering the data through all of thetemplates, we perform an additional round of clustering before we write these triggers to file and end
this stage. In this final round of clustering, we cluster triggers coming from diﬀerent templates when
they satisfy some requirements. As we have noted before, adjacent templates have a match close to
1 in order for a signal falling between templates to have a mismatch with the nearest template ofless than 3%. This results in templates near a particular template that produces triggers having agood chance of also crossing SNR threshold and producing triggers.
An algorithm called “trigscan” was devised that clusters these triggers across the template bankusing the metric calculated during the template bank generation [129]. In generating the templatebank, the portions of the metric that were used are those that ignored the time dimension ofparameter space, as mentioned in chapter 5. In contrast, trigscan employs the full 3D metric to do
the clustering. Trigscan works using the following algorithm:
1. sort the trigger list to be time ordered and start with the first trigger and use it to seed acluster of one trigger;
2. create an error ellipse of size �trigscan around the trigger using the metric at that point in(τ0 , τ3 , t) space;
3. for all triggers with a time t ∈ (ti − T, ti + T ), where ti is the time of the current trigger and
T is some time window, create error ellipses of size �trigscan at those points using their values
for the metric;
4. for any triggers whose error ellipse overlaps the current trigger’s, add them to the current
cluster;
5. repeat steps 2–4 for any trigger added to the cluster;
6. when no additional triggers can be added to the cluster, save the trigger with the largest SNR

71
in the clustered trigger list;
7. if there is a next trigger in time, use it to create the next cluster of one trigger and repeating
steps 2–5
We decided to use a trigscan metric scaling value of �trigscan = 0.06. Ignoring the time dimension
for a moment, let us assume we are in flat space with hexagonally spaced templates where thefurthest distance any point can be from a template is 0.03. In this case, the templates have aseparation of 2 × 0.03 cos(π/6) ≈ 0.052. If we are to create circles around each template and expand
them to radius 0.06, templates slightly more than twice the templates separation distance will havecircles that intersect. This can be seen pictorially in figure 6.4. Since the (τ0 , τ3 ) space in which wehave chosen to place are templates is almost flat, our chosen value of �trigscan = 0.06 allows triggers
occurring at the same time (i.e., ignoring the time dimension of the 3D space) to be clustered when
they are from next-to-nearest-neighbor templates in the template bank.
The trigscan-clustered trigger list is then saved to an xml file where each trigger contains information such as the GPS time it occurred, the SNR (ρ), the eﬀective distance (Deﬀ ), the phase at coalescence (φ0 ), the mass parameters (m1 , m2 , M, η, M), the chirp time parameters (τ0 , τ2 , τ3 , τ4 , τ5 ),the length of the template (ttmplt ), and the six metric components of the 3D metric at that point.
Additional information associated with the time analyzed that is also calculated and saved to file isthe horizon distance (equation (5.38)) for the detector using Sn for that 2048-second segment.

6.5

First-Stage Coincident Trigger Generation

A true GW signal from a CBC should aﬀect all detectors around the world, and indeed, to verify a
GW detection, it will have to be found by multiple means (i.e., associated electromagnetic, neutrino,
or GW observation). For this search, we are only considering data from the LIGO detectors, so werely on finding coincident GW triggers from diﬀerent detectors.
Additional information useful in evaluating the performance of the pipeline is the backgroundrate of triggers. The background is assumed to be due to accidental coincidences, which can be

72

Figure 6.4: Trigscan Example
Templates are plotted as blue dots, which are spaced such that the furthest distance any point canbe from a template is 0.03. We draw a black circle of radius 0.06 around a chosen template. Bluecircles are drawn around three example templates that can be trigscan clustered with the chosen
template (i.e., the blue circles intersect the black circle). A red circle is drawn around a single
template that is slightly too far to be clustered with the chosen template.
estimated by time shifts (described in section 6.1). This estimation fails for collocated detectors,
such as H1 and H2, due to time correlated background arising from the shared environment andcrosstalk between the detectors (e.g., stray optical signals or electric control signals). Because ofthis, we can not reliably estimate the background in this way, thus H1H2 coincident triggers are
excluded from this search. As these eﬀects should only occur between collocated detectors, there isno in-time correlation between the Hanford and Livingston detectors.
Below, section 6.5.1 describes how the coincidence algorithm is executed, while section 6.5.2
describes how the time shifts are performed.

6.5.1

Coincidence Algorithm

The coincidence algorithm that has been developed for finding coincident triggers between multiple
detectors is called The Hierarchical Inspiral Coicidence Algorithm (THINCA) and we commonly73
refer to this stage of the pipeline as the first THINCA stage. Here we take as input triggers from up
to 3600 seconds of time when a unique set of detectors was generating triggers. Potential coincident
triggers could lie across the boundary of two segments and thus not be recorded, however we find anegligible loss in eﬃciency due to this eﬀect. The coincidence test used for this search in THINCA
is similar to the trigscan clustering algorithm described above, called the “ellipsoidal” THINCA or“e-thinca” coincidence test. The main diﬀerence between the two is that trigscan builds up clustersof triggers from a single detector, while the e-thinca test checks whether triggers from two detectors
are coincident without agglomerating coincident triggers into clusters.
First we will describe the e-thinca coincidence test before moving on to the full coincidence
algorithm. The e-thinca test for two triggers from diﬀerent detectors is as follows
1. create an error ellipse of size �e−thinca around the triggers using their respective metrics;
2. if their error ellipses overlap, return yes, otherwise return no.
With the coincidence test defined, we now describe the full coincidence algorithm. It is as follows:
1. create a time sorted trigger list from all the triggers of the multiple detectors, and start withthe first trigger;
2. for all triggers with a time t ∈ (ti , ti + T ), where ti is the current trigger and T is some time
window. If the trigger is from a diﬀerent detector, perform the coincidence test between the
two triggers;
3. if coincidence test return yes, create a coincident trigger in a coincident triggers list and link
both single-detectors triggers to this coincident trigger;
4. repeat 2-3 for the next trigger in the list until end of list;
5. loop over coincident triggers, if all the single-detector triggers in the current coincidence arecoincident with an additional trigger from a diﬀerent detector, add that trigger to this coincidence;

74
6. loop over coincident triggers and remove coincident triggers that are subsets of higher-ordercoincident triggers.
After this procedure, the list of coincident triggers is saved to an xml file where the single-detector
triggers of a coincidence all have the same event id as well as all the information saved during thesingle-detector trigger generation.

6.5.2

Background Estimation: Time Shifts

As mentioned previously, time shifts are performed between diﬀerent detectors in order to estimate
the amount of background triggers we expect from the pipeline. These time shifts are chosen suchthat any coincidences coming from the time shifted data could not have come from a real GW
signal. To ensure this, we have chosen to create a vector that we multiply by a shift number in order
determine the amount of time to shift each detector. The vector contains a unit time shift for eachdetector, which must be unique.
For this search we have chosen the vectorv = {0, 10, 5} ,(6.1)

which we associate with the detectors H1, H2, and L1 respectively. This is used for 100 diﬀerent
time shifts ranging in shift number from -50 to 50, skipping 0. The maximum number of time shifts
we could do with these shift vectors is 360, however we have found 100 to be adequate to estimate
the false alarm rate for the loudest in-time triggers in the past. Additional time shifts could be
performed with smaller shift vectors, however there is a time correlation for triggers coming fromreal signals that must be avoided.
As we said before, THINCA takes in as input up to 3600 seconds worth of triggers from a unique
combination of multiple detectors. Using the above information, we find that when we shift triggers
from diﬀerent detectors against each other, we can end up with up to 2000 seconds where there are
no H1 triggers to search for coincidence with H2 triggers. In order to avoid this but keep with the

753600-second limit, we have chosen to wrap the triggers on a 3600-second ring defined by the original
start and end times of the coincident segment. In this approach, when a trigger’s shift time becomes
greater than the end time of the ring, we subtract the duration of the ring from the trigger’s time,
which moves it to the beginning of the ring. Conversely, when a trigger’s shifted time becomes lessthan the start time of the ring, we add the duration of the ring to the trigger’s time, which moves
it to the end of the ring.
Once we have shifted the triggers in time, we can apply the above coincidence algorithm toproduce coincident triggers for the individual time shifts. At the end of the algorithm, before wesave the coincident triggers to file, we unshift the triggers’ times so that they return to their originaltimes.
There a several nice features of this approach. The first is that it keeps the amount of analyzed
time in each time shift the same as the in time. In addition, this approach limits the maximum
time oﬀset between time-shifted coincident triggers to 3600 seconds. Since we expect the detectors’ behaviors to be roughly stationary over this time period, this more accurately estimates the
background for the in-time coincident triggers.
As discussed in section 6.1, there are noise correlations between H1 and H2 that bias the background estimation for H1H2 coincidences, which we find to also have an eﬀect H1H2L1 coincidences.
We have tested using a background estimation for H1H2L1 coincidences in which we shift H1 andH2 triggers by the same amount and found this to correct this bias. As this was done late in theanalysis, it was not folded into the final results, however it has become a standard procedure for theother searches we are performing in the S5 data.

6.5.3

Coincidence Tuning

The threshold in the algorithm that controls the maximum e-thinca value � we accept for coincidence
we refer to as the e-thinca threshold �∗ . Larger values of �∗ correspond to accepting larger error
ellipses for determining coincidence. We tuned the e-thinca threshold by looking at time-shifted
coincident triggers and triggers associated with software injections. Figure 6.5 shows the fraction

Fraction of Injections Found

76

0.615

5 × 10−4
0

0.614
0.613
0.612
0.611
0.0

0.2

0.4

0.6

0.8

E − Thinca Threshold

1.0

Figure 6.5: E-Thinca Tuning
The percentage of injections recovered (blue and left y-axis) and the expected percentage of
injections recovered due to background coincidences (red and right y-axis) as a function of thee-thinca threshold. We have oﬀset the blue and red curves for ease of comparison. The horizontal
red lines show the eﬃciency gain per additional found injection. We see that over the whole rangeof e-thinca thresholds, there is less than one expected injection due to the background.
of injections recovered as a function of the e-thinca threshold. There we plot both the fraction of
injections recovered and the expected fraction of injections whose time windows overlap time-shifted
coincidences Ffalse , which is given as

Ffalse =Ntime−shifts Twindow,
Ninj Ttime−shifts(6.2)

where Ninj is the number of injections made, Tinj = Ninj twindow is the total amount of time wewere searching for injections, twindow is the window used to determine whether or not an injection
is recovered, Ntime−shifts is the number of time-shift coincidences for a given e-thinca threshold, andTtime−shifts is the total amount of time analyzed for time-shifted coincidences.
Figure 6.5 shows that for injections into H1H2L1 time, we find that the last additional injection
is found at an �∗ = 0.5 and that there is less than one expected injection due to the background atthat value. Because of this, we decided to use �∗ = 0.5 as the threshold for this search.

77

6.6

Trigger to Template

In this stage of the pipeline, we create a new template bank from triggers found in coincidence
between detectors. This is done in the following manner:
1. use the coincident segments that a particular 2048-second single-detector segment participatedin;
2. keep triggers from only that detector;
3. search through triggers keeping only one trigger for each combination of mass parameters;
Once this is done, we save the remaining triggers to an xml file for use as the template bank inthe second-stage single-detector trigger generation.
One thing to note about this stage is that triggers from diﬀerent single-detector segments canwind up as templates for nearby single-detector segments. This is due to the fact that coincidence
segments can take in more than one single-detector segment from a given detector. But are wenot trying to reduce the number of templates for the second, computationally expensive stage ofmatched filtering whereas this eﬀect is causing an increase in the template bank? Even though thisis the case, the increase in the template bank due to this eﬀect is much smaller than the reduction
we are trying to accomplish.

6.7

Second-Stage Single-Detector Trigger Generation

Here we describe the second-stage single-detector trigger generation. As noted above, for this stage
we use the triggers found in coincidence with other detectors as the template bank for the second
round of matched filtering. This requirement significantly reduces the number of templates withwhich we perform this stage of matched filtering. This reduction is necessary because in addition
to performing the standard matched filtering we have described above, we also calculate two signalbased vetoes that test the consistency of the data to what we expect from a signal, one of whichis computationally intensive. This veto requires an additional p inverse Fourier transforms, which

78
causes this stage to be the time-limiting stage of the pipeline.
However, there are also drawbacks to this two-stage approach. Adding in the second stage of thepipeline adds significant complexity to the pipeline. In addition, coincident triggers found after the
first stage of coincidence may not be the same as after the second stage of coincidence. This adds
ambiguity when trying to track specific triggers through the pipeline.
Before the veto calculations begin, and after they are finished, the same algorithms are used asin the first-stage single-detector trigger generation described in section 6.4. Below, in section 6.7.1,
we describe how we tuned the χ2 and r2 signal-based vetoes (described in sections 5.2.3.1 and 5.2.3.2
respectively).

6.7.1

Signal-Based Veto Tuning

As in the e-thinca threshold, the veto parameters for both the χ2 veto and r2 veto are tuned byobserving how they respond to software injection triggers and how they respond to time-shifted
coincident triggers. The goal of the vetoes is to find an area of parameter space where there is a separation of time-shifted coincident triggers from software injection triggers that can be removed, thus
reducing the chance of a background trigger occurring while not aﬀecting the chance of recovering
a true signal.
The χ2 statistic, as described in section 5.2.3.1, has an expected behavior for software injection
triggers that grows proportional to ρ2 , where ρ is the SNR. This can be seen in figure 6.6. We aretrying to tune this search such that triggers associated with signals whose waveforms diﬀer slightly
from the templates are still recovered. The most extreme divergence from the template waveforms wehope to recover come from CBC waveforms whose component objects have significant spin angularmomentum, as discussed in section 6.2. By injecting these waveforms in as software injections, wefind that by loosening the χ2 veto parameters as compared to other injection waveforms, we are still
able recover these spinning signals. This is shown in figure 6.7. The parameters for the χ2 veto thatwere chosen for this search were a mismatch parameter of δ = 0.2 and a threshold of ξ 2 = 10.0.

79

Figure 6.6: χ2 Tuningχ values versus SNR for triggers associated with nonspinning software injections (red) andtime-shifted coincident triggers (black). We can clearly see the ρ2 dependence of the χ2 value forboth categories of triggers.
2

80

Figure 6.7: χ2 Thresholdχ versus SNR for software injections that include spin (black). The old χ2 threshold is plotted inred (tuned with nonspinning injections) along with the chosen χ2 threshold in blue (tuned withspinning injections).
2

81

Figure 6.8: r2 Tuningr veto duration versus SNR for software injection triggers (red) and time-shifted coincident
triggers (blue). The black line denotes the r2 veto that was chosen for this search.
2

The r2 veto (described in section 5.2.3.2) also is expected to have a dependence on ρ because it
derived from the χ2 statistic. Since there is a threshold below which the χ2 can fluctuate without22
contributing to the rduration
calculation, we expect the rduration
to be zero below a certain value of2
ρ = ρ∗ . Above that value of ρ, since we have not modeled the rduration
dependence on ρ, we choose2
to veto triggers with an rduration> A × ρB , where A and B are tuned in this search. For this search

we have chosen T = 6 sec. and r∗2 = 15.0. With these values fixed, we have tuned ρ∗ , A, and B
using figures such as figure 6.8. Using this information, we have chosen ρ∗ = 12, below which we
22veto triggers with rduration> 0, and above which we veto triggers with rduration> A × ρB . We havealso chosen A = 7.5 × 10−3 and B = 1.12764. All of these choices are loose enough such that theleave some space for signals to move due to noise fluctuations or waveform uncertainties without
being vetoed.
In figure 6.8 there are several interesting features that require explanation. The main featurethat jumps out is the bimodal distribution of software injection triggers. We have found a clean

82

Figure 6.9: Diﬀerent r2 Veto Populations
r veto duration versus SNR for software injections with total mass less than 20 M⊙ (left) andgreater than 20 M⊙ (right).
2

2
break between to two distributions where the distribution with lower rduration
is from software2
injections with M < 20M⊙ while the distribution with higher rduration
is from software injections

with M > 20M⊙ , where M is the total mass of the system, as can be seen in figure 6.9. The
total mass is used to define the frequency at which to stop the waveform calculation since thePN approximation breaks down at the frequency of the ISCO. From equation (2.38), the ending
frequency fISCO corresponding to M = 20M⊙ is ∼ 220Hz, which is roughly when the fISCO is in
the “sweet spot” of the detection band (i.e., where the noise is lowest).
Another feature in figure 6.8 is a handful of injections above the line denoting the veto. Closer
inspection of these triggers showed there was a glitch several seconds before the injections. An
example is seen in figure 6.10 showing the r2 time series for the same data stretch with and without
2
the injection present. This figure also shows how there is nonzero rduration
expected for injections

above a certain ρ.

83

Figure 6.10: r2 Sample Time Series
The r time series for the same example stretch of data with (left) and without (right) a software
injection present. The green line denotes the threshold value above which we count time samples
in calculating the r2 veto duration. Comparison of the two plots reveal features that are due to theinjection and present without the injection.
2

84

6.8

Second-Stage Coincident Trigger Generation

This is the final stage of the HIPE pipeline. We perform the second-stage coincident trigger generation using the same algorithm as in the first stage (described in section 6.5) with the addition ofseveral vetoes. These diﬀerent vetoes include vetoing times the instruments were not behaving nominally, what we call Data Quality (DQ) vetoes, as well as some signal-based vetoes using informationfrom coincident triggers.
In the following sections, we discuss the use of DQ vetoes, and the tuning of two signal-based
vetoes based on the collocation of the Hanford detectors, the H1H2 eﬀective distance cut and theH1H2 consistency check, described in section 5.2.3.3.

6.8.1

Application of Data Quality Vetoes

While operating the LIGO detectors, we monitor many auxiliary channels that measure the state ofthe detectors or their local environment, in addition to the GW channel. We use these channels todetermine whether or not we want to analyze data from the detectors at those times. This is basedon how severe the excursions are from nominal observation mode, what we call “science mode,” dueto known artifacts introduced into the GW channel data by instrumental or environmental eﬀects
coupling to the GW channel. DQ flags are then defined as time segments during which these eﬀects
are present in the instrumental or environmental channels. Based on several considerations, wedefine several diﬀerent categories of DQ flags for use as vetoes in this analysis.
Category 1 DQ flags are times when we know of a severe problem with the detector, bringing in toquestion whether the detector was actually in science mode. Some example cases are when calibration
for the data is not present, or when loud vibrations were caused in the detector environment in orderto test the response of the seismic isolation systems. We choose not to analyze these times because
the noise from those times is too high for the data to be useful and we do not want it to aﬀect our
median PSD for adjacent times.
We then analyze the data, producing single-detector triggers for use in the rest of the categorization of the DQ flags. In the categorization we take into account several considerations. The first is

85
the “eﬃciency ratio.” For this, we calculate the fraction of triggers from an individual detector that
occur during a particular proposed DQ flag and divide this by the “dead time,” which is the fractionof the time flagged by the DQ flag. When the eﬃciency ratio is much greater than one, this DQ
flag is highly correlated with false trigger production and is thus a good candidate for eliminatingexcess triggers. Another criteria we use in the categorization procedure is the “use percentage” of a
DQ flag. The use percentage is defined as the fraction time segments from a particular DQ flag that
contain single-detector triggers. This value is a measure of how strongly a particular channel couples
to the GW channel. Ideally, the use percentage for a good DQ flag veto would be one, meaning thatbackground triggers always occurred whenever that particular DQ flag was on. Finally, we measure
the “safety” of using these DQ flags as vetoes by comparing them with the times of hardware injections, which measure the response of the detector to a simulated GW signal. This helps to confirm
that the DQ flags will not be correlated with real signals. We require all of the DQ flags used as
vetoes to pass this safety check.
Category 2 DQ flags are defined as times when there is a reasonably well established coupling
between the GW channel and an auxiliary channel, the flag has a high eﬃciency ratio, particularly
with high SNR triggers, and when there is a use percentage of 50% or greater. An example is when
any of the data channels in the length-sensing and control servo reach their digital limit. These aretimes when we do not trust the data from the detector, however the data were not bad enough tomess up the PSD calculation. We allow these times to be included in our segment generation, butremove triggers that come from these times while calculating bounds on the rate of CBCs.
Category 3 DQ flags are defined as times when the coupling between the auxiliary channel and
the GW channel is less well established or when the use percentage is low, but we still find higheﬃciency ratio for the particular DQ flag. An example is when the winds near the detectors are over
30 Mph. We allow these times to be included in our segment generation, but, as with the Category
2 vetoes, remove triggers that come from these times while calculating bounds on the rate of CBCs.
Category 4 DQ flags are defined as times when the coupling between the auxiliary channel and
the GW channel is not well established, when the use percentage is low, when the overall dead time

86
is several percent or greater, or the eﬃciency ratio is not much greater than one. An example iswhen nearby aircraft pass overhead. We allow these times to be included in our segment generation
and also retain triggers that come from these times. This category is used in following up the loudest
events from the search.
Summarizing the uses stated above, we do not analyze data vetoed by category 1 DQ flags, we
veto but still follow up triggers in times defined by category 2 and 3 DQ flags, and we use thecategory 4 DQ flags as information when following up the loudest events from the search. These
veto categories significantly reduce the SNR of outlying triggers, as can be seen in figure 6.11.
One could argue that the procedure used to define these vetoes biases the search since it uses
in-time triggers. However, the triggers that are used are single-detector triggers before coincidence.
There are many more triggers at this stage than after coincidence so that triggers associated withreal GWs make up a small fraction and do not aﬀect the tuning procedure to a significant degree.
This procedure is in the process of being automated (as much as possible) to bias and to minimize
the chance that GWs will be vetoed.

87

Figure 6.11: Eﬀects of DQ Vetoes
Histogram of triggers for the H2 detector at the single-detector level before any coincidence requirement, clustered by the trigger with maximum SNR within 10 seconds, plotted after veto categories
1, 2, and 3. That is, category 2 vetoes remove the triggers in red, category 3 vetoes remove thetriggers in yellow, and triggers in the green histogram remain.

88

6.8.2

Multidetector Signal-Based Veto Tuning

As described in section 5.2.3.3, there are additional signal-based vetoes we can perform for thecollocated Hanford detectors. These vetoes, the eﬀective distance cut and the amplitude consistency
check, determine whether coincident triggers from two coaligned detectors satisfy certain amplitudeconditions given the detectors aligned antenna patterns.
The veto parameter for both the eﬀective distance cut and the amplitude consistency check is thesame κ∗ (cf., equation (5.37)), which is tuned by observing the fractional eﬀective distance diﬀerenceκ distribution for software injection triggers and time-shifted coincident triggers.
As described in section 5.2.3.3, κ should be near zero for triggers associated with signals (e.g.,
from software injections), while this is not necessarily the case for coincident noise triggers fromtime-shifted data. This can be seen in figure 6.12. From this figure, we see that the distribution of
the time-shifted coincident H1H2 triggers peaks around κ ∼ 0.7. We choose our threshold value tobe κ∗ = 0.6. At this value we veto the majority of the time-shifted coincident triggers, while vetoinga very small fraction of the software injections.
In figure 6.12 there is an interesting feature in κ plotted versus H1 SNR for the time-shifted
coincident triggers. Here we see strip of time-shifted coincident triggers that extends up to high
SNRs. We have found that this strip is consistent with what you would expect if you had found an
H1 trigger coincident with an H2 trigger with an SNR of 5.5, while H2 had a horizon distance of 16
Mpc and H1 had a horizon distance of 32 Mpc. This situation is shown with the thin black line.

89

Figure 6.12: H1H2 Eﬀective Distance Cut
The fractional eﬀective distance diﬀerence κ versus the H1 SNR (top left), the histogram of κ (topright), and the H2 eﬀective distance versus the H1 eﬀective distance (bottom) for time-shifted
coincident triggers (blue) and software injection triggers (red) in H1H2L1 time. The thick blacklines show the veto value of κ∗ = 0.6 and the thin black line shows a fit discussed in the text.

90

Chapter 7

False Alarm Rate Detection
Statistic
In this chapter, we describe the development of the false alarm rate (FAR) as the detection statistic
for the Search for Low Mass CBCs in the First Year of LIGO’s Fifth Science Run (S5) Data, referred
to as “this search.” This new detection statistic allows us to search over a large region of parameter
space without being limited by a high background FAR from a smaller subregion. In reference[21], coincident triggers were ranked by combined eﬀective SNR. However, for this search, we use
a statistic derived from the background FAR, thus enabling a much improved eﬃciency for a given
FAR.
The FAR is a statistic that can be defined for triggers when there is a measure of the backgroundtrigger rate. The time-shifted triggers provide an estimate of the FAR for each in-time coincident
trigger. By counting the number of time-shifted triggers with a combined eﬀective SNR greater thanor equal to the in-time coincident triggers’ combined eﬀective SNR, and dividing by the total amountof time we searched for time-shifted triggers, we calculate the FAR for each in-time coincident
trigger. This procedure is done separately for diﬀerent trigger categories, which have diﬀerent
background rates. The result is a varying mapping of combined eﬀective SNR to FAR for thediﬀerent trigger categories. The FARs from the diﬀerent categories can then be recombined into a
combined FAR (FARc) for the final combined set of triggers.
Section 7.1 discusses how diﬀerent trigger categories were chosen for this search, section 7.2
explains the FAR calculation, section 7.3 details how trigger categories are recombined, section

91
7.4 shows the comparison of diﬀerent detection statistics, section 7.5 presents comparisons of intime triggers to the many diﬀerent background trials, and section 9.3.2 describes implications formeasuring the background probability for the loudest triggers.

7.1

Definition of Trigger Categories

In this search, there are several aspects that influence the background rate. These include the typeof coincidence a trigger was found in (triple coincident triggers versus double coincident triggers),
the vetoes that were applied to the triggers, and the mass of the templates that correspond to thetriggers. All of these aspects combine in order to define the trigger categories.

7.1.1

Coincidence Types

The first aspect we will focus on is the coincidence type. In this search we are searching for coincident triggers in time and two mass parameters between three GW detectors, the Hanford 4 km
interferometer (H1), the Hanford 2 km interferometer (H2), and the Livingston 4 km interferometer(L1). A trigger is deemed a double coincidence if it passes the coincidence requirements between
two detectors. A trigger is promoted to a triple coincidence if there are three double coincident
triggers that are all found in coincidence with each other. Since triple coincident triggers have an
extra coincidence requirement, there is a much lower background rate for triple coincident triggers
than for double coincident triggers.
Say each detector has a probability of producing a background trigger for a particular mass
template at a particular time of �1 . Adding the double coincident requirement then reduces thatprobability by an additional factor of �1 , so that �2 = �1 �1 = �21 . Adding an additional coincident
requirement added another factor of �1 , so the probability of getting a triple coincidence trigger ata particular time is then �3 = �1 �2 = �31 .
This background rate reduction can be illustrated by looking at background H1H2L1 triggers
and H1L1 triggers produced during triple coincident time (figure 7.1). For this search the diﬀerent
trigger types we decided to use were H1H2L1, H1L1, and H2L1 triggers. We exclude H1H2 triggers

92

Figure 7.1: Eﬀect of Coincidence on Backround Trigger Rate
The cumulative rate versus combined eﬀective SNR for background H1L1 double coincident
triggers during triple coincident time (blue) and H1H2L1 triple coincident triggers during triple
coincident time (red). The additional coincidence requirement for H1H2L1 triggers is seen to
drastically reduce the background trigger rate.
from this analysis as their background cannot be eﬀectively estimated (section 6.5).

7.1.2

Veto Diﬀerences

Another aspect that aﬀects the background rate is the diﬀerences in vetoes that are applied to
diﬀerent triggers. One set of vetoes that are trigger dependent are the amplitude vetoes applied
between H1 and H2 triggers. These vetoes compare the eﬀective distances found between coincident
H1H2 triggers in H1H2L1 coincidences or, for H1L1 or H2L1 triggers, determine whether the other
Hanford detector should have also produced a trigger at that time. These vetoes can only be applied
when both Hanford detectors were in operation resulting in background rate diﬀerences between
H1L1 triggers from triple coincident time versus H1L1 triggers from double coincident time, andsimilarly for H2L1 triggers. This eﬀect is shown for H2L1 trigger in figure 7.2, which clearly shows
broader background distributions for H2L1 time where the veto in not applied even though there is

93

Figure 7.2: Eﬀect of Amplitude Veto on Backround Trigger Rate
The cumulative rate versus combined eﬀective SNR for background H2L1 double coincident
triggers during double coincident time (blue) and H2L1 double coincident triggers during triple
coincident time (red). The Hanford amplitude veto drastically reduces the brackground rateduring triple coincident times, even though there is an order of magnitude more triple coincident
time than double coincident time.
an order of magnitude more triple coincident time. This eﬀect is also stronger for H2L1 triggers in
triple coincident time than for H1L1 triggers in triple coincident time (figure 7.3). This is due tothe fact that the H2 detector is roughly half as sensitive as the H1 detector, which means that if
there is a trigger in H2, most of the time H1 should also have produced a trigger at that time.

94

Figure 7.3: Eﬀects of Diﬀerent Amplitude Vetoes for H1 and H2
The cumulative rate versus combined eﬀective SNR for background H1L1 double coincident
triggers during triple coincident time (blue) and H2L1 double coincident triggers during triple
coincident time (red). The Hanford amplitude veto is seen to be more stringent for H2L1 triggers
than for H1L1 triggers due to the sentivitiy asymmetry between H1 and H2.

95

7.1.3

Mass Diﬀerences

This search is over a much larger parameter space than a single search has covered before. Specifically, the mass space we cover with templates goes from a minimum total mass M = m1 + m2 of 2
M⊙ to a maximum total mass of 35 M⊙ with a minimum component mass of 1 M⊙ . In order to cover
this large space, we use around 6000 templates whose time durations vary from over 44 seconds tounder 0.34 seconds.
Due to the large variation in the length of the templates, we find that there is a larger variationin the templates’ responses to noise glitches in the detectors where high mass templates tend to pickup the noise glitches as triggers with higher eﬀective SNR. This can be seen in figure 7.4 where
we plot the cumulative distributions of background H1L1 double coincident triggers during triple
coincident time versus combined eﬀective SNR separated by small bands in the mean chirp mass
M = M η 3/5 , where η is the symmetric mass ratio η = m1 m2 /M 2 . Three populations are foundwith diﬀerent slopes for the distributions signifying diﬀerent background rates, thus, we divide themass space into three regions, for simplicity. The chirp mass divisions for these regions occur at
M = 3.48 M⊙ and M = 7.40 M⊙ .
After dividing the triggers of a given type into these three diﬀerent populations, we find that,indeed, the higher chirp mass templates have a broader combined eﬀective SNR distribution forbackground triggers, signifying diﬀerent background rates for the diﬀerent chirp mass regions (figure
7.5).

96

Figure 7.4: Combined ρ2eﬀ Distributions for Background Triggers for Diﬀerent Chirp Mass Bands
The cumulative rate versus combined eﬀective SNR for background H1L1 double coincident
triggers during triple coincident for eleven separate bands of “mean chirp mass” of equal widths.
Lines colored toward the blue end of the spectrum are lower mass bands while line colored toward
the purple end of the spectrum are high mass bands. Three populations are found with diﬀerent
slopes for the distributions signifying diﬀerent background rates. Similar results are obtained bylooking at plots for diﬀerent coincident trigger types.

97

Figure 7.5: Combined ρ2eﬀ Distributions for Background Triggers for Diﬀerent Chirp Mass Bins
The cumulative rate versus combined eﬀective SNR for background H1L1 double coincident
triggers during triple coincident time for the low chirp mass bin (blue), the medium chirp mass bin(red), and the high chirp mass bin (green). Higher chirp mass templates have a shorter time
duration, which leads to a larger response to noise glitches, and thus a broader combined eﬀective
SNR distribution for background triggers.

98

7.1.4

Trigger Categories

Due to all of the above eﬀects, we end up dividing the triggers from this search into 15 diﬀerent
categories during our three diﬀerent observation times. During triple coincident times we have 3
coincidence types subdivided into 3 chirp mass bins, resulting in 9 trigger categories. During each of
the double coincident times we have 1 coincidence type subdivided into 3 chirp mass bins, resulting
in 3 trigger categories for each double coincident time.

7.2

FAR: False Alarm Rate

The FAR can be calculated using any intermediate statistic ranking the “loudness” of triggers (suchas the combined eﬀective SNR for this search), where a higher intermediate statistic corresponds toa larger excursion for background noise behaviour. The FAR for any trigger is calculated usingN
FAR = � ,
Ti

(7.1)

i

where N is the total number of background triggers with intermediate statistic greater than or equal
to the one in question, and Ti is the amount of analyzed time in background trial i.
For individual trigger categories, the maximum FAR is found by

Nj
FARmax,j = � ,
Ti

(7.2)

i

where Nj is the total number of background triggers in the jth trigger category. This can beconverted into a minimum inverse false alarm rate (IFAR) trivially by taking the inverse.
The expected number of triggers below a particular FAR due to background is FAR × T0 where
T0 is the foreground time analyzed. Another way to say this is the expected number of triggers
above a particular IFAR due to background is T0 /IFAR where T0 is the foreground time analyzed,
as seen in figure 7.6.

99

Figure 7.6: Background IFAR Plots
The cumulative number of triggers versus IFAR divided by the analysis time for the nine
categories of background triggers in triple coincident time. The top left, top right, and bottom
plots show the three diﬀerent coincident trigger types for the low, medium, and high mass bins
respectively. Because they are background triggers, they should, and do, fall exactly on the
background line. The diﬀerent category triggers in each plot all start at the lower right corner ofthe plot, but since they overlap, only the top color is shown.

100

7.3

FARc: Combined False Alarm Rate

When combining categories, we think of the diﬀerent categories (coincidence type, trigger mass bin)as diﬀerent trials. If there are m trials, then we expect that of those trials there will be m triggers
with an IFAR ≥ T0 . A simple way to normalize the FAR to bring the expected number of triggers
with an IFAR of T0 back to one is
FAR� = m × FAR .

(7.3)

This can be seen in figure 7.7.
The simple way of combining shows that, moving from right to the left, there is a kink in the
background triggers whenever you reach a minimum IFAR for a trigger category. The reason thesekinks occur is because once we reach a minimum IFAR for a trigger category, there is then one
less “trial” for us to combine. A better thing to do at that point is to normalize the FAR by theremaining number of categories rather than the total. In analytic form, the FARc becomes


p
p
�
�
FARc = Θ (FARmax,j − FAR) × FAR +[Θ (FAR − FARmax,j ) FARmax,j ] ,
j=1

(7.4)

j=1

where p is the number of categories of triggers and Θ (x) is the Heaviside function. Calculating the
FARc in this way can be seen in figure 7.8.

101

Figure 7.7: Combined Background IFAR Plot
The cumulative number of background triggers versus IFAR divided by the analysis time for nine
categories of background triggers combined without correcting for the minimum IFAR of eachcategory. As in the previous plot, because they are background triggers, they should, and do, fallexactly on the background line. The vertical, red, dotted lines denote the minimum IFAR for eightdiﬀerent trigger categories, which are the locations of the kinks in the background line. The ninth
is not shown as that is where the background line terminates.

102

Figure 7.8: Background IFARc Plot
The cumulative number versus the combined IFAR (IFARc) divided by the analysis time for nine
categories of background triggers. Again, because they are background triggers, they should, anddo, fall exactly on the background line.

103

7.4

Comparison of Detection Statistics

Now that we have a method for calculating the FARc detection statistic, let us compare this statistic
against the combined SNR and combined eﬀective SNR in recovering software injections signals. Thisis done using a common comparison known as the Receiver Operating Characteristic (ROC) curve,
shown in figure 7.9.
The FARc statistic out performs both the combined SNR and combined eﬀective SNR for all
mass regions. It is interesting to note that the largest improvement is seen for the low mass bin, andthe least improvement in the high mass bin. This lends support to our assumption that the highmass templates have a higher false alarm rate.

104

Figure 7.9: Comparison of Detection Statistics
We compare the eﬃciency of recovering injected signals versus false alarm probability (FAP) for
three diﬀerent detection statistics in the low (top left), medium (top right), and high (bottom)
mass regions. The vertical line in each plot shows where the loudest in-time trigger is expected tooccur when it is consistent with background. The injected signals were spaced uniformily inlog(D). The plot has been made by reweighting each injection such that they represent a
population uniform in the volume (D3 ).

105

7.5

Plotting Background Trials: Lightning Bolt Plots

In addition to plotting the IFARc of foreground triggers, we can also plot the IFARc of each of theindividual background trials (i.e., time shifts). In this case, the IFARc has been normalized by

IFARci = IFARc ×

T0
,
Ti

(7.5)

where i corresponds the diﬀerent background runs. An example plot of this is shown in figure 7.10.
The plots give 100 examples of the expected distribution of background-only triggers to better judge
the significance (inconsistency with background) of the in-time coincident triggers.

106

Figure 7.10: In-Time IFARc Plot
The cumulative of coincident triggers with IFARc greater than the value of IFARc on the x-axis
for nine categories of foreground triggers from H1H2L1 observation time in this search. The
dashed line shows the background distribution estimate averaged over all the time slides. The greyline give the distribution for each of the 100 background trial taken separately. The yellow bandshows the N 1/2 one-standard deviation errors on the average background estimate.

107

7.6

False Alarm Probability

Using the FARc as the detection statistic allows a very easy calculation of the false alarm probability(FAP) and the background probability P0 [130, 131]. Assuming a Poisson distribution, the FAP is
the probability of getting any (zero) background triggers louder than the loudest zero-lag trigger(i.e., with a FARc lower than the FARc of the loudest zero-lag trigger since low FARcs are moresignificant) and the background probability is the probability of getting zero background triggers
louder than the loudest zero-lag trigger, given by

FAP = 1 − ex(7.6a)

P0 = ex ,(7.6b)

where x = −FAR × T , FAR is the FARc of the loudest zero-lag trigger, and T is the total analyzed
time searching for zero-lag triggers. In-time coincident triggers with low FAP are GW detection
candidates, and are subject to follow-up, as described in section 10.1.

108

Chapter 8

Ranking of Detection Candidates
In this chapter we describe the postprocessing pipeline that has been developed as an extension
to the standard pipeline (discussed in chapter 6) to analyze data from GW detectors, such as the
LIGO detectors, for GW signals from CBCs. In addition, we will also describe the tuning that wasdone for the Search for Low Mass CBCs in the First Year of LIGO’s Fifth Science Run (S5) Data,described in this thesis, and referred to as “this search.”
At the end of the standard pipeline, we have distilled the data down from one data stream
sampled at 16,384 Hz for each GW detector into files containing coincident triggers that have passed
our thresholds and vetoes. These files are organized by coincident time segments defined by times
when diﬀerent combinations of detectors’ data were being analyzed.
Recalling section 6.8, during these coincident time segments triggers are discarded if they occur
during times flagged by the category 1, 2, or 3 DQ flag veto files. As mentioned in section 7.1,
these vetoes aﬀect the application of signal-based vetoes, which, in turn, significantly aﬀect thebackground distribution of triggers. Because of this, we separate the coincident times (and triggers
within those times) according to the detectors that were operating without DQ vetoes active. We
then recalculate the amount of analyzed time in each type of coincidence time for the in-time data
as well as for each time shift, accounting for the dead time incurred from the application of DQ
flag vetoes. Next we cluster coincident triggers in order to reduce the final number of triggers to
a manageable level, and separate them by categories (mass regions and coincident trigger types).
Finally, we calculate the FAR for each of the in-time, software injection, and time-shift coincident

109

Figure 8.1: Postprocessing Pipeline
The stages of the postprocessing producing the final candidate triggers.
triggers, utilizing the in-time and time shifts’ amounts of analyzed time. The postprocessing pipeline
is represented in figure 8.1.
The diﬀerent stages of the postprocessing pipeline used in this search are discussed below:• separate triggers by veto times (detailed in section 8.1);• calculate nonvetoed analyzed times (detailed in section 8.2);• cluster triggers and separate them by categories (detailed in section 8.3);• calculate the FAR for the triggers (detailed in section 8.4).

110

8.1

Separating Triggers by Data Quality Flag Vetoes

The application of DQ flag vetoes causes a loss of detector live time that was not accounted forin previous searches (it was instead counted as a loss of eﬃciency). Accounting for it as a loss oflive time requires us to reclassify the coincident category of coincident triggers in the in-time andtime-shifted data.
The first thing that is done in the postprocessing pipeline is to separate triggers according towhich detectors were producing triggers that did not have DQ flag vetoes active. As mentioned
above, because the DQ flag vetoes aﬀect which signal-based vetoes we apply, we want to group
double coincident triggers occurring during triple coincident time where the third detector had an
active DQ flag veto with double coincident triggers coming from true double coincident time.
We work with single coincident science segments in which we separate triggers in a manner
appropriate for both time-shift and in-time triggers. The output is a number of SEPTIME files foreach type of time when two detectors were not being vetoed (for H1H2L1 times, the output files
can be H1H2 times, H1L1 times, H2L1 times and H1H2L1 times; for H1L1 times, the output is only
H1L1 time). In the case of in-time data, it is also dependent on whether or not there was a nonzero
amount of time of that type.
Figure 8.2 shows a visual example of extracting the correct time-shift triggers from a single triple
coincident segment. The vector used for the time shifts here is v = {0, 2, 1} for H1, H2, and L1
respectively. The top black bar at y = 10 shows the segment we are working with is [0 − 10), thenext bars show the times when each of H1 (red), H2 (blue), and L1 (green) are not vetoed (H1 notvetoed [1 − 10), H2 not vetoed [1 − 10), L1 not vetoed [1 − 6)). Just below those lines show thein-time triggers that are shifted when computing coincidence in time-shifted data.
The middle lines and triggers from y = 5 to y = −5 show the time shifts numbered 5,4, . . . ,1,1, . . . -5. The colored bars again denote the times the detectors were not vetoed, shifted by theappropriate amount for that time shift and detector. The triggers under the time-shift bars are thein-time triggers that were coincident in that time shift, shifted to their appropriate location. If you
pick a coincident trigger in a specific time shift, you can look just above it to tell which detectors

111
were producing nonvetoed triggers.
The final four lines at the bottom denote the four diﬀerent times triggers could be generated(i.e., H1H2L1, H1H2, H1L1, and H2L1 times from top to bottom). The unshifted trigger times,
which can be compared to the input trigger times at the top of the plot, are shown on each of thosefour lines. The septime program outputs those four lines are four separate trigger files.
At the end of this stage, for each coincident segment, we are left with files whose contents arecoincident in-time or time-shifted triggers that occurred during coincident times defined as times
when those detectors were analyzed and not vetoed by DQ flag vetoes.

112

Figure 8.2: septime Example
An example of septime extracting time-shift triggers from the appropriate coincidence times once
DQ flag vetoes are taken into account.

113

8.2

Calculate Analyzed Times

Now that we have sorted the triggers according to times when detectors were analyzed and notvetoed by DQ vetoes, what we will now refer to as “analyzed times,” we also want to calculate theamount of analyzed times per in-time and time-shifted data. We do this using a program called
antime that uses the science segment files and appropriate veto files. When detectors are vetoed,that time is removed from that type of time and added to the type of time defined by the detectors
not being vetoed (e.g., H1H2L1 time where H2 is vetoed is added to H1L1 time).
Summarizing the in-time analyzed times, before the application of category 2 ans 3 DQ flagvetoes, the in-time analyzed times were 13,237,711 s for H1H2L1 times, 1,179,439 s for H1L1 times,
1,487,038 s for H2L1 times, and 7,337,970 s for H1H2 times. Once the category 2 and 3 vetoes aretaken into account we have the following times from triple coincident times: 11,696,445 s for H1H2L1
times (88.4% of the original), 76,446 s for H1L1 times (0.6%), 470,764 s for H2L1 times (3.6%), and790,885 s for H1H2 times (6.0%). In addition, from the double coincident times we have: 615,994 s
for H1L1 times (52.3%), 949,899 s for H2L1 times (63.9%), and 6,951,276 s for H1H2 times (94.7%).
In total, we have 11,696,445 s of H1H2L1 time, 692,440 s of H1L1 time, 1,420,663 s of H2L1
time, and 7,742,161 s of H1H2 time.
Figure 8.3 shows the in-time and time-shifted analyzed times after category 2 and 3 DQ flagvetoes have been removed. The strong dependence on time-shift number seen in H1L1 and H2L1
time is due to correlations between H1 and H2 DQ vetoes. Since these vetoes are positively correlated
due to their shared environment, as H1 and H2 times are shifted against each other during H1H2L1
times, time that was initially being vetoed by both H1 and H2 DQ vetoes in the in-time turns into
time where only one of the Hanford detectors is being vetoed for the time shifts. Since the H1H2L1
time initially has an order of magnitude more time than H1L1 or H2L1 time, the small eﬀect ofsliding correlated H1 and H2 vetoes against each other becomes a more significant amount of timetransfered to H1L1 and H2L1 analyzed times in comparison to their initial amount of analyzed
times. This eﬀect is more pronounced as one moves away from the nonshifted in-time data because
the relative shift between H1 and H2 (10 s) is comparable to the typical duration of DQ flags.

114

Figure 8.3: Background Analyzed TimeThe analyzed times as a function of time-shift number can be seen here after the category 2, and 3
DQ flag vetoes are taken into account for H1H2L1 (top left), H1H2 (top right), H1L1 (bottom
left), and H2L1 (bottom right) times. The red bar in each plot denotes the in-time analyzed time.

115

8.3

Trigger Clustering and Separation by Category

In the next stage of the postprocessing pipeline, we prepare the triggers for the FAR calculation.
In order to reduce the total number of triggers from the search to a reasonable size, we cluster the
SEPTIME files using a sliding cluster window of 10 s, using the combined eﬀective SNR as theintermediate detection statistic. The output of this thep is one COIRE file per SEPTIME file. This
incurs an insignificant loss of signal recovery eﬃciency since, generally, when there is a glitch in thedata, the whole template bank produces triggers, masking possible signals during that time. The
clustering we do focuses the glitch into the loudest trigger during that time, which generally has thelargest combined eﬀective SNR in the high mass region of the template bank. This eﬀect could be
the cause of the louder background combined eﬀective SNR distribution for the high mass region,
described in chapter 7.1.
Next we separate triggers into diﬀerent categories due to the diﬀerences between the backgroundtrigger combined eﬀective SNR distributions (chapter 7.1). In that chapter we found that thebackground combined eﬀective SNR distributions can be aﬀected by multiple factors including thecoincidence level (i.e., double versus triple coincidence), the signal-based vetoes applied, and the massof the templates producing the triggers. The first two of these three factors leads us to separate
triggers by the diﬀerent observation times and the diﬀerent coincidence types (i.e., H1H2L1, H1H2,
H1L1, and H2L1 triggers in triple coincident time and H1H2, H1L1, and H2L1 triggers in double
coincident times). In addition, we separate these triggers at the chirp mass divisions of M = 3.48
M⊙ and M = 7.40 M⊙ .
These trigger categories are implemented by taking the above clustered COIRE files, grouping
them together by observation time, separating the triggers into the three mass regions, and subdividing the triggers in triple times further into coincidence types (H1H2L1 triggers, H1H2 triggers,
H1L1 triggers, and H2L1 triggers).

116

8.4

FAR Calculation

The final stage of the postprocessing pipeline is the FAR calculation. The FAR is calculated using
the time shifts and analyzed times as described in section 7.2. This is done separately for eachtrigger category using lalapps corse (the COincident Rate Statistic Estimator). Once the FAR
is calculated, the categories for a particular observation time are recombined as separate trials andtheir FARc is calculated (described in section 7.3). In the end, the IFARc is used as our detection
statistic where triggers with larger IFARc values are louder.

117

Chapter 9

Rate (Upper Limit) Calculation
Whether the search results in a small number of detection candidates or more, we can bound the rateof CBCs in the (nearby) universe. Depending on the significance of the loudest event, the confidenceband we establish may or may not exclude zero. In either case, we can establish bounds on the rate
for diﬀerent source mass ranges.
In this discussion, we follow the loudest event formalism described in [132, 133, 130, 131]. One
advantage of calculating bounds on a rate using the loudest event formalism is that for searches
that have steep background distributions, the expected upper limit you could get is always more
stringent than in fixed threshold formalisms. Another advantage is that one does not need to fix athreshold before looking at the in-time triggers. It should be noted that in the end, the formalism
chosen does not aﬀect the significance of the in-time triggers and is only necessary for calculating
rate bounds.
In previous CBC searches [17, 19, 18], we calculated rates of CBCs in terms of Milky Way
Equivalent Galaxies (MWEGs). However, as the horizon distances of our detectors expanded toinclude more galaxies, we found it was prudent [134] to report rates in terms of L10 s, where 1 L10 ≡1010 L⊙,B , L⊙,B is the solar blue-light luminosity, and 1 MWEG ≈ 1.7 L10 . We have used theseunits for our rate calculations in [21, 22], and we shall use it here as well. The reason for this change
is because L10 track the rate of star formation in a given galaxy and we think CBCs are proportional
to this [135]. One drawback of this is that old elliptical galaxies are ignored since they do not havesignificant blue-light luminosity [136].

118
In order to calculate an upper limit on a rate of gravitational wave signals due to a particularsearch, there are several things we need to know. These are related to how sensitive the search is
to signals from the universe, the background of the search, and the loudest event from the search.
Below we detail the calculations performed in arriving at a rate (figure 9.1).

119

Figure 9.1: Upper Limit Calculation Pipeline
The stages of the upper limit calculation, in which we information from the top of the figure toestablish bounds on the CBC rate in the nearby universe in terms of events yr−1 L−1
10 .

120

9.1

Posterior and Upper Limit Calculation

Calculating an upper limit on a rate of coalescences in the loudest-event formalism requires knowledgeof the cumulative luminosity to which the search is sensitive and a measure of the likelihood thatthe loudest event was due to the observed background, described in section 9.3. We combine thesewith the time analyzed to calculate the posterior on the rate for the search. Assuming a prior on
the rate p0 (µ), the posterior is given by [130]p (µ|CL , T, Λ) = p0 (µ)

CL T(1 + µCL T Λ) e−µCL T ,
1+Λ

(9.1)

where µ is the rate in events yr−1 L−1
10 , CL is the cumulative luminosity in L10 , T is the analyzed timein years, and Λ is a measure of the likelihood of detecting a single event with loudness parameter x(which in this analysis is FARc, cf., sections 7.3 and 8.4) versus such an event occurring due to theexperimental background, given by [130]Λ (x) =

�

−1 dCL
CL dx

��

1 dP0
P0 dx

�−1

,

(9.2)

and P0 the background probability given by equation (7.6b). In this search, P0 = dP0 /dx = ex , thus
the term involving these quantities disappears, which we will discuss in section 9.3.2.
In the limit that Λ → 0 (i.e., when the loudest event is absolutely due to background), the 90%
confidence limit we set on the rate µ90% is given by90% =� µ90%
0p(µ)dµ ⇒ µ90% =

2.303
,
CL T

(9.3)

and in the limit that Λ → ∞ (i.e., when the loudest event is absolutely due to signal), the 90%
confidence limit we set on the rate µ90% is given byµ90% =

3.9
.
CL T

(9.4)

121

9.2

Incorporating Systematic Errors

The posterior (9.1) assumes a known value of CL associated with the search. In reality, CL has
associated with it systematic uncertainties, which we model as unknown multiplicative factors ζ(fractional errors), each log-normally distributed about 1 with errors described in section 9.4. In
order to marginalize over the eﬀects of the systematic errors on the cumulative luminosity, we first
calculate a probability distribution for the cumulative luminosity for the search, pd (CL ). This iscalculated using the log-normal error distributions of the cumulative luminosity multiplier pd,error (ζ)
given bypd,error (ζ) =1√

ζσerror 2π− ln(ζ)2e 2σerror,(9.5)

where σerror is the fractional error of the cumulative luminosity calculated in section 9.4. The
pd,error (ζ) distributions are then convolved with each other to obtain a combined distribution
pd,comb (ζ). The distribution of the cumulative luminosity is then given as

pd (CL ) =�δ (CL − ζCL,m ) pd,comb (ζ)dζ ,(9.6)

where CL,m is the measured cumulative luminosity. To obtain a marginalized posterior, we integrate
the above posterior (equation (9.1)) times the distribution of cumulative luminosities (equation (9.6))
over the cumulative luminosity [130]p (µ|T, ζ) =�

pd (CL ) p (µ|CL , T, ζ)dCL .

(9.7)

The results of several experiments (e.g., diﬀerent types of S5 observing time and previous runs
such as S3 and S4) can be combined by taking the product of their likelihood functions; in the caseof uniform priors, this is equivalent to taking the product of their posteriors, allowing us to define
the rate upper limit µ at a confidence level α by numerically solving

α=

� µ�
0

i

pi (µ� )dµ� ,

(9.8)

122
where the pi (µ� ) are the marginalized posteriors from diﬀerent experiments calculated using a uniform prior on the rate. This is equivalent to using the posterior from previous observations as theprior for subsequent ones, sequentially.

9.3

Cumulative Luminosity and Background Probability

The standard quantities we need to have before we can calculate an upper limit on the rate ofcoalescences are the sensitivity of the search (i.e., the cumulative luminosity), and the probability
that no background triggers are louder than the loudest event.

9.3.1

Cumulative Luminosity and its Derivative

We measure the cumulative luminosity that the search is sensitive to in units of L10 . The cumulative
luminosity quantifies the potential sources of observable CBC, as measured by blue-light luminosityof the galaxies containing CBCs, which can be detected by our search. This can be calculated byintegrating the eﬃciency �(D) and physical luminosity L(D) both as a functions of distance overdistanceCL ≡�

�(D)L(D)dD .

(9.9)

In practice, we do this as a discrete sum

CL =

�

�i Li ,

(9.10)iwhere �i is the eﬃciency in distance bin i, and Li is the luminosity in distance bin i.
We calculate the eﬃciency �i in distance bin i as

�i =Nfound,i,Ninjected,i(9.11)

where Ninjected,i is the number of injections injected with distance D ∈ (Di,min , Di,max ), Di,min is theminimum distance in distance bin i, Di,max is the maximum distance in distance bin i, and Nfound,i123
is the number of found injections in that distance bin. Distance bins in which there are no injections
are taken to have �i = 0. An injection is counted as found if there is a trigger associated with theinjection with a detection statistic louder than the detection statistic of the loudest in-time trigger.
The physical luminosity is calculated by creating an “injection” file using a galaxy catalog. The“injections” are distributed according to the locations and luminosities of galaxies in the catalog.
More luminous galaxies end up have more “injections” associated with them. The luminosity weightW each “injection” carries is determined by the most luminous galaxy in the catalog, galaxy j. W
is given by
W = Nj /Lj ,(9.12)

where Nj is the number of “injections” in galaxy j, and Lj is the luminosity of galaxy j. The
physical luminosity Li in distance bin i is then given as

Li = ni W ,(9.13)

where ni is the number of “injections” from all galaxies with distance D ∈ (Di,min , Di,max ).
The derivative of the cumulative luminosity is computed numerically by recalculating CL using
several values of the detection statistic x around the largest value from the in-time triggers.

9.3.2

Background Probability

The background probability takes a simple form with the choice of the FARc for the detection statistic.
Assuming a Poisson distribution, the probability of getting zero background triggers louder than theloudest in-time trigger (i.e., with a FAR lower than the FAR of the loudest in-time trigger since low
FAR are more significant) is given by equation (7.6b) where, as before, x = −FAR × T , FAR is the
FAR of the loudest in-time trigger, and T is the total analyzed time searching for in-time triggers.
Since more significant triggers are given by lower FAR, triggers with higher values of x are
louder, so the derivative with respect to the detection statistic in the direction of louder triggers isin the positive x direction (d/dx). The derivative of the background probability with respect to the124detection statistic in the direction of louder triggers dPb /dx is

dPb= ex .
dx(9.14)

Combining equations (7.6b) and (9.14), we find that (1/P0 )(dP0 /dx) = 1, and this term in equation(9.2) goes away.
We therefore can compute all of the quantities in equation (9.1) and 9.2 needed to determine theposterior on the rate µ, before marginalizing over systematic errors.

9.4

Systematic Error Calculation

Systematic errors associated with CBC searches for GW signals include errors associated with detector calibrations, simulation waveforms, Monte Carlo statistics, and galaxy catalog distances andmagnitudes. Calculating these errors in terms of the cumulative luminosity is described below [131].
We can marginalize over these uncertainties using equation (9.7), to obtain the final posterior on therate, and use equation (9.8) to obtain the upper limits. This analysis can be repeated as functionof source mass.

9.4.1

Monte Carlo Errors

We refer to statistical errors associated with the eﬃciency calculation as Monte Carlo errors. Since
we calculate the eﬃciency as a function of distance, we calculate the error for a particular distancebin i using the binomial formula, which gives an error of zero when the eﬃciency is zero or one, or
when there are no injections in that bin:

�Monte Carlo,i =�Nfound,i (Ninjected,i − Nfound,i ),3
Ninjected,i+ 10−5(9.15)

where �Monte Carlo,i is the error in the eﬃciency, Ninjected,i is the number of injections injected into
distance bin i, and Nfound,i is the number of injections found in distance bin i. Finally, the associated125
fractional luminosity error σMonte Carlo is given in terms of �Monte Carlo,i , the physical luminosity ineach distance bin Li , and the measured cumulative luminosity CL,m as
σMonte Carlo =

9.4.2

1
CL,m

�

(9.16)

�Monte Carlo,i Li .

i

Calibration Errors

Calibration errors in the detectors are errors due to uncertainties in the absolute calibration of
the detector response to diﬀerential displacements caused by GWs [137]. These errors aﬀect theamplitude, and in turn the distance, at which we made injections to calculate the eﬃciency of oursearch, since the injections we made assuming a specific value of the noise floor. During S5, the onesigma uncertainty in the amplitude (and thus the distance) associated with the calibration was 8.1%
for H1, 7.2% for H2, and 6.0% for L1 [137]. To calculate the luminosity error due to a calibration
error ∆calibration , all of the injections are moved further in that detector’s distance by ∆calibrationD� = D (1 + ∆calibration ) ,(9.17)and the eﬃciency is recalculated using equation (9.11) where Ninjected,i is the number of injections
injected with distance D� ∈ (Di,min , Di,max ), Di,min is the minimum distance in distance bin i,
Di,max is the maximum distance in distance bin i, and Nfound,i is the number of found injections
with distance D� ∈ (Di,min , Di,max ). This gives us �calibration,i where the i denotes a particular bin
in distance. The standard eﬃciency �i is then subtracted from �calibration,i on a bin-by-bin basisd�calibration,i = �calibration,i − �i .

(9.18)

Finally, the associated fractional luminosity error σcalibration is given in terms of d�calibration,i , Li ,and CL,m as
σcalibration =

1
CL,m

�
�
i

d�calibration,i Li�.

(9.19)

126

9.4.3

Waveform Errors

Waveform errors are associated with how diﬀerent the true signals are from the signals we use tomeasure the eﬃciency of our pipeline (i.e., the mismatch between the true signals and our injections).
This error eﬀectively reduces the distances in our eﬃciency calculation since we do not recover allof the power available in the signal due to the mismatch between the signal and our injections. To
calculate this, all of the injections are moved further by the error ∆waveformD� = D (1 + ∆waveform ) ,(9.20)and the waveform fractional error is calculated as above (section 9.4.2). We calculate the waveform error in units of luminosity assuming a waveform mismatch of 10%. This value is obtainedby a detailed comparison of a variety of PN waveform approximants [8, 138], and it represents a
conservative estimate of this theoretical error on the true waveform.

9.4.4

Galaxy Errors

Galaxy errors are errors associated with our galaxy catalog [134] used to construct the physicalluminosity. Galaxy errors come in two types: distance errors and magnitude errors.
Future searches will be normalized per unit volume (Mpc3 ) to obtain a rate density, rather thanper unit L10 , in order to avoid these errors (which are not intrinsic to this search method). In thelarge volume limit (i.e., above a physical distance of 30 Mpc), there is a 0.02 L10 Mpc−3 conversion
that can be used to compare the two rate units.

9.4.4.1

Distance Errors

To calculate the error on the luminosity due to distance errors, the physical luminosity calculationis changed such that the distance to each of the “injections” in a galaxy is increasedD� = D (1 + ∆D ) ,(9.21)

127
where ∆D is the fractional distance error for a given galaxy. Also, the “injections” in a particulargalaxy’s weighting is changed to be2

W � = W (1 + ∆D ) ,(9.22)

where W is the original weighting for the galaxy. This form of weighting is chosen because thegalaxy’s luminosity is only known in terms of its magnitude and distance. If there is an error in the2

distance, the luminosity must change by (1 + ∆D ) in order to maintain the same magnitude. The
luminosity Ldistance,i is calculated with weighting W � and with “injections” at distances D� as�

Ldistance,i =n�(9.23)W� ,i=1

where n� is the number of “injections” from all galaxies with distance D� ∈ (Di,min , Di,max ).
Ldistance,i is then multiplied by the standard eﬃciency �i as a function of distance and integrated
over distance to obtain the cumulative luminosity CL,distance at 1σ error in distanceCL,distance =���i Ldistance,i

i

�

.

(9.24)

The associated fractional luminosity error σdistance in terms of CL,m as
σdistance =9.4.4.2

1
CL,m(CL,distance − 1) .

(9.25)

Magnitude Errors

To calculate the error on the luminosity due to magnitude errors, the physical luminosity calculationis changed such that the “injections” in a particular galaxy are weighted by

W � = W 10∆M /2.5 ,(9.26)

128
where ∆M is the (absolute) magnitude error for a given galaxy. This follows from the formula used
to calculate the luminosity of a galaxy

LB,i = LB,⊙

�

Di
10pc

�2

10(MB,⊙ −mB,i )/2.5 ,(9.27)

where LB,⊙ is the blue solar luminosity, Di is the distance to the galaxy, MB,⊙ is the (absolute) blue
solar magnitude, and mB,i is the (measured) apparent blue magnitude of the galaxy. The luminosity
Lmagnitude,i is calculated with weighting W � as

Lmagnitude,i =

n
�

W� ,

(9.28)

i=1

where n is the number of “injections” from all galaxies with distance D ∈ (Di,min , Di,max ). The
magnitude fractional error is calculated as above (section 9.4.4.1).

129

Chapter 10

Results of the Search for Low Mass
Compact Binary Coalescences in
the First Year of Data from
LIGO’s Fifth Science Run
In this chapter we describe the results from the Search for Low Mass CBCs in the First Year of
LIGO’s Fifth Science Run (S5) Data. In section 10.1 we discuss the follow-up procedures we use tobuild confidence in a possible GW event. We discuss the loudest events from the search in section
10.2. Finally, we present upper limits on the rate of coalescences in section 10.3.
During the analysis, and prior to unblinding the nonplayground data, we discovered an error
in the computation of the template metric. This metric is used in the placement of the bank andthe coincidence test. The error caused the metric distance between templates to be overestimatedfor the higher mass signals. This has the eﬀect of causing the template placement algorithm toover-cover the higher mass region (i.e., to produce a bank with less than the requested 3% loss
in SNR). This increased the computational cost of the search, but did not significantly reduce
the sensitivity. However, this error also aﬀected the coincidence algorithm by overestimating the
distance between triggers for high mass signals. Since the coincidence window was empirically tuned
on software injections and time-shifted coincidences the impact on the sensitivity of the search
was not significant. Consequently, the decision was taken to unblind the data using the original,
suboptimal analysis in order to begin studying any possible detection candidates and to use this130
result to compute the upper limit (in the absence of a detection). The decision was also taken to
perform a complete reanalysis of the data with the corrected metric to verify the (non)detection
statement from the original search by checking for any change in the list of the loudest in-timecoincident triggers (the detection candidates).

10.1

Follow-up Procedure for Coincident Triggers

The distribution of IFARc in-time coincident triggers is examined to identify the events with thelargest IFARc (the “loudest” events). Events with a FARc much smaller than 1 per observation
time are deemed significant. However, this takes into account only the GW channel data filteredby the search pipeline. Information from other channels are only used coarsely in the definitionof the category 2 and 3 DQ flags (section 6.8.1). Much more information is available in both the
GW channel and the auxiliary channels, but it is not yet automatically folded into the detection
statistic (this is a work in progress, which will be discussed more in section 11.1). We therefore havedeveloped a rather detailed follow-up procedure to examine this additional information in detail,albeit, in a necessarily subjective manner. Automating this procedure for time-shifts and simulated
GW signals (i.e., injections) to provide an improved, unbiased detection statistic is in progress.
In the meantime, we check our loudest coincident triggers with a checklist of tests designed to see if
a statistically significant trigger is believable as a detection candidate. The methods employed in thischecklist are tested against injection and time-shift coincidence triggers. Reference [139] describes
the tests that we perform on the coincident triggers and the data surrounding them. At present,
our standard tests include the following: We check the integrity of the data for corruption. We also
check the status of the detectors and the presence of any data quality flags in the surrounding data.
We assess whether there could have been environmental or instrumental causes found in auxiliary
channels at the time of the trigger. We check the appearance of the data at the time of the trigger
in the form of the strain, SNR, and χ2 time-series, and time-frequency spectrograms (an example isshown in figure 10.1).
In addition, for any statistically significant candidate that survives the tests listed above, we plan

131

Figure 10.1: Comparison of Triggers from Loudest Event to an Injected Signal
A “Qscan” of the triggers in the loudest event (top) in this search after category 2 and 3 DQ flagvetoes have been applied and of an signal injected into the data (bottom). The Qscans show
time-frequency content in the form of power in a Sine-Gaussian basis.
to do the following: Assess the coherence between the signals recorded by each individual detector
operating at the time of the event. Verify the robustness of the trigger against small changes in thepipeline (i.e., changes in the adjacent Fourier transform boundaries or changes in the calibration of
the data). Check the robustness across pipelines by employing other search techniques to analyze
the same data (i.e., CBC pipelines using diﬀerent templates or pipelines designed to search for
unmodeled bursts). Finally, we will check for coincidence with external searches for gamma-ray
bursts, optical transients, or neutrino events. (This last test is for information only, as a genuine
GW event might or might not be accompanied by other signals.)
We examine the distribution of in-time triggers after each category of DQ flag veto is applied. Incase there is a statistically significant outlier before category 2 or categories 2 and 3 DQ flag vetoes

132
are applied, we carry out a follow-up exercise to see if the veto that eventually rejected the event
was rightfully applied. There are two reasons that this could be important. First, a very strong GWfrom within the Milky Way could cause an instrumental saturation of the sort that we use as a veto;this kind of problem would be easy to diagnose if it were to occur, since the signal would be strong
enough for us to see in the moments leading up to the signal-induced saturation. Second, we want toguard against false dismissal of a candidate by other kinds of vetoes, which can have nonnegligible
dead-time associated with them. Some of our vetoes are associated with recognizable forms of false
signals; we check to be sure that a vetoed loud event looks like that kind of false signal, and not
like a genuine coalescence signal. In this search, there was a single statistically significant outlier inthe distribution of events before the application of categories 2 and 3 DQ flag vetoes. The follow-up
exercise confirmed that the category 3 DQ flag that vetoed that event was correctly applied, as
it occurred during a time that was especially glitchy and was recommended to be removed fromanalysis as the data was taken.

10.2

Loudest Events

At the end of our pipeline we are left with a set of coincident triggers that are potential detection
candidates. The cumulative distribution of events above an IFARc versus IFARc for the diﬀerent
observation times is shown in figure 10.2; the H1H2L1 distribution is also shown in figure 7.10.
This figure shows that the loudest candidates in all three observation times were consistent withthe estimated background and thus were likely accidental coincidences. Thus, the search yielded no
detection candidates. The results of the reanalysis were consistent with the original analysis anddid not produce any plausible GW signals. We report an upper limit on the rate of coalescences in
section 10.3.
As an exercise to prepare for future detections, we carried the loudest several events (such as
the three loudest events that appear in each of the histograms in figure 10.2) through the follow-up
procedure.
Even though we know our background is underestimated for H1H2 coincident triggers, we re-133viewed the two loudest H1H2 candidates using the detection checklist. In both of those cases, thewaveforms from the two interferometers failed to match each other in detail, thus ruling them out
as GW events.

134

Figure 10.2: IFAR Loudest Events
The cumulative distribution of events above a threshold IFAR, for in-time coincident events, shownas blue triangles, from all coincidence categories for the observation times H1H2L1 (top left), H1L1
(top right), and H2L1 (bottom right). The expected background (by definition) is shown as adashed black line. The 100 experimental trials that make up our background are also plotted
individually as the solid grey lines. The yellow band shows the N 1/2 one-standard deviation errors
on the average background estimate.

135
Table 10.1: Detailed Results of the BNS Upper Limit Calculation
BNS
Dhorizon = 30 Mpc
Coincidence Time
H1H2L1 H1L1 H2L1
Cumulative Luminosity (L10 )
250
230
120
Calibration Error
21%
3.9%
16%
Monte Carlo Error
5.4%
16%
13%
Waveform Error
26%
11%
20%
Galaxy Distance Error
14%
13%
6.1%
Galaxy Magnitude Error
17%
17%
16%
Λ [Eq. 9.2] �
0.30
0.41
0.72
�
−2
−1 −1
3.9 × 10
Marginalized Upper Limit yr L10
Summary of the search for BNS systems. The horizon distance, as defined in Eq. 5.38, for the 4
km detectors and averaged over the search. The cumulative luminosity is rounded to two
significant figures. The errors in this table are listed as logarithmic errors in the luminosity
multiplier based on the cited sources of error.
Table 10.2: Detailed Results of the BBH Upper Limit Calculation
BBH
Dhorizon = 80 Mpc
Coincidence Time
H1H2L1 H1L1 H2L1
Cumulative Luminosity (L10 )
4900
4900
1800
Calibration Error
23%
36%
25%
Monte Carlo Error
3.2%
11%
9.9%
Waveform Error
27%
34%
27%
Galaxy Distance Error
18%
18%
20%
Galaxy Magnitude Error
16%
16%
17%
Λ [See (9.2)] �
0.59
1.2
1.4
�
−3
2.5
×
10
Marginalized Upper Limit yr−1 L−1
10
Summary of the search for BBH systems. The rows of this table are defined as in Table 10.1.

10.3

Coalescence Rate Upper Limits

In the absence of detection, we set upper limits on the rate of CBCs per unit L10 , for several
canonical binary masses, as a function of the total mass of the compact binary system, and as afunction of the black hole mass for BHNS systems.
For each mass range of interest, we calculate the rate upper limit at 90% confidence level (CL)using the loudest event formalism [133, 130], described in chapter 9. We derive a Bayesian posterior
distribution for the rate as described in chapter 9 and reference [130]. The cumulative luminosity
for this search, calculated as described in chapter 9, can be found in Table 10.1, 10.2, and 10.3 forthe BNS, BBH, and BHNS canonical mass CBC systems, respectively.

136

Table 10.3: Detailed Results of the BHNS Upper Limit Calculation
BHNS
Dhorizon = 50 Mpc
Coincidence Time
H1H2L1 H1L1 H2L1
Cumulative Luminosity (L10 )
990
980
390
Calibration Error
22%
47%
25%
Monte Carlo Error
3.6%
9.7%
11%
Waveform Error
26%
40%
27%
Galaxy Distance Error
17%
16%
17%
Galaxy Magnitude Error
17%
18%
18%
Λ [See (9.2)] �
0.45
1.0
1.1
�
−2
1.1
×
10
Marginalized Upper Limit yr−1 L−1
10
Summary of the search for BHNS systems. The rows of this table are defined as in Table 10.1.

137

10.3.1

Upper Limits Neglecting Spin

We apply the upper limit calculation to three canonical binary masses as well as calculating the upper
limit as two functions of mass. Our three canonical binary masses are BNS (m1 = m2 = (1.35 ± 0.04)
M⊙ ), BBH (m1 = m2 = (5 ± 1) M⊙ ), and BHNS (m1 = (5 ± 1) M⊙ , m2 = (1.35 ± 0.04) M⊙ ). We
use Gaussian distributions in component mass centered on these masses, with standard deviations
following the ± symbols. The neutron star mass distribution was chosen according to an analysisof 50 known pulsars in binary systems at the time (reference [140]). However, the mass distribution
for black holes is still uncertain, therefore we have chosen a canonical value to represent stellar mass
black holes, which allows us to compare results across diﬀerent science runs.
We combine the results of this search from the three diﬀerent observation times in a Bayesianmanner, described in section 9.1, and the results from previous science runs [21, 22] are incorporated
in a similar way.
Assuming that spin is not important in these systems, we calculate upper limits on the rateof binary coalescences using our injection families that neglect spin (section 6.2). There are anumber of uncertainties that aﬀect the upper limit, including systematic errors associated with
detector calibration, simulation waveforms, Monte Carlo statistics, and galaxy catalog distances andmagnitudes [131]. We marginalize over these as described in section 9.4 and obtain upper limits on
the rate of binary coalescences of

R90%,BNS = 3.9 × 10−2 yr−1 L10 −1 ,(10.1a)R90%,BBH = 2.5 × 10−3 yr−1 L10 −1 ,(10.1b)R90%,BHNS = 1.1 × 10−2 yr−1 L10 −1 .

(10.1c)We also calculate upper limits for two additional cases: as a function of the total mass of thebinary, with a uniform distribution in the mass ratio q = m1/m2, and as a function of the massof the black hole in a BHNS system, holding fixed the mass of the neutron star at mNS = 1.35
M⊙ (figure 10.3).

138

Figure 10.3: Upper Limits versus Mass
Upper limits on the binary coalescence rate per year and per L10 as a function of total mass of thebinary system with a uniform distribution in the mass ratio (left) and as a function of the mass ofa black hole in a BHNS system with a neutron star mass of mNS = 1.35 M⊙ (right). The darker
area shows the excluded region after marginalization over the estimated systematic errors. The
lighter area shows the additional region that would have been excluded if the systematic errors of
our most sensitive observation time had been ignored.

139

10.3.2Upper Limits Including Spin

In the previous section, we reported upper limits on the rate of mergers for diﬀerent classes of objects
using injection waveforms generated assuming nonspinning objects. We can also evaluate the upper
limits using injection waveforms that take into account the eﬀects of spinning bodies.
Since the maximum possible rotational angular momentum S for a black hole of mass m is Gm2 /c,
it is useful to describe the spin of a compact object in terms of the dimensionless spin parameter��â = (cS) / Gm2 . The distribution of black hole spin magnitudes within the range 0 ≤ â ≤ 1,
as well as their orientations relative to binary orbits, is not well constrained by observations. To
illustrate the possible eﬀects of BH spins on our sensitivity to BBH and BHNS signals, we provide
an example calculation using a set of injections using the SpinTaylor approximant, which calculates
waveforms for systems whose component objects can include significant spin angular momentum.
The spin magnitude â for this set of injections is uniformly distributed between 0 and 1 (section 6.2).
However, assuming a canonical mass and uniform density for a component object, astrophysicalobservations of neutron stars show typical angular momenta corresponding to â � 1 [141]. In
addition, the spin eﬀects are found to be weak for the frequency range of interest for LIGO [16] so
the BNS upper limits in section 10.3.1 are valid even though we have ignored the eﬀects of spin.
Using the above injections, we obtain marginalized upper limits on the rate of binary coalescencesof

R90%,BBH = 3.2 × 10−3 yr−1 L10 −1 ,(10.2)

R90%,BHNS = 1.4 × 10−2 yr−1 L10 −1 .

(10.3)

We can see that the eﬀect of spin on the relatively low mass systems targeted in this search is rather
small, modifying the average distance we can see by a factor of order ∼5%–10%, and thus our upper
limit by a factor that is three times larger (since the volume we are sensitive to is proportional tothe cube of the distance) of order ∼15%–30%.

140

10.4

Summary of the Search

We can compare the above results to a prediction made for the upper limit on the rate of BNS
coalescences from reference [131]. The example upper limit in that reference is computed assuming
a horizon distance of 35 Mpc for the detector, an observation time of 1 yr, and a loudest eventfrom the analysis with a combined SNR of 10. The result they predict, after marginalization oversystemic errors is 7.3 × 10−3 yr−1 L10 −1 . If we adjust the parameters above to match the sensitivity(average horizon distance of 30 Mpc) and duration (0.37 yr of analyzed data) of this search, thepredicted result becomes 3.1 × 10−2 yr−1 L10 −1 , which is very close to the upper limit we obtain inequation (10.1a).
Aside from possibly detecting GWs, this is what Initial LIGO was proposed to achieve andwhat it was built for. S5 was the first science run since the LIGO detectors reached their design
sensitivity and also the first long science run LIGO has undertaken. The search for BNS GW signals,
the flagship search, concluded with a result as sensitive as expected, a great achievement for LIGO
and the LIGO Scientific Collaboration (LSC).
We can also extend this prediction to BBH systems by adjusting the predicted result by thefactor (MBNS /MBBH )5/2 , which comes from the mass dependence of the amplitude of a CBC GW inequation (2.43). This results in a predicted result for BBH systems of 1.2 × 10−3 yr−1 L10 −1 . Thisis further from the obtained result in equation (10.1b) than the BNS comparison, however this is
expected due to the higher FAR for higher mass templates. This shows that the IFARc detection
statistic prevents the contamination of the lower FAR for lower mass templates from the higher FAR
from the higher mass templates.

141

Chapter 11

Future Improvements for Compact
Binary Coalescence Searches
In previous chapters we focused on the procedures used for the Search for Low Mass CBCs in the FirstYear of LIGO’s Fifth Science Run (S5) Data, referred to as “this search,” including data analysis
improvements that were made based on the collective experience of analyzing previous science runs’
data. In this chapter, we focus on such improvements that are yet to be made, or in the processof being implemented within the LSC CBC Group that will further benefit future searches. These
include improvements in separating signals from noise with better detection statistics (section 11.1),
better estimating the background of a particular search (section 11.2), coherently combining thedata from multiple detectors (section 11.3), cleaning the data before analysis to remove detector
glitches (section 11.4), and lowering the latency of the searches to improve the scientific output(section 11.5).

11.1

Detection Statistic

The goal of developing new detection statistics is to get to the ideal detector response to Gaussian
noise. We are already close to this ideal for the low mass search, and in particular the low mass
portion of the low mass search. This is because the low mass waveforms a more broadband than thehigh mass waveforms, which makes the standard SNR calculation relatively insensitive to detector
glitches. As was seen in this search’s waveforms, including information about the background noise

142of the detectors allowed us to develop a detection statistic that prevents the contamination of the
low mass portion of parameter space by the inclusion of the high mass portion.
The addition of the background information can be seen as a step toward defining a likelihood
ratio for the detection statistic. The likelihood ratio Λ is the ratio of the probability that a particulartrigger c was associated with a true signal P (c|h) and the probability that a particular trigger was
due to background noise alone P (c|0) given byΛ=

P (c|h).
P (c|0)(11.1)

As we saw in section 9.3.2, the FARc detection statistic can be converted to a FAP (i.e., theprobability of getting any background triggers louder than that trigger) using equation (7.6a).
The final step toward a likelihood ratio detection statistic is calculating P (c|h) using the injectedsignals. This procedure is being demonstrated to gain additional sensitivity in the triggered search
for CBC signals associated with short GRBs [142]. The reason for this is that including the P (c|h)
term introduces more information that we can use to distinguish noise triggers from signal triggers.
However, there is a downside of introducing this term. As in calculating P (c|0) where included
information about how the pipeline responded to noise in the GW channel, for P (c|h) we include
information about how the pipeline responds to signals in the GW channel. However, there is adiﬀerence between the two calculations. For the former, there was a well defined way to characterize
the noise in the detectors (i.e., time shifting the data before analysis). Whereas for the latter,
since we do not know the astrophysical distribution of signals we end up making a choice for the
distributions of signals. The resulting likelihood ratio detection statistic then depends more strongly
on those assumptions of the signal distributions than the previous detection statistic.
A diﬀerent approach toward an improved detection statistic is to use a multivariate classifier
to rank order the triggers. Such a classifier develops a classification algorithm by taking in alarge number of parameters (such as the many recovered parameters of the coincident triggers werecord and information from auxiliary channels if the detectors). It does this for both background143
triggers and injection triggers and then looks for how the two classes of triggers separate in themultidimensional space. Diﬀerent classifications of this sort are currently being tested using triggers
from the analysis of this search. This technique can be powerful because it allows easy inclusion
even more information beyond that included in the present detection statistic, ρeﬀ (equation (5.40)),and the associated background and signal probability factors P (c|0) and P (c|h).

11.2

Background Estimation

One of the main questions we ask at the end of a search is “what is the probability that a particulartrigger came from the background noise?” As described above, the FARc detection statistic has asimple answer to this question, however the current procedure of calculating the FARc is limited in
the smallest FAP it can estimate, when there is only one background trigger with significance equal
to or greater than the trigger in question. Naı̈vely, since we perform 100 time-shifted analyses toestimate the background, we would expect to be able to estimate a minimum FAP of 1%. Unfortunately this was under the assumption that we had a single experiment, while in fact the FARc
calculation combines the results of nine diﬀerent trigger categories, which are eﬀectively diﬀerent
experiments. We then have to multiply the minimum FAP by this trials factor giving us a minimum
FAP that we can estimate of 9%.
One way to get around this problem and lower the minimum FAP is to increase the number oftime-shifts we perform. However there are two problems associated with this. First of all, we want toestimate the background for triggers with roughly the same combination of detector noise conditions.
This means that we cannot shift the data by too large an amount, otherwise the noise characteristics
from multiple detectors will no longer be overlapping. Since there are time correlations in the datafrom a single detector, we need to shift the data of multiple detectors by a minimum amount so we
do not include those correlations in our background estimation. The combination of the two eﬀects
gives us a maximum number of time-shifts we can perform. Second, on taking into account practicalconsiderations, as we increase the number of shifts, the amount of time it takes to analyze the data
increases roughly by the same factor.

144
A diﬀerent way to estimate the background is to build a model of the background that we
could then use to extrapolate to low FAPs. One way to do this is to produce a joint probability
distribution function (PDF) for the background starting from the PDFs from the original singledetector triggers. As stated in chapter 5, for each detector the PDF of the SNR squared can bemodelled as a χ2 distribution with two degrees of freedom for the Gaussian noise component. Since
real data also contains a nonGaussian component that dominates in the tail of the distribution,
we would add a Poisson distribution to capture these tails. To factor in the eﬀect of signal-based
vetoes, we would need to run those portions of the pipeline with and without the vetoes to see
their eﬀect in each detector. We would then multiply the PDF of each detector by the ratio of
the two to get the PDF including vetoes. Finally, since we would want to rank order the triggers
according to eﬀective SNR, we would construct the theoretical distribution for the value of the χ2veto (i.e., χ2 distribution with 2p − 2 degrees of freedom and noncentrality parameter given by thatfor glitches), and combine it with the SNR PDF to obtain an eﬀective SNR PDF for single detectors.
The joint PDF for coincident triggers could then be taken as the outer product of the single-detector
distributions, normalized such that its integral gives a rate that matches the number of triggers in
the time-shifted analyses divided by the sum of the time-shifted analyses observation times.

11.3

Coherent Analysis

Some of the more powerful vetoes we use in analyzing GW-detector data are the amplitude consistency vetoes between interferometers. These vetoes ensure that the strength of a trigger is consistent
between diﬀerent detectors. However, there is also phase information that can be compared. Thisadditional information can be included by matched filtering a coherent data stream, as described in
section 3.6 and references [143, 144].

145

11.4

Data Cleaning and Noise Regression

Currently, all of the loudest triggers produced by analyzing LIGO data are the result of glitches inthe detectors. As described in section 6.8.1, we remove times during which we expect the detectors
to be glitching. Currently, as the signal processing we use for templated CBC searches introduces
artifacts that can last for a few seconds around a glitch, there are wide windows we must add tothese times to ensure that triggers resulting from such glitches are removed. These windows canmake up a significant portion of each time we wish to veto, so trying to reduce or eliminate thesewindows could help save additional analysis time. One way to do this is to switch over from vetoing
such times to eliminating the glitches from the data processing. This could happen if instead of
vetoing the required times, plus their windows, the required times’ data was excised and replacedwith the appropriately colored Gaussian noise such that the transients from the excision boundaries
were minimized. These times could then be flagged as not involving real data but would allow theanalysis of adjacent times.
Another option, which would be much better when possible, would be the regression of noiseusing information from the auxiliary data channels. For situations where the transfer function
between the auxiliary data channel and the strain data channel are known, this would allow excess
noise from these sources to be virtually eliminated from strain data and improving sensitivity. For
channels with a time-varying coupling to the strain data that is recorded, adaptive filtering could
be used to eliminate both excess noise as well as transients caused by an increase and then decreaseof the coupling.

11.5

Low Latency Searches

Improving the latency of CBC searches for GWs would have a significant scientific impact. Currently,
there is an inherent latency in the trigger generation of roughly 30 minutes, since the analysis
proceeds in batch mode analyzing 2048 seconds of data at a time (chapter 6). This is a longenough time that any electromagnetic counterparts could have faded into the background by the146time telescopes received an announcement of a likely GW. There are several techniques that are
being utilized to reduce this latency. References [145, 146] take advantage of the chirp-like natureof CBC GW signals to break the signal up into diﬀerent bands, which can be combined to generatea trigger with an SNR that is roughly the same as the current full band analysis. These techniques
can reduce the latency of the search to a few seconds. Additionally, for BNS signals in Adv. LIGO,
if only the lowest frequency bands are used, a trigger for a GW could be produced a few minutes
before the arrival of the signal from the neutron stars’ collision, alerting the LIGO detectors toremain online for the next few minutes, potentially even changing the response of the detectors to
the high frequency portion of the inspiral by dynamically detuning the signal extraction opticalcavity, and alerting electromagnetic telescopes where and when to look, resulting in multimessenger
observations of the same event.

147

Chapter 12

Testing Eﬀects of a Massive
Graviton
The General Theory of Relativity has proven to be very predictive over the past century. As esotericas it may seem, its eﬀects can impact our every day lives, such as the General Relativistic corrections
needed for GPS navigation to work. However, we know the predictions of General Relativity must
break down at some small length scale since quantum eﬀects will become important at the Planckscale.
General Relativity could also have flaws at large scales. It has also been observed that stars
within galaxies and galaxies within galactic clusters orbit around the center of mass faster than theyshould based on the amount of visible matter [147, 148]. There are one of two conclusions that canbe drawn from this, either there is more matter there than we can see, or General Relativity does
not describe gravity at large distances.
Many people currently believe the former of these two options, positing the existence of dark
matter to explain the missing mass. The latter option requires a modification to General Relativity
at large distances, which could be accounted for in the form of a massive graviton.
In this chapter we will compare the analysis of Will [149] using 1.5PN waveforms to a new
analysis using complete inspiral-merger-ringdown (IMR) waveforms. The rest of this chapter is asfollows: we briefly summarize the eﬀects of a massive graviton (section 12.1), give a dispersion
relation for GWs in general (section 12.2) and specific to the analysis of GWs from CBC signals(section 12.3 and 12.4), use this dispersion relation to predict theoretical bounds for the mass of

148
the graviton from future CBC GW observations (section 12.5), and then compare these bounds tocurrent observational bounds and theoretical bounds from other techniques (section 12.7).

12.1

Eﬀects of Massive Gravitons

Endowing the graviton with mass would have many aﬀects on diﬀerent physical processes. These
include changing the dispersion relation for gravitational waves such that they no longer travel at
the speed of light [149], increasing the energy loss due to gravitational waves of a given frequency[150], and causing the Newtonian gravitational potential to have a Yukawa form instead of being
1/r which can aﬀect the orbits of planets [149] and gravitational lensing observations [151]. The restof this chapter will focus on the first of these eﬀects expect the final section (section 12.7) where wecompare these bounds to the bounds from other methods.

12.2

Dispersion of Gravitational Waves

In reference [149], Will derived the propagation eﬀects for two GWs traveling over cosmologicaldistances, which we repeat below. Let us distinguish these GWs using the presence or absence ofa “� ” on the properties of the two waves. Let us also denote properties of either wave at the timeof emission with an e subscript and at the time of arrival or measurement with an a subscript.
The emitted energies of these GWs is then given as Ee and Ee� and the emission time diﬀerence as∆te ≡ t�e − te . For clarity, we also recover the G and c terms for formulae in this chapter.
Let us assume the GWs propagate in a flat Friedman-Robertson-Walker homogeneous andisotropic spacetime, with background metric given as����
ds2 = −c2 dt2 + a2 (t) dr2 + r2 dθa + sin2 θdφ2 ,(12.1)

where a(t) is the expansion factor of the universe.
Consider a graviton moving radially with mass m, 4-velocity uα = dxα /dτ , and 4-momentum

149pα = muα . Let it be emitted at a distance r = re and detected at r = 0. The energy-momentum
relation is given as
m2 c4 = −pα pβ gαβ = E 2 − a−2 p2r c2 ,(12.2)

where E = p0 c is the energy of the particle, m is the mass of the particle, and pr is the particles
radial momentum. Combining this with the ratio pr c/E = pr /p0 = dr/dt, after eliminating E we
find
dr
c
=−
dt
a

�

a2 m2g c2
1+
p2r

�−1/2

.

(12.3)

Since the graviton will be following a geodesic, its 4-velocity will satisfy the geodesic equation
uα ;µ uµ = 0. Using the property of the metric that g µν ;α = 0, we also know that uα;µ uµ = 0. Expanding the covariant derivative using a Christoﬀel symbol, this is then given by uα,µ uµ − uν uµ Γν αµ = 0.
Expanding the Christoﬀel symbol and using its symmetries, we can write this as uα,µ uµ −uµ uν gµν,α =0. Specializing this to ur and using the fact that only the time and radial components of the 4velocity uα are nonzero, and both g0µ and grµ are independent of r, we then find uα,µ uµ = 0, whichshows ur to be conserved. This then implies that pr is a conserved component of the 4-momentum.
��
We choose to evaluate it at the moment of emission of the GW, p2r c2 = a2 (te ) Ee2 − m2g c4 .

Using the relation from quantum mechanics between a wave’s frequency f and its energy E = hf ,where h is Planck’s constant, and the definition of the Compton wavelength of a particle λ ≡ h/m,
we find that mg c2 /Ee = c/(λg fe ). Let us assume that the lowest frequency we are interested inis such that Ee � mg c2 . Then the dispersion relation above can be expanded to second order inc/(λg fe ) as
c
dr
=−
dt
a

�

c2 a2
1− 2
2a (te )λ2g fe2

�,(12.4a)a(t)dt .

(12.4b)which after integration gives

re =

� ta
te

c
c3
dt − 2a(t)2a (te )λ2g fe2

� ta
te

Let us apply this to the two GW emitted with energies Ee and Ee� with a time diﬀerence of150∆te ≡ t�e − te . Application to the first wave yields
re =

� ta
te

c
c3
dt − 2a(t)2a (te )λ2g fe2

� ta

c
c3
dt − 2a(t)2a (te )λ2g fe�2

� ta +∆ta

a(t)dt ,(12.5a)tewhile application to the second yields

re =� ta +∆tate +∆tea(t)dt .

(12.5b)te +∆teCombining these two results gives us
� ta +∆ta
ta

c
dt −
a(t)

� te +∆te
te

c
dt
a(t)

=

�

� ta +∆ta1
a(t)dta2 (te + ∆te )fe�2 te +∆te
�
� ta
1
− 2
a(t)dt.
a (te )fe2 te
c3
2λ2g(12.6)

This can be simplified even further in the case where the characteristic time it takes a to change islarge compared to the emission time diﬀerence, characterized by ∆te � a/ȧ, and where the time
diﬀerence between te and ta is large to the form
�
�
��
Dc
1
1
∆ta = (1 + Z) ∆te + 2,− 2
2λg fe�2fe(12.7)

where Z ≡ a0 /a(te ) − 1 is the cosmological redshift, and D is a distance given by(1 + Z)c
D≡
a0

� ta

a(t)dt ,(12.8)

tewhere a0 ≡ a(ta ) is the present value of the expansion factor. D is diﬀerent from standard cosmological distances, such as the luminosity distance DL ≡ a0 (1 + Z)χe . For a matter-dominatedspatially-flat universe, D and DL are given by

D=

�
�
2c
(1 + Z) 1 − (1 + Z)−5/2 ,5H0(12.9a)

151
�
�
2c
DL =
(1 + Z) 1 − (1 + Z)−1/2 ,H0(12.9b)where H0 is the Hubble constant given by c/H0 � 42Gpc. Both cases reduce to the familiar Hubble
law cZ/H0 in the limit Z � 1. It should also be noted that the frequency upon arrival to thedetector will have undergone a redshift such that fa = fe /(1 + Z).

12.3

Dispersion of Compact Binary Coalescence Gravitational
Waveforms

This dispersion relation above, equation (12.4), tells us that GW of diﬀerent frequencies will travel
at diﬀerent speeds. This can change the morphology of GW signals emitted from a single sourcethat originally were emitted with a well-defined frequency evolution.
With this in mind, let us analyze what happens to the inspiral waveform we calculated in chapter
2, where the orbital phase evolution Φ(t) is given in equation (2.36b). Taking the time derivative of
the phase evolution 2Φ(t), we find the frequency evolution πf (t) = dΦ/dt to be

f (t) =

�GM
c3 π 2

�1/2



256G3 M 3 η

�

5c9
405GM
16c3 η + tISCO − t

3/8

�

,

(12.10)where M = m1 + m2 is the total mass of the system, η = m1 m2 /M 2 is the symmetric mass ratio,
and we have defined tISCO to be the time when the frequency is that corresponding to the ISCO
fISCO = c3 /(G63/2 πM ).
Let us consider equation (12.7) for two waves emitted at t = tISCO − ∆t and t� = tISCO and findthe distance D at which they will arrive at the same time (i.e., ∆ta = 0). Plugging these two times
into equation (12.10), using the frequencies in equation (12.7), and solving for D where ∆ta = 0, we
find
D=

π2
53/4 M�

256ηM 3�

2λ2g ∆t405M
16η + ∆t��3/4.

(12.11)− 216π 2 M 2Taking the limit of this as ∆t → 0, we find the distance at which the end of the waveform starts to152

Figure 12.1: Waveform Inversion
Plots of inspiral waveforms at D � Dinv (left) and D � Dinv (right). The upper plots show
spectrograms of the waveforms below. This makes it easier to see what frequency information
arrives at what time.
become inverted (Dinv ) to be
Dinv =

5λ2g.
16π 2 ηM

(12.12)

Diﬀerent waveforms are plotted in figure 12.1 showing how the morphology of the GW wave
changes as it propagates due to dispersion eﬀects of a massive graviton. The end of the waveform
starts becoming inverted at D ∼ Dinv . At D � Dinv , the information at the end of the waveform(i.e., the higher frequency portion of the waveform) can be seen to arrive before information thatwas emitted earlier (i.e., the lower frequency portion of the waveform).

153

12.4

Recovery of Dispersed Compact Binary Coalescence Gravitational Waveforms

In chapter 2 we found the inspiral portion of CBC waveforms as a function of time. Converted to
the form found in [149], these are given as

h(t) = A(t)e−iΦ(t) ,(12.13a)�tΦ(t) = Φ0 + 2π tc f (t)dt ,M
A(t) = D
L

�

5M
tc −t

�1/4

(12.13b)(12.13c).

When searching for these signals in GW detector data, as discussed in chapter 2, it is convenient
to convert these waveforms to the frequency domain using the SPA (equation (2.41)). Converting
the SPA phase (equation (2.45)) to the form found in [149], Ψ(f ) is given as

Ψ(f ) = 2π

� fcf(tc − t)df � + 2πf tc + φ0 ,(12.14)

where the subscript c denotes a value at the time of coalescence. These waveforms are used to searchfor signals in the GW detector data with the technique known as matched filtering (discussed inchapter 5) where an inner product is defined as(h1 |h2 ) = 4�

��

0

∞

h̃∗1 h̃2df
Sn (f )

�,(12.15)where �[· · · ] denotes the real part of [· · · ], h̃∗ denotes the complex conjugate of h̃, and Sn (f ) is thepower spectral density of the particular detector.
Since these are the waveforms that are used for the detection of signals, all of the f ’s above
are measured frequencies, which as we have stated above, are related to the emitted frequencies by
f = fe /(1 + Z). We are now in a position to see how the massive graviton propagation eﬀects change
a signal of the form given in equation (2.41). The term which is modified by these eﬀects is the154
integral term in the phase Ψ(f ) (equation (12.14)). Applying equation (12.7) to this term we findit changes to be� fc
f(tc − t)df �

=
=

� fec
fe
� fec
fe

(tec − te )dfe� +� fec(tec − te )dfe� −

D
D
D
+ 2
− 2 2 fe .
2
2λg feλg fec2λg fecfeDdf � −
2λ2g fe�2 e

� fec
fe

D
dfe�
2
2λ2g fec(12.16)

These additional terms can be seen to have diﬀerent eﬀects on the waveform. The last term has
an fe1 dependence, which is the same dependence as the time oﬀset term in equation (12.14) andcan be interpreted as a change in the arrival time. The second-to-last term has an fe0 dependence,
which is the same dependence as the phase oﬀset term in equation (12.14) and can be interpreted
as a change in the phase at coalescence. The other additional term (the second term) has an fe−1
dependence, which is the same as the 1.0PN phase term. This is the term that will change themorphology of the waveform. Its eﬀect on equation (2.45) is to add an additional −βf −1 term,
πD
where β ≡ λ2 (1+Z). To 1.5PN order, Ψ(f ) is then given as
g

�
�
�
3
3715 55−5/3−1
Ψ(f ) = 2πf tc + φ0 −(πM f )++ η (πM f )
128η
756
9
�
−2/3
−16π (πM f )− βf −1

12.5(12.17)

Using the Dispersion of Gravitational Waves to Boundthe Mass of the Graviton

In reference [149], Will used the above formalism with restricted inspiral waveforms to 1.5PN orderin the phase to predict bounds on the mass of the graviton from the detection of signals fromAdv. LIGO and LISA. Will’s analysis used the rms error of associated massive graviton phase termfrom parameter estimation techniques. Specifically, the rms error of a parameter θα is given by

∆θ ≡
α

��

2

(θα − �θα �)

�

=

√

Σαα ,

(12.18)

155
where Σ is the inverse of the Fischer Information matrix Λ given by

Λαβ ≡

�

∂h ∂h
|
∂θα ∂θβ�.

(12.19)For this analysis, we will compare the results of Will’s analysis [149] using 1.5PN waveforms toa similar analysis using new IMR waveforms. The IMR waveforms we use are frequency domain
phenomenological waveforms whose phenomenological parameters have been tuned using Numerical
Relativity waveforms with symmetric mass ratios 0.25 ≥ η ≥ 0.16 [152, 153]. These waveforms are
given byh̃(f ) ≡ Aeﬀ (f )eiΨeff (f ) ,(12.20)where the eﬀective amplitude Aeﬀ (f ) and phase Ψeﬀ are given by

Aeﬀ (f ) ≡ C






(f /f1 )−7/6









(f /f1 )−2/3

if f < f1
if f1 ≤ f < f2wL(f, f2 , σ) if f2 ≤ f < f3









0
if f ≥ f3 ,Ψeﬀ (f ) ≡ 2πf t0 + φ0 +6�ψk f (k−5)/3 ,(12.21a)(12.21b)

k=0

�
−7/6
M 5/6 f1
5η
where C = DL π2/3
24 is the amplitude parameter for optimally oriented and located sources,� �−2/3f2
1
σw = πσis a normalization such that Aeﬀ is continuous through f2 , L(f, f2 , σ) = 2π2f1
(f −f2 )2 +σ 2 /4
is a Lorentzian function centered on f2 with width σ, t0 is the arrival time, and φ0 is the phaseoﬀset. For these waveforms, f1 corresponds to the where the PN amplitude evolution ends (i.e.,
where the inspiral ends and the merger begins), f2 corresponds to where the merger ends and the
ringdown begins, and f3 is the frequency at which we cut oﬀ the waveform. The ringdown portionof the waveform is modelled as a Lorentzian which agrees with the quasi-normal mode ringing of a
perturbed black hole from black hole perturbation theory.

156
The phenomenological amplitude parameters αk = {f1 , f2 , f3 , σ} are given byαk =ak η 2 + bk η + ck 3/5
η,πM(12.22a)where ak , bk , and ck are given in Table 12.1. The phenomenological phase parameters ψk are givenbyψk =xk η 2 + yk η + zk −k/5η,(πM)(5−k)/3(12.22b)where xk , yk , and zk are given in Table 12.1. The phenomenological phase parameters do not matchthe PN terms because of the choice of the authors in references [152, 153]. They allowed these terms
to vary while trying to maximize the “faithfulness” and “eﬀectualness” of the phenomenological
waveforms to recovering hybrid PN-Numerical Relativity waveforms.
Table 12.1: Phenomenological CoeﬃcientsParameter
ak
bk
ck
f1
6.6389 ×10−1-1.0321 ×10−1
1.0979 ×10−1
f2
1.3278-2.0642 ×10−1
2.1957 ×10−1−1
σ1.1383-1.7700 ×10
4.6834 ×10−2−1
f3
1.7086-2.6592 ×10
2.8236 ×10−1
Parameter
xk
yk
zk
ψ0
-1.5829 ×10−1
8.7016 ×10−2-3.3382 ×10−2
ψ2
3.2967 ×101
-1.9000 ×101
2.13452
ψ3-3.0849 ×10
1.8211 ×102
-2.1727 ×101
32
ψ4
1.1525 ×10-7.1477 ×10
9.9692 ×101
32
ψ6
1.2057 ×10-8.4233 ×10
1.8046 ×102
Unitless coeﬃcients describing the amplitude and phase of the phenomenological waveforms fromreference [152] to be used in equations (12.22). There is no ψ1 term as the 0.5PN phasing term is0.

The amplitude C given above is the amplitude one would find for an optimally oriented andoptimally located source. Averaging over random angles for the two sky location angles, the polarization angle, and the inclination angle, there is an overall factor of 2/5 we should apply to this. For
this reason we redefine C =2M 5/6 f1
5DL π 2/3

−7/6

from here on.

Applying equation (12.16) to this waveform results in a measured phase evolution ofΨeﬀ (f ) ≡ 2πf t0 + φ0 +

6
�

k=0

ψk f (k−5)/3 − βf −1 ,(12.23)

157
where t0 is the measured arrival time which has absorbed the massive graviton term with the f 1
dependence and φ0 is the measured phase oﬀset which has absorbed the massive graviton term withthe f 0 dependence.
In producing this result, we have made an assumption which we must check. This is that, atthe order to which we the phase evolution of the binary, energy loss due to massive graviton eﬀects
are negligible. Since we use phenomenological terms up to the equivalent 3.0PN order and ignore
massive graviton eﬀects of order (r/λg )2 when computing the energy loss at the source, this implies−6
we require r2 λ−2� 1. Combining this with v 2 ∼ M/r and v = 2πrf , we find the frequencies
g v
−3/5

this holds for are f > λg

M −2/5 /(2π), which is f > 2.6 × 10−3 (M/M⊙ )−2/5 Hz assuming thecurrent bound on λg > 1012 km. For the systems we consider, this is always satisfied.
The parameters we use in the Fischer Matrix calculation are ln C, φ0 , t0 , ln M, ln η, and β. The
partial derivatives of h̃ with respect to these are given by∂ h̃(f )= h̃(f ) ,∂ ln C(12.24a)∂ h̃(f )= −ih̃(f ) ,∂φ0(12.24b)∂ h̃(f )= 2πih̃(f ) ,∂t0(12.24c)∂ h̃(f )∂ ln M∂ h̃(f )∂ ln η

158

� �


 7 ∂ ln f1 f −7/6
6 ∂ ln M f1





� �−2/3


ln f1
f

 23 ∂∂ ln
M f1= CeiΨeff

if f < f1
if f < f1


�
�


∂ ln L(f,f2 ,σ)∂ ln w

wL(f,
f
,
σ)
+
if f2 ≤ f < f32
∂ ln M∂ ln M








0
if f ≥ f3




−7/6

if f < f1(f /f1 )






−2/3

�
�
6
if f1 ≤ f < f2�∂ ln ψk (k−5)/3 (f /f1 )
iΨeff
+ iCe
ψk
f

∂ ln M
k=0
wL(f, f2 , σ) if f2 ≤ f < f3









0
if f ≥ f3 ,= CeiΨeff� �−7/6f
7 ∂ ln f1


6
∂
ln
η
f1





� �−2/3
 23 ∂∂lnlnfη1 ff1

if f < f1
if f < f1


�
�


∂ ln L(f,f2 ,σ)∂ ln w

wL(f,
f
,
σ)
+

2
∂ ln η
∂ ln η








0

+ iCeiΨeff(12.24d)(f /f1 )−7/6








�
(f /f1 )−2/3

if f2 ≤ f < f3
if f ≥ f3
if f < f1

6 �
if f1 ≤ f < f2�∂ ln ψk (k−5)/3ψkf∂ ln η

k=0
wL(f, f2 , σ) if f2 ≤ f < f3









0
if f ≥ f3 ,∂ h̃(f )= −if −1 h̃(f ) .
∂β(12.24e)(12.24f)

When taking the inner products in calculating Λ, we find it is helpful to introduce the definition

159
ofof an integral Iab (p, q) given by

Iab (p, q) ≡ 4

� bf −p/3df ,q
2
2
a ((f − f2 ) + σ /4) gα (f )(12.25)where gα (f ) is the frequency dependence of the detectors PSD. The PSD can then be written as
Sn (f ) = S0 gα (f ). In this notation, the SNR can be written asρ ≡ (h|h)1/2�
�1/2
Cw 2 σ 2 f3
7/3 f1
4/3 f2= √f I (7, 0) + f1 If1 (4, 0) +I (0, 2).
4π 2 f2
S0 1 0(12.26)In the calculation of the inverse of the Fischer Matrix, we find all of the terms to be proportional
to√√S0 /C. Let us then define ∆ as ∆β ≡ ∆1/2 S0 /C. In setting an upper limit on the mass of thegraviton, we view ∆β as the upper bound on β which results in a lower bound on λg given by

λg >

�

D
(1 + Z)DL

�1/2 �

π 2 CM
1/2

S0 ∆1/2

�1/2.

(12.27)In Table 12.2 we summarize the bounds on λg we obtain from the complete IMR waveforms andcompare them to the results from Will [149]. This table uses two diﬀerent noise models, one forsignals which Adv. LIGO will be sensitive to, and one for signals LISA will be sensitive to. The
noise model for Adv. LIGO is given by S0 = 3 × 10−48 and gα (f ) = ((f /fn )−4 + 2 + 2(f /fn )2 )/5
where fn = 70Hz. Since Adv. LIGO is sensitive down to f0 = 40Hz, this is where we start the
integration instead of at 0Hz. The noise model for LISA is given by S0 = 4.2 × 10−41 and gα (f ) =(101/2 (f /fn )−14/3 + 1 + (f /fn )2 /1000 + 313.5(f /fn )−(6.398+3.518 log10 x) ) where fn = 10−3 Hz. For
LISA, we start the integration at the lower frequency corresponding to a T = 1yr integration time,such that f0 = (5MT −1 /256)3/8 /(πM).
We see from table 12.2 and figure 12.2 that for some signals, we can obtain better bounds
from the PN inspiral-only waveform than from the complete IMR waveform. If we repeat theabove analysis for the phenomenological waveform keeping only the inspiral portion with the f −7/6

160
Table 12.2: Massive Graviton Bounds
m1
m2
Distance
InspiralComplete Waveform
(M⊙ )
(M⊙ )
Detector(Mpc)
Bound (km)Bound (km)1210
10Adv. LIGO
1500
6.0 × 10
4.8 × 1012
12
50
50Adv. LIGO
3000
2.8 × 10
4.8 × 1013
4
4
15
10
10
LISA
3000
7.4 × 10
4.3 × 1015
5
5
16
10
10
LISA
3000
2.3 × 10
1.1 × 1016
6
6
16
10
10
LISA
3000
5.4 × 10
2.1 × 1017
7
7
16
10
10
LISA
3000
6.9 × 10
5.1 × 1017
8
8
16
10
10
LISA
3000
5.8 × 10
9.7 × 1017
9
9
16
10
10
LISA
3000
2.6 × 10
1.1 × 1018
The theoretical bounds on λg that could be placed based on a GW detection from diﬀerent mass
CBC signals. The column “Inspiral Waveform” refers to the bounds one could get from justlooking at the inspiral portion of the waveform up to the ISCO frequency. Where available, these
agree quantitatively with Will [149]. The column “Full Waveform” refers to the bounds one could
get from looking at the complete IMR waveform which has been characterized by Ajith [152] andAjith et al. [153].

amplitude frequency dependence and setting the cutoﬀ frequency to be fISCO , we find that the PN
waveform always beats the phenomenological waveform. This is due to the diﬀerences between thephenomenological phasing parameters and the PN phasing parameters Ajith et al. noted in [153].
The above observation implies the parameters of the PN waveform are less degenerate than theparameters of the phenomenological waveform.

161

Figure 12.2: Massive Graviton Bounds
Plots of theoretical bounds on λg that could be placed based on a GW detection from Adv. LIGO
mass systems (left) and LISA mass systems (right). For each mass of the Adv. LIGO systems, thebound was calculated using lesser of the distance corresponding to an SNR of 10 for the 1.5PN
waveform and a distance of 3 Gpc and a minimum frequency of 10 Hz. For each mass of the LISA
systems, the bound was calculated using a distance of 3 Gpc and a bandwidth corresponding to 1
year of integration time. The 1.5PN curves correspond to the bounds one could get from justlooking at the 1.5PN waveform up to the ISCO frequency, as was done in Will [149]. The IMR
curves correspond to the bounds one could get from looking at the complete IMR waveform which
has been characterized by Ajith [152] and Ajith et al. [153].

162

12.6

Beyond the Fischer Matrix Approach

The approach that has been taken by Will [149] and here to constrain the mass of the graviton
is the Fischer information matrix approach. This approximation is only valid in the large SNR
limit and does not include known issues involved in real data analysis pipelines, which search for
inspiral signals, such as recovery of a signal with incorrect parameters and separating the amplitude
correlations between distance, mass, and sky location. A full simulation is required to accurately
determine the achievable bounds on the mass of the graviton with one, or even several, observed
signals.

12.7

Other Bounds on the Graviton’s Mass

Let us compare these results to the results available from other physical eﬀects of a massive graviton.
The first eﬀect we shall look at is changes in planetary motion due to a Yukawa gravitational
potential. If this were the case, Kepler’s third law would be violated since the gravitational force
would no longer follow an inverse-square law. This leads to a bound on the mass of the graviton in
1−a2

the form of λg > ( 6ηpp )1/2 where ap is the semimajor axis of planet p, ηp is planet p’s bound on
the η parameter given by Talmadge et al. [154]. Both λg and ap are given in astronomical units.
The most stringent bound of this form comes from the orbit of Mars which limits the mass of thegraviton such that λg > 2.8 × 1012 km [149].
The next bound on the mass of the graviton is from binary pulsar observations. If the graviton
had mass, the orbits of binary pulsars would decay at a slightly faster rate than predicted by General
Relativity, due to additional energy loss from the leading order massive graviton terms in the powerradiated. Combining the observations of PSR B1913+16 and PSR B1534+12, the 90% confidencebound on the mass of the graviton is λg > 1.6 × 1010 km [150].
The final bound on the mass of the graviton we shall look at can be obtained by studying
gravitational lensing data. Assuming a specific distribution of dark matter, lensing data can bestudied looking for an eﬀect which would change the distribution of the variance of the power

163
spectrum γ 2 (θ) smoothed over a filter of radius θ. When comparing the predictions of a 1/r potential
versus the Yukawa potential, Choudury et al. [151] find the bound on the mass of the graviton to be
λg > 3 × 1021 km [151]. However, this upper limit is less robust than the others above as one must
assume a specific distribution of dark matter before calculating the bound.

164

Chapter 13

Toward Understanding Black Hole
Merger Dynamics
The stunning breakthroughs in Numerical Relativity over the past few years, starting with thoseby Pretorius [155] in spring 2005, have provided greater insight into the realm of strongly curved
spacetime and dynamical gravity. Non-spinning binary black hole merger simulations are quickly
becoming the norm, and numerical relativists have started looking at more interesting situations,
including the orbit and merger of spinning binary black holes. Surprising results from this sort ofsituation are exemplified in the extreme-kick configuration in which two identical, spinning black
holes are initially in a (quasi-)circular orbit, with oppositely directed spins lying in the orbital plane(figure 13.1).
As Campanelli, Lousto, Zlochower, and Merritt [156, 157] (CLZM) discovered and Healy et
al. [158] helped flesh out, of all equal mass, quasi-circular initial configurations, this one has the

m
z

SA

x

y

SB

m
Figure 13.1: Extreme-kick Configuration
Extreme-kick configuration for a black-hole binary: Identical holes, A and B with masses
m = M/2 move in a circular orbit with their spin angular momenta SA and SB antialigned andlying in the orbital plane.

165
largest kick speed for the final black hole. Not only that, it also exhibits intriguing orbital motions.
During the inspiral phase, as the holes circle each other, they bob up and down (in the z directionof figure 13.1), sinusoidally and synchronously. After merger the combined hole gets kicked up or
down with a final speed that depends on the orbital phase at merger (relative to the spin directions).
This bobbing-then-kick, as deduced by CLZM from numerical simulations, is graphed quantitatively
in figure 13.2.
Qualitatively, Pretorius [159] has oﬀered a lovely physical explanation for the holes’ bobbing(figure 13.2) in this configuration: In figure 13.3, taken from his paper, we see snapshots of the holes
at four phases in their orbital motion. In each snapshot, each hole’s spin drags space into motion(drags inertial frames) in the direction depicted by gray, semicircular arrows. In phase B, hole 1
drags space and thence hole 2 into the sheet of paper (or computer screen); and hole 2 drags space
and thence hole 1 also inward. In phase D each hole drags the other outward. This picture agrees inphasing and semiquantitatively in amplitude with the bobbing observed in the simulations (figure13.2).
However, momentum conservation dictates that, when the holes are moving upward together
with momentum pzA + pzB , there must be some equal and opposite downward momentum in theirgravitational field (in the curved spacetime surrounding them); and when the holes are moving downward, there must be an equal and opposite upward field momentum. Our aim is to quantitatively
answer the following questions: How is this field momentum distributed? What are the details of the
momentum flow between field and holes? And to what extent are other momentum-flow processes
responsible for the motion shown in figure 13.2?
The outline of this chapter is as follows: section 13.1 gives a brief overview of the Landau-Lifshitz
formalism, which can be used to calculate the energy and momentum content of the gravitational
field, section 13.2 sets up the calculation of momentum consveration by splitting the total momentum
into its respective pieces, section 13.3 calculates each of the pieces of the momentum both for thegeneral binary and in the extreme-kick configuration, and section 13.4 combines these results toshow that the momentum of the objects balanced by the field momentum.

166

0.5
SP2
SPA
SPB
SPC
SPD
SPE

0.3

z/M

0.1

−0.1

−0.3

−0.5

0

100

200

300

t/M
Figure 13.2: Bobbing Motion in Extreme-kick Configuration
Bobbing and kick of binary black holes in the extreme-kick configuration of figure 13.1, as
simulated by Campanelli, Lousto, Zlochower, and Merritt (CLZM) [157]. Plotted vertically (as afunction of time horizontally) is the identical height z of the two black holes, and thentransitioning through merger (presumably at t/M ∼ 170), the height of the merged hole, above theinitial orbital plane. This height versus time is shown for six diﬀerent initial configurations, eachleading to a diﬀerent orbital phase at merger. In all six configurations, the initial holes’ spins arehalf the maximum allowed, a/m = 0.5. The height z and time t are those of the “punctures” that
represent the holes’ centers in the CLZM computations, as defined in their computational
coordinate system, which becomes Lorentz at large radii. These z and t are measured in units ofthe system’s total mass M � 2m.

Figure 13.3: Pretorius’ Explanation
Pretorius’ physical explanation for the holes’ bobbing in the extreme-kick configuration, in which
the eﬀect is attributed qualitatively to frame-dragging.

167

13.1

The Landau-Lifshitz Formalism in Brief

Here we give a detailed analysis of momentum flow in generic compact binary systems. We begin
in this section with a brief review of the Landau-Lifshitz (LL) formulation of general relativity as anonlinear field theory in flat spacetime [160].
The Landau-Lifshitz formulation [160] starts by choosing an arbitrary coordinate system thatis asymptotically Lorentz, to this an auxiliary flat spacetime is added by asserting that the chosen(“preferred”) coordinates are globally Lorentz in the auxiliary spacetime (so in them the auxiliary
metric has components diag(−1, 1, 1, 1)). In this formulation, gravity is described by the physicalmetric density
gµν =

√

−gg µν ,(13.1)

where g is the determinant of the covariant components of the physical metric, and g µν are thecontravariant components of the physical metric.
The EFEs are then reformulated as a nonlinear field theory in the space of the flat, auxiliary
metric. Making use of the superpotential

H µανβ ≡ gµν gαβ − gµα gνβ ,(13.2)

H µανβ ,αβ = 16πτ µν ,(13.3)

the EFEs are given by

where τ µν = (−g)(T µν + tµνLL ) is the total eﬀective stress-energy tensor, indices after the comma
denote partial derivatives (covariant derivatives with respect to the flat auxiliary metric), and theLandau-Lifshitz pseudotensor tµν
LL (actually a real tensor in the auxiliary flat spacetime) is given
by equation (100.7) of LL [160] or equivalently equation (20.22) of MTW [1]. By virtue of thesymmetries of the superpotential (which are the same as those of the Riemann tensor), the field

168
equations in the form (13.3) imply the diﬀerential conservation law for 4-momentumτ µν ,ν = 0 ,(13.4)

which is equivalent to T µν ;ν = 0 (where the semicolon denotes a covariant derivative with respect
to the physical metric).
It is shown in LL and in MTW that the total 4-momentum of any isolated system (as measured
gravitationally in the asymptotically flat region far from the system) is

pµtot =

116π

�H µα0j ,α dΣj ,(13.5)

S

where dΣj is the surface-area element defined using the flat auxiliary metric, and the integral is overan arbitrarily large closed surface S surrounding the system. Diﬀerentiating this, we finddpµtot
1
=
dt
16π

�

H µα0j ,α0 dΣj .

(13.6)S

Let us look at H µα0j ,α0 = H µανj ,αν − H µαkj ,αk . The first term is −16πτ µj by virtue of the field
equations (13.3) and the antisymmetry of the superpotential on its last two indices (13.2). Thatsame antisymmetry on the second term −H µαkj ,αk permits us to write it as the curl of a 3-vector
field, whose surface integral vanishes by virtue of Stokes’ theorem. Combining these two results wefind the total 4-momentum satisfies the standard conservation law
dpµtot=−dt�τ µj dΣj .

(13.7)S

This proof has also been given in LL and MTW where it relied on an assumption that the interior
of S be simply connected, i.e., that it not contain any black holes.

169

S
E
∂B

∂A
A

B

Figure 13.4: Division of SpaceThe regions of space around and inside a compact binary system.

13.2

Momentum Conservation for a Fully Nonlinear Compact Binary

We now apply this LL formalism to a binary system made of black holes and/or neutron stars; seefigure 13.4. We denote the binary’s two bodies by the letters A and B, and the regions of spaceinside them by these same letters, and their surfaces by ∂A and ∂B. For a black hole, ∂A could be
the hole’s absolute event horizon or its apparent horizon, whichever one wishes. For a neutron star,∂A will be the star’s physical surface. We denote by E the region outside both bodies, but inside
the arbitrarily large surface S where the system’s total momentum is computed.
By applying Gauss’s theorem to equation (13.5) for the binary’s total 4-momentum and using theEFE (13.3), we obtain an expression for the binary’s total 4-momentum as a sum over contributions
from each of the bodies and from the gravitational field in the region E outside them:pµtot = pµA + pµB + pµfield .

(13.8a)

170
Here
pµA ≡

1
16π

�

H µα0j ,α dΣj(13.8b)∂A

is the 4-momentum of body A and similarly for body B, and�

pµfield ≡τ 0µ d3 x(13.8c)Eis the gravitational field’s 4-momentum in the surrounding space.
If either of the bodies has a simply connected interior (is a star rather than a black hole), then we
can use Gauss’s theorem and the EFEs (13.3) to convert the surface integral (13.8b) for the body’s
4-momentum into a volume integral over the body’s interior:

pµA =

�

τ 0µ d3 x .

(13.8d)A

By an obvious extension of the argument we used to derive equation (13.7) for the rate of changeof the binary’s total 4-momentum, we can deduce from equation (13.8b) the corresponding equation
for the rate of change of the 4-momentum of body A:
dpµA
=−
dt

�

∂A

k
(τ µk − τ µ0 vA)dΣk .

(13.9)

Here the second term arises from the motion of the boundary of body A with coordinate velocitykvA= dxkA cm /dt. Equation (13.9) describes the flow of field 4-momentum into and out of body A.

We shall use equations (13.8), (13.7), and (13.9), specialized to linear momentum (index µ madespatial) as foundations for our study of momentum flow in compact binaries.
The actual values of the body and field 4-momenta, computed in the above ways, will depend on the arbitrary coordinate system that we chose, in which to make the auxiliary metric
be diag(−1, 1, 1, 1) and in which to perform the above computations. This is a “gauge dependence,”
which we will fix by our choice of coordinates. In the remainder of this chapter we shall choose

171
Harmonic coordinates, so the gravitational field satisfies the Harmonic gauge condition

gαβ ,β = 0 ,(13.10)and we shall specialize the above equations to the 1.5 post-Newtonian approximation and use themto study momentum flow during the inspiral phase of generic compact binaries.

13.3

Post-Newtonian Momentum Flow in Generic Compact
Binaries

13.3.1

Field Momentum Outside the Bodies

In Harmonic gauge at leading post-Newtonian order, the Landau-Lifshitz formalism gives for thedensity of field momentum
τ 0j ej = −

g×H
3
+
U̇N g
4π
4π

(13.11)(equation (4.1a) of [161]). Here, to the accuracy we need, g is the Newtonian gravitational acceleration field (the gravitoelectric field ), H is the gravitational analog of the magnetic field (thegravitomagnetic field ), UN is the Newtonian potential and the dot denotes diﬀerentiation with respect to time, and ej is the jth basis vector of the flat-spacetime field theory that we are using.
This field momentum can be split into several pieces, only one of which will flow back and forth
between the field and the bobbing holes. Each of the components in equation (13.11) (i.e., g, H,
and UN ) have contributions from body A and body B, where

gA
UN,A

m
= − 2 nA ,
rA
m
=
,
rA

(13.12a)(13.12b)(equations (2.5) and (6.1) of [161]). In addition, H can be split into terms arising from either

172
velocity (H velo ) or spin (H spin ) contributions of a particular body:

velo
HA
spin
HA

4mA (nA × vA ),2
rA
(3nA · SA )nA − SA
= −2
,
3
rA
=

(13.12c)(13.12d)(equations (2.5) and (6.1) of [161]), with nA the unit radial vector pointing from the center of body
A to the field point, rA the distance from the center of body A to the field point, vA the vectorial
velocity of body A, and SA the vectorial angular momentum of body A. The fields for body B are
the same as equations (13.12), but with each subscript A replaced by a B.
The full field momentum density is written most concisely as

0j0j
τ 0j = τspin+ τvelo,(13.13a)0j
0jwhere τspinand τvelo
are the terms that depend on the spins and the velocities, respectively. These

terms are given by

0j
τspin
ej

=

mB
3 r 2 [3(SA · nA )(nA × nB ) − (SA × nB )]
2πrA
B
1 mA
−
5 (SA × nA ) + (A ↔ B),
2π rA(13.13b)

and

j0
τvelo
ej

=

�
mA
mB [4(nB · vA )nA − 4(nA · nB )vA ]
2
2
4πrA
rB
�
3mB (nA · vA )nBmA [(nA · vA )nA − 4vA ]
−
+
2
2
rB
rA
+(A ↔ B),

(13.13c)where (A ↔ B) means the same expression with labels A and B interchanged.
We are only interested in that portion of the field momentum that is induced by the holes’ spins,
since this is the portion that must flow back and forth between the field and the bobbing holes in173
order to conserve total momentum. As can be seen from equations (13.12), this portion arises fromone hole’s gravitoelectric field g coupling to the spin-induced part of the other hole’s gravitomagnetic
field H spinδτ 0j ej = −

spin
spin
gA × HB
gB × HA−.
4π
4π(13.14)

Combining equations (13.14) and (13.12) we obtain for the binary’s density of field momentum(that portion which must flow during bobbing)δτ 0j ej =mB
3 r 2 [3(SA · nA )(nA × nB ) − (SA × nB )] + (A ↔ B) .
2πrA
B

(13.15)In order to integrate this over the region E, we find it convenient to rewrite the bodies’ gravitoelectric (13.12a) and gravitomagnetic fields (13.12d) as

j
gK
= mK

�

1
rK

�

j
i
HK
= −2SK

,
,j

�

1
rK

�

,

(13.16),ij

where K is A or B and where, as before, rK is the (flat-space) distance of the field point from thecenter of mass of body K. Inserting equations (13.16) into expression (13.14) and manipulating thederivatives, we obtain the following expression for the field momentum density:

δτ

0j

�
� � � � �111q
= − �jpl SA mB+ (A ↔ B) .
2πrA ,q rB ,l(13.17),p

Notice that this expression for the momentum density is the curl of a vector field; or, equally well,it can be viewed as the divergence of a tensor field.
The total spin-induced, flowing field momentum is the integral of expression (13.17) over theexterior region E (cf., figure 13.4). Using Gauss’s law, that volume integral can be converted into
the following integral over the boundary of E
1q
δpjfield = − �jpl SA
mB
2π

�

∂E

�

1
rA

� �
,q

1
rB

�

,l

dΣp + (A ↔ B) .

(13.18)

174
The boundary of E has three components: the surface S far from the binary on which we computethe binary’s total momentum, and the surfaces ∂A and ∂B of bodies A and B. The integral over S
vanishes because the integrand is ∝ 1/r4 and the surface area is ∝ r2 and S is arbitrarily far fromthe binary, r → ∞. When integrating over the bodies’ surfaces, we shall flip the direction of thevectorial surface element so it points out of the bodies (into E), thereby picking up a minus sign andbringing equation (13.18) into the form

δpjfield=�� � � � �
1
1
1q�jpl SA mB
dΣp
2πrA ,q rB ,l∂A�� � � � �
11+dΣp + (A ↔ B) .
rA ,q rB ,l∂B(13.19)We presume (as is required by the PN approximation) that the bodies’ separation is large compared22
to their radii. Then on ∂A, we can write (1/rA ),q = −nqA /rA
and (1/rB ),l = nlAB /rAB
, where nA is

the unit vector pointing away from the center of mass of body A, nlAB is the unit vector pointing fromthe center of mass of body B toward the center of mass of body A, and rAB is the (flat-spacetime)
distance between the two bodies’ centers of mass. The first integral in equation (13.19) then becomes2
nAB /rB�∂A

2
nqA /rAdΣp . For simplicity we take the surface of integration to be a sphere immediately

above the physical surface of body A and ignore the tiny contribution from the region between
2 pthat sphere and the physical surface. On this sphere, we write dΣp = rAnA dΩA , where dΩA is thesolid angle element, and we then carry out the angular integral using the relation

�∂A

nqA npA dΩA =2(4π/3)δqp . Thereby we obtain for the first integral in (13.19) (4π/3)δqp nAB /rB
independent of theradius rA of the sphere of integration. (If the body is not spherical, the contribution from thetiny volume between our spherical integration surface and the physical surface will be negligible.)
Evaluating the second integral in equation (13.19) in the same way, and carrying out straightforward
manipulations, we obtain for the external field momentum

δpfield =22 (mB SA − mA SB ) × nAB .
3rAB(13.20)175
13.3.1.1

Field Momentum Specialized to the Extreme-Kick Configuration

For the extreme-kick configuration, which has mA = mB = m and SB = −SA , the field momentum(13.20) becomesδpfield =4 m
2 SA × nAB .
3 rAB(13.21)Figure 13.5 shows the z-component (perpendicular to the orbital plane) of the field-momentum
density δτ 0z , as measured in the orbital plane at four diﬀerent moments in the binary’s orbital
evolution. Only that part of the momentum that flows during bobbing (equation (13.15)) is pictured.
Red depicts momentum density flowing out of the paper (+z direction), and blue, into the paper.
The yellow arrows show the holes’ vectorial spins S, and the arrowed circle is the binary’s orbital
trajectory. In the top-left and bottom-right frames, the black holes are momentarily stationary
at the top and bottom of their bobbing. Nevertheless, the momentum density has a nontrivial
distribution. In the top-right and bottom-left frames, the black holes are moving downward andupward, respectively, with maximum speed. In both cases, the field-momentum density between
the two holes flows in the same direction as the bobbing, whereas the momentum surrounding thebinary is in the opposite direction and larger. This leads to net momentum conservation for thebinary, as discussed in section 13.4.
It is worth noting that the four figures, going counterclockwise from the top-left, are taken aquarter period apart in orbital phase. The first and third diﬀer by half an orbital period (as do thesecond and fourth); and, consequently, the momentum patterns of each pair are identical, but signs
are reversed (red exchanged with blue, as dictated by the symmetry of the configuration). Thisfeature is responsible for the sinusoidal bobbing.

13.3.2Centers of Mass and Equation of Motion for the Binary’s Compact
BodiesRestrict attention, temporarily, to a body that is a star rather than a black hole, and temporarily
omit the subscript K that identifies which body. Then, following the standard procedure in special

176

Figure 13.5: Field Momentum Density
The four pictures show the z-component of field-momentum density δτ 0z in the orbital plane atfour diﬀerent times, a quarter orbit apart. Red represents positive momentum density (coming outof the paper), and blue, negative (going into the paper). Only the piece of momentum density δτ 0z
that flows during bobbing (equation (13.15)) is depicted. The yellow arrows are the black holes’
vectorial spins; the large, black arrowed circle shows the orbital path of the two holes. In the top
left picture, one sees the density of momentum when the black holes are at the top of their bob(maximum z) and momentarily stationary. The gravitational-field momentum is zero, but the
momentum density itself shows rich structure. A quarter orbit later, in the top right, the holes are
moving downward (into the paper) at top speed. The momentum between the black holes (blueregion) flows into the paper with them, while surrounding momentum (red region) flows out of thepaper (+z-direction). A half orbit after the first picture, in the lower left, the holes are
momentarily at rest at the bottom of their bob (minimum z), the net field momentum is zero, andthe momentum distribution is opposite that in the first picture (as one would expect duringsinusoidal bobbing). Similarly, three quarters of the way through the orbit, in the lower right, theholes have reached their maximum upward speed, and the momentum distribution is identical to
the second picture, but with the opposite sign.

177
relativity (e.g., Box 5.6 of MTW [1]), we define the star’s center-of-mass world line to be that set of
events xµcm satisfying the covariant field-theory-in-flat-spacetime relationshipS αβ pβ = 0 .

Here pβ =�(13.22)τ 0β d3 x is the body’s 4-momentum andS αβ ≡

�β0[(xα − xα− (xβ − xβcm )τ α0 ]d3 xcm )τ

(13.23)is the body’s tensorial angular momentum. Here the integrals extend over the star’s interior, andbecause the star’s momentum is changing, we take the time component of xµcm to be the same as
the time at which the integral is performed, x0cm = x0 . (If the momentum were not changing, thisrestriction would be unnecessary; cf. Box 5.6 of MTW.)
In a reference frame where the body moves with ordinary velocity v j = pj /p0 , equation (13.22)
says S i0 = S ij vj . We wish to rewrite this in a more illuminating form, accurate to first order in thevelocity v. At that accuracy, we can evaluate S ij in the body’s rest frame, obtaining S ij = �ijk Sk
where Sk is the body’s spin angular momentum

Sk =

�

�klm (xl − xlcm )τ m0 d3 x .

(13.24)Using definition (13.23) of S i0 with x0cm = x0 , our definition (13.22) of the center of mass then takes
the concrete form
mxcm =
Here on the left side we have replaced p0 =��xτ 00 d3 x − v × S .

(13.25)

τ 00 d3 x by its value in the body’s rest frame, which is

the mass m, since the two diﬀer by amounts quadratic in v.
Notice that, when computed in the body’s rest frame so v = 0, the center of mass is mxcm =�xτ 00 d3 x, but when computed in any frame moving slowly with respect to the rest frame, thisexpression must be corrected by the term −v × S. This is called the physical spin supplementary

178
condition [162].
In our Harmonic coordinate system and at the 1.5PN order of our analysis, the dominant, timetime component of the EFEs (13.3) reduces to η µν g00 ,µν = 16πτ 00 . The type of analysis carriedout in section 19.1 of MTW [1] then reveals that in the star’s rest frame, the monopolar part of its
g00 is centered on the location xcm ; or, equivalently, when one expands the star’s g00 around xcm inits own rest frame, there is no dipolar 1/r2 term (no mass dipole moment). This well-known result(e.g., [163, 164]) can be used as an alternative definition of xcm —a definition that works for black
holes as well as for stars.
Using this monopolar-field-centered definition of xcm , Thorne and Hartle [164] have employedmatched asymptotic expansions (valid for black holes) to derive the equations of motion for a systemof compact bodies (e.g., a compact binary) (their equations (4.10) and (4.11)). For a compact binary,
the spin-induced contributions to these equations of motion at 1.5PN order are (equation (4.11c) of
Thorne and Hartle)

mA

dδvA
dt

=

mA
[6nAB (SB × nAB · vAB ) + 4SB × vAB
3
rAB−6(SB × nAB )(vAB · nAB )]mB+ 3 [6nAB (SA × nAB · vAB ) + 3SA × vAB
rAB−3(SA × nAB )(vAB · nAB )] .

(13.26)HerevA ≡dxcm A,dt

δvA =dδxcm A,dt

vAB = vA − vB

(13.27)are the velocity of (the center of mass of) body A, the spin-induced perturbation of that velocity,and the relative velocity of bodies A and B. The first two lines of equation (13.26) are due to frame
dragging by the other body (body B); the last two lines are a force due to the coupling of body A’s
spin to B’s spacetime curvature.

179
13.3.2.1

Equations of Motion Specialized to the Extreme-Kick Configuration

As in the previous section, we now specialize the above discussion to the extreme-kick configuration.
In this configuration we are trying to explain the bobbing motion of the bodies as they orbit eachother (i.e., the motion of the bodies perpendicular to the standard orbital plane). As noted above forthe equal mass case, the symmetries of this configuration are such that vA = −vB , SA = −SB , andvAB · nAB = 0. Since we are only interested in the portion of the equations of motion perpendicular
to the orbital plane, the necessary portions of the equations of motion we shall look at are:� 2
�d δxcm A
dt2
� 2
�FD
d δxcm A
dt2
SC

=
=

1
3 4SB × vAB ,rAB
1

3 3SA × vAB ,rAB(13.28a)(13.28b)where the subscript FD refers to the frame dragging piece, while the subscript SC refers to thespin-curvature coupling piece. The sum of these two eﬀects gives
� 2
�
d δxA
2
= − 3 SA × vA ,2
dt
rAB
spin eﬀects(13.29)where we have used the fact that SA = −SB and vAB = vA − vB = 2vA . We get the sameacceleration of hole B by replacing all subscript A’s by subscript B’s because SB = −SA , vB = −vA .
We can easily integrate this equation in time by noting that the spin precesses much more slowly
than the orbital motion so SA can be approximated as constant, and noting that vA rotates withangular velocity Ω =

�

3
M/rAB
=

�

3 , where M is the total mass and m is the mass of each
2m/rAB

hole. The result, after one integration, can be written asmmδvA = mδvB = − 2 SA × nAB .
rAB
These pieces of the total momentum we will call the kinetic momentum.

(13.30)180

13.3.3The Momenta of the Binary’s Bodies

Let us now turn our attention to the total momentum contained by the bodies. This is done fortwo cases, first for stars where we can integrate over the interior, and then for black holes where wemust convert the volume integrals to surface integrals. We shall find that both methods yield thesame result. We initially omit the star’s label A or B for ease of notation.

13.3.3.1

The Momentum of a Star

For a star we can derive an expression for the momentum pj =�τ 0j d3 x (with the integral overthe star’s interior) in terms of the star’s velocity v j = dxjcm /dt by diﬀerentiating the center-of-massequation (13.25) with respect to time. To allow for the possibility that the mass might change withtime, we set m =�τ 00 d3 x before doing the diﬀerentiation; i.e., we diﬀerentiate

xcm

�

A

τ d x=00 3

�A(13.31)xτ 00 d3 x − v × S .

�Using τ 00 ,0 = −τ 0k ,k and Gauss’s theorem, we bring the left side into the form v A τ00 d3 x −�xcm ∂A (τ 0j − τ 00 v j )dΣj . The last term arises from the motion of the surface of the star throughspace with velocity v. Manipulating the time derivative of the integral on the right side of equation(13.31) in this same way, we bring it into the form
�

∂A

�

A

τ 0j d3 x −�∂A

xj (τ 0k − τ 00 v k )dΣk = pj −xj (τ 0k − τ 00 v k )dΣk , where pj is the star’s momentum. Inserting these expressions for the left

side and the right-side integral into equation (13.31), noting that the star’s spin angular momentum
evolves (due to precession) far more slowly than its velocity, denoting the time derivative of its
velocity by dv/dt = a (acceleration), solving for p, and restoring subscript As, we obtain

pA = mA vA +

�

∂A

k
(x − xcm A )(τ 0k − τ 00 vA)dΣk + aA × SA .

(13.32)

Although we have derived this equation for a star, it must be true also for a black hole. The
reason is that all the quantities that appear in it are definable without any need for integrating overthe body’s interior, and all are expressible in terms of the binary’s masses and spins and its bodies’

181
vectorial separation, in manners that are insensitive to whether the bodies are stars or holes. To
illustrate this statement, in section 13.3.3.2 we deduce (13.32) for a black hole, restricting ourselvesto spin-induced portion of the momentum that is being exchanged with the field, δpA .
It is this δpA that interests us. Because the spin has no influence on τ 00 at the relevant order(which is δτ 0k ∼ gH and δτ 00 ∼ g 2 where g and H are the gravitoelectric and gravitomagnetic
fields), equation (13.32) implies thatδpA = mA δvA +�

∂A(x − xcm A )δτ 0k dΣk + aA × SA .

(13.33)

The acceleration aA of body A is, at the order needed, just the gravitoelectric field of body B at2
the location of A, aA = −(mB /rAB
)nAB . Performing the surface integral on a sphere just above
2 k
the body’s physical surface we can write x − xcm A = nA rA and dΣk = rAnA dΩA . Inserting these

into equation (13.33), we obtain

δpA = mA δvA +� �� �
kinetic �

�

term

mB
3
rA
δτ 0k nA nkA dΩA + 2 SA × nAB .
r∂A��� � AB ���surface term

(13.34)

SSC term

Here “SSC term” refers to the “spin supplementary condition” required to get the correct, physicalcenter of mass; see text following equation (13.25). In the surface term, the field momentum density
δτ 0k is given by equation (13.17). The second term (A ↔ B) is smaller than the first by M/rAB
3and thus is negligible. Inserting the first term into the integral, using (1/rA ),qp = (3nqA npA − δqp )/rA2
and (1/rB ),l = −nlB /rB
, and

δpA

�

njA nlA dΩA = 4π
3 δjl , we bring equation (13.34) into the form

2 mB
mB
= mA δvA − 2 SA × nAB + 2 SA × nAB� �� � 3 rAB
rAB
��
��
��
�
kinetic �term

surface term

1 mB
= mA δvA +2 SA × nAB .
3 rAB

SSC term(13.35)

182
13.3.3.2

The Momentum of a Black HoleIn section 13.3.3.1 we derived expression (13.32) for the momentum of a body in a binary, assuming
the body is a star so we could do volume integrals, and we then asserted that this expression is also
valid for black holes. The spin-induced portion of this expression that gets exchanged with the fieldas the body moves is given by equation (13.33), which reduces to (13.35). In this section we shallsketch a derivation of equation (13.35) directly from the surface-integral definition (13.8b) of a blackhole’s momentum,
δpjA =

116π

�δH jα0k ,α dΣk .

(13.36)∂A

To evaluate this surface integral up to desired 1.5PN-order accuracy turns out to require some 2.5PN
fields. Qualitatively, this can be anticipated because the superpotential we use in the surface integral
is sourced by the spin-orbit piece of field momentum, and therefore necessarily a nonleading PN term.
One can see this more clearly by expanding δH jα0k ,α in terms of the metric density and using thesymmetries of the superpotential H (which are the same as the Riemann tensor). In general, themomentum is given by
δpjA = −

1
16π

�

∂A

(gjk gα0 − gjα g0k ),α dΣk .

(13.37)In Harmonic gauge, however, gαβ ,β = 0, and the spatial metric is flat until 2PN order while thetime-space components are of 1.5PN order. As a result, the terms at lowest and next-to-lowest PN
order are contained within two terms,

δpjA =

116π

�(gj0 ,k + gjk ,0 )dΣk .

(13.38)∂A

In this expression, the momentum arises from linear terms involving the metric density, instead ofquadratic ones. As a result, one must keep pieces of the metric perturbation that are of higher PN
accuracy. (Note: if we were to evaluate the time derivative of δpA using the surface integral (13.9),
we would not face such a delicacy; the integrand there is quadratic and requires only 1.5PN fields
for its evaluation.)

183
To find the momentum in terms of the standard post-Newtonian parameters, we resort to a
standard way that the metric perturbations are written in recent post-Newtonian literature, e.g., by
Blanchet, Faye, and Ponsot [165]:g00= −1 + 2V − 2V 2 + 8X̂ ,(13.39a)gi0= −4Vi − 8R̂i ,(13.39b)gij= δij (1 + 2V + 2V 2 ) + 4Ŵij ,(13.39c)−g= 1 + 2V + 4V 2 + 2Ŵkk .

(13.39d)√For spinning systems, we adopt the notation of Tagoshi, Ohashi, and Owen [166] where O(m, n)means to order cm for nonspinning terms and χcn for terms involving a single spin χ. (Here χ =|S|/m2 is the body’s dimensionless spin.) In this notation, terms we are interested in are of theorder O(3, 6), while the above post-Newtonian potentials have been obtained up to the following
orders [166, 167]:

V = O(2, 5),

Vj = O(3, 4),Ŵjk = O(4, 5), R̂j = O(5, 6).

(13.40)In terms of these post-Newtonian potentials, V , Vi , R̂i , X̂ and Ŵij , the perturbed metric density is

g00
g0i
gij

�
�
= −1 − 4V − 2 Ŵkk + 4V 2 + O(6, 7),��= −4Vi − 8 R̂i + V Vi + O(6, 7),��1= δij − 4 Ŵij − δij Ŵkk + O(6, 7).
2(13.41a)(13.41b)(13.41c)

184
As a consequence, equation (13.38) is given by

δpjA

=

� ��
�
��
1
−4Vj − 8 R̂j(S) + V(m) Vj(S)16π ∂A,k�� �
1−4 Ŵjk(S) − δjk Ŵii(S)
dΣk ,
2
,0

(13.42)where a subscript (S) means keep only the parts of those potentials proportional to the spins of thebodies, and a subscript (M ) involves monopolar pieces of the potential without spins (proportional
to the masses of the bodies). Terms without a subscript have both pieces.
Tagoshi, Ohashi and Owen express the potentials V(M ) , Vj , R̂j(S) and Ŵjk(S) in terms of thebodies’ masses, vectorial velocities, vectorial spins, and vectorial separations, and distance to thefield-point location (their equations (A1a), (A1d), (A1f), and (A1g)). While the full equations arequite lengthy, the portions that generate momentum flow—those involving the coupling of the massof one body to the spin of the other—are somewhat simpler. For convenience, we give these portionsof the equations below, rewritten in our notation, with the typos noted by Faye, Blanchet, andBounanno [167] corrected:

V(M )

=

Vj

=

Ŵjk(S)R̂j(S)mA+ (A ↔ B) ,
rA
� �
�
�
j
mA vA
3mB
mB (nA · nAB )3mBkl+ �jkl SA
nlA − 2
−
+
n
AB
2
rA
2rA rAB
4rA rAB
4rA rAB

(13.43a)+(A ↔ B) ,(13.43b)���niA 1 �
l k
l j
l m
=
�ijl SAvA + �ikl SA
vA − δjk �ilm vA
SA + (A ↔ B) ,(13.43c)
2
rA 2
� �
�
�
�
mB
mB
mB
mB
mB
k
= �jkl SA
nlA − 2++ nlAB −+ 2+
2
2rA rABrAB s22rA rAB
2rAB rB
rA s2
�
��
mB
mB
l
+nB
+
rA s2
rAB s2
�
�
�
�
mB
2mB
j
l
i
k
k
i
k mB
+nA �ikl SA nA (nAB + nB )+ 3
− 2nAB nB 3
rA s2
s
s
�
�
��
mB
2mB
j
l
i k mBiik+nAB �ikl SA −2nA nb 3 + (nA + nB )nAB −− 3s
rAB s2s+(A ↔ B) .

(13.43d)185Here, as before, mA , vA , and SA are the mass, velocity, and spin angular-momentum of object
A; rA is the separation of body A from a point in space and rAB is the separation of the two
objects; and nA and nAB are unit vectors pointing along rA and rAB , respectively. A new quantity,s = rA + rB + rAB , has been introduced, in addition. Inserting these expressions into equation(13.42) gives us δH jα0k ,α , and the momentum of body A is then found by performing a surface
integral over A’s surface. The surface integrals are computed under the same assumptions as insection 13.3.1; namely the separation of the bodies is much larger than their radii, and each surface
of integration is a sphere immediately above a body’s surface. When they are computed, they givethe same result as equation (13.35),

δpA = mA δvA +1 mB
2 SA × nAB .
3 rAB(13.44)As before, the momentum for body B is given by exchanging A and B.
As a consistency check, we can evaluate the system’s total momentum by doing a surface integral
at infinity:δpjtot =116π�δH jα0k ,α dΣk .

(13.45)S

The quantity δH jα0k ,α is exactly the same as above, from which one can find

δptot = mA δvA +mB
2 SA × nAB + (A ↔ B) .
rAB

(13.46)

This, combined with the fact that

δptot = δpA + δpB + δpfield ,(13.47)

as well as equation (13.44), givesδpfield =2mB
2 SA × nAB + (A ↔ B),
3rAB(13.48)186
as found in section 13.3.1.

13.4

Momentum Conservation

The total spin-induced momentum perturbation, δptot = δpA + δpB + δpfield (equations (13.35) and(13.20)) isδptot = mA δvA + mB δvB +12 (mB SA − mA SB ) × nAB .
rAB(13.49)Momentum conservation requires that the time derivative of this δptot vanish. The time derivativeof the kinetic terms can be read oﬀ the equation of motion (13.26):

mA

dδvA
dδvB
+ mB
dt
dt

(13.50)= −(mB SA − mA SB ) × [vAB − 3(nAB · vAB )nAB ] .
By inserting nAB = (xcm A −xcm B )/rAB into the second term of equation (13.49) and diﬀerentiatingwith respect to time, we obtain the negative of expression (13.50). Therefore,dδptot /dt = 0 ;(13.51)

i.e., as the binary’s evolution drives spin-induced momentum back and forth between the bodies andthe field, the total momentum remains conserved, as it must.
Interestingly, during the summation of momentum terms, one finds that the surface terms inδpA + δpB have exactly cancelled the field momentum δpfield , leaving the total momentum as thesum of the bodies’ kinetic term and their SSC term.

13.4.1

Momentum Conservation for Extreme-Kick Configuration

Let us examine in detail how momentum conservation is achieved in the presence of the bodies’
bobbing. Our detailed analysis (above) breaks each object’s momentum perturbation δpA,B into
three terms, the kinetic momentum mδvA,B , a term due to the SSC, and a surface integral term

187
Table 13.1: Pieces of Total Momentum
Kinetic
Body FrameSpinSSC
Surface
Total
Dragging
CurvaturepA−4aB × SB
3aA × SAaA × SA − 23 aA × SA − 23 aA × SA
pB
−4aA × SA
3aB × SBaB × SB − 23 aB × SB − 23 aB × SB2
pfield3 (aA × SA + aB × SB )Spin-dependent, time-varying pieces of body and field momenta at 1.5PN order, for theextreme-kick binary (circular orbit with spins antialigned and in the orbital plane). The body
momenta are broken down into kinetic, SSC and surface terms and are expressed in terms of the2
bodies’ spins SA,B and Newtonian-order gravitational accelerations aA,B = −mB,A nAB,BA /rAB.
See equations ((13.26)) and ((13.35)) for a similar decomposition in a generic binary.
(see table 13.1). The total kinetic momentummδvA + mδvB = − (aA × SA + aB × SB ) �= 0 ,(13.52)

is not conserved because of the noncancellation between the frame-dragging and spin-curvature
coupling terms. The total body momentum δpA + δpB is not conserved either; it sums up to 2/3
the total kinetic momentum:2
δpA + δpB = − (aA × SA + aB × SB ) �= 0 ,3(13.53)

To achieve momentum conservation, there is a nonzero spin-dependent total field momentum distributed outside of the bodies, with

δpfield =2(aA × SA + aB × SB ) .
3(13.54)Note that this total external field momentum is only −2/3 the spin-dependent total kinetic momentum—
instead of the −1 that one might have naively expected.

188

13.5

Summary

We have studied generic binary compact objects and mometum flow that occurs between the bodies
and the surrounding spacetime using the Landau-Lifshitz formalism. In this analysis, we have foundthat the momentum of the absolute motion of the binary system is balanced by the momentum
contained in the surrounding gravitational field. These results were also specialized to the interestingcase of spinning black holes in the extreme-kick configuration. We think this formalism holds great
promise for understanding the dynamics of compact objects and anticipate additional insights andintuition gained from applying this sort of analysis to other interacting black hole environments.

189

Bibliography[1] C. W. Misner, K. S. Thorne, and J. A. Wheeler, Gravitation (W. H. Freeman and Company,
NY, 1973).
[2] K. S. Thorne, Reviews of Modern Physics 52(2), 299 (1980).
[3] L. Blanchet, B. R. Iyer, C. M. Will, and A. G. Wiseman, Classical and Quantum Gravity
13(4), 575 (1996).
[4] A. Buonanno and T. Damour, Physical Review D 59(8), 084006 (1999).
[5] A. Buonanno and T. Damour, Physical Review D 62(6), 064015 (2000).
[6] T. Damour, P. Jaranowski, and G. Schäfer, Physical Review D 62(8), 084011 (2000).
[7] T. Damour, Physical Review D 64(12), 124013 (2001).
[8] T. Damour, B. R. Iyer, and B. S. Sathyaprakash, Physical Review D 63(4), 044023 (2001).
[9] T. Damour, B. R. Iyer, and B. S. Sathyaprakash, Physical Review D 57(2), 885 (1998).
[10] A. Buonanno, Y. Chen, and M. Vallisneri, Physical Review D 67(10), 104025 (2003), erratumibid. 74 (2006) 029904(E).
[11] L. E. Kidder, Physical Review D 52(2), 821 (1995).
[12] L. Blanchet, T. Damour, B. R. Iyer, C. M. Will, and A. G. Wiseman, Physical Review Letters
74(18), 3515 (1995).
[13] L. Blanchet, Physical Review D 54(2), 1417 (1996).

190[14] L. Blanchet, G. Faye, B. R. Iyer, and B. Joguet, Physical Review D 65(6), 061501 (2002).
[15] L. E. Kidder, C. M. Will, and A. G. Wiseman, Physical Review D 47(10), R4183 (1993).
[16] T. A. Apostolatos, C. Cutler, G. J. Sussman, and K. S. Thorne, Physical Review D 49(12),6274 (1994).
[17] B. Abbott, R. Abbott, R. Adhikari, A. Ageev, B. Allen, R. Amin, S. B. Anderson, W. G.
Anderson, M. Araya, H. Armandula, et al. (LIGO Scientific Collaboration), Physical Review
D 69(12), 122001 (2004).
[18] B. Abbott, R. Abbott, R. Adhikari, A. Ageev, B. Allen, R. Amin, S. B. Anderson, W. G.
Anderson, M. Araya, H. Armandula, et al. (LIGO Scientific Collaboration), Physical Review
D 72(8), 082001, 082001 (2005).
[19] B. Abbott, R. Abbott, R. Adhikari, A. Ageev, B. Allen, R. Amin, S. B. Anderson, W. G.
Anderson, M. Araya, H. Armandula, et al. (LIGO Scientific Collaboration), Physical Review
D 72(8), 082002, 082002 (2005).
[20] B. Abbott, R. Abbott, R. Adhikari, A. Ageev, J. Agresti, P. Ajith, B. Allen, J. Allen, R. Amin,
S. B. Anderson, et al. (LIGO Scientific Collaboration), Physical Review D 73(6), 062001,
062001 (2006).
[21] B. Abbott, R. Abbott, R. Adhikari, J. Agresti, P. Ajith, B. Allen, R. Amin, S. B. Anderson,
W. G. Anderson, M. Arain, et al. (LIGO Scientific Collaboration), Physical Review D 77(6),062002, 062002 (2008).
[22] B. Abbott, R. Abbott, R. Adhikari, J. Agresti, P. Ajith, B. Allen, R. Amin, S. B. Anderson,
W. G. Anderson, M. Arain, et al. (LIGO Scientific Collaboration), Physical Review D 78(4),042002, 042002 (2008).
[23] B. Abbott, R. Abbott, R. Adhikari, J. Agresti, P. Ajith, B. Allen, R. Amin, S. B. Anderson, W. G. Anderson, M. Arain, et al. (LIGO Scientific Collaboration and K. Hurley), The
Astrophysical Journal 681(2), 1419 (2008).

191[24] A. J. Weinstein (for the LIGO Scientific Collaboration), Classical and Quantum Gravity 21(5),S677 (2004).
[25] B. Abbott, R. Abbott, R. Adhikari, A. Ageev, B. Allen, R. Amin, S. B. Anderson, W. G.
Anderson, M. Araya, H. Armandula, et al. (LIGO Scientific Collaboration), Physical Review
D 72(4), 042002 (2005).
[26] B. Abbott, R. Abbott, R. Adhikari, A. Ageev, J. Agresti, B. Allen, J. Allen, R. Amin, S. B.
Anderson, W. G. Anderson, et al. (LIGO Scientific Collaboration), Physical Review D 72(6),062001 (2005).
[27] B. Abbott, R. Abbott, R. Adhikari, A. Ageev, J. Agresti, P. Ajith, B. Allen, J. Allen, R. Amin,
S. B. Anderson, et al. (LIGO Scientific Collaboration and TAMA Collaboration), PhysicalReview D 72(12), 122004 (2005).
[28] B. Abbott, R. Abbott, R. Adhikari, J. Agresti, P. Ajith, B. Allen, J. Allen, R. Amin, S. B.
Anderson, W. G. Anderson, et al. (LIGO Scientific Collaboration), Classical and Quantum
Gravity 23(8), S29 (2006).
[29] L. Baggio, M. Bignotto, M. Bonaldi, M. Cerdonio, M. D. Rosa, P. Falferi, S. Fattori, P. Fortini,
G. Giusfredi, M. Inguscio, et al. (LIGO Scientific Collaboration and AURIGA Collaboration),
Classical and Quantum Gravity 25(9), 095004 (2008).
[30] B. Abbott, R. Abbott, R. Adhikari, J. Agresti, P. Ajith, B. Allen, R. Amin, S. B. Anderson,
W. G. Anderson, M. Arain, et al. (LIGO Scientific Collaboration), Physical Review D 76(6),062003 (2007).
[31] B. Abbott, R. Abbott, R. Adhikari, J. Agresti, P. Ajith, B. Allen, R. Amin, S. B. Anderson,
W. G. Anderson, M. Arain, et al. (LIGO Scientific Collaboration), Classical and Quantum
Gravity 24(22), 5343 (2007).

192[32] B. Abbott, R. Abbott, R. Adhikari, J. Agresti, P. Ajith, B. Allen, R. Amin, S. B. Anderson,
W. G. Anderson, M. Arain, et al. (LIGO Scientific Collaboration), Physical Review D 77(6),062004 (2008).
[33] B. Abbott, R. Abbott, R. Adhikari, P. Ajith, B. Allen, G. Allen, R. Amin, S. B. Anderson,
W. G. Anderson, M. A. Arain, et al. (LIGO Scientific Collaboration), Classical and Quantum
Gravity 25(24), 245008 (2008).
[34] B. Abbott, R. Abbott, R. Adhikari, P. Ajith, B. Allen, G. Allen, R. Amin, S. B. Anderson,
W. G. Anderson, M. A. Arain, et al. (LIGO Scientific Collaboration), Physical Review Letters
101(21), 211102 (2008).
[35] B. Abbott, R. Abbott, R. Adhikari, A. Ageev, B. Allen, R. Amin, S. B. Anderson, W. G.
Anderson, M. Araya, H. Armandula, et al. (LIGO Scientific Collaboration), Physical Review
D 69(8), 082004 (2004).
[36] B. Abbott, R. Abbott, R. Adhikari, A. Ageev, J. Agresti, B. Allen, J. Allen, R. Amin, S. B.
Anderson, W. G. Anderson, et al. (LIGO Scientific Collaboration), Physical Review D 72(10),102004 (2005).
[37] B. Abbott, R. Abbott, R. Adhikari, A. Ageev, B. Allen, R. Amin, S. B. Anderson, W. G.
Anderson, M. Araya, H. Armandula, et al. (LIGO Scientific Collaboration), Physical Review
Letters 94(18), 181103 (2005).
[38] B. Abbott, R. Abbott, R. Adhikari, J. Agresti, P. Ajith, B. Allen, R. Amin, S. B. Anderson,
W. G. Anderson, M. Arain, et al. (LIGO Scientific Collaboration), Physical Review D 76(8),082001 (2007).
[39] B. Abbott, R. Abbott, R. Adhikari, J. Agresti, P. Ajith, B. Allen, R. Amin, S. B. Anderson,
W. G. Anderson, M. Arain, et al. (LIGO Scientific Collaboration), Physical Review D 76(4),042001 (2007).

193[40] B. Abbott, R. Abbott, R. Adhikari, J. Agresti, P. Ajith, B. Allen, R. Amin, S. B. Anderson,
W. G. Anderson, M. Arain, et al. (LIGO Scientific Collaboration), Physical Review D 77(2),022001 (2008).
[41] B. Abbott, R. Abbott, R. Adhikari, P. Ajith, B. Allen, G. Allen, R. Amin, D. P. Anderson,
S. B. Anderson, W. G. Anderson, et al. (LIGO Scientific Collaboration), Physical Review D
79(2), 022001 (2009).
[42] B. Abbott, R. Abbott, R. Adhikari, P. Ajith, B. Allen, G. Allen, R. Amin, S. B. Anderson,
W. G. Anderson, M. A. Arain, et al. (LIGO Scientific Collaboration), The AstrophysicalJournal Letters 683(1), L45 (2008).
[43] B. P. Abbott, R. Abbott, R. Adhikari, P. Ajith, B. Allen, G. Allen, R. S. Amin, S. B. Anderson,
W. G. Anderson, M. A. Arain, et al. (LIGO Scientific Collaboration), Physical Review Letters
102(11), 111102 (2009).
[44] B. Abbott, R. Abbott, R. Adhikari, A. Ageev, B. Allen, R. Amin, S. B. Anderson, W. G.
Anderson, M. Araya, H. Armandula, et al. (LIGO Scientific Collaboration), Physical Review
D 69(12), 122004 (2004).
[45] B. Abbott, R. Abbott, R. Adhikari, J. Agresti, P. Ajith, B. Allen, J. Allen, R. Amin, S. B.
Anderson, W. G. Anderson, et al. (LIGO Scientific Collaboration), Physical Review Letters
95(22), 221101 (2005).
[46] B. Abbott, R. Abbott, R. Adhikari, J. Agresti, P. Ajith, B. Allen, R. Amin, S. B. Anderson,
W. G. Anderson, M. Araya, et al. (LIGO Scientific Collaboration), The Astrophysical Journal659(2), 918 (2007).
[47] B. Abbott, R. Abbott, R. Adhikari, J. Agresti, P. Ajith, B. Allen, R. Amin, S. B. Anderson,
W. G. Anderson, M. Arain, et al. (LIGO Scientific Collaboration), Physical Review D 76(8),082003 (2007).

194[48] B. Abbott, R. Abbott, R. Adhikari, J. Agresti, P. Ajith, B. Allen, R. Amin, S. B. Anderson,
W. G. Anderson, M. Arain, et al. (LIGO Scientific Collaboration and ALLEGRO Collaboration), Physical Review D 76(2), 022001 (2007).
[49] H. Collins, Gravity’s Shadow: The Search for Gravitational Waves (The University of Chicago
Press, Chicago, 2004).
[50] J. Weber, Physical Review Letters 17(24), 1228 (1966).
[51] J. Weber, Physical Review Letters 22(24), 1320 (1969).
[52] J. Weber, Physical Review Letters 25(3), 180 (1970).
[53] B. L. Brown, A. P. Mills, and J. A. Tyson, Physical Review D 26(6), 1209 (1982).
[54] E. Mauceli, Z. K. Geng, W. O. Hamilton, W. W. Johnson, S. Merkowitz, A. Morse, B. Price,
and N. Solomonson, Physical Review D 54(2), 1264 (1996).
[55] G. A. Prodi, L. Conti, R. Mezzena, S. Vitale, L. Taﬀarello, J.-P. Zendri, L. Baggio, M. Cerdonio, A. Colombo, V. Crivelli Visconti, et al., in E. Coccia, G. Veneziano, and P. G., eds.,
GRAVITATIONAL WAVES - Proceedings of the Second Edoardo Amaldi Conference (World
Scientific, Singapore, 1999), pp. 148–158.
[56] P. Astone, M. Bassan, P. Bonifazi, P. Carelli, M. G. Castellano, G. Cavallari, E. Coccia,
C. Cosmelli, V. Fafone, S. Frasca, et al., Physical Review D 47(2), 362 (1993).
[57] P. Astone, M. Bassan, P. Bonifazi, P. Carelli, E. Coccia, C. Cosmelli, V. Fafone, S. Frasca,
A. Marini, G. Mazzitelli, et al., Astroparticle Physics 7(3), 231 (1997).
[58] D. G. Blair, E. N. Ivanov, M. E. Tobar, P. J. Turner, F. van Kann, and I. S. Heng, PhysicalReview Letters 74(11), 1908 (1995).
[59] R. L. Forward, General Relativity and Gravitation 2(2), 149 (1971).

195[60] A. de Waard, Y. Benzaim, G. Frossati, L. Gottardi, H. van der Mark, J. Flokstra, M. Podt,
M. Bassan, Y. Minenkov, A. Moleti, et al., Classical and Quantum Gravity 22(10), S215
(2005).
[61] O. D. Aguiar, L. A. Andrade, J. J. Barroso, F. Bortoli, L. A. Carneiro, P. J. Castro, C. A.
Costa, K. M. F. Costa, J. C. N. de Araujo, A. U. de Lucena, et al., Classical and Quantum
Gravity 23(8), S239 (2006).
[62] R. L. Forward, Physical Review D 17(2), 379 (1978).
[63] R. Weiss, in Lincoln Research Laboratory of Electronics Quarterly Progress Report (1972), no.
105, pp. 54–76.
[64] K. Tsubono (TAMA Collaboration), in K. Tsubono, M.-K. Fujimoto, and K. Kuroda, eds.,
Proceedings of the TAMA International Workshop on Gravitational Wave Detection (Saitama,
1996) (Universal Academic Press, Tokyo, 1997), p. 183.
[65] K. Danzmann, in E. Coccia, G. Pizzella, and F. Ronga, eds., First Edoardo Amaldi Conferenceon Gravitational Wave Experiments (World Scientific, Singapore, 1995), p. 100.
[66] B. Caron, A. Dominjon, C. Drezen, R. Flaminio, X. Grave, F. Marion, L. Massonnet,
C. Mehmel, R. Morand, B. Mours, et al., Classical and Quantum Gravity 14(6), 1461 (1997).
[67] A. Abramovici, W. E. Althouse, R. W. P. Drever, Y. Gursel, S. Kawamura, F. J. Raab,
D. Shoemaker, L. Sievers, R. E. Spero, K. S. Thorne, et al., Science 256(5055), 325 (1992).
[68] LISA homepage, accessed April 1, 2009, http://lisa.nasa.gov/.
[69] S. Kawamura, M. Ando, T. Nakamura, K. Tsubono, T. Tanaka, I. Funaki, N. Seto, K. Numata,S. Sato, K. Ioka, et al., Journal of Physics: Conference Series 122, 012006 (2008).
[70] M. V. Sazhin, Soviet Astronomy 22, 36 (1978).
[71] S. Detweiler, The Astrophysical Journal 234(3), 1100 (1979).
[72] The Planck Collaboration, ArXiv Astrophysics e-prints (2006), arXiv:astro-ph/0604069.

196[73] M. R. Nolta, J. Dunkley, R. S. Hill, G. Hinshaw, E. Komatsu, D. Larson, L. Page, D. N.
Spergel, C. L. Bennett, B. Gold, et al., The Astrophysical Journal Supplement Series 180(2),296 (2009).
[74] B. G. Keating, P. A. R. Ade, J. J. Bock, E. Hivon, W. L. Holzapfel, A. E. Lange, H. Nguyen,
and K. W. Yoon, in S. Fineschi, ed., Society of Photo-Optical Instrumentation Engineers(SPIE) Conference Series (2003), vol. 4843 of Presented at the Society of Photo-Optical Instrumentation Engineers (SPIE) Conference, pp. 284–295.
[75] C. Pryke, P. Ade, J. Bock, M. Bowden, M. L. Brown, G. Cahill, P. G. Castro, S. Church,
T. Culverhouse, R. Friedman, et al., The Astrophysical Journal 692(2), 1247 (2009).
[76] A. C. Taylor, New Astronomy Reviews 50(11-12), 993 (2006).
[77] C. Bischoﬀ, L. Hyatt, J. J. McMahon, G. W. Nixon, D. Samtleben, K. M. Smith, K. Vanderlinde, D. Barkats, P. Farese, T. Gaier, et al., The Astrophysical Journal 684(2), 771 (2008).
[78] D. Samtleben (for the QUIET collaboration), ArXiv e-prints (2008), 0806.4334.
[79] P. Oxley, P. A. Ade, C. Baccigalupi, P. deBernardis, H.-M. Cho, M. J. Devlin, S. Hanany,
B. R. Johnson, T. Jones, A. T. Lee, et al., in M. Strojnik, ed., Society of Photo-Optical
Instrumentation Engineers (SPIE) Conference Series (2004), vol. 5543 of Presented at theSociety of Photo-Optical Instrumentation Engineers (SPIE) Conference, pp. 320–331.
[80] B. P. Crill, P. A. R. Ade, E. S. Battistelli, S. Benton, R. Bihary, J. J. Bock, J. R. Bond,
J. Brevik, S. Bryan, C. R. Contaldi, et al., in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series (2008), vol. 7010 of Presented at the Society of Photo-Optical
Instrumentation Engineers (SPIE) Conference.
[81] D. Baumann, M. G. Jackson, P. Adshead, A. Amblard, A. Ashoorioon, N. Bartolo, R. Bean,
M. Beltran, F. de Bernardis, S. Bird, et al., ArXiv e-prints (2008), 0811.3919.
[82] J.

Hough

and

S.

Rowan,

Livinghttp://www.livingreviews.org/lrr-2000-3.

Reviews

inRelativity

3(3)(2000),

197[83] P. R. Saulson, Physical Review D 42(8), 2437 (1990).
[84] D. Sigg (for the LIGO Scientific Collaboration), Classical and Quantum Gravity 25(11), 114041(2008).
[85] X. Siemens, B. Allen, J. Creighton, M. Hewitson, and M. Landry, Classical and Quantum
Gravity 21(20), S1723 (2004).
[86] Y. Gürsel and M. Tinto, Physical Review D 40(12), 3884 (1989).
[87] O. Krause, S. M. Birkmann, T. Usuda, T. Hattori, M. Goto, G. H. Rieke, and K. A. Misselt,
Science 320(5880), 1195 (2008).
[88] J. L. Blum and J. M. Miller, in Bulletin of the American Astronomical Society (2009), vol. 41,
p. 384.
[89] R. A. Hulse and J. H. Taylor, Astrophysical Journal 195, L51 (1975).
[90] CBC Group, Rate predictions for compact binary coalescence events (2009), http://www.lscgroup.phys.uwm.edu/ligovirgo/cbc/protected/papers/rates/rates.pdf.
[91] V. Kalogera, C. Kim, D. R. Lorimer, M. Burgay, N. D’Amico, A. Possenti, R. N. Manchester,A. G. Lyne, B. C. Joshi, M. A. McLaughlin, et al., The Astrophysical Journal Letters 601(2),L179 (2004), http://stacks.iop.org/1538-4357/601/L179.
[92] V. Kalogera, C. Kim, D. R. Lorimer, M. Burgay, N. D’Amico, A. Possenti, R. N. Manchester,A. G. Lyne, B. C. Joshi, M. A. McLaughlin, et al., The Astrophysical Journal Letters 614(2),L137 (2004).
[93] C. Kim, V. Kalogera, D. R. Lorimer, , and T. White, The Astrophysical Journal 616(2), 1109
(2004).
[94] C. Kim, V. Kalogera, and D. R. Lorimer, in A life with stars (2006).
[95] R. O’Shaughnessy, C. Kim, V. Kalogera, , and K. Belczynski, The Astrophysical Journal672(1), 479 (2008).

198[96] R. O’Shaughnessy, C. Kim, T. Fragos, V. Kalogera, , and K. Belczynski, The AstrophysicalJournal 633(2), 1076 (2005).
[97] R. M. O’Leary, F. A. Rasio, J. M. Fregeau, N. Ivanova, , and R. O’Shaughnessy, The Astrophysical Journal 637(2), 937 (2006).
[98] R. M. O’Leary, R. O’Shaughnessy, and F. A. Rasio, Physical Review D 76(6), 061504 (2007).
[99] R. M. O’Leary, B. Kocsis, and A. Loeb, ArXiv e-prints (2008), 0807.2638.
[100] M. C. Miller and V. M. Lauburg, The Astrophysical Journal 692(1), 917 (2009).
[101] M. Freitag, M. A. Gürkan, and F. A. Rasio, Monthly Notices of the Royal Astronomical Society
368, 141 (2006).
[102] E. Salpeter, The Astrophysical Journal 121, 161 (1955).
[103] M. J. Benacquista, Living Reviews in Relativity 9(2) (2006), http://www.livingreviews.org/lrr2006-2.
[104] S. Sigurdsson and E. S. Phinney, The Astrophysical Journal 415, 631 (1993).
[105] S. Sigurdsson and L. Hernquist, Nature 364, 423 (1993).
[106] S. R. Kulkarni, P. Hut, and S. McMillan, Nature 364, 421 (1993).
[107] D. Guetta and L. Stella, ArXiv e-prints (2008), 0811.0684.
[108] J. Goodman, The Astrophysical Journal 308, L47 (1986).
[109] B. Paczynski, The Astrophysical Journal 308, L43 (1986).
[110] D. Eichler, M. Livio, T. Piran, and D. N. Schramm, Nature 340, 126 (1989).
[111] R. Narayan, B. Paczynski, and T. Piran, The Astrophysical Journal 395, L83 (1992).
[112] J. S. Bloom, S. Sigurdsson, and O. R. Pols, Monthly Notices of the Royal Astronomical Society
305, 763 (1999).

199[113] K. Belczynski, V. Kalogera, , and T. Bulik, The Astrophysical Journal 572(1), 407 (2002).
[114] K. Belczynski, R. Perna, T. Bulik, V. Kalogera, N. Ivanova, and D. Lamb, in 36th COSPAR
Scientific Assembly (2006), vol. 36 of COSPAR, Plenary Meeting, p. 3470.
[115] D. J. Champion, D. R. Lorimer, M. A. McLaughlin, J. M. Cordes, Z. Arzoumanian, J. M.
Weisberg, and J. H. Taylor, Monthly Notices of the Royal Astronomical Society 350, L61(2004).
[116] K. Belczynski, R. E. Taam, V. Kalogera, F. A. Rasio, , and T. Bulik, The AstrophysicalJournal 662(1), 504 (2007).
[117] C. Hopman, D. Guetta, E. Waxman, , and S. P. Zwart, The Astrophysical Journal Letters643(2), L91 (2006).
[118] D. B. Fox, D. A. Frail, P. A. Price, S. R. Kulkarni, E. Berger, T. Piran, A. M. Soderberg, S. B.
Cenko, P. B. Cameron, A. Gal-Yam, et al., Nature 437, 845 (2005).
[119] A. M. Soderberg, E. Berger, M. Kasliwal, D. A. Frail, P. A. Price, B. P. Schmidt, S. R.
Kulkarni, D. B. Fox, S. B. Cenko, A. Gal-Yam, et al., The Astrophysical Journal 650(1), 261(2006).
[120] N. Weiner, Extrapolation, interpolation, and smoothing of stationary time series with engineering applications (M.I.T. Technology Press and John Wiley and Sons, Inc., New York,1950).
[121] B. J. Owen, Physical Review D 53(12), 6749 (1996).
[122] B. J. Owen and B. S. Sathyaprakash, Physical Review D 60(2), 022002 (1999).
[123] T. Cokelaer, Physical Review D 76(10), 102004 (2007).
[124] B. S. Sathyaprakash and S. V. Dhurandhar, Physical Review D 44(12), 3819 (1991).
[125] B. S. Sathyaprakash, Physical Review D 50(12), R7111 (1994).

200[126] LSC Algorithm Library software packages lal, lalwrapper, and lalapps, http://www.lscgroup.phys.uwm.edu/lal.
[127] K. S. Thorne, in S. W. Hawking and W. Israel, eds., Three hundred years of gravitation(Cambridge University Press, Cambridge, 1987), chap. 9, pp. 330–458.
[128] S. Droz, D. J. Knapp, E. Poisson, and B. J. Owen, Physical Review D 59(12), 124016 (1999).
[129] A. S. Sengupta, J. A. Gupchup, and C. A. K. Robinson, In Preparation (2008).
[130] R. Biswas, P. R. Brady, J. D. E. Creighton, and S. Fairhurst, ArXiv e-prints (2007), 0710.0465.
[131] P. R. Brady and S. Fairhurst, Classical and Quantum Gravity 25(10), 1050002 (2008).
[132] B. Allen, J. K. Blackburn, P. R. Brady, J. D. E. Creighton, T. Creighton, S. Droz, A. D.
Gillespie, S. A. Hughes, S. Kawamura, T. T. Lyons, et al., Physical Review Letters 83(8),1498 (1999).
[133] P. R. Brady, J. D. E. Creighton, and A. G. Wiseman, Classical and Quantum Gravity 21(20),
S1775 (2004).
[134] R. K. Kopparapu, C. Hanna, V. Kalogera, R. O’Shaughnessy, G. Gonzalez, P. R. Brady, andS. Fairhurst, The Astrophysical Journal 675(2), 1459 (2008).
[135] E. S. Phinney, The Astrophysical Journal 380, L17 (1991).
[136] R. O’Shaughnessy, K. Belczynski, , and V. Kalogera, The Astrophysical Journal 675(1), 566(2008).
[137] LIGO Scientific Collaboration, Calibration of the LIGO detectors for S5, Tech. Rep. In Preparation, LIGO Project.
[138] Y. Pan, A. Buonanno, J. G. Baker, J. Centrella, B. J. Kelly, S. T. McWilliams, F. Pretorius,
and J. R. van Meter, Physical Review D (Particles, Fields, Gravitation, and Cosmology) 77(2),024014 (2008).

201[139] R. Gouaty (for the LIGO Scientific Collaboration), Classical and Quantum Gravity 25(18),184006 (2008).
[140] S. E. Thorsett and D. Chakrabarty, The Astrophysical Journal 512(1), 288 (1999).
[141] R. N. Manchester, G. B. Hobbs, A. Teoh, and M. Hobbs, The Astronomical Journal 129(4),1993 (2005), http://www.atnf.csiro.au/research/pulsar/psrcat/.
[142] LIGO Scientific Collaboration and Virgo Collaboration, In Preparation (2009).
[143] S. Bose, A. Pai, and S. Dhurandhar, International Journal of Modern Physics D 9(3), 325(2000).
[144] A. Pai, S. Dhurandhar, and S. Bose, Physical Review D 64(4), 042004 (2001).
[145] F. Marion (on behalf of the Virgo Collaboration), in J. Dumarchez and J. T. T. Vân, eds.,
Gravitational waves and experimental gravity—Proceedings of the XXXVIIIth RECONTRESDE MORIOND (The Gioi Publishers, 2004), pp. 145–150.
[146] K. Cannon, A. Chapman, C. Hanna, D. Keppel, A. C. Searle, and A. J. Weinstein, In Preparation (2009).
[147] F. Zwicky, The Astrophysical Journal 86, 217 (1937).
[148] K. C. Freeman, The Astrophysical Journal 160, 811 (1970).
[149] C. M. Will, Physical Review D 57(4), 2061 (1998).
[150] L. S. Finn and P. J. Sutton, Physical Review D 65(4), 044022 (2002).
[151] S. R. Choudhury, G. C. Joshi, S. Mahajan, and B. H. J. McKellar, Astroparticle Physics 21(5),559 (2004).
[152] P. Ajith, Classical and Quantum Gravity 25(11), 114033 (2008).
[153] P. Ajith, S. Babak, Y. Chen, M. Hewitson, B. Krishnan, A. M. Sintes, J. T. Whelan,
B. Brügmann, P. Diener, N. Dorband, et al., Physical Review D 77(10), 104017 (2008).

202[154] C. Talmadge, J. P. Berthias, R. W. Hellings, and E. M. Standish, Physical Review Letters61(10), 1159 (1988).
[155] F. Pretorius, Physical Review Letters 95(12), 121101 (2005).
[156] M. Campanelli, C. Lousto, Y. Zlochower, , and D. Merritt, The Astrophysical Journal Letters659(1), L5 (2007).
[157] M. Campanelli, C. O. Lousto, Y. Zlochower, and D. Merritt, Physical Review Letters 98(23),231102 (2007).
[158] J. Healy, F. Herrmann, I. Hinder, D. M. Shoemaker, P. Laguna, and R. A. Matzner, PhysicalReview Letters 102(4), 041101 (2009).
[159] F. Pretorius, in M. Colpi, P. Casella, V. Gorini, U. Moschella, and A. Possenti, eds., Physics
of Relativistic Objects in Compact Binaries: From Birth to Coalescence (Springer Verlag andCanopus Publishing Limited, 2009), vol. 359, pp. 305–369.
[160] L. D. Landau and E. M. Lifshitz, Classical Theory of Fields (Addison Wesley, Redding Mass.,
1962), p. Sec. 100.
[161] J. D. Kaplan, D. A. Nichols, and K. S. Thorne, submitted to Physical Review D (2009).
[162] L. E. Kidder, Physical Review D 52(2), 821 (1995).
[163] K. S. Thorne, Reviews of Modern Physics 52(2), 299 (1980).
[164] K. S. Thorne and J. B. Hartle, Physical Review D 31(8), 1815 (1985).
[165] L. Blanchet, G. Faye, and B. Ponsot, Physical Review D 58(12), 124002 (1998).
[166] H. Tagoshi, A. Ohashi, and B. J. Owen, Physical Review D 63(4), 044006 (2001).
[167] G. Faye, L. Blanchet, and A. Buonanno, Physical Review D 74(10), 104033 (2006).