USING MACHINE LEARNING AND THE HOUGH TRANSFORM TO SEARCH FORGRAVITATIONAL WAVES DUE TO R-MODE EMISSION BY ISOLATED NEUTRON STARS

By
ANDREW L. MILLERA DISSERTATION PRESENTED TO THE GRADUATE SCHOOLOF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT
OF THE REQUIREMENTS FOR THE DEGREE OF
DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
2019

c 2019 Andrew L. Miller

ACKNOWLEDGMENTS
I have way too many people to thank for making this joint degree work.
First I thank my supervisors, Bernard Whiting and Pia Astone, for the tons of work they
did to make this joint degree happen. I am very grateful to have met both of them. In detail,I thank Bernard Whiting, who has been a great adviser and friend. He was always been open
to new ideas, and I greatly enjoyed these years working with him. And of course, I liked all the
Italian, German and Syrian food and drinks we had together. I also thank Pia Astone, whoalways made time to help me. No matter the problem I had, she had the solution, and she wasalways so willing to help. It didn’t matter how stupid my problem was, she always listened. I
will greatly miss talking to her and seeing her all the time.
I thank Cristiano Palomba, who was always there to bounce ideas off of and for following
my work so closely, and Paola Leaci, who helped me greatly when I tried to first use computing
clusters such as CNAF to submit jobs and for always having my back. Of course I thank Sergio
Frasca, for knowing the answer to every question and for the Italian snacks he has given me
over the past five years. Also I thank Guido Mueller, for being the best graduate student
coordinator, approving this joint degree, having a great sense of humor, and suggesting thesolution to support my last semester at University of Florida.
I thank Guenakh Mitselmakher, for approving the funds to support me during my lastsemester at the University of Florida, and Ettore Majorana, for organizing my exchangewith the Institute of Cosmic Ray Research in Tokyo and always finding funds to support myattendance at Virgo meetings and conferences, and Sabrina D’Antonio, for debugging mycodes and always being so welcoming of me in Rome. I also thank Anna de Grossi, for helping
Pia and me navigate the bureaucracy of an Italian university. I am very thankful for all of thesecretaries and administrators who approved this agreement at every level. Finally I thank myparents and my aunt, whose continued support and almost daily Skype calls got me throughmy Ph.D.

3

TABLE OF CONTENTS
page
ACKNOWLEDGMENTS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3

LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

8

ABSTRACT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

10

CHAPTER1

INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

12

2

THE ORIGIN OF GRAVITATIONAL WAVES . . . . . . . . . . . . . . . . . . . . .

14

2.1
2.2
2.3

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Einstein Field Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Gravitational Wave Solutions in the Linearized Regime . . . . . . . . . . . . .
2.3.1 Gravitational Wave Polarizations . . . . . . . . . . . . . . . . . . . . .
2.3.2 Effect of Gravitational Waves on Point Particles . . . . . . . . . . . . .
Sources of Gravitational Waves . . . . . . . . . . . . . . . . . . . . . . . . .
2.4.1 Compact Binaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.4.2 Core Collapse Supernovae Explosions . . . . . . . . . . . . . . . . . .
2.4.3 Stochastic Background . . . . . . . . . . . . . . . . . . . . . . . . . .
2.4.4 Neutron Stars . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

14
14
16
20
21
22
22
24
24
25

THE ASTROPHYSICS OF NEUTRON STARS . . . . . . . . . . . . . . . . . . . .

26

3.1
3.2

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Neutron Stars . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.1 Magnetic Field Strength . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.2 Break-up Frequency . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.3 Indirectly Measuring an Upper Limit Magnetic Field of a Neutron Star .
Gravitational Waves from Deformed, Rotating Neutron Stars . . . . . . . . . .
Continuous Gravitational Waves . . . . . . . . . . . . . . . . . . . . . . . . .
Transient Continuous Gravitational Waves . . . . . . . . . . . . . . . . . . .
Gravitational Waves from R-modes . . . . . . . . . . . . . . . . . . . . . . .
3.6.1 Velocity and Density Perturbations . . . . . . . . . . . . . . . . . . .
3.6.2 Frequency Evolution of R-modes . . . . . . . . . . . . . . . . . . . . .
3.6.3 Current Quadrupole Radiation . . . . . . . . . . . . . . . . . . . . . .
3.6.4 CFS Instability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Cowling Approximation . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Equation of State . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

26
26
27
27
28
29
32
34
35
36
37
39
39
41
41

THE DAWN OF MULTI-MESSENGER ASTRONOMY . . . . . . . . . . . . . . . .

42

4.1
4.2

42
42

2.4

3

3.3
3.4
3.5
3.6

3.7
3.8
4

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
GW150914 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

4.3

5

GW170817 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3.1 Constraints on Difference between the Speeds of Light and Gravity . .
4.3.2 Constraints on Equation of State . . . . . . . . . . . . . . . . . . . .

44
44
45

GRAVITATIONAL WAVE DETECTORS . . . . . . . . . . . . . . . . . . . . . . .

46

5.1
5.2
5.3

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Detection Principle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Noise Sources . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.3.1 Seismic Noise . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.3.2 Thermal Noise . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.3.3 Shot Noise . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.3.4 Radiation Pressure Noise . . . . . . . . . . . . . . . . . . . . . . . . .
Detector Design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.4.1 Simple Michelson Interferometer . . . . . . . . . . . . . . . . . . . . .
5.4.2 The Fabry-Perot Interferometer . . . . . . . . . . . . . . . . . . . . .

46
46
48
48
50
51
56
57
57
58

PROSPECTS OF DETECTING NEUTRON STARS: RATES AND IMPLICATIONS

59

6.1
6.2
6.3
6.4
6.5

59
59
61
61

5.4

6

7

8

9

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Merger Rates: What is Known . . . . . . . . . . . . . . . . . . . . . . . . .
Optimizing Sensitivity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Implications of a Detection . . . . . . . . . . . . . . . . . . . . . . . . . . .
Rates and Possible Detections of Newly Born Neutron Stars in Current andFuture Detectors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

63

DATA ANALYSIS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

68

7.1
7.2
7.3

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Ideal Detection Technique: Matched Filtering . . . . . . . . . . . . . . . . . .
Time Delays . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.3.1 Roemer Delay . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.3.2 Shapiro Delay . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.3.3 Einstein Delay . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

68
68
70
71
71
73

THE ORIGINAL FREQUENCY HOUGH . . . . . . . . . . . . . . . . . . . . . . .

75

8.1
8.2
8.3
8.4
8.5

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Construction of Short Fast Fourier Transform Databases . . . . . . . . . . . .
Peakmap Creation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
The Frequency Hough Transform . . . . . . . . . . . . . . . . . . . . . . . .
Theoretical Sensitivity Estimate . . . . . . . . . . . . . . . . . . . . . . . . .

75
75
78
81
83

THE GENERALIZED FREQUENCY HOUGH . . . . . . . . . . . . . . . . . . . . .

86

9.1
9.2
9.3

86
86
88

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
The Coordinate Transformation . . . . . . . . . . . . . . . . . . . . . . . . .
Detection Statistic and Candidates . . . . . . . . . . . . . . . . . . . . . . .

5

9.4
9.5

9.6
9.7
10

Follow-up . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
Grids in the Parameter Space . . . . . . . . . . . . . . . . . . . . . . . . . . 93
9.5.1 Grid on x0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
9.5.2 Grid on k . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
9.5.3 Grid on n . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
Theoretical Estimate of Sensitivity . . . . . . . . . . . . . . . . . . . . . . . . 97
Empirical Estimate of Sensitivity . . . . . . . . . . . . . . . . . . . . . . . . 101

SEARCHING FOR A REMNANT OF GW170817 . . . . . . . . . . . . . . . . . . . 106
10.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
10.2 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
10.3 Parameter Space Explored . . . . . . . . . . . . . . . . . . . . . . . . . . . .
10.4 Candidate Selection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
10.5 Follow-up Procedure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
10.6 Upper Limits Procedure . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
10.7 Prospects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

11

TWO-DIMENSIONAL FAST FOURIER TRANSFORM LOW PASS FILTER . . . . . 118
11.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
11.2 Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
11.3 Preliminary Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
11.4 Outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

12

118
118
119
121

MACHINE LEARNING ALGORITHMS . . . . . . . . . . . . . . . . . . . . . . . . 124
12.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
12.2 Artificial Neural Networks . . . . . . . . . . . . . . . . . . . . . . . . . . . .
12.3 Convolution Neural Networks . . . . . . . . . . . . . . . . . . . . . . . . . .
12.4 Random Forests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
12.4.1 Constructing a Decision Tree . . . . . . . . . . . . . . . . . . . . . . .
12.4.2 Constructing a Forest with the Random Subspaces Method . . . . . . .
12.4.3 The Splitting Function: Central Axis Projection . . . . . . . . . . . . .
12.4.4 Voting and Predicting . . . . . . . . . . . . . . . . . . . . . . . . . .

13

106
106
109
112
113
114
115

124
125
128
130
132
132
133
134

CHARACTERIZATION OF CONVOLUTIONAL NEURAL NETWORKS AND RANDOM FORESTS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135
13.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
13.2 Construction of Input Images: Reduced Resolution Time/Frequency Maps . . .
13.3 Training Efficiency . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
13.4 Using a Forest of Neural Networks . . . . . . . . . . . . . . . . . . . . . . .
13.5 Robustness Towards Signals Not Trained on . . . . . . . . . . . . . . . . . .
13.6 Impact of Threshold Choice on False Alarm Probability . . . . . . . . . . . . .
13.7 Error Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
13.8 Comparison to Generalized Frequency Hough . . . . . . . . . . . . . . . . . .

6

135
136
136
137
138
140
141
142

13.9 Preliminary Characterization of Random Forests . . . . . . . . . . . . . . . . 144
14

MACHINE LEARNING-BASED SEARCH FOR A REMNANT OF GW170817 . . . . 146
14.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
14.2 Search Design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
14.3 Parameter Space Explored . . . . . . . . . . . . . . . . . . . . . . . . . . . .
14.4 Description of Search Results . . . . . . . . . . . . . . . . . . . . . . . . . .
14.5 Upper Limits . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

15

ALL-SKY SEARCH FOR LONG GRAVITATIONAL WAVE TRANSIENTS IN O3DATA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152
15.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
15.2 Search Setup and Parameter Space Explored . . . . . . . . . . . . . . . . . .
15.3 Search Progress and Preliminary Results . . . . . . . . . . . . . . . . . . . .
15.4 Upper Limits . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
15.5 Prospects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

16

146
146
148
148
150

152
153
153
154
155

CONCLUSIONS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 157

BIOGRAPHICAL SKETCH . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 164

7

LIST OF FIGURES
Figurepage

2-1 Gravitational wave polarizations . . . . . . . . . . . . . . . . . . . . . . . . . . . .

23

4-1 Equation of state constraints from GW170817 . . . . . . . . . . . . . . . . . . . .

45

5-1 Interferometer schematic. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

47

6-1 Changing TF F T and skip time in parameter space . . . . . . . . . . . . . . . . . .

62

6-2 Changing TF F T and skip time in as a function of braking index . . . . . . . . . . .

62

6-3 Strain curves current and future detectors . . . . . . . . . . . . . . . . . . . . . .

66

6-4 Effective reach of a search . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

67

7-1 Shapiro delay diagram . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

74

8-1 Scheme of the Frequency Hough hierarchical method. . . . . . . . . . . . . . . . .

76

8-2 Auto-regressive spectrum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

78

8-3 Peakmap in all-sky Frequency Hough . . . . . . . . . . . . . . . . . . . . . . . . .

81

8-4 Hough map . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

82

9-1 Peakmap before and after transformation . . . . . . . . . . . . . . . . . . . . . . .

87

9-2 Power-law Hough map . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

88

9-3 How to increase TF F T during the follow-up . . . . . . . . . . . . . . . . . . . . . .

91

9-4 Peakmaps before and after phase evolution correction in the follow-up procedure . .

92

9-5 Display of enhancement of critical ratio in the follow-up . . . . . . . . . . . . . . .

93

9-6 Characterizing step in the grid over the spindown-related parameter k . . . . . . . .

95

9-7 Number of points in k grid as a function of n . . . . . . . . . . . . . . . . . . . .

96

9-8 Number of points in grid on n as a function of n . . . . . . . . . . . . . . . . . . .

98

9-9 Number of points in the grid on n as a function of TF F T

99. . . . . . . . . . . . . .

9-10 Theoretical minimum amplitude signal detectable as a function of n . . . . . . . . . 102
9-11 Minimum detectable amplitude as a function of frequency . . . . . . . . . . . . . . 103
9-12 Efficiency vs. amplitude for n = 5, Tobs = 5000 s . . . . . . . . . . . . . . . . . . . 104
9-13 Generalized FrequencyHough summary . . . . . . . . . . . . . . . . . . . . . . . . 105

8

10-1 Possible remnants of the merger . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
10-2 Impact of n on spindown behavior of power-law signals . . . . . . . . . . . . . . . . 110
10-3 Parameter space of GW170817 remnant search using Generalized Frequency Hough

111

10-4 Parameter space to which we are most sensitive at merger time in GW170817 search 112
10-5 Upper limits on a remnant of GW170817 in each search configuration . . . . . . . . 116
10-6 Best distance reach upper limits from four pipelines on a remnant of GW170817 . . 117
11-1 Time/frequency map and transformed version in Fourier space . . . . . . . . . . . . 120
11-2 Time/frequency map before and after 2D FFT filter . . . . . . . . . . . . . . . . . 121
11-3 Quantifying enhancement with 2D FFT filter . . . . . . . . . . . . . . . . . . . . . 122
11-4 2D FFT filter applied to a signal and glitch . . . . . . . . . . . . . . . . . . . . . . 122
12-1 Example of a time/frequency map with injected signal . . . . . . . . . . . . . . . . 126
12-2 Artificial neural network diagram . . . . . . . . . . . . . . . . . . . . . . . . . . . 128
12-3 Random forest and decision tree example . . . . . . . . . . . . . . . . . . . . . . . 131
12-4 Central axis projection example . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134
13-1 Resolution reduction example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137
13-2 CNN efficiency versus number of training injections. . . . . . . . . . . . . . . . . 138

13-3 False alarm probability study as a function of number of neural networks used

. . . 139

13-4 Robustness sensitivity curves for CNNs . . . . . . . . . . . . . . . . . . . . . . . . 140
13-5 False alarm probability as a function of detection statistic . . . . . . . . . . . . . . 141
13-6 Spread in error on false alarm probability and efficiency obtained experimentally . . . 142
13-7 Comparison of neural networks with the Generalized Frequency Hough . . . . . . . 143
13-8 Random forest efficiency versus number of trees in the forest . . . . . . . . . . . . 145
14-1 Machine learning search parameter space . . . . . . . . . . . . . . . . . . . . . . . 149
14-2 Distance upper limits machine learning search . . . . . . . . . . . . . . . . . . . . 151
15-1 Candidates in O3 C00 Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155
15-2 Upper limits in O3 C00 analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156

9

Abstract of Dissertation Presented to the Graduate Schoolof the University of Florida in Partial Fulfillment of theRequirements for the Degree of Doctor of Philosophy
USING MACHINE LEARNING AND THE HOUGH TRANSFORM TO SEARCH FORGRAVITATIONAL WAVES DUE TO R-MODE EMISSION BY ISOLATED NEUTRON STARS
By
Andrew L. Miller
December 2019
Chair: Bernard F. Whiting
Major: Physics
We present two new methods to search for gravitational waves (GWs) from isolated
neutron stars lasting O(hours − days), and the results of three searches of real LIGO/Virgo
data. One method is a generalization of a method to detect continuous waves (CWs) thatlooks for signals that follow a specific model; the other uses machine learning to detectsignals that deviate from this model. In two of our searches, we tried to find a remnant of
the first-ever binary neutron star (BNS) merger, GW170817; in the third search, we looked
everywhere in the sky for any kind of long duration transient signal in LIGO’s third observing
run. For the first method developed, the Generalized Frequency Hough, we present empiricaland theoretical estimates of its sensitivity, an extensive follow-up procedure of GW candidates,
the parameter space to which it is sensitive, and computational constraints when using it. For
the second method, an application of convolutional neural networks, we demonstrate thatits sensitivity is comparable to that of first method but is much faster to apply and can findsignals that the first method cannot. Estimates of rates and sensitivities in future detectors, aswell as implications of a GW detection from a neutron star, are also discussed. Additionally we
explore the prospects of using a two-dimensional Fast Fourier Transform to improve the qualityof time/frequency images to be used by either method.
Presentiamo due nuovi metodi per la ricerca di onde gravitazionali da stelle di neutroni
isolate, di durata dell’ ordine di ore o giorni.Uno dei due metodi è adatto alle ricerca di segnaliche seguano una legge di potenza, l’altro, basato su procedure di machine learning, può invece

10cercare segnali non modellizzati. Le sensibilità dei due metodi vengono qui confrontate e si
mostrerà come siano entrambi adatti alla ricerca degli stessi segnali. Entrambi i metodi sono
stati usati in una ricerca reale, qui presentata, fatta al tempo dell’ evento GW170817, percercare un possibile segnale dovuto all’ oggetto rimasto dopo la fusione delle due stelle dineutroni. Notiamo che questi metodi possono entrambi essere utilizzati anche per la ricerca disegnali prodotti da esplosione di supernova. In questo lavoro presentiamo anche stime del tasso
di fusione di stelle di neutroni e ne discutiamo le prospettive per i rivelatori gravitazionali diprossima generazione. Discuteremo inoltre le implicazioni della scoperta di onde gravitazionali
da stelle di neutroni. Infine presenteremo una procedura, basata su trasformate di Fourier
a due dimensioni, che abbiamo studiato nell’ottica di migliorare la qualità delle immagini
tempo/frequenza, che sono utilizzate come dato di ingresso da entrambi i metodi

11

CHAPTER 1
INTRODUCTION
In 1915 Einstein proposed a new formulation of gravity, General relativity (GR), to address
issues that Newtonian gravity could not explain, e.g. the perihelion of Mercury, and to makethe theory of gravity compatible with relativity. In GR, there is no longer a force due to gravity;
rather, mass curves the space and time around it, and this curvature tells other matter how to
move. Einstein’s insights have provided the basis for the field of GW physics, as well as for agreat understanding of how the universe works. GR has surpassed all tests in the weak-fieldand now parts of the strong-field regime after observations of binary black hole (BBH) and
BNS mergers by LIGO/Virgo (1 , 2 ). GWs are a direct consequence of Einstein’s theory ofgravity, and can be thought of as very tiny ripples in space-time that propagate from theirsource in all directions. This dissertation focused on developing methods to pull out thesesignals deeply embedded in the noise of our detector.
In Chapter 2, we detail the theory behind GWs; in Chapter 3, we explain the physics ofneutron stars and begin to apply the theory of GWs to these objects. The past few years haveseen significant progress in GW physics and astronomy, so we discuss some implications of thetwo most popular discoveries in Chapter 4. We describe our detectors and the various typesof noise we encounter in Chapter 5. Chapter 6 deals with possible implications of a detectionof an isolated neutron star for the equation of state and the rates of neutron star mergers.
Afterwards in Chapter 7, techniques in data analysis and corrections to the data are elucidated.
We explain a method to detect continuous gravitational wave signals from isolated neutronstars in Chapter 8 , and show its generalization to power-law signals in Chapter 9. We search
for a remnant of binary neutron star merger GW170817 in Chapter 10 using the Generalized
Frequency Hough. A technique used to enhance the quality of images is shown in Chapter 11.
We then go over in detail different types of machine learning algorithms in Chapter 12 andcharacterize their sensitivities and robustness in Chapter 13. In Chapter 14 we use the machinelearning algorithms in a proof-of-concept search for the same remnant of GW170817, but we12

cover a different parameter space than our previous search. Afterwards in Chapter 15, we
present a preliminary analysis on the data from LIGO’s third observing run, the C00 calibration.
Finally we draw some conclusions and make plans for future work in Chapter 16.

13

CHAPTER 2
THE ORIGIN OF GRAVITATIONAL WAVES
2.1

Introduction

This chapter focuses on developing the framework with which we analyze and think
about gravitational waves from relativistic objects in the universe. We discuss in Section 2.2
the development of Einstein’s field equations (EFEs). Then, in Section 2.3 we describe how
GWs result naturally from linearizing the EFEs, as well as the energetics, polarizations, anddetectability of GWs. Finally in Section 2.4, we categorize the four major sources of GWs. This
chapter is heavily based on the information presented in (3 ).
2.2

Einstein Field Equations

To formulate a theory of gravity, one needs to have an equation of motion, that is, amodel for the dynamics of massive and massless particles. These equations of motion are also
called field equations, when treating gravity as a field. The key insight that Einstein madewas that the curvature of space-time was related to the mass/density, meaning the equations
should be of this form (3 ):Kµν = κTµν

(2–1)

where κ = 8πG(fixed from requiring that EFE reproduce Newton’s equations for smallc4velocities), G is the gravitational constant, c is the speed of light, Tµν is the energy-momentum
tensor, and Kµν is the tensor that is some nonlinear combination of the space and time
coordinates themselves and their derivatives. We describe space and time in the so-called
metric coefficients gµν , and the derivatives of the metric are expressed in terms of the Ricci
tensor Rµν , and the Ricci scalar R. The metric describes the bending of space and time due
to a particular mass distribution. We define the Ricci scalar and tensor as derivatives of themetric coefficients:

14Rµν = ∂α Γαµν − ∂ν Γαµα + Γαµν Γβαβ − Γβµα Γανβ ,

(2–2)

1
Γσµν = g σα (∂µ gνα + ∂ν gαµ − ∂α gµν )
2

(2–3)

where:

and:
∂
.
∂xα

(2–4)

R = g µν Rµν .

(2–5)

∂α =
The Ricci scalar is:

Tµν has these components:• T00 : energy/volume• T0i : c × momentum/volume of source particles in the xi -direction• Tij : internal stresses of source, related to pressure
The form of the equation has certain requirements:
1. The slow-velocity limit (Newtonian) implies that the derivatives of gµν should have nohigher than linear terms in second order derivatives of gµν .
2. Tµν is symmetric, so Kµν should be as well.
Under these conditions, the field equations take the form (3 ):Kµν = aRµν + bRgµν + λgµν

(2–6)

λ can be chosen to be 0 or not depending on the cosmological model. When we allow λ = 0,
we fix that Kµν must be linear in second-order derivatives of gµν .
In Newtonian physics, we have conservation of energy and momentum for a perfect fluid.
The same is true in general relativity as long as we can neglect the gravitational field. When15

we cannot, as in this case when deriving the EFE, the “conservation” equation represents theequation of motion of the matter in the presence of a gravitational field. This conservation
condition is expressed by:∇µ T µν = 0(2–7)

Equation 2–7 places tight constraints on the form of EFE. By applying this equation andsome mathematical identities that are coordinate-independent, we can determine that b = 1/2.
When taking the limit as velocities become small (where Newton’s laws apply), it can be seenthat a = −1, so we have:
18πG− Kµν = Rµν − gµν R = − 4 Tµν .
2
c
2.3

(2–8)

Gravitational Wave Solutions in the Linearized Regime

To see where GWs appear in GR, we begin with the following two claims: (1) assume weare very far from the source such that the overall metric can be described by the metric ofspecial relativity ηµν , and (2) a GW will cause a small perturbation to this metric hµν . Thiscan be written as.

gµν = ηµν + hµν(2–9)

where |hµν |  1.
Substituting this metric into the EFE involves a lot of algebra, the highlights of which are
the following:
1. All terms O(h2µν ) and higher are neglected.
2. A trace-reversed perturbation h̄µν = hµν − 21 ηµν h is actually used to simply the resulting
algebra.

16

3. A gauge condition resulting from the invariance of hµν under an infinitesimal transformation of the coordinates (xµ → xµ + ξ µ ) is applied to simply the final equation, in an
analogous way as is done in electromagnetism with the vector potential Aµ .
• Note: gauge conditions are allowed because the transformed perturbation h̄µν
and hµν both satisfy the EFE, so nothing about the dynamics of the system haschanged. It’s similar to adding a constant to an equation for a particle’s position inclassical mechanics: it doesn’t affect the acceleration terms.
At this point, we arrive at an equation for h̄µν that looks similar to a wave equation:

2 h̄µν = −16πG µν
T .
c4(2–10)

where 2 = ∂µ ∂ µ . This equation holds so long as the gauge condition is satisfied:∂µ h̄µν = 0(2–11)

There is however still some residual gauge freedom, so we also fix:
2 ξ µ = 0.

(2–12)

Solving Equation 2–10 in the absence of Tµν gives wave solutions:ρh̄µν = Aµν eikρ x ,(2–13)

where Aµν is a constant/complex/symmetric tensor. In order for Equation 2–10 to be satisfied
by Equation 2–13, kρ must be null, so kρ k ρ = 0 implying the waves travel at the speed of light.
If we write kρ = ( ωc , ~k), we see immediately that ω = ck, the identical dispersion relation forlight.
One can also solve Equation 2–10 with T µν by using a Green’s function, and it can beshown in this way that h̄µν is related to the volume integral over T µν at the retarded time:174Gh̄ (ct, ~x) = 4ckl

ZT kl (ct − |~x − ~y |, ~y ) 3
d ~y .
|~x − ~y |(2–14)

Here, ~y is the position of the source, and ~x is the field point at which one measures h̄µν .
We have been able to eliminate the T 00 and T 0l terms because of the conservation ofenergy/momentum tensor.
Because we are far from the source of GWs, a Taylor expansion of the |~x − ~y | piece inEquation 2–14 can be done:|~x − ~y | = r − ~y · n̂ + ...

(2–15)

where r is the magnitude of the distance between us and the star, and n̂ is the normal vector
pointing from us to the source. We then take only the leading term:
4G 1h̄ (ct, ~x) = 4
c r
kl

Z

T kl (ct − r + ~y · n̂, ~y )d3 ~y .

(2–16)

Moreover, we can also expand T kl at r = ∞ where each term in the expansion corresponds
to a certain multipole moment. The expansion of T kl corresponds to an expansion in powersof v/c. By neglecting higher order moments, we are implicitly assuming that the speed of thesource is small, and that the metric is not at all influenced by the speed of source.
We can see this more clearly if we take the Fourier Transform of T kl :

T

kl

Z
=

d4 k kl ~ −iω(t−r/c+~y·n̂/c)+i~k·~yT (ω, k)e.
(2π)4(2–17)

Upon expanding the exponential piece, we have increasing powers of v/c:−iω(t−r/c+~
y ·n̂/c)+i~k·~
y

e

=e

−iω(t−r/c)
ω i i 1  ω 2 i j i jy y n n + ...
1−i y n +−ic2cThis expansion in the Fourier space is equivalent to expanding T kl in position space:18(2–18)

T kl = T kl + y i ni ∂ 0 T kl +1 i j i j 0 2 kl
y y n n (∂ ) T + ...
2c2(2–19)In this way, different spatial integrals over the various terms in Equation 2–18 represent
different moments. By expanding T kl , we arrive at:

4Gh̄ (ct, ~x) = 4
c
kl




1 kl 1kl,m
S + nm Ṡ
+ ... ,rc

(2–20)

where:

Zij

S (t) =
S ij,k (t) =Z

d3 xT ij (t, ~x)(2–21)

d3 xT ij (t, ~x)xk .

(2–22)

It is better to write the T ij components of the energy momentum tensor in terms of the energy
density T 00 and the momenta densities T 0i . Defining the energies:Z1
M = 2 d3 xT 00 (t, ~x)cZ1i
M = 2 d3 T 00 (t, ~x)xi
c
Z
1
ij
M
= 2 d3 T 00 (t, ~x)xi xjc(2–23)(2–24)(2–25)...
and the momenta:

P

i

P i,j

Z
1
=
d3 xT 0i (t, ~x)cZ1=d3 xT 0i (t, ~x)xj
c

...

19

(2–26)(2–27)

and using energy/momentum conservation laws, we can find:

1 ij
M̈
2

1 ... ijk 1  i,jk
j,ik
k,ij
+
=
P̈
+
P̈
−
2
P̈
.
M
6
3

S ij =(2–28)

Ṡ ij,k(2–29)

Equation 2–28 is the time-derivative of the quadrupole moment; the first and second terms inEquation 2–29 are the octupole moment and the current quadrupole moment, respectively.
The leading order term to the emission of gravitational radiation is the the second time
derivative of the quadrupole term. Neglecting all higher order moments, the equation for theamplitude tensor of GW emission reads:
2Gh̄ij (ct, ~x) = − 6 Q̈ij (ct − r),cr(2–30)

where Qij = M ij − δ ij Mkk is the quadrupole moment tensor. We will analyze these moments ina later chapter when we discuss GWs from isolated neutron stars and from r-modes.
Note that this approximation breaks down when the source’s motion is relativistic because
curvature and velocity are no longer independent; however, this is not so much a problem
for GWs that come from an isolated neutron star, since we don’t use post-Newtonian theory,
an expansion of the energy due to gravitational waves in powers of velocity of the source, toconstruct the waveform.
2.3.1

Gravitational Wave Polarizations

According to Equation 2–13, the Aµν term has 16 terms, however with an appropriate
choice of gauge, there will only be two independent terms. Aµν is symmetric (Aµν = Aνµ )because T µν is, which means 16 independent terms become 10 independent ones. But, Aµν
also satisfies the gauge condition: Aµν kµ = 0, so this reduces the number of degrees of
freedom to 6.
From here, it is best to work with a plane wave propagating in the z direction, sokρ = (ck, 0, 0, k). Solving Equation 2–11 for ξ ν with plane wave solutions, and doing some20

algebra, all terms vanish except for 2 independent terms in Aµν : A11 , A21 , A12 , and A22 . To
be traceless, A11 = −A22 = A+ and symmetric, A12 = A21 = Ax . We are now in thetransverse-traceless (TT) gauge, a consequence of removing gauge degrees of freedom fromAµν . Of course, the transformation to the TT gauge works for a wave propagating in anarbitrary direction as well.
2.3.2

Effect of Gravitational Waves on Point Particles

To understand how a GW with two polarizations will affect the motion of a particle,
we must consider a ring of particles through which it passes. We begin with the geodesic
equation:d 2 xµdxν dxρ+ Γµνρ
= 0,
dτ
dτ dτ

(2–31)

where xµ is the position four-vector and τ is the proper time (or in general any parameter bywhich one can parameterize the path of a particle).
Equation 2–31 is derived in a metric-independent way by using the variational principle,
which states that freely falling test particles will follow the path in a curved space-time suchthat the action is extremized. In classical mechanics, this corresponds to the fact that a
particle will always follow the shortest distance between two points.
If we only considered one point particle, we would find that the particle stays at rest
regardless of the passing GW, because in the TT gauge, the coordinates are attached toindividual particles. Therefore, we need to understand the difference in distance between two or
more particles as the GW passes.
Let us imagine that we have two particles on the x-axis, at coordinate distance L/2
from the origin. If we calculate their distance as a GW perturbation comes through from thez-direction, kρ = (ω, 0, 0, ω) with a single polarization:21

ds2 = gµν dxµ dxν

(2–32)

L20 = (η11 + h11 )L2(2–33)L0= (1 − h+ cos ωt)L2

1≈ L 1 − h+ cos ωt .
2

(2–34)(2–35)

Here we have simply used Equation 2–13 to eliminate the A+ term in favor of h+ and taken
the real part of 2–13.
If we consider two particles on the y-axis separated by L, we similarly get:
1L0 ≈ L 1 + h+ cos ωt .
2

(2–36)The difference in sign of the h+ term results in an out-of-phase oscillation on the y-axis withrespect to the x-axis.
h× is related to h+ by a rotation of π/4 radians, so the motion of the test particles will
also differ by this amount.
Based on Equation 2–35, we can see that the proper distance between two particles
changes by an extremely small small distance 12 h+ and oscillates harmonically as the GW
passes through. We show the effect of a GW on a ring of point particles in Figure 2-1.
2.4

Sources of Gravitational Waves

We will describe here the four major sources of gravitational radiation to which groundbased interferometers are sensitive.
2.4.1

Compact Binaries

The only sources of GWs that have been detected so far are compact binaries, composedof either two black holes, two neutron stars or one black hole and one neutron star. In a binary
system, the two objects orbit each other at a frequency Ω, and due to the emission of GWs,
they lose orbital energy and angular momentum, causing them to get closer to each other22

Figure 2-1: The effect of a plus/cross polarized GW on a ring of point particles. The wave istraveling perpendicular to the page.
and orbit more quickly. Eventually, the two objects will collide in a huge explosion of GWs,
forming a remnant whose properties depend on the masses, spins, tidal deformabilities, etc. ofthe two original objects. The GW emitted from a compact binary is modeled in three distinct
phases: (1) the inspiral, (2) the merger and (3) the ringdown. In the inspiral phase, the phaseevolution of the binary is modeled using post-Newtonian expansion, which expresses the phaseas a function of powers of v/c.
Next, the merger phase is modeled by numerically solving the equations of generalrelativity for a binary, since the two-body problem is not analytically solvable.
For black holes only, the ringdown is modeled using perturbation theory, and has ananalytic solution: quasi-normal modes. The black hole essentially “rings” like a bell as it emits
the remainder of its GW energy as it settles down into a stable state and is geometricallysymmetric, so it no longer emits GWs.
Binaries are the most understood source of GWs because their waveforms can be exactly
described in these three phases, though their rates were uncertain. When we know the23waveform, we can use the most sensitive detection filtering technique to pull out this signal
from the data called matched filtering, discussed later.
2.4.2

Core Collapse Supernovae Explosions

It is possible that non-axisymmetric explosions of stars within our galaxy will emit
detectable GWs. The rate of such events around 1-10 per century (4 ), and the last known
event has occurred about 100 years ago (5 ), so there is anticipation that another will occur
soon.
Signals from supernovae are much more complex than from binaries, since neutrino
interactions, hydrodynamics, general relativity, shocks, and fluid instabilities all play a role indetermining the time/frequency behavior of a GW. In fact, the GWs do not come just from an
non-spherical structure, but from an asymmetric emission of neutrinos, from a non-sphericalstructure, or forces on the fluid caused by the pressure of the neutrinos. Many models of theseinteractions exist (6 , 7 ), though the searches for supernovae typically do not assume a signal
morphology (8 , 9 ). To claim a detection, these searches assign a loudness criteria to extremely
short duration events, and use all available detectors to look for coincident events.
As an example of a mechanism to produce anisotropic neutrino emission, convective
motion occurs due to the fact that negative lepton and entropy gradients are being formedin the outer parts of the proto-neutron star (6 ). This situation may give rise to dynamical
turnovers of material heated by neutrinos, which would cause large deviations from aspherically symmetric explosion (10 ). Specifically near the neutrino-sphere, neutrinos can beconvectively transported, leading to anisotropic emission, caused by fluctuations in temperatureand density in this region.
2.4.3

Stochastic BackgroundAn incoherent superposition of all GW events in the universe will give rise to a backgroundof events, characterized mostly by BBH mergers in the LIGO frequency bands. We also expect
a stochastic GW background from the waves that were generated in the early universe as well,analogous to the Cosmic Microwave Background we observe electromagnetically. Possible

24sources, detectable with future generation instruments, include waves from inflation, phasetransitions, and cosmic strings. Each of these sources provides a unique contribution to thestochastic GW background (11 ).
2.4.4

Neutron Stars

Low mass x-ray binaries (LMXBs) are binary systems where one object is a neutron staror black hole that accretes matter from its companion, e.g. Scorpius X-1 (12 ). LMXBs canemit continuous (monochromatic) gravitational waves because this accretion should have spunup the rate of rotation of the neutron star to be much higher (> 1000 Hz) (13 ) than what
we currently observe, which is O(100) Hz (14 ). Something is limiting the spin-up of thesesystems: it is therefore possible that the system radiates away angular momentum in the formof gravitational waves that balances the accretion torque (15 , 16 ). The amplitude of the GW
that balances these torques increases with the amount of observed X-ray flux independent of
the distance the source is from us, so these kinds of binaries are very promising sources of GWs(17 ).
Isolated neutron stars are expected to emit GWs due to a small deformation on thesurface of the star, or due to some fundamental oscillation mode (r-modes, for example,
discussed later). A detection of a GW from a neutron star will give us direct insight into thebulk motion of the star, something that is impossible with electromagnetic radiation frompulsars. Moreover, this deformation or oscillation mode is related to equation of state (EoS)
parameters, which would allow us to at least rule out some EoS. The rest of this thesis will be
devoted to detecting GWs from isolated neutron stars.

25

CHAPTER 3
THE ASTROPHYSICS OF NEUTRON STARS
3.1

Introduction

Neutron stars are formed in stellar collapses when the progenitor star’s mass is between
8-25 Msun (18 ). Despite our ability to observe neutron stars, we still do not know a lot about
the physics that governs their behavior and their emissions of GWs. We describe some defining
physical aspects of neutron stars, including their magnetic fields and fast rotation rates, inSection 3.2. Then we analyze the types of GWs we expect from deformed neutron stars in
Section 3.3. In Section 3.4, we describe a likely source of GWs from isolated, old neutronstars: continuous waves (CWs). Next, we describe a type of GW that has recently received
a lot of attention in the LIGO/Virgo communities in Section 3.5: transient continuous waves(tCWs). We analyze additional possible GWs due to perturbations throughout the neutron star,so-called r-modes, in Section 3.6. We conclude this chapter by discussing the implications of a
GW detection for the neutron star EoS in Section 3.8.
3.2

Neutron Stars

There have been only about 2000 neutron stars detected using electromagnetic radiationin our galaxy, but O(108 ) are expected to exist (19 ). The neutron stars we have detected arecommonly called pulsars, because of how regularly they shine light at us. They can be usedas very precise clocks, and are in fact used as “arms” in GW experiments (20 ). The radiusof a neutron star is O(10 km), but its mass is on the order of a solar mass, implying a verylarge density. The fastest spinning neutron star we have seen rotates at 716 Hz (21 ), which,
together with high densities, make these compact objects a relativistic laboratory. They also
have magnetic fields on O(1010 ) T.
The neutron star is unsurprisingly composed mostly of neutrons: during stellar collapse
the force of gravity is so strong that the outward pressure of electrons and protons repelling off
each other is not enough to sustain the star (as it would be in the case of white dwarfs). The
neutron star is supported by the pressure between all the neutrons that have been created due

26to gravity forcing electrons and protons into such close proximity that they combine and form
neutrons (22 ). The pressure that the neutrons exert on each other is enough to stop further
collapse of the star to a black hole.
In the following subsections, we describe how the magnetic field strength and the
maximum rotational frequency of a neutron star can be estimated.
3.2.1

Magnetic Field Strength

The strength of the magnetic field in a neutron star can be understood intuitively inthe following way: astronomers can measure stellar magnetic fields by looking at the spectrathat come from a star. Typically, the spectra show dips in energy at certain frequencies
that correspond to certain elements being present in the star’s atmosphere. If the star has amagnetic field, an absorption line is split into evenly spaced lines, which is due to the Zeeman
effect: a magnetic field splits the energy levels in an atom. By measuring the spacing in the
spectra, one can determine the magnetic field of a star. These fields are typically weak, O(0.1)
T maximum.
We can use this information to determine the magnetic field of a neutron star formed afterstellar collapse. If we imagine a star that has the radius of our sun Rs = 7 × 105 km, a massR
~ · dA~of 10 Msun , and a magnetic field of 0.1 T, we can calculate the magnetic flux φ = B
through a ring that has the radius of the star before and after the collapse, and equate the two(since energy is conserved):
Bns = Bs
3.2.2

Rs
Rns

2

≈ 109 T.

(3–1)

Break-up Frequency

It is straightforward to compute the rotation frequency at which a neutron star is nolonger stable, the so-called break-up frequency fbu . If we assume that a point particle isrotating on the surface of the neutron star, then the only two forces acting on the particle (inthe rotating frame) are gravity and the centrifugal force. Equating the two gives:

27

1
fbu =
2π

r

GM
≈ 1833 Hz
R3



M
Msun

1/2 

R
10 km

−3/2
,

(3–2)

where M is the mass of the star and R is its radius. Moreover, we know that a neutron star
has a maximum mass of around 3Msun (otherwise, the neutron degeneracy pressure cannot
sustain the inward force due to gravity and collapse to a black hole would occur), so we canplace an upper bound on the maximum rotation frequency of the neutron star, ∼ 3000 Hz.
This is a very simple calculation: clearly the material of which a neutron star is composed
will also affect fbu , meaning that the EoS will play a role. However, the calculation is meant
to give an idea of what kind of GWs to expect from a rotating neutron star: recalling that theGW frequency is twice the rotational frequency, we expect GWs with frequencies no higher
than ∼ 6000 Hz.
This calculation raises another question: if the break-up frequency is ∼ 3000 Hz, why
have we not seen any pulsars that rotate at this frequency? The answer is not known, butsome combination of the emission of GWs and electromagnetic radiation may be responsible
for this restriction on the maximum spin (23 ).
Additionally, if the neutron star is accreting matter from a companion star, the rate at
which it is accumulating mass will impact the maximum frequency allowed, since there will bean additional torque applied to the neutron star, spinning it up (24 , 25 ).
3.2.3

Indirectly Measuring an Upper Limit Magnetic Field of a Neutron Star

We can calculate the power radiated due to an arbitrary object in terms of the second
time derivative of the magnetic dipole moment m̈, which is (26 ):Pmag =

dEµ0 |m̈|2
=
,
dt
6πc3

(3–3)

where µ0 is the permeability constant. Now, we assume that a neutron star acts like amagnetic dipole, and that this dipole is inclined at an angle α with respect to the central axisof the star. The dipole rotates with angular velocity Ω, meaning m = m0 eiΩt . Taking two

28

time derivatives, and recalling that the magnetic moment of a uniformly magnetized sphere ism = BR3 /µ0 :Pmag =Ω4 B 2 R6 sin2 α
.
6πc3

(3–4)

We also know the time derivative of the rotational energy of the star is given by:dProt =
dt



1 2
IΩ
2


= IΩΩ̇.

(3–5)

We can now attempt to bound the magnetic field of the neutron star by assuming that all
the rotational energy is converted into magnetic dipole radiation, so we equate Pmag and Prot ,and write Ω = 1/P and Ω̇ = − PṖ2 we see that:sB=

6πc3 I Ω̇
=R6 sin2 α Ω3

r

6πIc3
P Ṗ .
R6 sin2 α(3–6)

The magnetic field is proportional to the spin period P and its derivative Ṗ . So all one has
to do is measure how quickly the pulsar is rotating and spinning down, which is done throughmeasurements of how often we receive light pulses over time.
3.3

Gravitational Waves from Deformed, Rotating Neutron Stars

Typically a neutron star is modeled as an ellipsoid, a rotating rigid body with threesemi-axes a, b, c of different lengths. To determine the strength of GWs, we need to calculateh̄ij in Equation 2–30, and for this, we must determine the moment of inertia tensor Qij . Since
we do not expect the neutron stars to be moving quickly (since their break-up frequency ismildly relativistic), we can write
T 00 ≈ ρc2 ,(3–7)

where ρ is the density of the neutron star, which we assume to be uniform over the star volume
in these calculations (this is not true in practice). Now we need to compute the principle
moments of inertia of an ellipsoid. From classical mechanics, we know the moments of inertia
of a rigid body are:

29

ij

Z

Q =

d3 xρ(~x)(r2 δ ij − xi xj ),

(3–8)

V

where we have removed the trace for use in the quadrupole formula 2–30, and that acoordinate system can be chosen such that the moment of inertia tensor is diagonal. The
components are:Z

M 2
(b + c2 )5ZV
M 2(a + c2 )I2 =d3 xρ(x2 + z 2 ) =5ZVM(a2 + b2 )I3 =d3 xρ(x2 + y 2 ) =
5
V
I1 =

d3 xρ(y 2 + z 2 ) =(3–9)(3–10)(3–11)Let us now allow the rigid body to rotate, which we will describe the rotation as about thez-axis by a rate ω with the rotation matrix R: cos ωt sin ωt 0 



R=
 − sin ωt cos ωt 0  ,


0
0
1
From now on primed coordinates refer to the rotating frame x’,y’,z’, while unprimed refer tothe inertial frame. So Q’ij is time-independent, while Qij is time-dependent. They are related
through tensor products:

Q = RT Q’R

30(3–12)

so
Q11 = I1 cos2 ωt + I2 sin2 ωt = I1 + I2 +Q12 =I1 − I2
cos 2ωt2

I1 − I2
sin 2ωt2(3–13)(3–14)

Q22 = I1 sin2 ωt + I2 cos2 ωt = I1 + I2 −I1 − I2
cos 2ωt2(3–15)

Q33 = I3

(3–16)

Q13 = Q23 = 0.

(3–17)

From the moments of inertia in the inertial frame, we see that GWs are radiated at twice
the rotational frequency of the neutron star, ωgw = 2ω. Now, we simply need to take thetime-derivative of the appropriate components of Qij .
2
Using Equation 2–30, seeing that the moments of inertia have the same pre-factor I1 −I,
2and exchanging angular GW frequency for frequency (ω = 2πf ), we can deduce that thepre-factor is:

h0 =4π 2 G I3 
4π 2 G (I1 − I2 )
2
f
(t)
=
f (t)2 ,
4
4
c
d
c
d

(3–18)

2
where  = I1I−I
is the level of asymmetry on the star and d is the star’s distance from us.
3

The amplitude h0 is modulated depending on the inclination angle of the observer withrespect to the source of the radiation. This modulation affects h+ and h× differently.
We can now derive the power due to GW emission, which will be useful in Section 3.5
when we parameterize our analysis methods and search in terms of braking indices:... ...
dE
G= −Lgw = − 5 < Qij Qij >dt
5c(3–19)
 cos 2ωt sin 2ωt 0 
I1 − I2  sin 2ωt − cos 2ωt 0  .
Qij =

2 


0
0
1

31

(3–20)

Taking three time derivatives and substituting in I3 and : sin 2ωt − cos 2ωt 0 
...



Qij = 4ω 3 I3  
 − cos 2ωt − sin 2ωt 0  .


0
0
1

(3–21)

Plugging this into Equation 3–19 and time-averaging over the resulting sum of sin2 ωt andcos2 ωt terms:
dE
32G= Ė = − 5 ω 6 I32 2 .
dt
5c

(3–22)

Equating 3–22 with the rotational energy Erot = 12 I3 ω 2 of a neutron star, we can obtain thespindown as a function of the frequency:ω̇ = −

32G 2 5I3  ω .
5c5(3–23)

This is the power law behavior that we have developed a method to search for.
3.4

Continuous Gravitational Waves

CWs are quasi-infinite, quasi-monochromatic signals that we expect from small deformations on the surface of a quickly rotating neutron star. One can imagine a perfect sphere witha small “mountain” or “crater”: it is this deformation that causes a time-varying quadrupole
moment and hence the emission of CWs. Because the deformation is small, the signal istypically even smaller than those we expect from coalescing binaries. However, the signal
is quasi-infinite, so the greatest limitation on sensitivity is how much data we have. This isbecause by integrating for longer times, the signal-to-noise ratio increases with the square rootof the observation time. For an unchanging frequency f0 , we can see this clearly by looking
at Equation 7–13 for a continuous wave signal, but the same result applies even when thesource has a spindown or is Doppler modulated due to the earth/source relative motion. The
amplitude h(t) of a CW signal is a sinusoid h(t) ∼ cos(2πf0 t). The Fourier Transform of asinusoid is a delta function, meaning the ideal matched filter is also a delta function.

32

In CW searches, one year of data are analyzed in small chunks of time, typically 10008000 s, because of the intrinsic spindown of the source and the modulation of the signal dueto the Doppler shift. To elaborate further, a Fast Fourier Transform (FFT) length of theentire run is ideal for a sinusoidal signal, since all the power can be confined to the smallest
frequency bin possible. However, even after correcting for the Doppler modulation of the earth,
the source can still have an intrinsic spindown (as it must, because the decrease in rotationalenergy of the neutron star is what powers the GW emission). Therefore, unless one knows theintrinsic spindown a priori, one has to choose FFT lengths such that the signal will not spinout of one frequency bin during TF F T , the time of the FFT. Therefore, CW searches use FFTlengths of O(hours − days) depending on the computation power available for the search. The
all-sky Frequency Hough (FH) pipeline can only afford TF F T = 8192 s up to 128 Hz, because
the number of points in the sky increases with the square of the frequency and because thevariation due to the Doppler motion at these frequencies must remain within one frequency bin;however other searches such as Einstein-at-Home (27 ) can take day-long FFTs because theysource the huge computation load to individuals’ computers (28 ).
The model for a CW signal is:f¨
f = f0 + f˙(t − t0 ) + (t − t0 )2 + . . .
2(3–24)

where f0 is the frequency at time t0 , f˙ is the spindown and f¨ is the second-order spindown.
Since the frequency does not change very much, the amplitude of CW signals is effectively
constant, as shown in Equation 3–18. However, earth’s rotation induces a Doppler modulation,and there are source-independent relativistic effects that must be corrected for– more
information is in Section 7.3. At the detector, the signal is a time-varying function:h(t) = F+ h+ (t) + F× h× (t),where F+,× are the beam pattern functions given by (29 ):33(3–25)

F+ = sin ζ[a(t) cos 2ψ + b(t) sin 2ψ](3–26)

F× = sin ζ[b(t) cos 2ψ − a(t) sin 2ψ],

(3–27)

where ψ is the polarization angle of the GW, ζ is the angle between the interferometers’ arms,
and a(t) and b(t) are written explicitly in (29 ) but depend on the detector site location, source
position, the frequency of earth’s rotation, and a reference phase at t = 0. Additionally,h+ (t) = h0,+ cos(2πf0 t)(3–28)h× (t) = h0,× sin(2πf0 t),(3–29)

and:1 + cos2 ι
2(3–30)

h0,× = h0 cos ι.

(3–31)

h0,+ = h0

where ι is the inclination angle between us and the neutron star, and h0 is given by Equation
3–18.
3.5

Transient Continuous Gravitational Waves

Transient continuous gravitational waves (tCWs), or long duration transients, are signalsthat come from isolated neutron stars and last O(hours-days). They are “long” with respect
to burst and CBC signals, but they are “short” with respect to CWs. This portion of theparameter space only began recently due to a lack of promising sources (30 –33 ). For thesesignals we assume a power-law model for the time/frequency evolution of the GW signal:

f˙ = −kf n34(3–32)

where n is the braking index, and k is a proportionality constant that contains the physics ofthe emitted radiation, defined to be positive. This radiation is modeled as coming from therotational energy of the neutron star; hence as the energy is emitted, the neutron star spins
down. Integrating Equation 3–32, we get:f0

f (t) =

1 + k(n − 1)f0n−1 (t − t0 )

1
 n−1

(3–33)For all braking indices except n = 7 (r-modes), we model the amplitude change as a functionof time and frequency as in Equation 3–18 (29 ). This amplitude evolution as a function oftime explicitly is (34 ): 2
4π 2 GIzz  2t 1−nh(t) =
f 1+
c4
d 0
τ

(3–34)

where τ = [(n−1)kf0n−1 ]−1 is the spindown time scale,  is the ellipticity (degree of deformation)of the neutron star, and Izz is the z-component of the star’s moment of inertia.
For r-modes (n = 7), we use the following relation (35 , 36 ):rh(t) =

2048 GJM R3
αf (t)3375
c5 d(3–35)

where α is the saturation amplitude, the amplitude at which non-linear effects stop the growthof the r-mode and J = 1.635 × 10−2 is a dimensionless angular momentum for polytropic EoSmodels. R-modes will be discussed in Section 3.6.
3.6

Gravitational Waves from R-modes

A promising source of GWs from neutron stars comes from toroidal perturbations
throughout the star, so-called r-modes. These modes are analogous to Rossby waves on earth,and have a Coriolis restoring force. In Subsection 3.6.1 we delve into the physics of r-modes:
how r-modes come about from a rotating object, and the kinds of GWs they emit.

35

3.6.1

Velocity and Density Perturbations

r-mode oscillations are caused by small changes in the velocity and density of fluid matterin the star. To compute the form of these perturbations, we write the multipoles of Equation
2–20 as (37 , 38 ):∞m=l

1X X llm+ Jlm (Tµν )lm ]hµν =∂0 [Dlm Tµνr l=2 m=−l(3–36)

where Dlm is the mass quadrupole moment and Jlm is the current quadrupole moment.
We can then use this expression for hµν , coupled with an expression for the power of theGW, to write:∞∂0 Egw =

m=l

1 X X
< |∂0l+1 Dlm |2 + |∂0l+1 Jlm |2 >,32π l=2 m=−l(3–37)and knowing that we are dealing with waves of the form eiωt , it follows that:∂0 Egw = −ωrX

Nl ωi2l+1 (|δDlm |2 + |δJlm |2 ),(3–38)

l≥2

where ωr = ωi + mΩ. Here ωr is the rotational frequency of the star in the rotating frame, andωi is the rotational frequency in the inertial frame. They are simply off by a factor of mΩ, but
it is this factor that makes all the difference, as will be discussed in a following section 3.6.4.
Additionally, we can express the mass and current quadruples in terms of these velocity anddensity perturbations δv and δρ respectively (37 ):Z∗δρrl Ylm
dV
r
Z
2
l
∗B
=
rl (ρδv + δρΩ)Ylm
dV,
c l+1

δDlm =(3–39)

δJlm(3–40)

B
where the Ylm are the standard spherical harmonics and Ylm
are the magnetic multipoles,

defined in (39 ).
Writing explicitly the r-mode eigenfunction ξ:

36



ξ X
Dlm
B iωr t=0,∂φ , −Dlm ∂θ Ylm
e ,
r
sin
θ
l,m

(3–41)

the velocity (δv = ∂0 ξ) and density perturbations equations can be derived (38 , 40 ): r lYllB eiωr t
R "#r
 r l+12l∂ρlδρ = αR2 Ω2 ρ+ δψ Yl+1,l eiωt
∂P 2l + 1 l + 1 R

δv = αRΩ

3.6.2(3–42)(3–43)

Frequency Evolution of R-modes

To begin, we treat the neutron star as a rotating fluid using Euler’s equation:∂~v~ × ~v ) + ρΩ
~ × (Ω~ × ~r) = −∇P − ρ∇Φ,
ρ+ ~v · (∇~v ) + 2ρ(Ω∂t

(3–44)

where v is the velocity of the fluid, P is the pressure of the fluid and Φ is the gravitational
potential.
We now consider small toroidal displacements (no perturbation in the radial direction) of
the fluid star, as well as perturbations in velocity, density, pressure, and gravitational potentialof the star. And we decompose the small toroidal displacement vector ξ~ into harmonics of theazimuthal angle, which can be done so long as the object is azimuthally symmetric.
This vector can be written as ξ~0 ei(mφ−ωr t) , where we have decomposed this perturbation
into harmonics m with respect to the azimuthal angle, and assumed periodic time dependence.
We define the displacement vector ξ in terms of the velocity perturbation:~
δ~v = ∂0 ξ,(3–45)

where δv is the perturbation in fluid velocity. In this section δX refers to the perturbation ofquantity X. In this case we want to see how the inertial toroidal frequencies in the inertial androtating frames are related, which can be accomplished through the derivative transformation:

37

dξ~i
= ∂t ξr + ~v · ∇ξrdt1−ωi = −ωr + RΩ ∂φ ξ~rR
ωi = ωr − mΩ

(3–46)(3–47)(3–48)

Applying these perturbations and linearizing Equation 3–44, we obtain the following (37 ):δρ1∂0 δ~v + 2Ω × δ~v = −ωr2 ξ~ + 2iωr Ω × ξ~ = 2 ∇p − ∇δp − ∇δΦ.
ρ
ρ

(3–49)Taking the curl of Equation 3–49, analyzing the radial part, neglecting terms of O(Ω2 )and knowing that the radial velocity is much smaller than the horizontal velocity, we get:~ × ~v )r + 2(~vh · ∇h )Ω~ r + 2Ωr · (∇h~vh ) = 0.
∂t (∇(3–50)

Other necessary equation include the Poisson equation for the perturbed gravitational potential,
an EoS, and the conservation of mass equation.
Considering only linear order terms in Ω, so dropping ~vh · ∇h ∝ Ω3 , and assuming no
spindown of the star (∂t Ω = 0), we get in the rotating frame:id h~∇ × ~v + 2Ω = 0.
dt

(3–51)

Writing v = iωr ξ, and using toroidal eigenvectors written in terms of the quadrupole moment:
ξ
=
r



1
m
m
0, Dlm
∂φ Yl , −Dlm ∂θ Yleiωt ,sin θ(3–52)

we get:∇ × ξ = Tlml(l + 1) m
Yl ,r2

where the Tlm are function of r. We can then obtain (37 ):38(3–53)

ωr−mode =2mΩ.
l(l + 1)

(3–54)

Converting this frequency to the inertial frame, we get:ωi = ωr − mΩ =

2mΩ
− mΩ
l(l + 1)

(3–55)

Now, requiring l = m = 2, we obtain that ωi = − 43 Ω, which slightly differs from the GW
frequency being twice the rotational frequency from a small deformation derived above.
3.6.3

Current Quadrupole Radiation

It is not the mass quadrupole that radiates most strongly for r-modes, but the current
quadrupole moment, the second term in equation 2–29. The angular frequency of an r-mode,to O(Ω3 ), is given by not neglecting the scalar product ~v · ∇ term in Equation 3–50:2mΩωr−mode =l(l + 1)



R3 Ω2
1−η
,
GM

(3–56)

where the constant η ∼ 0.7. Because the Ω2 factor is weighed down by M and η, the Ω term
actually drives the frequency evolution more quickly than the Ω2 term. This implies that thespindown of a neutron star due to r-mode emission will be due to the current-quadrupole
principally.
3.6.4

CFS Instability

We study GWs from r-modes because of their ability to be driven unstable by GW
radiation. Due to the difference in frames, a mode that is prograde in the rotating frame
may be retrograde in the inertial frame, and as GWs are radiated away, angular momentum is
actually subtracted from the already negative angular momentum of the star, thus causing a‘runaway’ effect whereby the amplitude of the mode grows until some saturation mechanism
terminates the growth. The instability condition is derived in (37 ).
To better understand this instability condition, we start with the linearized Euler Equation
3–49. In the lowest order approximation, the potential can be expressed in terms of the density,

39

the pressure perturbation and the potential perturbation. Using this, we can integrate theenergy due to the perturbations (41 ):
1
E=
2

Z

3



2

d xρ|δv| +δP− δφ δρρ

(3–57)

Using expressions for the perturbations found in (35 ) and taking a time derivative, we find:

dE
=−
dt

Z

[2ηδσ ab σab + ζ|δσ|2 ]d3 x − ωi (ωi + mΩ)X

Nl ωi2l (|δDl,m |2 + |δJl,m |2 )(3–58)

l>=2

where η is the shear viscosity of the medium through which the r-modes propagate and ζ
is the bulk viscosity.
The first term tends to be negative since viscosity damps (or decreases) the energy of ther-mode over time. However, the second term may change sign, depending on if piece outside
the sum is positive or negative. So the instability condition, that is, the way to ensure that the
perturbation energy increases over time, is:ωi (ωi + mΩ) < 0(3–59)

We assume that only l = m r-modes exist in neutron stars. Using Equation 3–55, we
obtain:ωi = −(l − 1)(l + 2)Ω
l+1

(3–60)

So the instability condition Equation 3–59 is:ωi (ωi + mΩ) = −2

(l − 1)(l + 2) 2Ω,
(l + 1)2(3–61)

which is negative definite for l >= 2, hence forcing the second term in Equation 3–58 tocontribute positive energy over time, causing a growth of the r-mode.

40

3.7

Cowling Approximation

The Cowling Approximation (42 ) has been hinted at in the above r-mode formalism, but
specifically refers to neglecting the perturbations of the metric due to a source’s relativistic
velocity. Many studies have explored the impacts of this approximation in the context
of different oscillation modes of a star (43 , 44 ). Essentially, using enough nodes in the
eigenfunctions of whichever mode one considers improves the approximation. Calculations ofgravitational waves due to r-modes have also been done within this framework (35 ). While alinear relationship between the r-mode frequency and rotational frequency of the star is a good
approximation of the r-mode frequency, the coefficient of proportionality depends strongly on
the relativistic parameter GM/(c2 R), which is O(0.1) for neutron stars (45 ).
For the work presented here, the frequency of the r-mode is assumed to be produced
only from l = m = 2 mode, where higher order modes are neglected. This is certainly anover-simplification of the gravitational wave emission, especially if the r-mode grows overtime and hydrodynamic effects become important (35 ). Our methods can therefore not probe
deeply the causes of r-mode evolution, but can discover the dominant behavior.
3.8

Equation of StateA central problem of neutron star physics is determining the EoS. General relativity does
not provide this for us: we can solve Einstein’s equations to get the Tolman-OppenheimerVolkoff (TOV) equations (46 ), and then only assume an EoS p = p(ρ, ...) to solve the TOV
equations. Therefore, we need to rely on observations and assumptions about the neutron star
structure when trying to derive the EoS. GW observations of an isolated neutron star will allow
us to constrain its mass and radius, thereby giving constraints on the tidal deformability of
the neutron star, the moment of inertia, and the ellipticity. This is because we can measure
the rotational frequency and spindown of the star, and in the case of signals predominately
considered here, the braking index. Along with an estimation of the amplitude of the GW
signal, and an electromagnetic counterpart to constrain the distance, we can deduce separately
what the moment of inertia and the ellipticity of the neutron star are (47 )41

CHAPTER 4
THE DAWN OF MULTI-MESSENGER ASTRONOMY
4.1

Introduction

The era of GW astronomy began about 4 years ago with the detection of the first BBH
merger (2 ). From subsequent detections in LIGO/Virgo’s first and second observing runsof more BBHs and one BNS (48 ), physicists have been finally able to use GW to probe
interesting questions about the universe in the strong-field regime. In Section 4.2, we describe
the first ever detection of GWs by LIGO, as well as how this detection constrained speed of thegraviton. Afterwards we discuss in Section 4.3 the first BNS observed by LIGO/Virgo (1 ) thatallowed for one of the most extensive observation campaigns in scientific history (49 ).
4.2

GW150914

The first ever GW detection was made on 14 September 2015 by the two LIGO detectors(2 ). A BBH merger spanned the frequency range 35 to 250 Hz, and peaked in gravitationalwave strain at 1 × 10−21 , an extremely powerful event. From this merger, the LIGO/Virgocollaborations showed that general relativity works in the strong field regime, and that BBHs
do exist in nature, though alternative theories to explain the specific nature of this object are
possible (50 ). The implications of this merger also include higher-mass stellar black holes thanpreviously thought, meaning that models of stellar evolution need to be changed. Moreover,
we could establish some confidence in merger rates of BBHs. We discuss one interesting
implication of GW150914: the limit we can place on the speed of gravity.
In theories about massive gravity, the speed of the graviton is less than the speed of light,
meaning different modes (frequencies) travel at different speeds, given by:(~ω)2 = (~ck)2 + (mg c2 )2

42

(4–1)

where mg is the mass of the graviton and the term (mg c2 )2 is very small compared to theothers. If the waves propagate at the speed of light, we have:ω 2 = c2 k 2(4–2)

ω dω = c2 k dk

(4–3)

ω dω
= c2 .
k dk(4–4)

The speed of massive gravity, specifically the speed of a wave packet (the group velocity), is:
"

2 #
dω
c2
1 mg c2= vg =
≈c 1−
.
dk
ω/k
2
~ω

(4–5)

The duration of the signal observed in GW150914 was about ∆t = 0.1 s, from fmin = 30
Hz to fmax = 150 Hz. We can use this information to place an upper bound on the massof gravity by assuming that the two modes fmin and fmax were delayed by the length of thesignal ∆t. The time of arrival of a mode is:
"

2 #
D
1 mg c2
D
t=
≈
1+
.
vg
c
2
~ω

(4–6)

where D is the distance to the source. Writing an equation using fmin and fmax and subtracting them, defining ∆t = tmax − tmin (where the subscripts min and max refer to the frequencies),
we arrive at an upper bound on the mass of the graviton:

m2g c42 2≤ 8π ~



c∆t
D



mg c2 ≤ 1 × 10−22 eV.

2
2
fmax
fmin
2
2
fmin
− fmax


(4–7)
(4–8)

This is a much more stringent constraint than what previous experiments have given by about
two orders of magnitude (51 ).

43

4.3

GW170817

On 17 August 2017, LIGO and Virgo saw a binary neutron star merger 40 Mpc away (1 ).
What followed was one of the most extensive astronomical campaigns in history to see thismerger in every lens of the electromagnetic spectrum. We were able to see this merger in every
frequency band, which gave us insight into the remnant of the merger, as well as the physics ofneutron stars. Searches for a possible remnant of this merger will be presented in chapters 10
and 14.
4.3.1

Constraints on Difference between the Speeds of Light and Gravity

The GWs and the gamma ray burst within ∆t = 1.7 s of each other (52 ), which allows us
to constrain how much faster the GW could have been travelling than light. So we can placeconstraints on how much the speeds of light and gravity differ with the following logic:∆v = vgw − vemDD−
tgw temtem − tgw
= D
tem tgw ∆t∆v = vem
.
tgw
=

(4–9)
(4–10)(4–11)(4–12)

To zeroth order in Equation 4–6, assuming a massless gravity, tgw ≈ D/vem ; therefore:∆t∆v≈ vem .
vemD

(4–13)In the paper (53 ), the authors first assume that the GWs and EM waves were emitted at thesame time, and attribute ∆t = 1.74 ± 0.05 s entirely to faster GWs. This provides an upper,
positive bound to the difference in the two speeds, where the authors use ∆t = 1.79 s andD = 26 Mpc to be conservative. Then, the authors consider that the two waves were emitted
at times differing by more than 1.74 ± 0.05 s, allowing then the EM waves to make up some partof the travel. They consider ∆t = 10 s, and arrive at:44

− 3.74 × 10−15 ≤
4.3.2

∆v
≤ 6.67 × 10−16 .
vem(4–14)

Constraints on Equation of State

LIGO/Virgo observations of the masses of the neutron stars can begin to constrain thepossible EoS of the neutron star (1 ). In a typical matched filtering search, the phase evolutionof the signal is required to be exactly known, and this EoS information enters as one of thehigher post-Newtonian corrections to the phase evolution, in a parameter called the tidal
deformability Λ. Λ is higher for less compact (so higher radius) stars, and smaller for morecompact stars.
Assuming different priors for spins of the neutron star, LIGO/Virgo placed different
constraints on the tidal deformability, and ruled out many equations of states that produced
less compact neutron stars- see Figure 4-1. The reason that such priors are necessary is
because LIGO/Virgo is not as sensitive to the individual masses as it is to the chirp mass, since
there are degeneracies in the phase evolution of the signal.

Figure 4-1: This figure, taken from (1 ) shows the constraints on the tidal deformability parameter Λ for low and high spin priors. The color is the probability density. Predictions forΛ1,2 from different equations of state are drawn on the plot, and some equations of state thatpredict less compact starts are excluded at 50 and 90% confidence.

45

CHAPTER 5GRAVITATIONAL WAVE DETECTORS
5.1

Introduction

No GW detection would have been possible without the ingenious designs for an
instrument to measure extremely tiny perturbations in space-time. In this chapter we discuss
the basics of the detector used in GW physics. Section 5.2 demonstrates the key principles that
motivated the design of current GW detectors. Afterwards in Section 5.3, we describe manyof the sources of noise which experimentalists needed to account for when designing such asensitive instrument. Finally in Section 5.4, we explain some of the key components of the GW
detectors.
5.2

Detection Principle

LIGO and Virgo employ interferometry to detect GW signals. An interferometer is a
device that interferes two beams of light sent in perpendicular directions. A schematic ofan interferometer is shown in Figure 5-1. A (quite powerful) laser serves as input to theexperiment, which is then split into 2 perpendicular beams by a beam splitter. A simplified
scheme of the working principle is that the beams travel through a cavity formed by thetwo mirrors, bounce off the mirrors, rejoin at the beam splitter, and are forwarded to thephotodetector. What is measured is an interference pattern: if the lengths of the arms areexactly equal, and there is no noise, we expect that the two rays will constructively interfere. If
the light rays are out of phase with each other completely, we expect destructive interference
and a completely dark interference pattern.
GWs alter the position of the mirrors as discussed in Section 2.3.2. This movement isa function of the source properties that cause the GW. In essence, every black hole merger,
neutron star, supernova, etc. will change the distance between the two mirrors in a slightly
different way depending on the masses, spins, deformation, level of asymmetry, etc. A GW
is therefore like a fingerprint, allowing us to confidently identify the source of such a mirror
movement, if we are sure that noise does not interfere.

46

The GW therefore changes the arrival time of the light rays at the photodetector, and thiscauses different kinds of interference patterns, because the light rays are now out of phase.
Figure 5-1 is a much more simplified version of what LIGO/Virgo actually use because thenoise and the precision needed in this experiment required the development of many advanced
techniques. We describe many of the noise sources in LIGO/Virgo in this chapter. More detailsof the detector design are discussed in the Subsection 5.4.

Figure 5-1: This figure shows a simplified version of a LIGO/Virgo interferometer. A laser
beam is split into each of the arms, reflected off the mirrors, goes through O(100) roundtrips
in the Fabry-Perot (FP) cavity, and recombines at the beam splitter. The output is read by aphotodiode, where the phase shift between the light rays in each arm is measured.

47

5.3
5.3.1

Noise Sources

Seismic Noise

Seismic noise is due to the shaking of the earth at many frequencies, and is a very strong
source of noise. We need orders of magnitude of attenuation to have a shot of seeing GWs of
any kind. To achieve a large attenuation, LIGO/Virgo uses a system of pendula that, for eachpendulum in the cascade, reduces Sn by a factor of fr2 /f 2 , where fr is the resonance frequencyof the pendulum. To show this, we consider a simple pendulum with mass m hanging from a
rod of length l whose support point can move in time, with a known function (following thederivation in (54 )):~R(t)= X(t)~i + Y (t)~j.

(5–1)

~
where R(t)is the overall displacement of the system, and X(t) and Y (t) describe thehorizontal direction ~i and vertical direction ~j time-dependent displacements, respectively. The
position of the mass is:~ + l(sin θ~i + cos θ~j).
~r = R(5–2)

where l is the length of the pendulum and the angle θ is from the vertical to the pendulum.
These two quantities completely describe the system.
We need to write the Lagrangian and derive the equation of motion for the mass msuspended by the pendulum, which we will then use to determine the impact of different forceson the mass. We write the kinetic energy T and potential energy V as:

T =
V

1 2 1
~˙ 2 + l2 θ̇2 + 2lθ̇(Ẋ cos θ − Ẏ sin θ)
mv = m(|R|22= −m~g · ~r = −mgl cos θ − mgY.

(5–3)
(5–4)

where dots denote time derivatives. Now we can write the Lagrangian and then the equation of
motion:48

ml2 θ̈ + ml(g − Ÿ ) sin θ + mlẌ cos θ = 0.

(5–5)

If there is no displacement (X = Y = 0), we recover the equation of motion for a simple
pendulum.
We will focus on a small, periodic displacement in the x-direction: X = X0 cos(ωt), whereω is the driving frequency (which could be the GW frequency or a seismic disturbance). The
equation of motion takes the form, for small angles and Y (t) = 0:lθ̈ + gθ + Ẍ = 0(5–6)

lθ̈ + gθ = ω 2 cos ωt.

(5–7)

Because the oscillation is driven, the pendulum will, after some transients in the beginning,
oscillate at the same frequency as the support point:θ = θ0 cos ωt.

(5–8)

Plugging θ into the equation of motion, we arrive at:
θ0 =X0 ω 2,l ω02 − ω 2(5–9)

where ω02 = g/l is the oscillation frequency of the unperturbed simple pendulum.
The new position of the mass in the x-direction is: x = X + lθ, which comes from Equation
5–2 in the small-angle approximation. We now would like to compute the transfer function F ,which relates the movement of the suspension point of the pendulum to the movement of the
mass:

x
lθ
ω2
=1+
=1+ 2
X
X
ω0 − ω 2
ω02
fr2
F (ω) =
=
F
(f
)
=
.
ω02 − ω 2
fr2 − f 2

F (ω) =

49

(5–10)(5–11)

This transfer function represents an attenuation in displacement: if something attempts tomove the suspension point, so long as the rate at which it tries to move the suspension point
is not near or below the resonance frequency, the motion of the mass will be attenuated: if
f >> fr , we see that F (f ) ∼ fr2 /f 2 . If we keep the resonance frequency of the pendulum
significantly lower that the GW frequency we would like to find, we can apply N pendulums toget (fr2 /f 2 )N levels of attenuation in the motion of the test mass.
5.3.2

Thermal Noise

Thermal noise is caused by Brownian/random motion of the atoms that compose thedetector. It affects both the mirrors and the suspensions. We can express the power spectral
density of the noise using the fluctuation-dissipation theorem:

S(ω) = 4kB T Re Z(ω),(5–12)

where Z is the impedance, T is the temperature and kB is the Boltzmann constant.
This theorem states that as the force varies in time/frequency (fluctuations), there will bea loss of mechanical energy of the system (dissipation). It is expressed above in the frequency
domain.
In LIGO/Virgo, the thermal noise is calculated for the suspensions/mirrors by modelingthese systems as damped harmonic oscillators with a dampening term proportional to thevelocity of the atoms in motion:ẍ + ω02 x + γ ẋ =

F (t)
,
m

(5–13)

where γ is the damping factor, and F (t) is the driving force. The displacement spectral density
is given in terms of the fluctuation-dissipation theorem (55 ):X(ω) =

1(4kT Re Z(ω))1/2 .
ω|Z(ω)|

This term is found by casting the equation of motion in the form:

50(5–14)

F (ω) = −iωZ(ω)x(ω),(5–15)For a damped harmonic oscillator, allowing γ(ω) → φ(ω) (for generality), we can write byinspection:im 2
ω − ω02 + iω 2 φ(ω)ω

(5–16)

ktω02 φ(ω)X(ω) =mω[(ω02 − ω 2 )2 + ω04 φ2 (ω)]

(5–17)

Z(ω) =The displacement noise then is:2

where φ(ω) is the loss angle of the oscillator whose fundamental mode is ω0 , and takes
different forms for thermal noise at the mirror, suspensions, etc.
5.3.3

Shot Noise

Shot noise results from the fact that we are counting photons bouncing off each mirror
to measure the phase shift of the light, and there is an uncertainty in the number of photons
that arrive at the photodiode. To calculate the spectral density of this noise, we do thefollowing: let Nγ be the number of photons that we would like to count at the photodetectorin observation time T . The angular frequency of the laser light is ωL . Then the power due tothese photons at the photodetector is:P =1
Nγ ~ωL ,T(5–18)

where ~ is the reduced Plank’s constant. When we count a number of independent and discrete events, the probability of such photons hitting the mirrors follows a Poisson distribution,and in the limit that we count many (we count around O(109 )), the uncertainty ∆Nγ is:

∆Nγ =

p
Nγ .

51

(5–19)

So the fluctuation is the observed power of the laser beam is (replacing Nγ by ∆Nγ in theabove equation):

1p
Nγ ~ωL
T

1/2
~ω
∆P =
P
.
T

∆P =

(5–20)(5–21)

To determine P , we must write equations for the photon’s electric field as it travels througheach cavity and account for how it is affected by the passing of a GW. The electric field Ex,yalong each arm will have the form:
1
Ex,y = ∓ E0 e−iωL t+2ikL Lx,y ,2(5–22)and the total field Etot = E1 + E2 is what we will see at the photodetector. Its square modulus
is|Etot |2 = E02 sin2 [kL (Ly − Lx )] = E02 sin2 φ0(5–23)

Knowing just this and that P ∝ |Etot |2 and E0 ∼ P0 we can compute the change in the powerof the laser just due to noise:(∆P )shot =

~ωL
P0
T

1/2
| sin φ0 |(5–24)The next step is to determine (∆P )GW , which will affect the phase of Ex and Ey thatrecombine at the photodetector. To do this, we need to compute the arrival times of thephotons from each arm when affected by the GW (in the transverse-traceless gauge).
Based on the analysis in Sections 2.3.1 and 2.3.2, the perturbation hµν has only h+ andh× components, so the metric is written as:52

ds2 = gµν xµ xν = −c2 dt2 + (1 + h+ )dx2 + (1 − h+ )dy 2 + 2h× dxdy + dz 2

(5–25)

where gµν = ηµν + hµν .
Consider a GW with + polarization traveling in the z-direction:h+ (t) = h0 cos ωgw t(5–26)For a photon, ds2 = 0, and when it travels only in the x-direction, dy 2 = dz 2 = 0. The photon
is affected by the passing GW, and we can express the distance traveled from the beam splitterto the mirror Lx as (the forward trip of the light):c
Lx = c(t1 − t0 ) −2

Z t1

dt0 h+ (t0 )(5–27)

t0

where t0 is the time that the light is emitted and t1 is the time that the light reaches the endof the arm of length Lx . For the return trip, the expression is the same, just with t1 → t2 andt0 → t1 . t2 is then the time that the photon returns back to its origin.
By summing the forward and backward parts:
2Lx 1+
t2 − t0 =
c
2

Z t2

dt0 h+ (t0 ).

(5–28)

t0

In the case of no GWs, we just get the standard expression: t2 − t0 = 2Lc x . Since the correction
h+ is small, we can let t2 → t0 + 2Lc x .
2Lx 1+
t2 − t0 =c2

Z t0 +2Lx /ct0

We integrate this expression and get:53

dt0 h0 cos ωgw t.

(5–29)

2Lx h0 Lx+sinc(ωgw Lx /c) cos(ωgw (t0 + Lx /c))cc
2Lx h0 Lx=+
sinc(ωgw Lx /c)h+ (t0 + Lx /c)cct2 − t0 =(5–30)(5–31)

It is more useful to define t2 = t as the time which we observe the light at the photodetector, and then compute the corresponding t0 . This means that h+ (t0 + Lx /c) →h+ (t − 2Lx /c + Lx /c) based on Equation 5–28 in which we neglect the h+ integral:

tx0 = t −2Lx Lx−h+ (t − Lx /c) sinc(ωgw Lx /c)cc(5–32)and x → y for the other arm (with a sign flip of the third term due to the minus sign present
in the dy 2 term in the metric).
Now, we can compute the phase shift in the electric field due to the GW. The free propagation of the light in the cavity does not change the phase, and the reflection/transmission of
the light only give additional 1/2 factors in the amplitude, so:∆φgw,x = ωL tx0 =ωL Lx
h0 cos ωgw t sinc(ωgw Lx /c).
c

(5–33)

where this phase shift is between the time observed t and the time that the beam left theorigin t0 . The analysis for the y-direction is the same, except that h+ → −h+ because of themetric, so:∆φgw,x = −∆φgw,y .

(5–34)At the photodetector, we see the difference of the phase shifts in x and y:∆φmich = ∆φgw,x − ∆φgw,y = 2∆φgw,x .
The total phase shift is:54(5–35)

φtot = φ0 + ∆φmich(5–36)and the power P ∼ |Etot |2 is (in the spirit of Equation 5–23):

P = P02 sin2 (φ0 + φgw,x )=P0(1 − cos[2φ0 + ∆φmich (t)]) .
2(5–37)(5–38)

Now, we can calculate the change in power due to a GW signal by subtracting from P thecontribution due to only noise Pnoise = P0 sin2 φ0 , and by allowing cos φmich ∼ 1 andsin φmich ∼ φmich :P − Pnoise = (∆P )gw =P0∆φmich | sin 2φ0 |2(5–39)

From here, we can compute the signal-to-noise ratio of the GW signal:


S
N



∆Pgw
∆Pshot
!

1/2
P0
=
| cos φ0 |∆φmich T 1/2 .
~ωL=(5–40)(5–41)

Comparing Equation 5–41 with Equation 7–13, we can read off Sn :
Sn
=
h0



~ωL
P0

1/2

1
.
| cos φ0 |∆φmich(5–42)√
To get a better idea of how Sn behaves, let us choose φ0 = π/4 so cos φ0 = 1/ 2, and let us
analyze the simple case of ωgw Lx /c  1, so we can Taylor expand Equation 5–35 and obtain:
2LωLc(5–43)~c2.
2P0 ωL(5–44)∆φmich = h01
Sn1/2 =Ls

55

A key finding is that the shot noise gets better as more laser power is used. While true,
another noise source, radiation pressure noise, increases with increasing laser power, as we will
discuss in Subsection 5.3.4. But in Fabry-Perot (FP) cavities (see Section 5.4.2) , the types ofcavities actually used in LIGO/Virgo, the shot noise is frequency dependent. We will discuss
FP cavities in a later section.
5.3.4

Radiation Pressure Noise

Radiation pressure noise is due to the fact that a different number of discrete photons
are hitting the mirrors at any given time, causing a time-varying, stochastic force to push themirrors. The number of photons hitting the mirrors fluctuates as a function of time because
the counting statistics are Poisson-distributed- see Equation 5–19). Therefore, we need to
calculate the impact that this force F has on the sensitivity of the detector.
The photons are reflected, resulting in a doubling in the momentum transferred. Therefore
the power that is incident on the mirror is: P = 2F c. The noise results because the force on
the mirrors is always changing slightly, meaning that we will use Equation 5–24 to describe thisfluctuation:
r
∆P
~ωL P
∆F = 2=2c
c2 T(5–45)

where we have set φ0 = π/2.
We would now like to compute the power spectral density of the force. We begin bycalculating the auto-correlation of the fluctuations in the force, ∆F . The fluctuations are
independent since we are describing a random process, so we can write:< ∆F (t)∆F (t) >=

1
SF ,
2T

(5–46)

where SF is the spectral density of the force. Plugging in Equation 5–45 into Equation 5–46,
we get:56

r
1/2
SF = 22~ωL .Pc2(5–47)To relate the power spectral density of the force to the power spectral density of the displacement (or the strain), we need to specify the transfer function. Let’s assume a simple system
where this stochastic force is the only one acting on the mirror of mass M : F = M ẍ. The
transfer function will relate the power spectral density of the force to that of the displacement.
We can evaluate the transfer function by solving for the relationship between F and x in
Fourier space:F̃ (f ) = −4π 2 M f 2 x̃(f ) = T(F →x) x̃(f ).

(5–48)

since derivatives in time become multiplications in frequency domain, and T(F →x) (f ) is thetransfer function. Here f is the frequency at which the photons (or the GW) cause the mirror
to oscillate. Then:

Sx1/2 (f ) =2M (2πf )2r

2~ωL P.
c2(5–49)To relate the power spectral density of the displacement of the mirrors to that of the strain, we
can use ∆L = Lh = Tx→h h for a simple interferometer, which means the transfer function isjust L:4
Sn1/2 =
M L(2πf )2
5.4
5.4.1

r

2~ωL P.
c2(5–50)

Detector Design

Simple Michelson Interferometer

We can calculate the optimal length of a simple interferometer by maximizing the phaseshift due to GWs 5–35. Fully expressing the sinc in Equation 5–33, we find:57

φmic =2ωLcos(ωgw t) sin(ωgw Lx /c).
ωgw

(5–51)

Maximizing with respect to Lx and writing ωgw = 2πf , we arrive t:
c1= 750 km
Lx =
4f
5.4.2



100 Hzf.

(5–52)The Fabry-Perot Interferometer

Because it is not feasible to build a detector of 750 km on earth, the collaborations havethought of a clever way to effectively increase the arm length of the detector using FP cavities.
Essentially, light is reflected back and forth between the beam splitter and the end mirror manytimes before being read into the photodetector.
The difficulty in this approach lies in ensuring that the mirrors in the cavity are as close
to perfect reflectors as possible. Moreover, it is best to keep the cavities in resonance, since
this ensures that the electric field that is shot into the cavity (and therefore the number ofphotons in the cavity) remains large. This means that various beams within the cavity are all
constructively interfering with each other.
The sensitivity (the amount of phase variation we are able to see) is enhanced in a
FP cavity by a factor π2 F relative to a Michelson interferometer, where F is the Finesse, aquantity that is primarily related to how well the mirrors reflect light. This factor can becomputed by using the reflected electric field amplitude at both mirrors, originally coming fromworking through equations that relate incoming, transmitted and reflected electric fields toeach other for a certain number n of roundtrips (55 ).

58

CHAPTER 6PROSPECTS OF DETECTING NEUTRON STARS: RATES AND IMPLICATIONS
6.1

Introduction

Isolated neutron stars are formed as a result of BNS mergers or supernovae, so it isinteresting to discuss the rates at which these events are expected to happen. Astrophysically,
we can gain some insight into the rates of mergers by observing how often we see gamma ray
bursts (GRBs), which could be indicative of mergers, as in the case of GW170817, or stellar
collapse. However, in this chapter we will focus on what GW observations alone can tell us.
In Section 6.2, we explain how event rates are calculated theoretically and empirically basedon GW data. We then describe ways to optimize our sensitivity to long duration transient
GW signals in Section 6.3. In Section 6.4, we discuss what we can learn from a hypotheticaldetection of a power-law GW signal. Finally, we calculate constraints that can be placed on therates of BNS mergers accessible to current and future detectors in Section 6.5.
6.2

Merger Rates: What is Known

In a 12 month period (20 May 2018- 19 May 2019), an average of 1.25 GRBs per day
above threshold, and 0.25 GRBs per day subthreshold (56 ), were observed by Fermi-GBM and
Swift-BAT (57 , 58 ). It is possible that some of these GRBs are due to mergers of neutronstars.
There are different ways to estimate the rates of neutron star coalescences. Using one
gravitational wave observation (GW170817) and an observation time equal to the duration−3of LIGO’s second observing run, LIGO/Virgo inferred a merger rate of 1540+3200
yr−1−1200 Gpc

(1 ) by applying a Bayesian method (59 ). The method is as follows: many injections are done
uniformly throughout the parameter space, and if the injection is recovered (above a certain
SNR) with a false alarm rate lower than a given threshold, then a signal in this portion of theparameter space is detected. By spreading out injections across the source parameter space
and at different redshifts, the space-time volume < V T > to which the detector is sensitive isestimated. Then, given a certain confidence level c, the rate Rc can be estimated (59 ):

59

Rc =

−ln(1 − c)
<VT >

(6–1)

More precisely, this space-time volume is (60 , 61 ):Z zmax< V T >= T0

dVc 1
Nf oundp(θ)
dθdzdz 1 + z
Ninj

(6–2)

c
where dVis the differential co-moving volume of injections in the astrophysical parameterdz

space, p(θ) is the distribution function for the astrophysical population of neutron stars fromwhich injection source parameters are drawn , Nf ound is the number of found injections, Ninj
is the total number of injections done, z is the redshift, T is the observation time and θ are all
Rz
c 1
dz is in fact the astrophysical volume in whichsource parameters. The quantity 0 max dV
dz 1+zthe injections were done.
Other ways to estimate the rate of BNS mergers include using just GRBs as indicators ofmergers, or assuming a proportionality between merger rates and stellar birth rates in spiralgalaxies nearby. The latter birth rate is estimated from the measured blue luminosity, which
can be done because of the binaries’ short lifetimes in Milky Way-like galaxies, and slow rate ofchange in galaxies’ star formations rates (62 , 63 ). Moreover, observations of pulsars with thesame spin period allow a calculation of the scale factor (64 , 65 ). The scale factor relates thenumber of pulsars identifiable in a galaxy to those potentially found by pulsar surveys (whose
search volume is smaller than the galaxy), meaning it is a measure of how many similar pulsars
to ones we have already found in BNS systems exist in a galaxy. By dividing the scale factor
by the estimated lifetimes of binaries, the merger rate is obtained, taking into account all
detected BNS systems. Corrections to the scale factor can be made because (1) the size of thedetectable emission cone of pulsars is uncertain, (2) the pulsars are just too faint to be seen by
any pulsar survey, or (3) the orbital periods are so short that the pulses are not distinguishable
due to high orbital frequency (66 ).

60

6.3

Optimizing Sensitivity

Since our interest will be in transient CW signals, we will have to deal with very highspindowns, meaning that they progress to lower frequencies (and hence improved detector
sensitivity) very quickly. However, higher frequency signals have higher amplitudes, as seenin Equation 3–34, but if a signal has a lower spindown, longer TF F T can be used, improving
sensitivity (see Equation 9.6). This means that for different combinations of f, f˙, Tobs and
TF F T (parameters defined in 3.5), the sensitivity will change. We therefore compute a figureof merit for the sensitivity for various combinations of these parameters, showing in which
portions of the parameter space it is better to analyze short and stronger signals versus longer
and weaker signals.
The figure of merit is:S(f ) ∝

(TF F T Tobs )1/4 2
√
f ,
Sn

(6–3)

which is computed for every frequency that the signal passes through during Tobs .
For braking index n = 5 and Tobs = 105 s (∼ 1 day) and for a range of initial frequencies
and spindowns, we plot the optimal TF F T and the amount of time skipped if the GW begins
at t = 0, using the noise power spectral density from O2 Livingston data. Signals at higher
frequencies/spindowns tend to require shorter TF F T , which means that it may actually bebetter to exclude some data so that we can use longer TF F T in order to optimize sensitivity forthis analysis. This is because we can use longer TF F T for signals with lower spindowns. This isshown in Figure 6-1.
The optimal TF F T and amount of time to skip are also both sensitive to the braking index.
We display this in Figure 6-2.
6.4

Implications of a Detection

If a fixed power law signal is detected by our methods, described in chapters 9, 12 and 13,
we will have an estimation of f0 , k and n. From these parameters, and from an estimation of
the distance d due to an electromagnetic counterpart, we can break the degeneracy between I

61

Figure 6-1: For a wide range of frequencies and spindowns, we plot the optimal TF F T in seconds (left) and how much time to skip in the analysis after a NS is born (right), with thecriteria that the sensitivity in Equation 6–3 is maximized.

Figure 6-2: We vary the braking index for a fixed f0 = 2000 Hz, f˙0 = −10−1 Hz/s, Tobs = 105
s. Left: optimal TF F T ; right: amount of time skipped, both as a function of braking index.
and  in equations 3–23 and 3–33 if we assume that all the energy emitted by the pulsar goes
into GWs. This section follows closely the analysis presented in (47 ).
For the case n = 5, we can arrive at:=qt
5c 1 + τ f˙0
4π 2

h0 d

f03

(6–4)

and the z-component of the moment of inertia is:I3 =c3 h20 d2 f0
5G (1 + τt ) f˙0

62

(6–5)

Moreover, we can also deduce an ellipticity-independent quantity from equations 3–22 and
3–23:
Ė= I3 = 4πω ω̇Zρ(r)r4 dr(6–6)

Ė can be estimated by fitting waveforms to the time-series data with different amplitudes h0
for given recovered parameters of f0 ,f˙0 , and n (under the assumption that all energy emitted
goes into GWs) (12 , 67 ):Ė =2π 2 c3 2 2 2
d f h05G(6–7)

In terms of the GW frequency and spindown f and f˙:Ė= 4π 2 I3 = 16π 3˙ffZρ(r)r4 dr(6–8)

With these parameters, we can begin to deduce what the moment of inertia of theneutron star is, giving us insight into its equation of state.
6.5

Rates and Possible Detections of Newly Born Neutron Stars in Current andFuture Detectors

The existing gravitational wave detectors have been able to see BBHs and BNSs out toO(Gpc) and O(Mpc) distances, respectively. However, due to longer arms and a triangulardesign (68 , 69 ), future detectors, specifically the US-based Cosmic Explorer (CE) and Europebased Einstein Telescope (ET) will be able to see black hole binaries out to redshifts of z ∼ 20.
Specifically the event rates inferred by GW detections in advanced LIGO/Virgo imply hundredsof events per year for future detectors.
It is more instructive now to consider a theoretically calculated rate, co-moving volume,
etc. as opposed to one based on injections described in Section 6.2. Essentially, the rate
(number of detections per time) Ṅ is related to the source density N˙ and the probability that
the SNR ρ of the signal is above a threshold (70 –72 ):63Ṅ = 4π N˙Z ∞drr2 P (ρ2 > ρ2det )

(6–9)

0

We need to relate the probability to the flux of gravitational waves (energy per unit area
per unit time) from a single source, which is (3 ):

FGW =

c3
< ḣ2+ (t) + ḣ2× (t) >16πG(6–10)

When assuming isotropic emission, we can find EGW over the entire time of GW emission:

EGW = 4πr2 T FGW
Defining h2rss =R∞−∞(6–11)

df |h̃+ (f )|2 + |h̃× (f )|2 , computing a time average over the duration T of

the GW, and taking the Fourier transform of Equation 6–10, we arrive at:

EGW =π 2 c3 2 2 2r f0 hrssG(6–12)

Now that we have a relationship between the strain hrss and EGW , we can directly relate the
SNR ρ to EGW . For a narrowband signal, ρ is simply related to hrss as (70 ):Z ∞|F+ h̃+ (f ) + F× h̃× (f )|2
ρ =2
df
S(f )
−∞
2

(6–13)We expand the square, and drop second order terms in h because for elliptical polarizations,h+ and h× are orthogonal, and for linear polarizations, h× = 0. Since we are considering
a narrowband signal for a short duration, the beam pattern functions can be seen as being
constant in time. Knowing this, we arrive at:

ρ2 = Θ2

h2rss
S(f )

(6–14)

where Θ2 is the combination of the detector beam pattern functions and source polarizationand inclination angles. Now, we can use Equation 6–13 to relate ρ and EGW :

64

ρ2 =

Θ2 G
EGW
α π 2 c3 S(f )r2 f 2(6–15)

where α = 1, 8/15, 2/5 for isotropic, linear or circularly polarized emission.
Using Equation 6–15, we can rewrite the probability term in Equation 6–9 as:P (ρ2 > ρ2det ) = P (Θ2 > r2 /r02 )(6–16)

where:

r02 =

EGW
G
2
3
απ c S(f )f 2 ρ2det(6–17)

Thus the detection rate is:
 Z ∞

4 3 ˙
2
2
2
dxx P (Θ > x )Ṅdet = πr0 N 3
3
0(6–18)

Defining the effective detection range Ref f as the radius of a sphere with volume V such thatthe detection rate is N˙ V :
 Z ∞
1/3
2
2
2
Ref f = r0 3
dxx P (Θ > x )(6–19)

0and evaluating the integral numerically (there is little difference among the value of thisintegral for different polarizations), we arrive at:
Ref f ≈

Egw
G
2
3
2π c S(f )f 2 ρ2det

1/2(6–20)and the rate is:4 3˙Ṅ = πReffN
3(6–21)We plot the strain curves of current and future detectors in Figure 6-3 (we have taken thenoise spectral densities from (73 )) and the effective distance we can see, for ρdet = 8 and BNS

65

masses of 1.4Msun in Figure 6-4. Ṅ is determined by how many signals we actually see (and in
the case of no detections, how many injections in the space-time volume < V T > we recover).

Figure 6-3: Strain curves for advanced LIGO at design sensitivity, advanced Virgo with squeezing employed and Einstein Telescope at design sensitivity. Einstein Telescope will benefit froman order of magnitude increase in sensitivity relative to LIGO and Virgo.

66

Figure 6-4: The effective reach of a search for binary neutron stars is plotted as a function of detection frequency for different detectors. Parameters: m1 = m2 = 1.4Msun ,EGW = 10−2 Msun c2 ,ρdet = 867

CHAPTER 7DATA ANALYSIS
7.1

Introduction

In this chapter we will discuss various data analysis techniques used in GW searches. We
begin in Section 7.2 by describing the most sensitive technique, matched filtering, and why it
is difficult to apply it in CW searches. In Section 7.3, we discuss the ways in which GWs are
influenced as they propagate from their source to the earth, which result in many time shiftsthat must be accounted for when running a real search.
7.2

Ideal Detection Technique: Matched Filtering

Let d(t) = h(t) + n(t), where d(t) is the data that we take, h(t) is the signal we would liketo find, and n(t) is the noise we do not want. We would like to filter the data, i.e. multiply it,by some optimally chosen function K(t). We define the output of the data multiplied by somefilter function K(t) over time as:

dˆ =

Z +∞
dtd(t)K(t).

(7–1)

−∞

The signal power S is the expected value of dˆ if a signal is present in d(t). We calculate this
as:

Z +∞
S =

dt < d(t) > K(t)(7–2)

dt < h(t) > K(t)(7–3)

df < h̃(f ) > K(t),

(7–4)

−∞
Z +∞

=
−∞
Z +∞

=
−∞

68

where < .. > is an ensemble average, and then we determine the noise to be:
N 2 = σ 2 =< d(t)2 > − < d(t) >2= < d(t)2 >Z +∞dtdt0 K(t)K(t0 ) < n(t)n(t0 ) >
=
−∞
Z +∞
0 0
=
df df 0 K̃(f 0 )K̃(f )e2πi(f t−f t ) < ñ(f )ñ(f 0 ) > .

(7–5)
(7–6)
(7–7)
(7–8)

−∞

< n(t) = 0 > because the noise is assumed to be random. Understanding that different Fouriercomponents of the noise are uncorrelated, we only have a contribution to the integral whenf = f 0:
Z +∞

2

N =

1
df Sn (f )|K̃(f )|2 ,
2
−∞

(7–9)

where Sn (f ) is the spectral power density of the noise, units of Hz −1 . Sn (f ) is measured anddepends on the different types of noise sources that are present at each frequency.
To find this optimal filter function, we can rewrite S/N into inner products between thenoise power spectral density and the signal. When viewed in this way, the max S/N is whenthe “signal vector” is parallel to the “noise vector”. Essentially, this means that the ideal filter
function is the h(t) signal itself:


S
N

2

Z ∞
=4
0

|h̃(f )|2df.
Sn (f )

(7–10)For a CW signal, taking the Fourier transform of Equation 3–25 (55 ):1h̃(f ) = δ(f − f0 )[F+ h0,+ + F× h0,× ],2(7–11)

Accounting for the fact that the observation time is fixed at T , the Dirac Delta function must
be replaced by a regularized Dirac Delta function (55 ):

69



S
N

2

2

Z ∞

= |F+ h0,+ + F× h0,× |df
0= = |F+ h0,+ + F× h20,× |2

δ(f − f0 )δ(0)Sn (f )T
Sn (f0 )(7–12)(7–13)We see that for CW signals, the SNR increases with the square root of the observation time.
Here, we remind that h0,+ and h0,× are functions of just f0 and time, so they can be pulled
out of the integral in Equation 7–13. The beam pattern functions are also narrowly peakedaround the sidereal rotational frequency of the earth, and can therefore be pulled out of theintegral as well.
Matched filtering is very computationally expensive because one needs to try many
combinations of parameters to obtain the best match for a potential signal. Moreover, it relieson the fact that one can fully trust the model. Unlike in BBH and BNS systems, where thephase evolution is known, the physics governing the GW emission of isolated neutron stars is
unknown.
7.3

Time Delays

The GW that arrives at our detectors has been altered by a variety of general relativistic
factors that must be accounted for. In CW data analysis, what we wish to do is transform thedata we receive on earth (a non-inertial reference frame with many time delays) to the solar
system barycenter (SSB) frame, which is assumed inertial relative to the source of GWs.
We can express the overall time delay in terms of three separate delays (55 ):tb = t + ∆R + ∆S + ∆E ,(7–14)

where tb is the time in the SSB frame, t is the time we measure at the detector, ∆R is the
Roemer delay, ∆E is the Einstein delay, and ∆S is the Shapiro delay. We explain these delaysin the following subsections.

70

7.3.1

Roemer Delay

The Roemer delay (74 ) is a correction that at first order results because the earth/sun
distance changes, which affects our transformation from t to tb (55 ):∆R =ros · n̂,c

(7–15)

where ros = roe + res + rsb is the changing observer-SSB distance due to changes in distance
between observer/earth due to earth’s rotation about its axis (roe ), between earth/sun due
to earth’s revolution around the sun (res ), and between sun/SSB due to the sun’s revolution
around the SSB (rsb ). n̂ is a unit vector pointing towards the source from the SSB.
7.3.2

Shapiro Delay

The Shapiro delay (75 ) refers to the fact that the gravitational field of the solar systemcauses light or GWs to travel along its geodesics, affecting the arrival times (55 ). This is aweak-field effect, so we can start with the Newtonian metric:

ds2 = −[1 − 2φ(~x)]c2 dt2 + (1 − 2φ(~x))d~x2(7–16)

where φ = −GM/(rc2 ) is the Newtonian potential, and we assume that GWs propagate at the
speed of light, so ds2 = 0. Keeping only terms O(φ), we obtain:s

1 − 2φ(~x)
1 + 2φ(~x)

1 − φ(~x)
1 − φ(~x)≈ ±×
1 + φ(~x)
1 − φ(~x)

cdt = ±

cdt ≈ ±[1 − 2φ(~x)](7–17)(7–18)(7–19)Let us define te as the time of emission, and tobs as the time we observe the GWs, and ~rp as
the position of the pulsar and ~robs as the position of the observer at tobs :71

Z ~rp(1 − 2φ(~x))d~xc(tobs − te ) =(7–20)

~
robs

Z ~rp
= |~rp − ~robs | − 2

φ(~x)d~x(7–21)

~
robs

Z ~rp
= |~rp − ~rb | + (~rb − ~robs ) · ~n − 2

φ(~x)d~x(7–22)~robs

We have made a clever transformation of the coordinate distance, which is valid because
the distance between the SSB and the pulsar is much greater than the distance between thebarycenter and observer. The only reason for this transformation is to express the delay interms of the Roemer delay and identify the Shapiro delay piece:



Z
1
1
2 ~rp
tobs ≈ te + |~rp − ~rb | + (~rb − ~robs ) · ~n −φ(~x)d~xccc ~robs(7–23)

We can identify the piece in (...) as tSSB , the time that the signal would arrive if there were no
general relativistic effects. The second term represents the Roemer delay, and the third term is
the Shapiro delay, which depends on the mass and distance/geometry to the source:2∆S = −cZ ~rpφ(~x)d~x(7–24)~robs

The biggest contribution to ∆S is the sun, so we can write:
2GM∆S = − 3
c

Z D
0

dρ
r

(7–25)

where D is the distance between the sun of mass M and the source, and we are integrating over the earth/source distance ρ. This is shown schematically in Figure 7-1. By the
Pythagorean theorem, we can write:

r2 = (dES + ρ cos θ)2 + ρ sin θ272(7–26)

where θ is the angle between earth and the source, respectively, and dES is the earth/sun
distance, on average 1 AU.
By combining equations 7–25 and 7–26, we find the Shapiro delay to be:


 
2GM
2D
− ln(1 + cos θ) > 0
∆S =
ln
c3
des

(7–27)

The first term is a constant offset in the arrival time of the signal, dependent on the earth andthe sun distances/mass, and the distance of the sun from the source. The second term is only
based on the relative source/earth location. The Shapiro offset is always a delay, hence why∆S > 07.3.3

Einstein Delay

The Einstein delay (76 ) is caused by the fact that the clock we use on earth to measure
the arrival time of the GWs measures its own proper time, and needs to be converted tothe coordinate time of the SSB. Starting from the same metric Equation 7–16, writing
ds2 = c2 dτ 2 , and dividing everything by c2 dt2 , we see:


1/2
v2
dτ
= (1 + 2φ) − (1 − 2φ) 2dtc(7–28)

where φ = φ(~xobs ) and v = d~xobs /dt is the velocity of the earth around the sun (and its
rotational velocity, but this is a small correction). Taylor expanding and keeping terms O(v)and O(φ):
dτv2=1+φ− 2dt
2c

(7–29)

φ is the gravitational redshift; the v 2 /(2c2 ) term is the transverse Doppler shift. Integrating
and rewriting so that we express everything in the coordinate time:
Z t
t≈τ+


v2
− φ dt2c2

73(7–30)

Figure 7-1: Diagram showing the geometry of the source relative to the earth and sun, used inour calculation of the Shapiro delay. All bold quantities in this figure refer to vectors.
where both φ and v are time-dependent, and the second term is the Einstein delay:
Z t
∆E =


v2
− φ dt2c2

74(7–31)

CHAPTER 8
THE ORIGINAL FREQUENCY HOUGH
8.1

Introduction

The FH is an example of a hierarchical method, and is shown in Figure 8-1. Hierarchicalmethods reduce computational time in exchange for a small loss in sensitivity, which canbe made up for by analyzing data for a long time (O(years)). The data are divided into
small time segments which are analyzed coherently, then combined incoherently. “Coherent”
analysis refers to performing FFTs on individual data segments and retaining the signal’s
phase information: FFTs transform the strain h(t) to the frequency plane of the detector, sothat a periodic signal would show up as a spike in the frequency plane. Each FFT is actuallyinterlaced by 50% to reduce the known effect of truncating an infinite time series. The data
are smoothed by applying a flat-cosine window function.
We discuss how we create a database of FFTs of different lengths in Section 8.2. In
Section 8.3, we describe the creation of thresholded time-frequency maps, which are theinput to our method, the FH Transform, examined in detail in Section 8.4. Finally we give atheoretical estimate of the FH’s sensitivity in Section 8.5.
8.2

Construction of Short Fast Fourier Transform Databases

The first step in the method involves constructing a so-called short Fast Fourier Transform
Database (SFDB), which is essentially the raw strain data Fourier transformed with one offour TF F T that cover slightly different frequency ranges. There are currently 4 databases: one
up to 128 Hz with TF F T = 8192 s, one up to 512 Hz with TF F T = 4096 s, one up to 1024
Hz with TF F T = 2048s and finally one up to 2048 Hz with TF F T = 1024 s. The databases
are constructed in this way to confine the frequency evolution due to earth’s rotation to one
frequency bin. We can estimate the maximum TF F T that can be used as a function of thefrequency with the following logic: we know the Doppler shift equation is:
f (t) = f0

~v · n̂
1+
c

75


(8–1)

Figure 8-1: Scheme of the Frequency Hough hierarchical method.
It is the change of the velocity of the detector with respect to the source due to earth’s
rotation and orbit around the sun that leads us to consider the maximum TF F T we can use.
We imagine that in time ∆t, the velocity changes by ∆v. Hence the change in frequency is:∆f = f0

∆v
c

(8–2)

If the earth turns an angle ∆θ in the integration time T , then we have ∆θ = ΩT . We
need now to relate this to ∆v, so we can imagine that the earth is moving at a certain speed
v1 = RΩ cos θ and that a nearby part of the earth is moving at a speed v2 = RΩ cos(θ + ∆θ),
where R is the average radius of the earth. Then:76

v2 − v1 = ∆v = RΩ(cos(θ + ∆θ) − cos θ)(8–3)

Assuming T  86400 s (one day) so that ∆θ  1, we can expand the first term with atrigonometric identity. Neglecting O(∆θ2 ) terms, we arrive at:∆v = RE Ω∆θ = RE Ω2 T

(8–4)

where RE is the distance to the rotation axis at the latitude θ (where the detector is at
the initial time). Writing Ω = 2π/TE , we arrive at the maximum coherence time Tmax :TE
Tmax =
2π

r

c
RE f0

s
≈ 1.1 × 104

100 Hzf0(8–5)

where Tmax is the maximum time one can analyze without having to do a Doppler (sky)correction and TE is the rotational period of Earth. We can see that lower frequencies permit
longer TF F T , and vice-versa.
We can repeat this order-of-magnitude analysis for the orbital motion of the earth aroundthe sun, and show that its effect is actually less important than that of the rotation around theearth.
For each FFT in the SFDB, the average noise is estimated based on an auto-regressive
mean that tracks the variation in the noise without being sensitive to peaks in the spectrum.
Here, auto-regressive means that the estimation for the noise in the spectrum is based not juston the frequency at which the estimate is done, but also the frequencies nearby. In this way,large “glitches” can be taken out of the data, as they will be literally seen as spikes. But a(transient) continuous gravitational wave signal will be weak, probably not so high above thenoise, so the contribution due to a signal will not be included in the spectrum estimation. See(77 ) for more details.
A visual representation of the calculation of the auto-regressive spectrum and selection ofpeaks is shown in Figure 8-2.

77

Figure 8-2: (a) The auto-regressive and power spectrum are shown and (b) This is a zoomedin version of a peakmap showing the selection of just local maxima in time and frequency.
8.3

Peakmap Creation

We now would like to create a time/frequency representation, but one that is useful forour adaptation of the Hough transform. To do this, we take each FFT in the database anddivide it by an estimation of the auto-regressive spectrum. This is called “whitening” the data.
At this point, we have essentially normalized power, and we apply a threshold of 2.5 to exclude
likely sources of noise (77 ).
In more details, we calculate the power spectral density Si (the periodogram) for each
FFT:

Si =| F F T |2

(8–6)

A
σ

(8–7)

FFT =

A is the power at a particular frequency and σ is the standard deviation of the noise.
Then the ratio between the power spectral density and the auto-regressive spectrum SiAR
is calculated:78

Ri (j) =

Si
SiAR

(8–8)

This ratio is computed across all frequency bins j for the ith FFT, and is usually around1 unless a peak at some frequency in the spectrum is present (see Figure 8-2 for a visual
representation of this). We compare the ratios to a threshold value defined based on thenumber of total candidates we want, therefore the threshold must be chosen carefully as it
affects the sensitivity of the search. A peak occurs is defined by a particular FFT’s start timei and frequency bin j. The amplitude of the peak does not factor into this analysis, however
we only select peaks that are local maxima: this reduces the computational time and decreases
sensitivity to glitches (huge spikes in the spectra due to noise) in exchange for a small loss insensitivity.
We create the Hough map from the peakmap; therefore, the selection of peaks is a criticalstep in our analysis. Even if a peak is due to noise and is above the threshold, it will appear in
the peakmap, and if a signal does not have a peak in the peakmap, it is lost. We can write theprobability that a peak is due to noise as p0 = P (θ; 0) and pλ = P (θ; λ), if a peak is due tosignal, where θ is the threshold and λ is the spectral amplitude of the signal, defined as:

λ=where h(f ) =R∞−∞

4|h̃(f )|2,
TF F T Sn(8–9)h(t)e−2πif t dt is the Fourier transform of h(t), the GW signal, and Sn is

the spectral noise density, which changes as a function of frequency. In Gaussian noise, theprobability distribution of the peakmap is exponential, so we can calculate p0 . The probabilityof having a peak in a frequency bin between two values in the spectrum, x and x + dx is:p(x) = e−x dx. The probability that this peak is a local maximum, that is, the two values
nearby are not peaks, is (1 − p(x))2 . Thus,
Z ∞
p0 =
θ

1
e−x (1 − e−x )2 dx = e−θ − e−2θ + e−3θ3

79

(8–10)

The probability of finding a peak due to a signal pλ is calculated based on knowing thatthe probability density of the spectrum is a non-central χ2 distribution with 2 degrees of
freedom and non-centrality parameter λ. This probability density p(x; λ) can be written as:√p(x; λ) = e(−xλ/2) I0 ( 2xλ)(8–11)

where I0 is the modified zeroth order Bessel function whose mean value x̄ = 1 + λ/2. Then,
applying the small signal approximation, λ  x, we can expand I0 and the exponential term to
get a simpler expression for p(x; λ):p(x; λ) ≈ e

−x



λ λ
1− + x
2 2(8–12)

Again we can compute the probability of selecting a local maximum:Z ∞Z xp(x; λ)pλ =p(x0 ; λ)dx0 dx(8–13)

0θ

Doing the integration, and keeping terms up to O(λ), we arrive at:λpλ = p0 + θ(e−θ − 2e−2θ + e−3θ )2(8–14)

Knowing these two probabilities, we can discuss now how we determine the threshold usedin the peakmap.
The threshold is chosen based on the maximization of the average value of our detection
statistic. The FH uses a simple detection statistic:CR =y − z̄
y − N p0=p.
σ
N p0 (1 − p0 )(8–15)The Hough map is binomial because each peak is either a noise peak or a signal peak, so we
can use the expressions for the mean and standard deviation of a binomial distribution. y is thenumber of peaks at in the Hough map, z̄ is the average number count of the noise, and σ isthe standard deviation of the noise. It is upon maximizing the expectation value of the critical80

ratio that we can arrive at determining an adequate threshold. Assuming a signal is present, we
can write that the average value of the number count will be y = N pλ . Thus (78 ):√N (pλ − p0 )= N ψ(θ; λ).
CR = pN p0 (1 − p0 )(8–16)

By plotting ψ(θ; λ) to see at which value ψ is maximal, we can determine the appropriate
threshold. Actually, plotting this function shows that the ideal threshold is θ = 2, but we useθ = 2.5 to reduce computational time at a small cost in sensitivity.

Figure 8-3: (a) Pre-Doppler corrected peakmap for a signal at f = 355.725 Hz. The signal
is visible at a sinusoidal curve, and a glitch is present at a higher frequency. (b) Doppler corrected peakmap. Signal is visible as horizontal line, glitch has had its power spread in multiple
frequency bins.
8.4

The Frequency Hough Transform

The FH is a clever implementation of the Hough transform, a pattern recognition
technique originally used to find tracks in particle physics experiments. This method maps thepoints (the peaks) in the time/frequency domain of the detector to the frequency/spindown
domain of the source. This method requires that the signal’s time/frequency evolution belinear, so a Taylor series expansion is used as shown in 3–24, where the second order spindownf¨ is typically neglected. The method works as follows: for each position in the sky, thepeakmap is corrected for the Doppler shift. We see how powerful this Doppler correction can

81be to distinguish a glitch and signal in Figure 8-3 Then, for each peak in the time/frequency
plane, a line is drawn for different choices of f˙ at different times t:

f˙ = −f0
f+.
(t − t0 ) (t − t0 )(8–17)

From this information, the values of f0 (the y-intercept) are calculated, and thenhistogrammed in the f0 /f˙ plane- see Figure 8-4. A candidate signal will have a large numberof peaks that are traced back to a particular f0 .
In practice, the FH employs an adaptive procedure that weights less noisy areas of thedetector and accounts for the differing detector sensitivity across the sky (78 ).

Figure 8-4: This is a zoomed in version of a Hough map created for a signal at f=
355.725Hz, f˙ = −2.69 × 10−18 , longitude = 178.0718◦ , and latitude = −26.1775◦ . The color
bar represents the number of FFTs that captured a peak at a particular frequency/spin-down.

82

8.5

Theoretical Sensitivity Estimate

Sensitivity is defined as the minimum signal amplitude that a method can detect at a
certain confidence level (typically, 50%, 90% or 95% is used) in many experiments (usually thepercentage of 100 signals at a given frequency/amplitude is sufficient). It can be calculated
independently of any search; hence, it provides a good benchmark when computing upper
limits in the real search-one’s sensitivity cannot be better than the theoretical estimation!
The crucial detail is that the probability distribution of the Hough map, though binomial,is assumed to be Gaussian:Pn (θthr ; λ) = √

1
2πσ 2

(n−µ)2e− 2σ2(8–18)

This assumption works well if the number of FFTs is large, and the threshold on the peakmap
is such that the probability of selecting a peak is not too close to 0 or 1 (see equations 8–10
and 8–14. Moreover, it is clear that this assumption can break down in frequency bands thathave many disturbances, which occur frequently in the detectors. The probability distributionof the Hough map is therefore not completely Gaussian.
Integrating this probability distribution for signals that have surpassed our threshold
number count nthr (which can be related to the critical ratio by Equation 8–16):
Z ∞Pn>nthr (λ) =Pn (θthr ; λ)dn = Γ(8–19)

nthr

where Γ is the confidence level at which one claims a detection. Then:!
nthr − N pλerfc p= 2Γ
2N pλ (1 − pλ )(8–20)

Now, we must also include the false alarm probability into our analysis. We perform asimilar integral, expect that λ = 0 so the probability distributions we use are p0 :Z ∞Pn (θthr ; 0)dn = Pf∗a

nthr

83

(8–21)

therefore:
nthr (N, θthr , Pf∗a ) = N p0 +p2N p0 (1 − p0 )erfc−1 (2Pf∗a )(8–22)

Now, we can write Pf∗a = Ncand /N , since we only select the candidates with the highest
number count, but it turns out that Pf a is so small that the term that is linked to it can beneglected, since it is compared to a term with p0 , the probability that there are noise peaks.
We combine equations 8–20 and 8–22, along with the expressions for p0 and p1 , equations
8–10 and 8–14, in the small signal approximation. Essentially, the equation can be put into aform in which it is solvable as a quadratic equation, and then of the two possible solutions, theone that gives a physically meaningful interpretation is chosen. In this case, the solution is theone in which λ increases as Γ > 0.5, since we expect stronger signals to be more easily andconsistently detectable. We arrive at:

λmin ≈

2
θthr

s


√
p0 (1 − p0 ) 
−1
CR
−
2erfc
2Γ
thr
N p21

(8–23)

where p1 = e−θthr − 2e−2θthr + e−3θthr .
Now, we must relate this quantity to h0,min , the minimum signal amplitude that we candetect. In order to do this, we first write an expression for the GW signal:h(t) = F+ h0+ cos (φ(t)) + F× h0× sin (φ(t))
1 + cos2 ι
h0+ = h02
h0× = h0 cos ι(8–24)

where φ(t) = 2πf0 t (since the spindown is so small and is confined to once frequency, we canneglect the spindown and assume a constant frequency signal).
Then, we perform the exact Fourier transform of a sinusoidal signal, and account for thefact that when we take the Fourier transform with a finite length (as we do in practice), weshould apply a window function. We obtain:84

h̃(f ) ≈ TF F T(F+ A+ − iF× A× ) sin(π(f − fk )TF F T )2π(f − fk )TF F T(8–25)

where A+ and A× are the detector response functions that depend on the angle between thestar rotation axis and the line of sight, given by:

A+ =a0 + a1c cos Ω⊕ t + a1s sin Ω⊕ t + a2c cos 2Ω⊕ t+(8–26)

a2s sin 2Ω⊕ t

A× =b1c cos Ω⊕ t + b1s sin Ω⊕ t + b2c cos 2Ω⊕ t+(8–27)

b2s sin 2Ω⊕ twhere Ω⊕ is the sidereal angular frequency of the earth, whose coefficients a0 , a1c , etc. dependon the relative source/detector positions, and are given in (79 ).
Since we want to know our sensitivity at any sky location, frequency, etc., we take theaverage over the two beam pattern functions and the two amplitude functions. This involves
averaging over the sky position, the inclination angle, the wave polarization angle, and thefrequency of the signal. After performing these averages, relying heavily on equations B15-B19
in (78 ) we arrive at:

h0,min ≈

4.02

s

1/2

N 1/4 θthr

Sn (f )
TF F T

p0 (1 − p0 )p1

1/4 q√CRthr − 2erfc−1 (2Γ)

(8–28)

This is the minimum amplitude at which the FrequencyHough can detect a signal under theassumption of Gaussian noise.

85

CHAPTER 9
THE GENERALIZED FREQUENCY HOUGH
9.1

Introduction

In order to search for tCWs (described in 3.5), we had to generalize the FH to be able tosee power-law signals, or develop an entirely new method. In the end we were able to adapt
the existing FH Transform to search for a different type of signal. This chapter describes themechanics of this new method, as well as the different problems that we encountered whenusing the Generalized Frequency Hough (GFH).
The main reason that this method works is shown in Section 9.2. We define how we canclaim a detection of a GW, i.e. what our detection statistic is, in Section 9.3. Afterwards in
Section 9.4 we discuss the general procedure to follow up interesting candidates that we may
receive in a search. The construction of the grids in the parameter space that we would use
to run a search, expressed in terms of braking indices, spindowns and initial frequencies, isdemonstrated in Section 9.5. We then estimate the theoretical and empirical sensitivity of our
method in Sections 9.6 and 9.7 respectively, ensuring that they are in reasonable agreement.
9.2

The Coordinate Transformation

The FH works well if the signal’s frequency evolution is linear; however, when we consider
young neutron stars formed after a merger or in a supernova, the frequency evolution is highly
nonlinear. We assume a power-law model for the frequency evolution of a GW signal as in
3–33. However, this equation is not linear and therefore we cannot apply the FH. But if wechange coordinates (31 ):

x=

1
f n−1

; x0 =

1
f0n−1

(9–1)

Under this transformation, Equation 3–33 becomes the equation of a line:x = x0 + (n − 1)k(t − t0 )Points in the (t − t0 , x) plane are mapped to lines in the (x0 , k) plane:

86

(9–2)

k=−

x0
x
+
(n − 1)(t − t0 ) (n − 1)(t − t0 )(9–3)

We now have an analogous equation to equations 3–24.
We show in the left panel of Figure 9-1 a peakmap created with an injected r-mode
signal. In the right panel of Figure 9-1, we see a peakmap in the transformed space. It is this(Doppler-corrected) map that is the input to the GFH. To create these r-mode signals, we
simply simulate their expected phase evolution basd on equation 3–33, because performing
numerical simulations and identifying the fluid motions as r-modes is very difficult (80 ).

Figure 9-1: A peakmap of an injected r-mode (n = 7) signal in white noise with h0 = 1.5 ×
10−22 , f0 = 993.25 Hz, f˙0 = −1.99 × 10−1 Hz/s before and after the change of coordinates.
The color represents normalized power.
Notice that the peaks are denser near the low values of x (which correspond to higher
frequencies). This effect is also seen in the FH map; it occurs because in this new space,
the noise is not uniform. The change of coordinates takes points that are equally spaced infrequency and concentrates them at higher frequencies (lower x values) and spreads them out
at lower frequencies (high x values). So when the GFH is performed, more points in t/x are
mapped to lower x0 values. The non-uniformity of the noise can affect the sensitivity; therefore,
we construct peakmaps in a larger frequency band than we plan to actually analyze (20-25%
larger). This is by far the simplest way to deal with the non-linearity of the noise in this new
space at the cost of a very small increase in computational cost.

87

In Figure 9-2, we show the output of the GFH: a histogram in x0 /k space. The number
count corresponds to the number of FFTs in which a signal with parameters x0 and k wasfound. The lines intersect at the true parameters of the injected signal, and the highest
number count refers to the most likely parameters of the signal.

Figure 9-2: Hough map clearly showing the injected signal recovered, with the parameters
given in the caption of Figure 9-1. Lower values of x0 correspond to higher frequencies, andvice-versa. The color corresponds to the number of lines that hit each x0 /k bin. The referencetime t0 is taken to be at the beginning of the time segment analyzed.
9.3

Detection Statistic and Candidates

In both the original and generalized FHs, we evaluate a detection statistic called thecritical ratio to determine if a candidate is significant:

CR =

y−µ
σ

88

(9–4)

where y is the number count in the Hough plane in a particular bin (a potential candidate), µand σ are the average and standard deviation of the number counts in the map due to noise.
In each of the maps, we select significant candidates. Each candidate will have an
x0 , k, n, longitude and latitude. From x0 and k, we can find f0 and f˙0 with the inverse
transformation of Equation 9–1. Candidates are chosen in each “square” of the Hough map,
whose size (in bins) is determined by the number of candidates we wish to select in eachmap. We employ this way of selecting candidates because we want to cover the parameter
space as uniformly as possible, and so that we are not blinded by major disturbances, nor
the non-uniformity of the noise. Because the noise is nonuniform, one will always see higher
number counts at lower values of x0 , even in pure Gaussian noise. So we cannot simply selectthe candidates with the highest number counts; we must evaluate the critical ratio in each(almost independent) subsection of the Hough map.
Once the analysis has been performed on all detectors and a set of candidates has been
produced form each one, we then perform coincidences, which simply find if candidates
selected in each detector are close to each other in the parameter space. We define “close” inthe parameter space as being within a certain number of x0 and k bins. For the candidates ineach detector, the distance between them in the parameter space is calculated in the following
way, as was done in (78 ):sd=x0,2 − x0,1
δx0

2


+

k2 − k1
δk

2
(9–5)

where δx0 is the size of the bins in x0 , and δk is the varying size of the bins in k (see 9.5).
Note that in general we should also include error in the sky location, but for now we are only
considering a search where the position of the source is fixed.
9.4

Follow-up

We now employ a follow-up procedure for significant candidates returned by the search.
For each candidate, we have an estimation of its f0 , f˙0 , longitude, latitude, and n. We cantherefore reconstruct the frequency evolution of the signal f (t), and more importantly the89phase evolution. We make a heterodyne correction h’(t) = h(t)ejφsd (t) , where the phase φsd is(81 ):Z Tobsφsd (t) = 2π ·t0

f0

1 dt’ + φt0
 n−1

(9–6)

1 + k(n − 1)f0n−1 (t’ − t0 )We are essentially removing all phase variation for a GW signal in the data; hence, even a
rapidly varying signal will be monochromatic. In practice, however, this correction will not be
perfect, so there will be a residual spin-up or a spindown.
This issue can be addressed in a follow-up search by increasing TF F T by a factor of
about 10, but is generally limited by the amount of follow-ups one can afford. Moreover, thefrequency variations that CW searches can handle are orders of magnitude smaller than fortCW searches. Therefore, we have to quantify how much we can increase TF F T as a functionof how wrong we are in our correction. Additionally, because of 3–33, errors in our returned f0
and f˙0 parameters will affect the new frequency evolution differently.
In general, we have found that if we would like to increase the TF F T by a factor of 10,
we would need to consider a grid around the candidate’s parameters over-resolved by a factorof 10 in each parameter, determined by requiring that when using the new FFT length, acandidate would not spin out of one frequency bin during the observation time. This would
mean that one candidate would require 100 separate follow-ups. See Figure 9-3.
To actually perform the follow-up, it is easiest to use the Band-Sampled Data (BSD)
framework (81 ). The framework allows us to use any TF F T to correct for the signal’s phaseevolution, and is very flexible and easy to use. BSDs allow for both frequency and time bands
to be extracted when doing injections, cleaning, etc. They are constructed from SFDBs by
taking the inverse FFT of the data and constructing an analytic representation of the signal.
An example: imagine that we run a search using the generalized FH and recover a
candidate at n = 7 (r-mode) with f0 = 795 Hz and f˙0 = −0.25 Hz/s using an initial
TF F T,0 = 2 s. The signal would appear in the peakmap (before the Generalized FH is applied)
as in the left-hand panel of Figure 9-4, which is an injection done within the BSD framework.

90

Figure 9-3: We chose the original TF F T,0 = 2 s, f˙0 = 0.25 Hz/s. We show on the y-axis thefactor by which we can increase TF F T during the follow-up as a function of frequency suchthat a signal would still be confined to one new frequency bin if we use an over-resolution
factor of 10 in both x0 and k. Blue, red, yellow → n = 3, 5, 7 (also top, middle, bottom)
Knowing these parameters, we can then correct the signal’s phase evolution and use longer
FFTs, since the signal’s residual spindown will be lower. In the right panel, by applying
Equation 9–6, we obtain a monochromatic signal, where we have used TF F T = 2TF F T,0 . Aftercorrecting for the phase, we then apply the original FH (78 ) because any residual spindown
remaining (when the phase correction is not perfect) could be small enough such that thetime/frequency behavior is linear. We have done tests by correcting for the wrong frequency,
spindown and/or braking index, and have found that if the parameters are off by around one
bin at the most, we can still recover the signal after coincidences for certain amplitudes.

91After performing the original FH, we select candidates in this Hough map, where acandidate is defined by having a frequency, spindown and sky location, and do coincidences
between the two detectors. If there is no coincidence, the candidate is vetoed; otherwise, we
use the candidates’ recovered parameters to correct the peakmap that was fed to the original
FH (the right panel of Figure 9-4). This correction requires a simple shifting of bins, which
is done for the Doppler correction. After this second correction, we project this peakmap
onto the frequency axis and calculate a new critical ratio. The left panel of Figure 9-5 shows,for one detector, the projection of the peakmap created with TF F T,0 ; the right panel shows
the projection of the peakmap made with TF F T = 2TF F T,0 , that has been created with thefollow-up procedure.
This process is done for each detector individually. We then weight the critical ratio by
each detector’s sensitivity around the candidate frequency, and compute an average criticalratio. If this weighted average critical ratio is greater than that computed on the original
1/4

TF F T,0 peakmaps by a factor of TF F T , then we have a candidate that requires more analysis;otherwise, the candidate is vetoed.

Figure 9-4: Left: peakmap of an injected r-mode (n = 7) signal in white noise with h0 =1 × 10−22 , f0
= 795 Hz, f˙0= −2.5 × 10−1 Hz/s. Right: peakmap after phaseevolution correction; a monochromatic signal is seen at f0 = 795 Hz. The peakmap on theleft was constructed with TF F T = 2 s; on the right, as is the case of the follow-up, we used
TF F T = 4 s, twice as large as in the original analysis. The signal was injected and its phase corrected using the BSD framework.

92

Figure 9-5: A signal was simulated with f0 = 1045 Hz,f˙0 = −0.03125 Hz/s, n = 7, andh0 = 8 × 10−23 in white noise. Left: we performed the Generalized FH, obtained candidates, andcorrected the phase evolution in original peakmap TF F T,0 = 4 s, then projected that peakmap
onto the frequency axis. Right: we then corrected for the phase evolution of this signal while
increasing the TF F T by a factor of 2, ran the original FH, and used the best candidate’s parameters to shift the bins in the peakmap as we do for the Doppler shift. In one detector,
the critical ratio is increased by about 20%; when using both detectors, the weighted average
1/4
critical ratio is increased by about 22%, which roughly equals TF F T
9.5

Grids in the Parameter Space

Each parameter we use in the GFH is discretized in some way. In the following sections,
we describe how we construct the grid in each of the parameters. None of the steps are
uniform.
9.5.1

Grid on x0

The grid in x0 is obtained by taking the derivative of Equation 9–1:
δfdx0 = (n − 1) n
f

(9–7)

The step dx0 changes as a function of the frequency analyzed, as well as the TF F T =1/δf chosen. However, having a nonuniform grid in x0 is computationally very heavy, especiallyat higher frequencies, since we cannot fill every bin simultaneously as we can do if the grid isuniform. Therefore, we over-resolve the grid by taking f = fmax . This results in more bins,
and a splitting of the signal across some x0 bins in some cases. However, it makes the analysis

93computationally feasible. The computation time increase is negligible since we fill all x0 bins
simultaneously when constructing the Hough map.
9.5.2

Grid on k

We must also construct a grid in k, since the Hough map is constructed across all times
for a fixed k. An expression for the step can be derived analytically if we take the log10 ofEquation 3–32:

f˙ = −kf n(9–8)

log |f˙| = n log f + log k(9–9)

Equation 3–32 forms lines in this space, where different braking indices correspond to lines withdifferent slopes. If we consider a transformation f → f + δf and k → k + dk, we can find dk
such that the spindown remains constant when moving one frequency bin δf = 1/TF F T . After
solving equations 9–8 and 9–9 simultaneously, we find:
dk = k1+

−nδff0,max!
−1(9–10)

Since δf  f0,max , we can Taylor expand to better understand the behavior of this grid:dk ≈ −nkδff0,max(9–11)The grid depends on n and k, so it is nonuniform and changes each time we perform theHough.
The minimum and maximum values for k are related to the maximum and minimum
spindowns we wish to analyze, respectively:94

kmin =
kmax =

f˙0,minn
fmax
f˙0,maxnfmin(9–12)(9–13)To understand the properties of the k grid, we plot in Figure 9-6 how the range in k andthe step size dk change as a function of braking index:

Figure 9-6: The ranges of the step dk for different TF F T . Red, blue, and black correspondto TF F T = 2, 4, 8 s (lowest, middle, highest curves) for frequency bands of [1000 2400] Hz,[600 1400] Hz, and [350 800] Hz respectively. Tobs = 5000 s and ∆t = 1 × 106 s. The
upper (lower) curves correspond to the maximum (minimum) values for dk for each choice of
TF F T /frequency band.

95

It is clear that dk decreases with braking index, since the spindowns are smaller. However,
the range (kmax − kmin )/dk decreases with increasing braking index, which means that weshould expect fewer points in the grid on k at higher frequencies. This is shown in Figure 9-7.

Figure 9-7: The number of points in the grid on k as a function of braking index. Less points
are required to capture the variation in spindown at higher braking indices. Red, blue, andblack lines correspond to TF F T = 2, 4, 8 s for frequency bands of [1000 2400] Hz, [600 1400]
Hz, and [350 800] Hz respectively.
9.5.3

Grid on n

The search is model-dependent, and yet the power law models for neutron stars are quite
uncertain. Therefore we search across different braking indices, which can correspond to either
a combination of emission mechanisms- electromagnetic and GW emission from a neutronstar – or a braking index that is very slowly varying in time. We construct a grid of different
braking indices between n = 2.5 and n = 7 with a step that is calculated so that the frequency
96

variation stepping from n1 = n to n2 = n + dn for the signal duration is confined to one
frequency bin:f1 (t) =f0
1 + k(n1 − 1)f0n1 −1 (t − t0 )

 n 1−1f0

f2 (t) =(9–14)1

 1
1 + k(n2 − 1)f0n2 −1 (t − t0 ) n2 −1
1∆f = f2 (t) − f1 (t) ≤ δf =TF F T(9–15)(9–16)

We empirically find the value of dn at each braking index n, k and f0 such that for theduration of the analysis Tobs , the frequency variation remains within one frequency bin. δf isfixed by our choice of TF F T , which depends on the maximum spindown we wish to analyze.
Therefore, the grid depends on the FFT length, the braking index, the spindown range andfrequency range.
We provide here some plots describing the behavior of the grid on n. In Figure 9-8, weshow how the number of points changes in the grid as a function of the braking index. At
higher braking indices, the magnitude of the spindown decreases quicker than at lower braking
indices, meaning that fewer points in the grid n are required to cover the parameter space.
Additionally, the number of points in this grid changes as a function of TF F T for two
reasons: (1) the frequency bin size becomes smaller with increasing TF F T and (2) the
frequency dependence as a function of time becomes much weaker as the spindowns weanalyze decrease (since f˙ = 1/TF2 F T and TF F T is increasing), see Figure 9-9.
9.6

Theoretical Estimate of Sensitivity

We derive an analytic expression for the sensitivity of a semi-coherent search for longtransient periodic signals, such as those emitted by a newborn magnetar, assuming the analysis
is done with the generalized FH transform. The expression is in fact a generalization of thesensitivity computed for standard continuous wave signals: 9–22.

97Figure 9-8: The plot shows the number of points in the grid on n as a function of the braking index. Parameters used in this plot: TF F T = 2 s, f0,max = 2000 Hz, f˙0,max = 1/TF2 F T
n = [2.5, 7], histogrammed into 100 n bins of size δn ∼ 4 × 10−2 .
As in the previous sensitivity estimation, we assume a period signal but that can vary(rapidly or not) in time. Recall that our semi-coherent analysis requires that a signal not spinout of a frequency bin during TF F T . So, in each TF F T , the frequency and therefore amplitudeof the signal is constant.
In order to compute an average sensitivity we first introduce the quantity
Λ=i=N
TF F T X
λ(f (ti ))2Tobs i=1

(9–17)

where λ(f (ti )) is:
λi = λ(f (ti )) =4|h̃(f (ti ))|2
TF F T Sn (f (ti ))98(9–18)

Figure 9-9: This plot shows the number of points in the grid on n as a function of TF F T for a
fixed observation time Tobs = 105 s and maximum frequency f0,max = 2000 Hz. Note that for
longer TF F T , the spindowns we are analyzing are smaller. We see a sharp decay in the numberof points in the grid because the frequency is not varying so strongly, and therefore Equation
3–33 reduces to Equation 3–24.
where h̃(f (ti )) is the Fourier transform of the GW signal, N = Tobs /TF F T is the number of
FFTs used, and Sn (f (ti )) is the detector noise power spectrum. This quantity is essentially thesame as that quoted in 8–9, except that instead of having one contribution due to the strengthof the signal in the observation time, we have many different contributions, each within one
TF F T .
The quantity Λ is an average of the λi over the observation window and satisfies thecondition that if λi = const then Λ = const, i.e. we are back to the standard situation of
CW signals in which the frequency and the amplitude do not significantly change over the99

observation time. One can think of Λ as the signal amplitude decreasing over time (since theamplitude is a function of frequency), but the “strength” of the signal is also proportional
to how long it lasts; hence, we have placed the summation in Equation 9–17. As the signal
amplitude h0 varies with frequency as f 2 we write
h0 (ti ) = Af 2 (ti ) = AFi(9–19)

where:
A=

4π 2 GIzz c4d(9–20), This A depends only the star ellipticity , moment of inertia and distance d. Combining
equationsΛ≈4

TF F T 2 X Fi2 2.4308
A
TF F T
2Tobs
Sn (fi ) 25πi(9–21)

Following the same procedure given in (78 ), we obtain:

Amin =

4.02
1/2

N 1/4 θthrr

N
TF F T

X F2i
Sn (fi )i!−1/2 p0 (1 − p0 )p21

1/4 rCRthr −√

2erfc (2Γ)−1

(9–22)
where θthr is the threshold for peak selection in the whitened spectra, Sn is the noise spectral
density of the detector, p0 is the probability of selecting a peak above the threshold θthr if thedata contains only noise , p1 = e−θthr - 2e−2θthr + e−3θthr , CRthr is the threshold we use to
select candidates in the final FH map, and Γ is the chosen confidence level.
The minimum detectable strain at a given confidence level can be obtained from Equation
9–22 using a suitable “frequency” (indeed h0,min = Amin · frequency2 ). We use the initial
P F2
frequency f0 . The sensitivity depends on the signal evolution through the ratio i Sn (fi i ) . To
calculate the sensitivity we must fix n, f0 , f˙0 and Tobs . That is, for a given detector and given
search parameters (TF F T , θthr , CRthr , Γ etc.) and for each specific signal model, we will
have a different sensitivity. If fi = const, which is a good approximation for a standard CW
P F2
Nf4
signal case, then Fi = f02 and as a consequence i Sn (fi i ) = Sn (f0 ) . Inserting this expression100

in Equation 9–22, we recover the standard sensitivity expression in (78 ). If the emission of
GWs is due to r-modes, the signal amplitude scales with f 3 . In this case Fi = f 3 (ti ) and
A = 1.11 · 10−9 Iα. The distance reach of the search is easily obtained inverting Equation 9–22:

TF F Tdmax = 5.72·10−9 I38 −3 √
Tobs

X F2
i
S
(f
n i)
i

!1/2 

p0 (1 − p0 )N p21

−1/4 sθ
√thr −1CRthr − 2erfc (2Γ)
(9–23)

with I38 the star moment of inertia (with respect to the rotation axis) in units of 1038 kg · m2 and−3 the star ellipticity in units of 10−3 . For r-mode emission Equation 9–23 must be modified
replacing the numerical coefficient by 1.11 · 10−12 and it should be noted that  is in fact themode amplitude α.
9.7

Empirical Estimate of Sensitivity

We empirically estimate the minimum amplitude detectable at Γ = 90% confidence
at some initial frequencies for fixed braking indices of n = 3, 5, 7. We use 100 injections in
Gaussian noise between f0 and f0 + 10 Hz and define a detection as when a candidate selected
in a Hough map is within a distance of 3 x0 /k bins from the injection. In a real search, there
is a coincidence step between candidates selected in different detectors, but for this empiricalsensitivity estimation, we only do coincidences between a candidate found in one detector
and the injection. The granularity (that is, the spacing in h0 ) is 10% of the theoretical h0 for
the corresponding f0 . The minimum detectable amplitudes are shown in Figure 9-11 where
the theoretical sensitivity curves for the same parameters are also plotted. In this plot, wehave used varying estimations of CRthr to reflect the fact that due to the flat noise curve
we assumed, signals with the same amplitude in higher frequency will generate candidates
with higher critical ratios because the amplitude of the GW signal is proportional to f 2 - see
Equation 3–34. The plot demonstrates a good agreement between experimental and theoreticalresults. Additionally, we show plots of detection efficiency as a function of amplitude for a

101

Figure 9-10: Theoretical h0,min for different initial frequencies f0 for n=3, 5, 7 (red, green,
black) where the noise distribution has been weighted with the O2 Livingston sensitivity curvefor f˙0 = 1/16 Hz/s. There appears to be a dependence on the details of Equation 9–22.
few initial frequencies in Figure 9-12. The false alarm probability is ∼ 0.01%, but in generalchanges in each Hough map.
In this case the typical distance reach is of the order of 0.5-1 Mpc, obtained from thetheoretical estimates shown in Figure 9-10 for O2 Livingston data, assuming a moment ofinertia and star ellipticity. The moment of inertia could be larger than the canonical 1038 kg · m2
for a neutron star born after a coalescence. The other parameter which may have some
impact is the threshold for candidate selection CRthr , which depends on how many candidate
follow-ups we can afford. The distance reach is typically smaller for shorter Tobs .

102

Figure 9-11: Minimum detectable amplitude h0,min as a function of initial frequency f0 forbraking indices n = 3, 5, 7 (red, green and black / bottom, middle, and top). The continuous
curves are the theoretical sensitivity estimates; the curves marked with circles are the sensitivity
estimates obtained from injections in white noise at a level consistent with real O2 Livingston
1
data about 1 hour after GW170817 (Sn = 7.94 × 10−24 √Hz); parameters: Tobs = 5000
2
˙s, TF F T = 4 s, θthr = 2.5, f0 = 1/TF F T , Γ = 0.9. To compute the theoretical sensitivity
estimates, we use a varying CRthr reflective of the average critical ratio we recover from injections at 90% confidence, which changes because signals with higher frequencies with the same
amplitudes as those at lower frequencies will be recovered with higher critical ratios.
It is worth noting that in Figure 9-10, the sensitivity depends heavily on the details ofEquation 9–22: specifically, the decaying frequency is a strong function of the braking index,
which can induce changes in the sensitivity curve at various frequencies. This is why at somefrequencies, signals n = 3 or n = 5 have lower values of minimum detectable amplitude.

103

Figure 9-12: Detection efficiency plotted as a function of h0 for n = 5 for Tobs = 5000 s and
TF F T = 4 s. Different colors correspond to different initial frequencies f0 . The curves seemto follow the sigmoid distribution, which is as expected. We theoretically estimate the false
alarm probability to be order of 0.01% using equation 61 of (78 ); however the grids on x0 andk change in each Hough map, so the false alarm probability also changes. Moving from left toright tends to correspond to decreasing f0 .
If we are able, thanks to observations in the EM band or to robust theoretical models, torestrict the possible range of the search parameters, f0 , k, n, we can in principle make a deeper
search. For instance, by making a first rough search over a limited range of braking indexes we
could make a second step searching over the braking index residuals (with respect to the initialrough grid) and this would allow us to increase TF F T and then gain in sensitivity.

104

(a) Time series strain O2 Livingston(b) Time/frequency map made with (c) Transformation of 9-13(b) usingdata
TF F T = 8 s.
equation 9–1(d) Hough map showing injected (e) Phase-corrected time/frequency(f) Hough map, phase-corrected sigsignal’s parameters
map using 9-13(d)nal’s parameters

Figure 9-13: These plots illustrate the Generalized Frequency Hough pipeline, showing how
we estimate the parameters of a GW signal from an isolated neutron star starting with thestrain time series. Note that in Figure 9-13(d), a reference time equal to midway through theobservation time (t0 ∼ 2000 s) was used, which corresponds to a particular.

105

CHAPTER 10
SEARCHING FOR A REMNANT OF GW170817
10.1

Introduction

GW170817 was an extremely exciting event due not only to the physics we learned
from the merger, but also because of the unknown nature of the product of such a merger.
Therefore the LIGO/Virgo collaborations wanted to determine the identity of the remnant
using four separate pipelines. In Section 10.2, we motivate why we should search for a remnant
of a BNS merger and what we could hope to learn from a GW detection. Afterwards in Section
10.3, we show the portion of the parameter space our method, the GFH, was sensitive to inthis search. The next three sections, Sections 10.4, 10.5, and 10.6 deal with practical aspectsof the search, that is, the selection of candidates after the GFH transform, the follow-up of
potentially interesting signals, and the upper limits we were able to place on GWs originatingfrom a remnant of GW170817 respectively. We finish the chapter with some discussion of
future searches in Section 10.7.
10.2

Motivation

GW170817 was the first BNS merger ever observed by LIGO/Virgo. It has created many
opportunities in astrophysics for deep study into the nature of neutron star mergers andneutron star physics. One of these aspects is answering the question, what remained after
the merger? The mass of the two neutron stars, and the mass of the remnant, have quite
wide uncertainties, so it is not clear what formed after the merger. There are essentially 4
options as to what could happen: (1) immediate collapse to a black hole, (2) formation of
hypermassive neutron star (HMNS) that is supported by differential rotation lasting for O(s),then uniform rotation, then collapsing to a black hole, (3) a supramassive neutron star (SMNS)
that exists is sustained by normal rotation for O(hours), then collapses to a black hole, and(4) a stable neutron star that would in theory, after allowing for some time for the remnant tostabilize, emit continuous GWs (82 ). In Figure 10-1 the possible remnant scenarios are shown106

schematically as a function of the initial mass of the binary, as well as their associated timescales.

Figure 10-1: Diagram taking from (82 ) showing the stages of a binary system’s evolution as
a function of the total initial mass of the system. The total mass itself, and its relation to the
EoS-dependent MT OV affects the product of the merger, and each remnant has an associatedlifetime.
It is possible that electromagnetic observations could give us information on the
post-merger remnant; however, in reality it is often ambiguous. For example, the kilanova(r-process) model cannot fully explain the spectra that we see, so hybrid models have been
considered: early EM emission is powered by radioactive decay of r-process elements, andemission at later times is powered by the spindown of a long-lived remnant (83 ). Additionally,it has been shown that some stiffer equations of state can produce magnetic fields consistent
with the electromagnetic observations (whereby measuring the period and rate of change ofpperiod of a pulsar, we can calculate the magnetic field since B ∝ P Ṗ - see Equation 3–6).

107

This relation comes from equating the rotational energy of the star with the power due to arotating magnetic dipole- see Section 3.2.3.
Based on the derivation in 3.2.3, we can also see immediately that:Ω̇ = −

B 2 R6 sin2 α 3Ω
6πc3 I(10–1)

A rotating magnetic dipole emission mechanism gives us a braking index of n = 3.
Additionally, if the remnant was a black hole, fall back accretion would occur and couldin theory power the electromagnetic emission. Fall back accretion happens when neutrino
emission causes ejecta matter to cool very quickly so that matter falls back into the remnant(the direction of matter is reversed by rarefaction wave), if the velocity of the matter is notmuch larger than the escape velocity (84 ). But this has been show to only explain emission at
early times, where a spinning-down neutron star can explain the full spectra: the magnetic fieldrequired in O(1012 ) G, but the ellipticity is high O(10−3 ).
No single electromagnetic emission model can explain the observed spectra in the
UV/Visible/ infrared bands (30 ). However, detecting a neutron star as a remnant of thismerger could provide many insights, despite the ambiguous electromagnetic evidence. We
know the distance just from GW observations: the mass ratio can be determined fromthe time/frequency evolution and the strain on the detector tells us about the strength of
the signal. The location is know somewhat well with just LIGO/Virgo but extremely wellwith electromagnetic observations. But what we still don’t know is: what is the EoS of theneutron star, the percentage of their energy that goes into GWs, how deformed they are,their moments of inertia, and if there is a difference between pre- and post-merger equationsof states. We can only understand these things through GW observations: if we are able tomeasure the mass of the neutron star remnant, we can then constrain the product Iz . We
describe our search procedure for a remnant of GW170817 that was done jointly with threeother pipelines (30 ).

108

10.3

Parameter Space Explored

Since we knew the location and the exact time at which to search, we were able to explore
a relatively large range of the physical parameter space of a neutron star. We search acrossbraking indices ranging from n = [2.5, 7], spindowns from [1/256, 1/4] Hz/s, initial frequencies
f0 = [100, 2000] Hz, in three different configurations (three different TF F T = 2, 4, 8 s startingon average 1, 2 and 5 hours after the merger). We cover different frequency and spindown
ranges for each configuration, but we always cover the same range of braking indices.
We start the search 1 hour after the merger for a few reasons: (1) we expect that
immediately after the merger, the physics will be very complicated and not well-describedby a power law frequency/time evolution, (2) the spindown of the source should be smaller,
which means that we can use longer TF F T and observe for a longer time, and (3) we canprobe signals with high initial spindowns and frequencies that may have moved into ourfrequency/spindown window within this one hour delay.
We construct which time/frequency bands to analyze for each value in our grid on thebraking index on the basis of the spindown and the starting frequency of a hypothetical signal.
Essentially, we assume that if there is a signal that starts at f0 = 500 Hz with a maximum
spindown of 1/4 Hz/s, we determine the amount of time it takes for the signal’s spindown
to reach a value of 1/42 = 1/16 Hz/s, which would be when we can change the TF F T to be
4 seconds. We record the amount of time this takes, and the frequency band that would becovered by the signal. This is illustrated in Figure 10-2, as well as how different braking indicesaffect the bands that we analyze.
Each time/frequency band for each initial frequency we consider at each braking index
is slightly different. This means that when we would like to analyze which portions of theparameter space we are sensitive to, we must account for the fact that we will have different
TF F T and durations for each class of signal we analyze.
Figure 10-3 shows which portions of the parameter space we are most sensitive to 1
hour after the merger, where “most sensitive” refers to a hypothetical signal starting with the109

Figure 10-2: The spindown of two signals with f0 = 500 Hz with different braking indices(black/top: n = 2.5; blue/bottom: n = 7) as a function of time. The horizontal lines correspond to spindowns f˙ = 1/TF2 F T , f˙ = 1/(2TF F T )2 , and f˙ = 1/(4TF F T )2 . When a signal
reaches a small enough spindown (intersects one of the horizontal lines), we can then analyze
it for a longer duration with a higher TF F T . Within each rectangle, the times and frequencies
that will be analyzed with a given TF F T , which will be different depending on n, f0 and f˙0 . We
explore all signals whose paths lie between the curves corresponding to n ≥ 2.5 and n ≤ 7.
specified frequency/spindown as being fully contained in the frequency band we are analyzing.
The white space (“the holes”) that are visible in Figure 10-3 are indicative of signals thatare for the most part contained in the frequency band we analyze, but spin out a little. We
are still sensitive to this signals, but not as sensitive as if we extended the frequency band
analyzed. But, we can still recover these signals, since they are contained within the map that
we analyze.

110

Figure 10-3: The frequency/spindown space that is explored by analyzing different portionsof frequency and time with different TF F T for n = 5. The durations Tobs of the analysis arecolored, and the three separate blocks correspond to TF F T = 2, 4, 8 s moving from right toleft.
Figure 10-4 shows the corresponding source parameter space, where we have “backtracked” the GW spindowns/frequencies in the Figure 10-4 to determine the true initialfrequencies/spindowns. It is clear that we can probe very high frequencies and spindowns thatare not normally accessible to us because of the noise present in the detector, and because of
the amount of sensitivity we would lose to try to make TF F T < 1 s. The so-called spindown
time scale is colored in this figure, defined as:

τ=

1
f0
=
n−1
kf0 (n − 1)f˙0 (n − 1)

111

(10–2)

Figure 10-4: The source frequencies and spindowns that we are most sensitive to in our search
scheme (Figure 10-2), with a choice to start analyzing the data one hour after a merger, forn = 5. The spindown timescale τ is colored. This plot includes signals that lie in the “holes”
of Figure 10-3.
To determine these initial parameters, we try many different initial frequency/spindown
combinations to determine which parameters will fall within one frequency bin of thefrequencies given in Figure 10-3 after allowing this simulated signal to decay for a certain
amount of time.
10.4

Candidate Selection

In each Hough map, we must select a certain number of candidates, chosen as acompromise between losing sensitivity, decreasing false alarm probability and decreasing
computational cost of the follow-up. A candidate is defined as having an f0 f˙0 , n, longitude
and latitude. In the Generalized Hough pipeline, we choose candidates with the highest number

112

counts/critical ratios in each box of the Hough map, where the size of each box is determinedby the number of candidates to be chosen. To give an idea, if we only select a small, O(20)candidates, from the map, the boxes are large; if we select O(200) candidates, the boxes are
smaller.
We employ this “box selection” method because of the non-uniformity of the noise. We
found that if we simply select candidates in the standard way (in rectangular strips), we are
predisposed to select only candidates in the highest frequency (low x value) space with low
spindown (low k). With the box selection method, we can select candidates with more uniform
f0 and k, which are both physical parameters. Additionally, we select the second strongest
candidate in the box if it is sufficiently far away in k. This is done so that we are not blinded
by noise artifacts.
10.5

Follow-up Procedure

Since we constructed the search in three configurations, we would first like to compare
candidates recovered in each one. Recall that each configuration is done at a different averagetime after the merger. What we can then do is do coincidences between candidates that
we recover in each configuration. We simply take the recovered parameters in, for example,
the TF F T = 2 s configuration, allow the signal to spin down for a certain amount of time(depending on how long the analysis window was), and then compare the signal’s parameters(x0 ,k and n) to the candidates’ parameters recovered in the second block (TF F T = 4 s). The
coincidence distance iss
d=x0,2 − x0,1
δx0

2


+

k2 − k1
δk

2


+

n2 − n1
δn2,(10–3)

which is the same as Equation 9–5, except that we now can consider a variation in the braking
index n. We chose to impose the constraints that n and k should not vary, by definition.
However, it is entirely possible to look for signals that do have variation in these parameters,
and if we want to explore signals with varying braking indices, we would need to account for
this.

113

After doing coincidences between TF F T = 2 s and TF F T = 4 s, we were left with somecandidates using a coincidence distance of d = 3 bins. Then, we did coincidences between thesurviving candidates and TF F T = 8 s, but no candidate survived. Therefore, we can conclude
that there was not a detectable signal present in all three portions of the parameter space weanalyzed.
Our search returned about 500 coincident candidates (too numerous to list) with anaverage critical ratio above 5. We chose a threshold of CR = 5 because it corresponds toa false alarm rate of < 0.001%, and this is the value typically used in CW searches. We
then looked at the critical ratio recovered in Hanford and Livingston, and rejected coincident
candidates whose critical ratio was more than 4 times different in one detector over the other.
We expect a true GW signal to have a difference in critical ratio that is proportional to (1) theorientation of the detector with respect to the source location and (2) the sensitivity difference
between the two detectors. For this search, most of the candidates were vetoed because thecritical ratio difference was too great. We only performed the follow-up on one candidate.
The procedure for the follow-up is the following: we correct the raw time-series data for thephase evolution of the candidate we are following up as described in 9.4 in each detector. We
then decide how much to increase the TF F T by (we have chosen a factor of 2, and would
have increased it more if the candidate survived this phase of the follow-up). Now, we create
another peakmap for each detector and perform the original FH on it. This is because we have
hopefully removed most of the frequency variation from the signal, therefore the signal shouldbe almost monochromatic. After performing the FH, we select candidates in the new maps,
and do coincidences. In the case of the search, the candidate did not survive; therefore, allcandidates were vetoed.
10.6

Upper Limits Procedure

When a search doesn’t find anything, we perform upper limits, that is, we ask thequestion, “what is the minimum amplitude of a signal we could have found?” Or, “how far
away could we have seen?” Now, there is a key difference between upper limits and sensitivity:

114

the ULs are computed using the critical ratios of candidates in each frequency band that
we received for a search. That is to say, we essentially perform the same procedure as withsensitivity, but we require that a given h0 , that the critical ratio of 90% of a detected signal
has to be greater than the maximum critical ratio in the frequency band in which the search
was performed. So we compute upper limits in the following way: we inject 100 signals every
50 Hz in each configuration (TF F T = 2, 4, 8 s) at 1, 2 and 5 hours, to be reflective of thesearch that was actually done. We only inject signals that each configuration analyzed: forexample, we did not use TF F T = 8 s to analyze frequencies above 1200 Hz 5 hours after themerger, so we do not injects at f = 1200 Hz in that configuration. We include in Figure 10-5
our upper limits on the amplitude h at 1, 2, and 5 hours after the merger. The maximum
possible spindown was used when injecting these signals, because our theoretical sensitivity
estimates show that the sensitivity is worse for signals with higher spindowns, so our ULs areconservative. We also randomized over cos ι, and injected signals with frequencies f to f + 50
Hz randomly.
Additionally we plot the upper limits obtained by all four pipelines, given in (30 ), in Figure10-6. Each method had its own way of doing injections (different signal durations, parameters
over which the analysis was randomized, etc.); however, each performs similarly. It is evident
that our sensitivity was not good enough to see a remnant of a neutron star merger: our reachis between 0.5-1 Mpc and GW170817 was at 40 Mpc. This analysis would have been more
successful if Einstein Telescope (68 ) or Cosmic Explorer (69 ) were online.
10.7

Prospects

No signal coming from a remnant of GW170817 was detected by our pipeline or others.
But we knew this would happen before doing the analysis: the point was to be ready for future
BNS mergers, especially if one happens closer to us. Moreover, third generation detectors are
predicted to have high enough sensitivity to detect signals at 40 Mpc away: in fact, if Einstein
Telescope was already built and on during GW170817, we may have been able to detect asignal from the remnant, if it was a neutron star.

115

Figure 10-5: Upper limits placed in each configuration. Each point represents 100 injections
done with the maximum possible spindown that an FFT can handle: f˙ = 1/TF2 F T between f
and f + 50 Hz. The durations of signals analyzed are also different, ranging from about 1 hourto 18 hours. Higher frequencies tend to have longer signals, since the spindown decays quicker
for lower frequencies. Note that the frequency on the x-axis is the frequency referenced to 1, 2,
or 5 hours after the merger (TSF T = 2, 4, 8 s respectively) The behavior of the upper limits is
consistent with the sensitivity estimation in white noise, given in Figure 9-11.

116

Figure 10-6: Upper limits placed by each pipeline that ran a search for a remnant of
GW170817. Our method is labelled FreqHough. The methods all perform comparably; however, the merger happened at a distance that is beyond the reach of current analyses. If future
detectors were online, we would have had a better chance of seeing a GW signal.

117

CHAPTER 11
TWO-DIMENSIONAL FAST FOURIER TRANSFORM LOW PASS FILTER
11.1

Introduction

If we think about the time-frequency map as an image, the r-mode waveform is a continuous, low frequency signal. The noise that surrounds the waveform is quite discontinuous,and can be expressed as a high frequency signal (many Fourier coefficients are needed tomodel the rapidly changing frequency of the signal). Hence, a low-pass filter would allow us
to reduce the effect of the noise by filtering out the high frequency noise and enhancing the
low-frequency r-mode signal. Our implementation of the low-pass filter is similar to matchedfiltering, which is a data analysis technique used to enhance the SNR of potential gravitational
wave signals. However, matched filtering requires O(106 ) templates of the signal because if theparameters of the filter are even slightly off, we miss the signal. Our filtering technique obtains
an enhancement of the quality of the image even if the signal we are trying to find does notmatch the signal with which we are filtering. In Section 11.2, we describe the two-dimensional
FFT image enhancement technique. Section 11.3 showcases some preliminary results we haveobtained, showing that this technique can be used to enhance images. We then describe thefuture work that needs to be done in Section 11.4.
11.2

MethodTo simulate the signals, we inputted amplitudes based on Equation 3–34 and frequencies
following 3–33. We then integrated Equation 3–33 calculate the phase evolution of the signal.
After this, we created a time/frequency spectrum by taking a one-dimensional FFT of theamplitudes- see the top panel of Figure 11-1.
The low-pass filter takes the form of a two-dimensional Fast Fourier Transform (2D FFT)of a simulated time/frequency behavior r-mode signal. The 2D FFT is defined as:A(u, v) = F (a) = ΣΣa(t, f )e−i(ut+vf )118(11–1)

where a(t, f ) are the amplitudes of the signal in the time/frequency plane (the colors in the
top panel of Figure 11-1), t is time, and f is frequency. u and v are the ’new frequencies’in Fourier space and don’t have any physical meaning. The notation F (a) denotes the twodimensional FFT of an image a, which is a matrix of amplitudes at each time/frequency.
When a low-frequency signal is transformed from the time-frequency plane to the Fourier
space by a 2D FFT, it becomes concentrated at the corners of the new image, and the highfrequency signals in the original image are spread throughout the new image. To visualize thistechnique, we center the information that is actually present in the corners of the image, see
the bottom panel of Figure 11-1).
To do the filtering, we use the following equation:

Y = F −1 (|F (X)| × |F (M )|)(11–2)

where the bars |..| refers to taking the modulus of the image, Y is the new image, X is theoriginal image that we want to enhance, and M is the filter (the ideal model of the signal).
We calculate the following factor before and after applying the filter:
An =

Σh2
M

(11–3)

where Σh2 represents summing the square of the amplitudes along the path of the signal intime-frequency space, and M is the median of the image. To calculate this factor after theenhancement, we normalize the image, apply a threshold of 0.9 to the image and then sum thesquare of the surviving time/frequency amplitudes. To calculate the enhancement, we thendivide the two factors.
11.3

Preliminary Results

We have found that signals spinning up and spinning down in the time/frequency plane
carve out an “ellipse” in the Fourier space, when looking at the centered image; hence, by
combining multiple signals, we can cover much of time/frequency parameter space with only afew signals- see Figure 11-1. This enhancement resembles a “matched filter”, in the sense that119

Figure 11-1: Top: time/frequency map of the ideal signals used as a filter. Bottom: the“ellipse” carved out in the new parameter space after performing a 2D FFT on the top image.
we are using our knowledge of the waveform of the signal; however, unlike matched filtering,
we do not need to use the exact signal to obtain an enhancement (85 ).
Figure 11-2 shows an example of what the filtering technique does to noisy images. We
see that the filter used (the right panel) does not match the signal we are searching for (theleft panel). Yet after the filter, the signal looks better than it did before (middle panel), and isenhanced based on our calculation of the enhancement ratio. Additionally, we have quantified
the level of enhancement that we receive using six filters that do not match the signal. See
Figure 11-3.
We applied a filter containing nine different signals to the case where both an r-mode
signal and a comparatively stronger glitch are present together. See Figure 11-4. We can see

120

clearly that the r-mode signal is enhanced relative to before, and the glitch remains at about
the same strength.

Figure 11-2: Left: r-mode signal before the filter is applied; middle: r-mode signal after thefilter is applied; right: filter applied to transform left image to middle image. We see in the
middle image that the r-mode signal is clearly enhanced, even though the filter we used did notmatch the signal we were looking for.
11.4

Outlook

The input SNR of the image determines how well the signal can be enhanced: if the signal
is too strong, the signal is not enhanced so much, but the noise is, lowering the SNR. If thesignal is too weak, then the filter cannot help us enhance the signal. This problem has beenfurther explored and solved in (86 ) through a different calculation of enhancement. However,there is a range of input SNRs for which the signal can be enhanced by a factor of 2-8 basedon how we calculate SN Rn . In the future, we would like to remove the threshold described in
11.2 and sum just along the line defined by only the signal to calculate the SN Rn . We havebegun developing“masks” to perform this summation, whereby we are able to extract just thesignal amplitudes corresponding to the actual signal frequencies from the image.
It is clear that a combination of filters enhances the image but the enhancement is low
compared to filtering with one signal. Hence we anticipate using a combination of filters
to obtain a coarse estimation of whether a signal is present in the data, then follow up thisprocedure with one-signal filters.
Additionally, filtering using Equation 11–2 is the simplest filter to use but not necessarily
the best one, so we plan to experiment with other types of filters to see if they further enhance121

Figure 11-3: Enhancement of f0 = 1400 Hz, α = 1 signal using different filters. The colors correspond to the different filters used. Red: f0 = 1500 Hz, one signal; orange: f0 = 1400 Hz,
1 signal; green: f0 = 1300 Hz, 1 signal; light blue: f0 = 515-1500 Hz, 4 signals; dark blue:
f0 = 500-1500 Hz, 9 signals; brown: f0 = 1000 Hz α = 0.1, 1 signal .

Figure 11-4: Left: a weak r-mode signal obscured by a glitch; middle: filtered image, showing
r-mode signal has been enhanced; right: filter used to transform left image to middle image.

122

the images. For example, we see in Figure 11-3 that the 1000 Hz signal does not reallyenhance the 1400 Hz signal, mostly because the shapes of the signals in the time/frequency
plane and the new Fourier plane after the 2D FFT are too different. However, some preliminary
tests indicate that a more robust filter can be created such that signals that appear sufficiently
different in time/frequency can still give enhancements. This work has been extended in (86 ).

123

CHAPTER 12
MACHINE LEARNING ALGORITHMS
12.1

Introduction

In the LIGO/Virgo collaborations, many methods have been developed to search fortransient signals that last O(s) and for continuous waves (CWs), quasi-infinite signals fromspinning neutron stars (87 –89 ). For intermediate duration signals lasting O(hours − days),
a cross-correlation-based method was developed some years ago (90 , 91 ). But now, twocategories of these types of signals exist, requiring two different analysis procedures: (1) CWlike signals with still slowly varying frequencies (92 , 93 ) and (2) signals with large frequency
variations. Only recently has the scientific community taken an interest in running a search forsignals in the second category (30 ), which required the development of new methods, bothmodeled (31 , 32 ) and unmodeled (33 ). We expect that recently born neutron stars (i.e. after
a BNS merger or supernova) will emit strong gravitational radiation on this timescale (82 ).
As the LIGO/Virgo detectors become more sensitive in the coming years, we should be
able to not only detect more binary neutron star mergers, but also CWs from neutron stars.
The all-sky and directed searches for CW signals are computationally expensive, since one
must search the entire data set in every position in the sky. Even for long duration transient
signals, signals that last for O(hours − days), modeled searches such as those based on the
Generalized Frequency Hough (GFH) (31 ) and the Adaptive Transient Hough (32 ) consume
weeks of computing time, especially at intermediate spindowns, and require that the signal
follow approximately a power law.
We therefore explore machine learning as a way to cut computation costs without losing
sensitivity. The key idea of machine learning is that one gives it many examples of somethingthat one wants to be characterized (hand-written digits, cats/dogs, gravitational-wave signals,etc.), and the algorithm learns how to distinguish between the different possibilities.
After sufficient training on different GW signals, MLAs can generalize and can perform
well relative to existing methods. We describe the type of methods we have employed to124supplement the GFH: Artificial neutral networks in Section 12.2, convolution neural networks in
Section 12.3 and random forests in Section 12.4.
12.2

Artificial Neural Networks

Artificial neural networks (ANNs) (94 ) are a biologically-inspired machine learningalgorithm (MLA) that takes as input an image and outputs a statistic that determines if theimage contains a pattern or not. It consists of an input layer, which is a collection of images as
in Figure 12-1 that is flattened, one or more hidden layer with a certain number of nodes thatallows nonlinear combinations of the pixels in the image, and an output layer. The hidden layer
essentially assigns weights to each pixel in an image, and the weights are constantly updated
throughout the training process as more examples are added. At this level, an activation
function (hyperbolic tangent, softmax, or logistic) is used to introduce non linearity into the
network. Typical activation functions include sigmoid (95 ):σ(aj ) =1,
1 + e−aj(12–1)

where aj is the output of a layer in the network, softmax:eajp(aj ) = P aj ,e

(12–2)

f (a) = max(0, a)(12–3)

and rectified linear unit:

The input data takes the form of a data matrix X, where each row corresponds to an
unrolled image (an A x B-sized image becomes a 1xAB vector), and AB corresponds to thenumber of dimensions in the time/frequency map.
In standard neural network architectures, there is first the “forward pass” in the network(96 ), where a guess for whether an image has a signal or not is computed. We show an

125Figure 12-1: Injected signal with f0 = 1040 Hz, h0 = 7.17 × 10−23 . This image (a matrix) is
flattened into a row vector and given as input to the Artificial Neural Networks. We will see
this image again in Figure 13-1.
(1)

example for one input image x1j . To begin, the weight matrix wkj is multiplied by the unrolled(1)

image and a bias node wk0 is added, and then this output is fed into Equation 12–1.
If we have N images with dimension d = AB, so an input matrix of N x (d+1) dimensions,(1)then the first weight matrix wij has size (nhidden + 1) x (d + 1). Matrix multiplication is
performed between the input matrix and the first weight matrix, and then an activation
function is applied to smooth the output. Each hidden layer in the network will have anassociated weight matrix.
In the simplest case of one hidden layer, we would now have a matrix with size N x(2)nhidden + 1. The final step of the network requires another weight matrix wij of size nhidden + 1

126

x nclass , where nclass = 2 (signal or noise). Equation 12–4 is repeated again for this weightmatrix, and the output of the network is a N x nclass matrix, where all reference to the hidden
layer has been removed. Now, this matrix is fed into the softmax activation function 12–8.
Let us consider one training example. In equations, the forward propagation (97 ) of the
networks looks like this:From input to hidden layer 1:dX

x01k = σ

!
(1)

(1)

wkj x1j + wk0(12–4)

j=1

and from hidden layer 1 to output:!

dX

p1k = p(2)

(2)

wkj x01j + wk0(12–5)

j=1

where p1k is the outputted probability after the forward propagation. The error is thencalculated between the returned and expected values for the image, and through backpropagation (98 ), the error is attributed in different amounts to each of the nodes in thehidden layer. The error per example δk is:δ (2) = ak − pk

(12–6)(1)

Then, this error is then converted into error on the weights wij :(1)δ (1) = wij δ (2)

∂
0
σ(z1k
)
∂aj

(12–7)

0where z1k
is given by the term inside the sigmoid function in Equation 12–4.

Changes to the weights are done according to which nodes were most or least responsible
for the error at the output layer. The goal is to find the optimal weights such that the error
on the validation set (a set of data used to avoid overfitting) is smaller than a certain value. A
basic neural network can be seen in Figure 12-2.

127

The training data are typically divided into smaller sets, called batches, and each batchis pushed through the network. One can choose the number of times that the network should
see the full dataset, which is called the number of epochs. If the network is trained too
much, overfitting occurs, which means that the network learns the training data too welland can therefore not generalize to new situations. To prevent this, an early stopping criteria
can be employed: if the sum of the mean squared error between the outputs of the forward
propagation within each epoch continues to decrease (until a threshold on minimum error isreached), we continue to train the network (47 ).

Figure 12-2: This diagram shows the layout of an artificial neural network. The left-hand panel
shows forward propagation through the network, while the right panel shows back propagation. The Θ1,2 are matrices of weights that transform the input images into a guess as towhich class an image belongs. g is the sigmoid function. The ∆s represent the errors that arepropagated back through the network after a guess is made in order to adjust the matrices ofweights. This was taken from (99 ).
12.3

Convolution Neural Networks

Convolutional neural networks (CNNs) are essentially ANNs but instead of using thetypical matrix dot product to pass between layers, the image is convolved with a weight matrix
called a “kernel”. The kernel moves over subsets of the image and returns a measure ofoverlap between the kernel and the image, before this output is fed into an activation function.

128

The principle of determining the ideal weights for ANNs at each layer remains the same
for CNNs, except now one can use many convolutional layers, different activation functions, arandom pruning of the network, etc. in order to create a better classifier.
The motivation to use CNNs in gravitational wave physics come from their quick andsuccessful implementation in other fields. For example, AlexNet (100 ) was able to quickly train(on a Graphics Processing Unit (GPU)) to classify images into 1000 separate classes with amuch smaller error rate than previous algorithms. This work revolutionized the field of artificial
intelligence, and made deep learning possible on a reasonable timescale. CNNs and other deeplearning algorithms can be found in almost every field now, from speech recognition (101 ) to
medical imaging (102 ).
We implement a CNN shown to work for both supernovae (103 ) and CW signals. The
key aspects of this CNN are 5 convolution blocks, described in detail below, each containing
zero-padding, convolutions, ReLU and max pooling steps, and 2 fully-connected and dropout
layers followed by a softmax layer at the end.
1. Input (N x 62 x 37 x 1)
2. 5 convolution blocks• ZeroPadding 2D• Convolution 2D (3x3x16)• ReLU• MaxPool(2x2)
3. Flatten
4. 2 blocks of:• Fully-Connected (16 nodes each)• Hyperbolic Tangent• Dropout
5. Fully-Connected (2 nodes)
6. Softmax

129

7. Output
The zero-padding layer ensures that the size of the image stays constant throughout
the convolution block by surrounding the images with zeros. This is necessary because themaximum pooling layer takes the maximum value of every non-overlapping 2x2 block of theimage, effectively cutting the size of the image in half each time.
The flatten layer unrolls an H x W x N kernel into a vector of length H · W · N .
In the fully connected layers, each output of the previous layer is used as in input to eachof the 16 nodes in this layer. These layers are essentially what each hidden layer in an ANN
is. The last fully-connected layer in our CNN outputs a vector whose size is just the numberof classes in which we are trying to classify: 2 in this paper (signals or noise). The values
correspond to the CNNs’ confidence in classification of the image.
The dropout layer is used to prevent overfitting of the network. It essentially zeroes theactivation of randomly chosen hidden nodes of the model at each iteration of the learning
phase, akin to pruning a decision tree.
To convert the “confidence” of the CNN into a probability distribution Pi , the final layer
is a softmax activation function:e xi
Pi = PN

n=1 e

xn

,

(12–8)

where N is the number of classes (signal or noise) and xi are the outputs of the fullyconnected layer (essentially the result of convolutions of the input image with the kernel).
12.4

Random Forests

In machine learning, decision trees are a predictive method for classifying data. They are
created to make a series of decisions based on the parameters of the input data. The trees
grow in size and complexity in order to fully separate data points that belong to one class oranother. However, individual decision trees are likely to overfit the data. Therefore, we employ
a “forest” of trees, where each tree is as independently generated as possible in order to obtain130

a monotonic increase in classification accuracy as more trees are added to the forest (104 ,105 ). We describe in Subsection 12.4.1 how we construct our random forest.
We can see in Figure 12-3, taken from (106 ), an example of a random forest composedof 5 decision trees. The two classes (signal and noise) are defined by the two colors light blueand black, and at each numbered node, a decision is made to split the training set, given in
the trunk of the tree, into different leaf nodes. For our purposes, each of these nodes will befurther split, until each leaf contains only blue dots or black dots.

Figure 12-3: A random forest composed of 5 trees (taken from (106 )). The training set, in thebark of the tree, is composed of two classes: light blue dots and black dots, representing noiseand signal. At each numbered node, a decision is made, and the data are split. In our random
forests, we will actually continue this figure until all leaf nodes contain only instances of one
class, so each leaf will have exclusively blue dots or black dots.

131

12.4.1

Constructing a Decision Tree

To make a decision tree, input data and labels are required, as in the case of ANNs/CNNs.
But one additional piece is required: a splitting function. The splitting function determines inwhat way to separate the training data into the two classes we would like to identify, and thishappens at each node of the tree.
To begin, we take the entire (down-projected) data matrix X and give it to the “rootnode” of the tree. The splitting function is then applied, and X is split into two groups
X1 and X2 and y is split into y1 and y2 based on some criterion (described in a subsequent
section). Then the two individual groups are again given to the splitting function, and morenodes are created until the classes are fully split, that is, until each leaf node contains elementsfrom only one class.
12.4.2

Constructing a Forest with the Random Subspaces Method

The random subspaces method (RSM) was proposed in (104 , 105 ) to generate independent decision trees that will have 100% accuracy on their corresponding training sets. Previous
implementations had relied on randomly selecting entire training examples and giving them todifferent trees; however, it is not guaranteed that each tree will classify the entire training set
correctly. In the random subspaces method, each feature vector is down-projected, meaning acertain number of features are randomly selected and the rest of the features are set to 0 in
the unselected dimensions. For example, if we want to create a subspace of 8 dimensions, andwe have a 16 dimensional feature vector:

f = 1 2 3 4 . . . 16
one tree will randomly select

f = 5 2 6 8 3 1 16 15132

and another tree in the forest will select

f = 4 3 1 5 11 9 15 14
All training examples (all feature vectors) are down-projected into the desired subspacedimension before being sent into the tree. Each tree will be the same number of dimensions,but will live in randomly different subspaces. It was shown in (105 ) that these trees will havedifferent generalization errors and be accurate on the entire (down-projected) training set.
12.4.3The Splitting Function: Central Axis Projection

Central axis projection (CAP) (104 ) is a fast and feature-independent algorithm that has
been used to make splits in decision trees. Many conventional splitting techniques randomly
select a single feature to split at each node, which does not allow us to take into accountany relationships that exist between the individual peaks. To begin, all training examples
are simultaneously given to the root node. The classes whose means are farthest apart by
Euclidean distance are found, and then a line is drawn between the coordinates of those means.
All the data are then projected onto this line, called the central axis, following the standard
equation (projecting ~a onto ~b):Pba =~a · ~b~b|b|2(12–9)

We start at the beginning of the line, and take discrete steps, placing a hyperplane ateach step that is perpendicular to the line. We then count the number of points on eachside of the hyperplane, and choose the hyperplane that minimizes the number of points thatare not on the side of the hyperplane where the majority of that class’s points fall. This
choice of hyperplane causes the split in the dataset, separating some points correctly andsome not. Each daughter node of the tree will then perform the same selection of hyperplane,
and the process will continue until the tree is fully split. For example, if the feature space133

is two-dimensional, a line is drawn perpendicular to the central axis at each step; in a three
dimensional space, a plane is drawn- see Figure 12-4 for an example.

Figure 12-4: Left: original data in a 2D subspace with central axis drawn between meansof the two classes. Right: all data in this 2D subspace have been projected onto the central
axis- its endpoints are shown in red. In a finite number of steps, we will move from the first
endpoint to the second one, determining first which side the majority of each class’s points fallon. Then we minimize the number of points on the wrong side of the hyperplane.
12.4.4

Voting and Predicting

Once many trees have been created in a forest, we need a way to combine their estimateson new testing data, that is, their classification if the new data have a gravitational wave signal
or not. Each testing feature vector is down-projected into the corresponding subspaces of eachtree in the forest, and since the trees are fully split, each feature vector will end up at a leaf
node corresponding to whether it is a signal or not. After the feature vector is sent down eachtree, we count the number of votes for signal and the number for noise. In general, the classwith the greatest number of votes wins, but in this case, it is simply that the class with the
majority of votes wins.

134

CHAPTER 13
CHARACTERIZATION OF CONVOLUTIONAL NEURAL NETWORKS AND RANDOMFORESTS
13.1

Introduction

MLAs have been shown to achieve extremely good sensitivities in both white and realnoise, comparable with matched filtering and some unmodeled methods (103 , 107 –111 ).
However, what is lacking from these studies is a quantitative analysis of (1) how to train, (2)how much to train, (3) robustness towards signals on which the networks were not trained,(4) differences in architectures’ effects on detection efficiency/false alarm probability, and(5) applying the networks to a real search. Moreover, these algorithms have only focused on
BBH signals or supernovae, and have not yet been expanded to search long duration transient
signals from isolated neutron stars. We address these points for ANNs, developed in (47 )and CNNs, developed in (103 ). We characterize the performance of the networks under avariety of conditions, and propose a possible design for a search using CNNs to provide triggers
that would be followed up by another method. We also compare the ANNs and CNNs to an
established method to search for long duration transient signals, the GFH. The work in thischapter is based heavily on (112 ).
In Section 13.2, we define what are the inputs to the algorithms. Then we quantify how
much we should train as a function of source parameters in Section 13.3. We experiment withusing a forest of ANNs in Section 13.4. Next, the robustness of the methods towards different
signals on which they were trained is quantified in Section 13.5. The detection statistic we use
to determine success of the algorithms or not, as well as the effects of non-stationary noise on
the algorithms’ performance, is described in Section 13.6. The error associated with efficiency
and false alarm probability is empirically determined in Section 13.7 for a set of O(1000)
models. A comparison is made between ANNs/CNNs and the GFH in Section 13.8. Then we
design a possible hierarchical search using CNNs as a way to quickly generate triggers to befollowed up by the GFH in Section 14.2. Finally we do a quick analysis of decision trees’ ability
to detect GW signals in Section 13.9.

135

13.2

Construction of Input Images: Reduced Resolution Time/Frequency Maps

The pattern on which we train the networks is in the time/frequency domain, so we first
construct a time/frequency map by fast Fourier transforming each chunk of data (interlacedby half and windowed). Then, we estimate the noise power spectral density by using an
auto-regressive spectrum to whiten the data (77 ). Afterwards we plot the result in thetime/frequency plane, as shown in the left panel of Figure 13-1. Ideally, the more informationthat we can give to the networks, the better they can learn to identify and distinguish signals
from noise. However, there are too many pixels on which to attempt to train the ANNs/CNNs,
so we reduce the resolution of each time/frequency map in an analogous way as done in (47 ).
We split the time/frequency map into different square blocks and take the maximum value thatappears in each 16 x 16 block. The size of this box was chosen based on an optimal resolution
reduction factor determined in (47 ). Additionally, we apply a threshold to each local maximum
we select: if it is less than a threshold of 2.5 (this choice of threshold is described in (78 )),we put a 0 at this point in the reduced time/frequency map. The value 2.5 is equivalent toselecting pixels that are at least 2.5 standard deviations away from the mean level of noisein the map. Selecting a local max and applying a threshold gives us the so-called peakmap,
constructed in a similar way as described in (78 ), but our method of resolution reduction
ensures that the peakmaps have the same size, which is necessary for the MLAs to work well.
In Figure 13-1 we show the before and after resolution peakmaps.
13.3

Training Efficiency

To run a real search, we need to train the CNNs to recognize both signals and noise.
However, when we train on noise maps, we have instructed the networks to see thesetime/frequency maps as not having signals. We would therefore have to throw out this data,
so we wish to minimize how much data we need to use for training. Thus we quantify at what
point additional training of CNNs no longer adds appreciably to the detection efficiency. Figure13-2 shows a CNN trained with different numbers of injections and noise maps on real data.
We apply a threshold on the probability pthr outputted by Equation 12–8: pthr = 0.9, chosen136

Figure 13-1: Left: injected signal with f0 = 1040 Hz, h0 = 7.17 × 10−23 . Right: reduced
resolution peakmap by a factor of 16 on each axis. This reduced peakmap is the input to themachine learning algorithms. When we reduce the resolution of the image, we pick the pixelin each 16x16 block that has the highest value, but if this value is less than a threshold (2.5,
essentially not higher than 2.5 standard deviations from the mean), we zero this pixel.
based on the false alarm probability analysis done in Section 13.6 . Note that pthr = 0.9
does not imply a detection efficiency of 90%; rather, pthr is our detection statistic. Even for
training on ∼ 20000 injections, we see that the detection efficiencies are comparable to trainingon ∼ 100000 injections. This appears to be independent of signal strength. However, whentraining on fewer injections, we are in fact training on fewer noise maps, which means thefalse alarm probability is higher. These results depend on which pthr we use, but they offer aguideline for training a single network on a variety of amplitudes.
13.4

Using a Forest of Neural Networks

A single ANN was able to achieve high detection efficiency at the expense of a high false
alarm probability (47 ). It is therefore necessary to consider ways to lower this false alarmprobability. One way is to train on more noise maps, but there are practical limitations to thissince the length of an observing run is finite, and we do not want to lose too much data. The
solution we propose is to use a forest of ANNs that vote on whether each map contains asignal or just noise. But naively training the entire data set on each ANN will not decrease thefalse alarm probability, so we must consider that each ANN should see only an independent
portion of the image to train on. Down-projecting the image into a subspace allows us to137

Figure 13-2: A CNN was trained on a variety of amplitudes in real noise with a different number of signal injections. It appears that the detection efficiency does not change much as thenumber of training injections is increased. Therefore, it is not necessary to do many injections
to achieve high detection efficiencies. For this plot, pthr=
0.9. Each of 32 amplitudes isrepresented by a line in the above figure.
preserve accuracy on the training set and to generate independent ANNs. While the detection
efficiency does not change that much, the false alarm probability drastically improves- see
Figure 13-3.
13.5

Robustness Towards Signals Not Trained on

Machine learning methods have become useful because they can be trained to identify
certain patterns, and can also generalize to find similar ones, but we must quantify if andhow the detection efficiencies will change. The signals we have described thus far have hadconstant braking indices, which means that the GFH was a strong method to search for them.

138

Figure 13-3: False alarm probability as a function of the number of neural networks we use in
our forest for different amplitudes in white noise. The false alarm probabilities were calculatedwith a probability threshold of 0.7. A forest of neural networks was trained for each signal amplitude. It is clear that even with the addition of 1 or 2 independently trained neural networks,
the false alarm rate greatly decreases. At this probability threshold, a single neural network
trained on the full images has a false alarm rate ∼10 times greater than a forest of neuralnetworks each trained on a randomly chosen portion of the image.
However, if the braking index varies too quickly in time, the GFH cannot recover the signal,
no matter how strong the signal is, because the signal deviates too much from the power law
model in Equation 3–33. Figure 13-4 shows two sets of sensitivity curves for a CNN trained on
signals with a braking index of n = 7. In the left panel, the curves are the result of testing on
injections with fixed braking indices (n = 5 or n = 7). The CNNs performed worse at higher
frequencies, which is consistent with the sensitivity curve of the LIGO detector. In the rightpanel, the CNNs were tested on injections with varying braking indices, with variation δn/δt

139

between [−10−4 , 10−4 ] /s. In total, this variation corresponds to a change ∆n = Tobs δn= 0.1:
δt
the GFH cannot detect signals with such large ∆n. Comparing both plots in Figure 13-4, we
see that despite being trained only n = 7 in the [600, 750] Hz band, the CNNs can detectsignals with or without large δn/δt and have comparable efficiencies across all frequency bands.

Figure 13-4: A CNN was trained on injections with n=
7 and tested on 16000 (500 injections/amplitude) with fixed braking indices (left) and varying braking indices (right)
in different frequency bands. The random variation in braking index is between δn/δt=−4−4[−10 , 10 ]Hz/s. This variation over the observation time is too large for the GFH to recover a signal. A threshold on probability pthr = 0.9 was used in both plots. The sensitivity
curves for the different frequency bands are consistent with the detector’s sensitivity curve(we are less sensitive at higher frequencies). The CNNs can obtain comparable sensitivities inboth the fixed and varying braking index cases. However, even after cleaning both known andunknown noise lines, the threshold on probability must be high for some frequency bands toachieve a false alarm probability of 1%- see Figure 13-5.
13.6

Impact of Threshold Choice on False Alarm Probability

The output of the CNN can be interpreted as a probability that a time/frequency map
contains a signal or not. If we wish to decrease the false alarm probability, we can use athreshold that is higher than pthr = 0.5. To determine an optimal threshold, we test our
network on 500 noise maps generated at random times. As shown in Figure 13-5, to guarantee
a false alarm probability of 1%, we should have pthr = [0.7, 0.95] depending on the frequency
band. We also varied the number of noise maps used for training, and found that ∼ 23 days ofdata (2000 noise maps) was enough to obtain a false alarm probability of ∼ 1% at pthr = 0.9.

140Figure 13-5: We plot the false alarm probability as a function of the threshold on the output ofa CNN (Equation 12–8) for 500 noise maps per frequency band, which correspond to ∼ 11.5
days of new data. The CNN was trained on 2000 noise maps (23 days of data). Different
frequency bands, which have different levels of noise, require different pthr to achieve the same
false alarm probabilities.
13.7

Error Analysis

The CNNs operate as a black box: we control the input, the structure of the networkand the output. Because of this, there is a good degree of randomness in the possible outputsof the network, meaning that each measurement we make has an associated error. Thisrandomness is caused by the dropout layer, which throws away half of the nodes in the
network, and the fact that the network sees randomly shuffled data in different batches in eachtraining session. We therefore train a network 1255 times with the same input data in orderto quantify the range of measurements of both detection efficiency and false alarm probability.

141

We discovered that the errors do not obey a binomial distribution for pthr = 0.9. For the [600,750] Hz frequency band, we train on 38400 injections and 1200 noise maps. The sensitivity
and false alarm probability curves are shown in Figure 13-6.

Figure 13-6: 1255 CNNs were trained on 38400 injections with n = 7 and 1200 noise maps andtested on 44800 (1400 injections/amplitude) and 1400 noise maps. The error bars represent
the 90% confidence level. Left: efficiency is plotted as a function of amplitude. Intermediate
strength signals seem to have the greatest spread in error. Right: false alarm probability isplotted as a function of threshold. The errors are largest around 0.4, and show that higher
thresholds are necessary to truly guarantee low (∼ 1%) false alarm probabilities.
13.8

Comparison to Generalized Frequency Hough

We have compared the detection efficiencies of our CNNs to those of the GFH (31 ), to
AlexNet CNNs and to ANNs, developed in (47 ). In Figure 13-7, we plot the sensitivity curves
for each of these four methods, where we have used pthr = 0.9 for the ANNs and both CNNs.
The ANN performs the worst, not even reaching 90% efficiency, while both architectures of
CNNs and the GFH achieve similar efficiencies. The CNNs are even slightly more efficient thanthe GFH at some amplitudes. However, the CNNs cannot perform parameter estimation andthe false alarm probabilities are about 100 times greater for both CNNs than for the GFH, see
Figure 13-5. A “detection” for the GFH is defined as when the parameters x0 = 1/f0n−1 and kof the loudest candidate in the map are within a distance of 3 bins of the injected signal. For
the CNNs, a detection occurs when the returned probability is greater than pthr = 0.9.

142

We also show in Figure 13-7 the ideal efficiency of a CNN (green curve). This is thebest-case scenario because we train and test the CNN at a fixed GPS time, meaning that the
CNN learns perfectly the noise structure and hence the false alarm probability is 0. The green
curve represents essentially an upper limit on the CNN’s sensitivity. In contrast, the red curve
corresponds to the realistic case, where the CNN is affected by the non-stationary noise of the
LIGO/Virgo detectors. This result shows that we must be aware of the changing noise in thedetector when determining how efficient a CNN will be in detecting power-law signals.

Figure 13-7: We show sensitivity curves for an ANN, the GFH, a CNN and AlexNet CNN
trained on 2000 injections/amplitude and tested on 1000 injections/amplitude in real noise.
pthr = 0.9. It appears that different networks obtain similar efficiencies to the GFH, though thefalse alarm probabilities vary: ∼ 1% for the ANN/CNNs, and ∼ 0.01% for the GFH. The green
curve, corresponding to a CNN trained and tested at a fixed GPS time, represents what theefficiency would be if the noise was the same throughout the run.

143

13.9

Preliminary Characterization of Random Forests

The previous work on MLAs (47 ) resulted in high false-alarm probabilities using neuralnetworks, support vector machine (SVM), and constrained subspace classifiers. Random forests
are advantageous here because the trees separate features from each training example until
no further separation is possible, growing in complexity as the data grow in complexity. It
was shown that a random forest generation package could be used to search for GWs from
GRBs, and that the forest could generalize well when signals of different morphologies areused in the testing set, with good false alarm rates (113 ). Moreover, random forests, andspecifically the random subspace method (explained later), allow us to use all the data in thecreation of each tree- no resolution reduction of the initial time/frequency maps is needed.
The number of subspaces available to down-project into grows with the dimensionality of thedata (2d , where d is the dimensionality of the data). As we increase the number of trees in
the forest, we are exploring all parts of the feature space in a random, non-redundant way.
Additionally, the random forest can be generalized to use other types of splitting functions, forexample, a full SVM at each node. The splitting function that minimizes computation time tocreate, store (and predict with) each tree, and that maximizes detection efficiency, could be
investigated. Random forests provide a convenient way to explore different types of MLAs as
splitting functions while preserving generalization accuracy.
We have trained and tested our random forest algorithms on injected signals in white
noise with the following parameters: durations about ∼ 1700 s, f0 = 777-945 Hz, α =10−2 -10−1 , d = 10−2 -10−1 Mpc, SN R ∼ 10-260, n = 5, size of training set= 1800, andsize of testing set= 200. Each tree was created in a different 10 dimensional subspace of thefeature space, whose resolution was 585 dimensions. Originally, the resolution of the image
was much higher, O(10,000) dimensions; however, the random forest has a hard time working
with such large data streams, since the data are all given to the root node of the tree at thesame time to be split. In Figure 13.9, we show the results of our analysis. It is clear that as thenumber of trees in the forest increases,the detection efficiency improves. But note that each144

tree only analyzed such a small (random) portion of the data, and still a detection efficiency of∼ 90% was achieved. This is primarily due to the fact that a portion of the testing sample wascomposed of signals that were too weak to be seen by any algorithm.

Figure 13-8: Detection efficiency of the forest as a function of the number of trees present.
The leveling off at 90% occurs because there were some signals analyzed by the random foreststhat were too weak to be seen by any algorithm.

145

CHAPTER 14
MACHINE LEARNING-BASED SEARCH FOR A REMNANT OF GW170817
14.1

Introduction

To test our procedure, we ran a search on one week of O2 Livingston and Hanford data
after GW170817, spliced in 2000 s time/frequency maps, following the search design outlinedin the previous section, to look for a remnant of this binary neutron star merger. Though we
did not expect to find a signal, we wanted to ensure that we are prepared for future binary
neutron star mergers.
Performing this search took about 30 minutes. The CNNs took about 2 minutes toanalyze one week of data from [100,1900] Hz and to generate triggers. Then, the follow-up of
those 50 triggers by the GFH took about 15 minutes, when running this follow-up on 14 cores.
Finally, the follow-up for the candidate remaining after the GFH took less than a minute.
In Section 14.2, we describe the design of our machine learning/ GFH search for aremnant of GW170817. Afterwards in Section 14.3, we display the portion of the parameter
space that our search covered. We discuss the results of the search in Section 14.4, and giventhat no candidate signals were found, we place upper limits on the possible GW emission froma remnant for both fixed and time-varying braking indices in Section 14.5. These upper limits
are consistent with those found in (30 ).
14.2

Search Design

We plan to use the CNNs in a real search in order not only to speed up the analysis,
but also to probe a portion of the parameter space to which the GFH is not sensitive. It
seems that the CNNs are of comparable efficiency to the GFH, so they can be used to quickly
analyze time/frequency maps to determine if a signal is present or not. CNNs can analyze atime/frequency map in a microsecond, while the FH and GFH could take ∼ 1 minute (31 , 78 ).
If the map contains a signal, we can then apply the GFH to estimate the parameters of thesignal. Moreover, the CNNs are robust towards varying braking indices and different frequency

146

bands with a limited amount of training data. This means that the CNNs can be trained on asmall frequency band/time range and then tested across many different frequencies/times.
Once a CNN is trained on ∼ 20000 injections (at different amplitudes) across 150 Hz and2000 s maps, for about 22 days of data for a false alarm probability of 1% (which is about
2000 noise maps), we can apply the following procedure for each detector:
1. Load short fast Fourier transform Databases (77 ).
2. Choose appropriate FFT length depending on maximum f˙0 , create complete time/frequency
map and peakmap (78 ).
3. Clean known noise lines and persistent lines.
4. Reduce resolution by a factor of 16 on each axis.
5. Run CNNs on each reduced time/frequency map.
6. Repeat for each detector’s data. Do coincidences between time/frequency maps that
produce triggers.
7. If a trigger remains after the coincidence step, run GFH on the peakmap and obtainparameters of the signal.
8. Perform coincidences between candidates returned by the GFH, follow-up, etc. using
methodology described in (31 ).
9. Repeat for each 150 Hz/2000 s map.
We must emphasize that if a coincidence occurs between time/frequency maps at the levelof the CNNs, the GFH can estimate the parameters of a signal only if it follows a power-law
behavior. Otherwise, we can only constrain that the variation in braking index is no more thana certain δn.
As stated in Section 13.3, we cannot test a CNN on the same noise maps on which wetrain, since the CNN has learned that those maps do not contain a signal. But, if we create
different CNNs that are trained and tested on different portions of the data, we can avoid thisproblem. The logic is the following: suppose we have data for 90 days, and that we need totrain a CNN on only 2000 noise maps (23 days of data) to guarantee reasonable false alarm147

probability/detection efficiency comparable to the GFH. At a minimum, we can train on thefirst 23 days, and then run an analysis on the last 67 days. Afterwards, we could train ondays 24-47, then analyze days 1-23 and 48-90. This would then require the same number of
injections, just a doubling of the training time, which takes ∼ 30 minutes. The more we repeat
this “shifting”, the more times that a different CNN will see the same data.
We have shown that the CNNs are robust towards signals in different frequency bands in
Figure 13.5; hence we have not used injections in all bands. Instead, we have just trained on
signal and noise maps in one band.
14.3

Parameter Space Explored

We display the parameter space to which we are most sensitive in our search, assuming
a power-law signal. Using the criteria that a signal must not spin out of the frequency band
we are analyzing with the CNNs and the GFH, we determine which signal parameters we havethe best change of observing. Figure 14-1 plots the initial frequencies, spindowns and spindown
timescales to which we are most sensitive for our search configuration, fixing n = 7. The
spindown timescale is colored.
14.4

Description of Search Results

Coincidences between the outputs of Livingston and Hanford resulted in 50 time/frequency
maps having a value of p > pthr = 0.9. A grid over n was constructed for each of these maps
from n = [2.5, 7], in each 150 Hz band between f = [100, 1900] Hz. We then applied the GFH
individually to the peakmaps generated at these times (not the maps that the network saw, but
the maps constructed by the procedure in (78 )). The GFH returned 431 coincident candidates,
where a coincidence occurs if the Euclidean distance between candidate parameters x0 /k ineach detector is less than 3 bins. Most candidates had a very low detection statistic, called thecritical ratio (CR), defined as:CR =y − ȳ
σ148(14–1)

Figure 14-1: Parameter space explored in our search for a remnant of GW170817, for n = 7,
invoking the criteria that a signal must be fully contained within a frequency band analyzed to
be detectable.
where y is the number count returned by the GFH, ȳ and σ are the mean and standard
deviation of the full Hough map, respectively.
We required CR > 3.5, corresponding to a false alarm probability of 0.02%. After
applying this threshold, 1 candidate remained (for TF F T = 4 s). We then use the full BandSampled-Data (BSD) (81 ) follow-up procedure described in (31 ). The procedure corrects
the peakmaps for the time/frequency evolution of the signal, resulting (theoretically) in a
monochromatic signal. TF F T is then increased to improve the signal-to-noise ratio and theoriginal Frequency Hough is applied. If the critical ratio of a candidate is less now than before
the follow-up, the candidate is vetoed. In each iteration, TF F T is doubled. The candidate

149

remained for a few iterations until we used a TF F T = 64 s, when it was finally vetoed. We
therefore conclude that there was no detectable gravitational wave emission from a remnant of
GW170817.
14.5

Upper Limits

We determine the weakest signals we can see by performing 250 injections, each lasting
2000 s, at different amplitudes with parameters uniformly distributed between n = [2.5, 7],
f0 = [250, 1900] Hz, f˙ = [1/64, 1/16] Hz/s, cos ι = [−1, 1], polarization angle ψ = [−90, 90]
degrees, throughout the entire week after the merger. Another set of injections was alsodone with δn/δt = [−10−4 , 10−4 ] /s, analyzed only by the CNNs. A detection is obtained if
a candidate returned by the GFH is within 3 bins of the injections and if its CR > 3.5, therequirement used in the follow-up in the search. Figure 14-2 reports upper limits on distanceon a possible power-law signal coming from a remnant of GW170817 at the 50% confidence
level. We have plotted five curves: one using the CNNs followed up by the GFH to estimate
the parameters of the signal (magenta), one curve using just the CNNs to detect time-varying
braking index signals (blue), and three other curves for comparison purposes. The red one isfrom the original search (30 ), adjusted to be at the 50% confidence level, and the other two
curves were obtained by just using the CNNs (green) and GFH (black) to calculate the upper
limits. It seems that the previous upper limits are better because the signals were simulated for2000-16000 s consistent with the search length at n = 5 (30 ). This should imply an increaseof at most 40% in distance reach at higher frequencies, and n = 5 corresponds to the best
sensitivity by about 20% relative to other braking indices at lower frequencies for the GFH,
determined by using Equation 9–22 and looking at Figure 9-10. Additionally the grid in kconstructed during the GFH is not optimal for signals of such a short duration, especially athigh frequencies, so we are currently investigating ways to improve this grid.
For each curve, we determined the distance reach by assuming a moment of inertia of3
Izz = 100Msun
G2 /c4 ≈ 4.34 × 1038 kg· m2 , and calculating the maximum allowed ellipticity  byenergy conservation, as was done in (30 ).

150

Figure 14-2: Upper limits on distance reach at 50% confidence for the search design presented
in Section 14.2. CNNs followed up by the GFH (magenta curve) and only CNNs trying todetect time-varying braking index signals (blue). For comparison curves obtained by using justthe CNNs (green), the GFH in this search (black) and the previous post-merger search (red)(30 ) are shown. The CNNs are sensitive to time-varying braking indices as well, even thoughthese upper limits are done on fixed braking indices. Reasons why the GFH performed betterin the first post-merger search include the use of a fixed n = 5, boosting sensitivity at low
frequencies, and longer duration signals at high frequencies, improving sensitivity there (31 ).
Additionally the GFH in this search performs worse at higher frequency because the grid in k
is much finer and for signal durations of 2000 s, the signal is split among neighboring k bins.
Alternative ways to optimally construct the grid in k are currently being investigated.

151

CHAPTER 15
ALL-SKY SEARCH FOR LONG GRAVITATIONAL WAVE TRANSIENTS IN O3 DATA
15.1

Motivation

As stated in Section 6.2, many GRBs are observed above and below threshold each day
by Fermi-GBM and Swift-BAT (57 , 58 ). We do not know how far away these sources of
GRBs are, nor do we know what caused them. Due to the fact that GW170817 generated analmost sub-threshold GRB event, it is possible that mergers of binary systems are producing
these GRBs. If this is the case, long transient signals similar to those we searched for after
GW170817 may be occurring frequently.
Because we expect that these signals last for less than O(days) and have large spindowns,
we do not need to be concerned about the Doppler modulation. This means that our searches
are by default all-sky searches, or at least require so few grid points in the sky that any
signal not “properly corrected for” would simply be shifted uniformly in time/frequency
space, retaining the characteristics of the power-law signal for a particular braking index.
Moreover, our frequency bins are quite large compared to CW searches, so any kind of Doppler
modulation is confined to one frequency bin.
To illustrate that our search is essentially independent of position in the sky , we begin
with the Doppler shift equation:
~vorb · n̂f (t) = f0 1 +c(15–1)

where vorb is the velocity of the earth around the sun and n̂ is a unit vector pointing from theearth to the source of GWs. The maximum Doppler shift occurs when the earth is moving
directly towards or away from the source; hence:∆fdopp,max = f0

vorb
≈ 10−4 f0c

We impose that this shift must be less than a frequency bin dfmax :152(15–2)

∆fdopp,max <= dfmax =1
TF F T,max(15–3)

Hence, for a signal of f0 = 2000 Hz, TF F T,max <= 5 s. As long as we restrict our search tofrequencies at 2000 Hz or below and use TF F T <= TF F T,max , we can perform an all-sky searchwithout having to worry about Doppler modulations shifting potential GW signals by more than
one frequency bin.
15.2

Search Setup and Parameter Space ExploredThough the CNNs can perform an inherently unmodeled search, they are clearly likely topick up signals similar to those on which they were trained. So, even though we only train on asubset of the parameter space we are searching for, the CNNs can generalize to find new ones,
as we showed in 13.5. Therefore we are sensitive to all power-law and power-law-like signals,
regardless of f0 , f˙0 and n. However, our sensitivity is actually limited by how much time wecan observe and which TF F T we use.
We search 150 Hz x 4000 s time/frequency maps created using TF F T = 4 s for signals
lasting at most 4000 s, whose initial frequencies could be between 100 Hz and 2000 Hz and˙ = 1/16 Hz/s.
whose spindowns could be no greater than f0,max
One CNN was trained on the first calibration of data, C00, which had not been cleaned,for each detector (Livingston and Hanford). For Livingston, we used 55756 signal maps andnoise maps from 09-Apr-2019 00:38:17 to 07-May-2019 09:35:19.For Hanford, we used 54172
signal maps and noise maps from 08-Apr-2019 20:56:39 to 06-May-2019 07:11:05. The
networks were trained only on science data from each detector.
The analysis of coincident science time data began on data taken after 07-May-2019
09:35:19 and ended on July 1, which spans almost 2 months of data.
15.3

Search Progress and Preliminary Results

Time/frequency maps covering approximately 2 months of Livingston/Hanford data wereanalyzed by the CNNs. Thirty-six of these maps were in coincidence. On these 36 maps, we
perform the GFH follow-up as described in 14.2, covering a range of braking indices from 2.5

153to 7. This results in 180 Hough transforms to do, which is in sharp contrast to the O(105 )Hough transforms needed in the three configurations covering only a day of data in which
the first GW170817 post-merger search was run (see Chapter 10 for more details). After
performing coincidences among the candidates returned by the GFH, imposing that the average
critical ratio of candidates must be greater than 5, removing candidates near known noiselines, and vetoing candidates whose critical ratio in Hanford and Livingston differs by greaterthan a factor of 4 (all these procedures were employed in (30 )), we are left with 53 coincident
candidates to follow up. Considering a small window around the parameters of each coincident
candidate, and all combinations of these parameters, we will follow up 837 candidates. The
parameters of these candidates are shown in Figure 15-1. These potential signals occur at
the following times: 19-May-2019 05:38:31, 20-May-2019 00:09:01, 03-Jun-2019 10:54:09,
06-Jun-2019 05:51:13, 07-Jun-2019 08:36:57, 16-Jun-2019 18:24:57, and 23-Jun-2019
07:03:17.
A coincidence distance equal to 2 bins will be used now, and a greater critical ratio
threshold will be employed, since the signals should be stronger and be closer in the parameter
space in both detectors at this stage. Unfortunately the data structures needed to perform thefull follow-up procedure are not ready at the time of writing; in the meantime, preparing forthe likely scenario of no detection, we place stringent upper limits on the expected all-sky GWemission.
15.4

Upper Limits

Given that we do not expect candidates to survive our follow-up procedure, we place
upper limits on our distance reach, at 50% confidence. We perform 100 injections across manydifferent amplitudes in each frequency band analyzed, whose parameters were uniformly chosen
in reference time t0 , n = [3, 7], f0 = [fmax − bw/4, fmax ], where bw = 150 Hz is the bandwidth
in the time/frequency map, f˙0 = [1/(2TF F T )2 , 1/TF2 F T ]. The injections were first passed to theCNNs, and then the coincident maps were given to the GFH for further analysis, as done in thereal search. The upper limits, seen in Figure 15-2, were calculated in the same way as done

154

Figure 15-1: Candidates returned by the CNNs and GFH in the O3 all-sky search analysis.
Around 250 Hz, the candidates are actually spread from 240 to 250 Hz. We believe that many
candidates survived the coincident step around this frequency because the LIGO interferometers are most sensitive to potential signals the [100,250] Hz band.
in the O2 machine learning post-merger remnant search (see Section 14.5 for more details).
These are more stringent upper limits than obtained after the O2 post-merger searches.
15.5

Prospects

This search was done for only one choice of TF F T , but this is certainly not optimal if theGW signal has a smaller spindown than 1/TF2 F T or lasts for longer than the time we analyzed.
We therefore plan to employ different configurations, each containing time/frequency maps
with different TF F T = 2, 4, 8, 16 s and durations to cover this parameter space with better
sensitivity.

155

Figure 15-2: Upper limits calculated for our all-sky machine learning-based search for long transient GW signals. Distance reach is at 50% confidence. We injected 100 signals per amplitude,across a range of braking indices, initial frequencies, spindowns, and GPS start times in linewith the search parameter space.

156

CHAPTER 16
CONCLUSIONS
We have described in detail the background of GW physics, beginning from the basics ofgeneral relativity, and neutron star astrophysics, demonstrating the problem of determining theequation of state. R-modes, and the GWs they emit from being excited in a neutron star, werealso discussed. Recent results from real GW detections by the LIGO/Virgo collaborations were
presented, and physical conclusions concerning the speed of the graviton and the equation ofstate were drawn. We also explained how the LIGO/Virgo detectors work at an introductory
level, and described various data analysis strategies used both in our procedures and within
the LIGO/Virgo collaborations. Moreover, implications of a detection of GWs from an isolated
neutron star, and extrapolations of our sensitivity in future detectors, were discussed. The
major outputs of this thesis are two new, extensively developed methods to search for GWs
from isolated neutrons lasting O(hours − days) and two searches of O2 data for a remnant of
GW170817, and one almost-finished all-sky search of O3 data for long transient signals. It isclear from our results that machine learning can be useful in GW physics but needs to be usedcarefully, and that the GFH is a robust parameter estimator for power-law signals. However,both methods can be improved and applied on future data when a supernova occurs or twoneutron stars merge closer than 40 Mpc.

157

Bibliography
1. B. P. Abbott et al., Physical Review Letters 119, 161101 (16 2017).
2. B. P. Abbott et al., Physical Review Letters 116, 061102 (2016).
3. M. P. Hobson, G. P. Efstathiou, A. N. Lasenby, General relativity: an Introduction for
Physicists (Cambridge University Press, 2006).
4. W. Li et al., Monthly Notices of the Royal Astronomical Society 412, 1441–1472(2011).
5. S. P. Reynolds et al., The Astrophysical Journal Letters 680, L41 (2008).
6. E. Müller, H.-T. Janka, Astronomy and Astrophysics 317, 140–163 (1997).
7. K. N. Yakunin et al., Physical Review D 92, 084040 (2015).
8. S. Klimenko, I Yakushin, A Mercer, G. Mitselmakher, Classical and Quantum Gravity 25,
114029 (2008).
9. P. J. Sutton et al., New Journal of Physics 12, 053034 (2010).
10. S. A. Colgate, M. Herant, W. Benz, Physics Reports 227, 157–174 (1993).
11. N. Christensen, Reports on Progress in Physics 82, 016903 (2018).
12. A. L. Watts, B. Krishnan, L. Bildsten, B. F. Schutz, Monthly Notices of the Royal
Astronomical Society 389, 839–868 (2008).
13. G. B. Cook, S. L. Shapiro, S. A. Teukolsky, The Astrophysical Journal 423, L117
(1994).
14. R. Wijnands, M. van der Klis, Nature 394, 344 (1998).
15. L. Bildsten, The Astrophysical Journal Letters 501, L89 (1998).
16. G. Ushomirsky, C. Cutler, L. Bildsten, Monthly Notices of the Royal Astronomical
Society 319, 902–932 (2000).
17. C Messenger et al., Physical Review D 92, 023006 (2015).
18. N. K. Glendenning, Compact stars: Nuclear Physics, Particle Physics and General
Relativity (Springer Science & Business Media, 2012).
19. M. Camenzind, Compact Objects in Astrophysics (Springer, New York, 2007).

158

20. R. S. Foster, D. Backer, The Astrophysical Journal 361, 300–308 (1990).
21. J. W. Hessels et al., Science 311, 1901–1904 (2006).
22. S. L. Shapiro, S. A. Teukolsky, Black holes, white dwarfs, and neutron stars: The physics
of Compact Objects (Wiley, New York, 1983).
23. P. Arras, presented at the Binary Radio Pulsars, vol. 328, p. 317.
24. G. B. Cook, S. L. Shapiro, S. A. Teukolsky, The Astrophysical Journal 424, 823–845(1994).
25. F. Gittins, N. Andersson, arXiv preprint arXiv:1811.00550 (2018).
26. J. D. Jackson, Classical electrodynamics, 1999.
27. B. P. Abbott et al., arXiv: 1707.02669 (gr-qc) (July 2017).
28. B Knispel et al., Science 329, 1305–1305 (2010).
29. P. Jaranowski, A. Krolak, B. F. Schutz, Physical Review D 58, 063001 (1998).
30. B. P. Abbott et al., The Astrophysical Journal 875, 160 (2019).
31. A. Miller et al., Phys. Rev. D 98, 102004 (10 2018).
32. M. Oliver, D. Keitel, A. M. Sintes, Physical Review D 99, 104067 (2019).
33. L. Sun, A. Melatos, Physical Review D 99, 123003 (2019).
34. N. Sarin, P. D. Lasky, L. Sammut, G. Ashton, arXiv preprint arXiv:1805.01481 (2018).
35. B. J. Owen et al., Physical Review D 58, 084020 (1998).
36. M. G. Alford, K. Schwenzer, Monthly Notices of the Royal Astronomical Society 446,
3631–3641 (2014).
37. N. Andersson, K. D. Kokkotas, International Journal of Modern Physics D 10, 381–441(2001).
38. B. J. Owen, Physical Review D 82, 104002 (2010).
39. K. S. Thorne, Reviews of Modern Physics 52, 299 (1980).
40. J Provost, G Berthomieu, A Rocca, Astronomy and Astrophysics 94, 126 (1981).
41. L. Lindblom, B. J. Owen, S. M. Morsink, Physical Review Letters 80, 4843 (1998).
42. T. G. Cowling, Monthly Notices of the Royal Astronomical Society 101, 367 (1941).

159

43. L. Lindblom, R. J. Splinter, The Astrophysical Journal 348, 198–202 (1990).
44. S. Yoshida, Y. Kojima, Monthly Notices of the Royal Astronomical Society 289, 117–
122 (1997).
45. S. Yoshida, U. Lee, The Astrophysical Journal 567, 1112 (2002).
46. J. R. Oppenheimer, H. Snyder, Physical Review 56, 455 (1939).
47. A. Mytidis, A. A. Panagopoulos, O. P. Panagopoulos, A. Miller, B. Whiting, PhysicalReview D 99, 024024 (2019).
48. B. Abbott et al., Physical Review X 9, 031040 (2019).
49. B. P. Abbott et al., Astrophys. J. Lett 848, L12 (2017).
50. P. Pani, C. F. Macedo, L. C. Crispino, V. Cardoso, Physical Review D 84, 087501(2011).
51. L. S. Finn, P. J. Sutton, Physical Review D 65, 044022 (2002).
52. A. G. A. von Kienlin, C. Meegan, the Fermi GBM, GCN 21520 (2017).
53. B. P. Abbott et al., The Astrophysical Journal Letters 848, L13 (2017).
54. G. Gonzales, Moving pendulum, http://www.phys.lsu.edu/faculty/gonzalez/
Teaching/Phys7221/PendulumWithMovingSupport.pdf.
55. M. Maggiore, Gravitational Waves: Volume 1: Theory and Experiments (Oxford
University Press, 2008), vol. 1.
56. P. Shawhan, LSC/Virgo, “Fermi-GBM and Swift GCN Notice Statistics”, LIGO
Document T1900297-v4, 2019.
57. C. Meegan et al., The Astrophysical Journal 702, 791 (2009).
58. H. A. Krimm et al., The Astrophysical Journal Supplement Series 209, 14 (2013).
59. B. P. Abbott et al., The Astrophysical Journal Letters 832, L21 (2016).
60. V. Tiwari, Classical and Quantum Gravity 35, 145009 (2018).
61. V. Tiwari, S. Fairhurst, M. Hannam, The Astrophysical Journal 868, 140 (2018).
62. J. Abadie et al., Classical and Quantum Gravity 27, 173001 (2010).
63. E. S. Phinney, The Astrophysical Journal 380, L17–L21 (1991).

160

64. R. Narayan, The Astrophysical Journal 319, 162–179 (1987).
65. C. Kim, V. Kalogera, D. Lorimer, The Astrophysical Journal 584, 985 (2003).
66. V. Kalogera, R Narayan, D. N. Spergel, J. Taylor, The Astrophysical Journal 556, 340(2001).
67. R. Bondarescu, S. A. Teukolsky, I. Wasserman, Physical Review D 79, 104003 (2009).
68. M Punturo et al., Classical and Quantum Gravity 27, 194002 (2010).
69. B. P. Abbott et al., Classical and Quantum Gravity 34, 044001 (2017).
70. P. J. Sutton, arXiv preprint arXiv:1304.0210 (2013).
71. H.-Y. Chen et al., arXiv preprint arXiv:1709.08079 (2017).
72. L. S. Finn, D. F. Chernoff, Physical Review D 47, 2198 (1993).
73. M. Evans et al., “Unofficial sensitivity curves (ASD) for aLIGO, Kagra, Virgo, Voyager,
Cosmic Explorer and ET”, https://dcc.ligo.org/LIGO-T1500293/public, 2019.
74. W. van Straten et al., Nature 412, 158 (2001).
75. I. I. Shapiro, Physical Review Letters 13, 789 (1964).
76. J. H. Taylor, J. M. Weisberg, The Astrophysical Journal 345, 434–450 (1989).
77. P Astone, S Frasca, C Palomba, Classical and Quantum Gravity 22, S1197 (2005).
78. P. Astone, A. Colla, S. D’Antonio, S. Frasca, C. Palomba, Physical Review D 90,
042002 (2014).
79. P. Astone, S. D’Antonio, S. Frasca, C. Palomba, Classical and Quantum Gravity 27,
194016 (2010).
80. S. Caride, R. Inta, B. J. Owen, B. Rajbhandari, Phys. Rev. D 100, 064013 (6 2019).
81. O. Piccinni et al., Classical and Quantum Gravity 36, 015008 (2018).
82. L. Baiotti, L. Rezzolla, Rept. Prog. Phys. 80, 096901 (2017).
83. Y.-W. Yu, L.-D. Liu, Z.-G. Dai, The Astrophysical Journal 861, 114 (2018).
84. S. A. Colgate, The Astrophysical Journal 163, 221 (1971).
85. A. Vajpeyi, A. Miller, P. Astone, S. Frasca, Bulletin of the American Physical Society 62(2017).

161

86. L. Pierini, Ph.D. Thesis, Sapienza (2019).
87. B. P. Abbott et al., Physical Review D 96, 062002 (2017).
88. B. Abbott et al., Physical Review D 100, 024004 (2019).
89. K. Riles, Modern Physics Letters A 32, 1730035 (2017).
90. E. Thrane et al., Physical Review D 83, 083004 (2011).
91. E. Thrane, V. Mandic, N. Christensen, Physical Review D 91, 104021 (2015).
92. R. Prix, S. Giampanis, C. Messenger, Physical Review D 84, 023007 (2011).
93. D. Keitel et al., arXiv preprint arXiv:1907.04717 (2019).
94. T. Hill, L. Marquez, M. O’Connor, W. Remus, International journal of forecasting 10,
5–15 (1994).
95. M. I. Jordan, T. M. Mitchell, Science 349, 255–260 (2015).
96. Y. LeCun, Y. Bengio, G. Hinton, Nature 521, 436 (2015).
97. H. Adeli, Kybernetes 28, 317–318 (1999).
98. R. Hecht-Nielsen, in Neural networks for perception (Elsevier, 1992), pp. 65–93.
99. A. Ng, in Machine learning.
100. A. Krizhevsky, I. Sutskever, G. E. Hinton, presented at the Advances in neural
information processing systems, pp. 1097–1105.
101. G. Hinton et al., IEEE Signal processing magazine 29 (2012).
102. H. Greenspan, B. Van Ginneken, R. M. Summers, IEEE Transactions on Medical
Imaging 35, 1153–1159 (2016).
103. P Astone et al., Physical Review D 98, 122002 (2018).
104. T. K. Ho, presented at the Proceedings of 3rd International Conference on Document
Analysis and Recognition, vol. 1, pp. 278–282.
105. T. K. Ho, IEEE Transactions on Pattern Analysis and Machine Intelligence 20, 832–844,
issn: 0162-8828 (Aug. 1998).
106. P. T. Baker et al., Physical Review D 91, 062004 (2015).
107. D. George, E. Huerta, Physical Review D 97, 044039 (2018).

162

108. D. George, E. Huerta, Physics Letters B 778, 64–70 (2018).
109. H. Gabbard, M. Williams, F. Hayes, C. Messenger, Physical Review Letters 120, 141103(2018).
110. C. Dreissigacker, R. Sharma, C. Messenger, R. Zhao, R. Prix, arXiv preprint
arXiv:1904.13291 (2019).
111. T. D. Gebhard, N. Kilbertus, I. Harry, B. Schölkopf, arXiv preprint arXiv:1904.08693(2019).
112. A. L. Miller et al., Physical Review D 100, 062005 (2019).
113. T. S. Adams et al., Physical Review D 88, 062006 (2013).

163

BIOGRAPHICAL SKETCH
Andrew Miller earned his Bachelor of Science in physics at The College of New Jersey
in 2015. He decided to immediately pursue his Ph.D. in physics at the University of Florida
and University of Rome, La Sapienza. His interests are in detecting gravitational waves fromisolated neutron stars, machine learning algorithms, pattern recognition techniques, future
generation gravitational wave detectors, and axion-like dark matter particle searches. He earned
his Ph.D. in physics from the University of Florida and University of Rome, La Sapienza jointly
in December 2019.

164