of 39
LIGO-P060010-06-Z
All-sky search for periodic gravitational waves in LIGO S4 data
B. Abbott,
14
R. Abbott,
14
R. Adhikari,
14
J. Agresti,
14
P. Ajith,
2
B. Allen,
2, 51
R. Amin,
18
S. B. Anderson,
14
W. G. Anderson,
51
M. Arain,
39
M. Araya,
14
H. Armandula,
14
M. Ashley,
4
S. Aston,
38
P. Aufmuth,
36
C. Aulbert,
1
S. Babak,
1
S. Ballmer,
14
H. Bantilan,
8
B. C. Barish,
14
C. Barker,
15
D. Barker,
15
B. Barr,
40
P. Barriga,
50
M. A. Barton,
40
K. Bayer,
17
K. Belczynski,
24
J. Betzwieser,
17
P. T. Beyersdorf,
27
B. Bhawal,
14
I. A. Bilenko,
21
G. Billingsley,
14
R. Biswas,
51
E. Black,
14
K. Blackburn,
14
L. Blackburn,
17
D. Blair,
50
B. Bland,
15
J. Bogenstahl,
40
L. Bogue,
16
R. Bork,
14
V. Boschi,
14
S. Bose,
52
P. R. Brady,
51
V. B. Braginsky,
21
J. E. Brau,
43
M. Brinkmann,
2
A. Brooks,
37
D. A. Brown,
14, 6
A. Bullington,
30
A. Bunkowski,
2
A. Buonanno,
41
O. Burmeister,
2
D. Busby,
14
R. L. Byer,
30
L. Cadonati,
17
G. Cagnoli,
40
J. B. Camp,
22
J. Cannizzo,
22
K. Cannon,
51
C. A. Cantley,
40
J. Cao,
17
L. Cardenas,
14
M. M. Casey,
40
G. Castaldi,
46
C. Cepeda,
14
E. Chalkey,
40
P. Charlton,
9
S. Chatterji,
14
S. Chelkowski,
2
Y. Chen,
1
F. Chiadini,
45
D. Chin,
42
E. Chin,
50
J. Chow,
4
N. Christensen,
8
J. Clark,
40
P. Cochrane,
2
T. Cokelaer,
7
C. N. Colacino,
38
R. Coldwell,
39
R. Conte,
45
D. Cook,
15
T. Corbitt,
17
D. Coward,
50
D. Coyne,
14
J. D. E. Creighton,
51
T. D. Creighton,
14
R. P. Croce,
46
D. R. M. Crooks,
40
A. M. Cruise,
38
A. Cumming,
40
J. Dalrymple,
31
E. D’Ambrosio,
14
K. Danzmann,
36, 2
G. Davies,
7
D. DeBra,
30
J. Degallaix,
50
M. Degree,
30
T. Demma,
46
V. Dergachev,
42
S. Desai,
32
R. DeSalvo,
14
S. Dhurandhar,
13
M. D ́ıaz,
33
J. Dickson,
4
A. Di Credico,
31
G. Diederichs,
36
A. Dietz,
7
E. E. Doomes,
29
R. W. P. Drever,
5
J.-C. Dumas,
50
R. J. Dupuis,
14
J. G. Dwyer,
10
P. Ehrens,
14
E. Espinoza,
14
T. Etzel,
14
M. Evans,
14
T. Evans,
16
S. Fairhurst,
7, 14
Y. Fan,
50
D. Fazi,
14
M. M. Fejer,
30
L. S. Finn,
32
V. Fiumara,
45
N. Fotopoulos,
51
A. Franzen,
36
K. Y. Franzen,
39
A. Freise,
38
R. Frey,
43
T. Fricke,
44
P. Fritschel,
17
V. V. Frolov,
16
M. Fyffe,
16
V. Galdi,
46
J. Garofoli,
15
I. Gholami,
1
J. A. Giaime,
16, 18
S. Giampanis,
44
K. D. Giardina,
16
K. Goda,
17
E. Goetz,
42
L. M. Goggin,
14
G. Gonz ́alez,
18
S. Gossler,
4
A. Grant,
40
S. Gras,
50
C. Gray,
15
M. Gray,
4
J. Greenhalgh,
26
A. M. Gretarsson,
11
R. Grosso,
33
H. Grote,
2
S. Grunewald,
1
M. Guenther,
15
R. Gustafson,
42
B. Hage,
36
D. Hammer,
51
C. Hanna,
18
J. Hanson,
16
J. Harms,
2
G. Harry,
17
E. Harstad,
43
T. Hayler,
26
J. Heefner,
14
I. S. Heng,
40
A. Heptonstall,
40
M. Heurs,
2
M. Hewitson,
2
S. Hild,
36
E. Hirose,
31
D. Hoak,
16
D. Hosken,
37
J. Hough,
40
E. Howell,
50
D. Hoyland,
38
S. H. Huttner,
40
D. Ingram,
15
E. Innerhofer,
17
M. Ito,
43
Y. Itoh,
51
A. Ivanov,
14
D. Jackrel,
30
B. Johnson,
15
W. W. Johnson,
18
D. I. Jones,
47
G. Jones,
7
R. Jones,
40
L. Ju,
50
P. Kalmus,
10
V. Kalogera,
24
D. Kasprzyk,
38
E. Katsavounidis,
17
K. Kawabe,
15
S. Kawamura,
23
F. Kawazoe,
23
W. Kells,
14
D. G. Keppel,
14
F. Ya. Khalili,
21
C. Kim,
24
P. King,
14
J. S. Kissel,
18
S. Klimenko,
39
K. Kokeyama,
23
V. Kondrashov,
14
R. K. Kopparapu,
18
D. Kozak,
14
B. Krishnan,
1
P. Kwee,
36
P. K. Lam,
4
M. Landry,
15
B. Lantz,
30
A. Lazzarini,
14
B. Lee,
50
M. Lei,
14
J. Leiner,
52
V. Leonhardt,
23
I. Leonor,
43
K. Libbrecht,
14
P. Lindquist,
14
N. A. Lockerbie,
48
M. Longo,
45
M. Lormand,
16
M. Lubinski,
15
H. L ̈uck,
36, 2
B. Machenschalk,
1
M. MacInnis,
17
M. Mageswaran,
14
K. Mailand,
14
M. Malec,
36
V. Mandic,
14
S. Marano,
45
S. M ́arka,
10
J. Markowitz,
17
E. Maros,
14
I. Martin,
40
J. N. Marx,
14
K. Mason,
17
L. Matone,
10
V. Matta,
45
N. Mavalvala,
17
R. McCarthy,
15
D. E. McClelland,
4
S. C. McGuire,
29
M. McHugh,
20
K. McKenzie,
4
J. W. C. McNabb,
32
S. McWilliams,
22
T. Meier,
36
A. Melissinos,
44
G. Mendell,
15
R. A. Mercer,
39
S. Meshkov,
14
E. Messaritaki,
14
C. J. Messenger,
40
D. Meyers,
14
E. Mikhailov,
17
S. Mitra,
13
V. P. Mitrofanov,
21
G. Mitselmakher,
39
R. Mittleman,
17
O. Miyakawa,
14
S. Mohanty,
33
G. Moreno,
15
K. Mossavi,
2
C. MowLowry,
4
A. Moylan,
4
D. Mudge,
37
G. Mueller,
39
S. Mukherjee,
33
H. M ̈uller-Ebhardt,
2
J. Munch,
37
P. Murray,
40
E. Myers,
15
J. Myers,
15
T. Nash,
14
G. Newton,
40
A. Nishizawa,
23
K. Numata,
22
B. O’Reilly,
16
R. O’Shaughnessy,
24
D. J. Ottaway,
17
H. Overmier,
16
B. J. Owen,
32
Y. Pan,
41
M. A. Papa,
1, 51
V. Parameshwaraiah,
15
P. Patel,
14
M. Pedraza,
14
S. Penn,
12
V. Pierro,
46
I. M. Pinto,
46
M. Pitkin,
40
H. Pletsch,
2
M. V. Plissi,
40
F. Postiglione,
45
R. Prix,
1
V. Quetschke,
39
F. Raab,
15
D. Rabeling,
4
H. Radkins,
15
R. Rahkola,
43
N. Rainer,
2
M. Rakhmanov,
32
M. Ramsunder,
32
K. Rawlins,
17
S. Ray-Majumder,
51
V. Re,
38
H. Rehbein,
2
S. Reid,
40
D. H. Reitze,
39
L. Ribichini,
2
R. Riesen,
16
K. Riles,
42
B. Rivera,
15
N. A. Robertson,
14, 40
C. Robinson,
7
E. L. Robinson,
38
S. Roddy,
16
A. Rodriguez,
18
A. M. Rogan,
52
J. Rollins,
10
J. D. Romano,
7
J. Romie,
16
R. Route,
30
S. Rowan,
40
A. R ̈udiger,
2
L. Ruet,
17
P. Russell,
14
K. Ryan,
15
S. Sakata,
23
M. Samidi,
14
L. Sancho de la Jordana,
35
V. Sandberg,
15
V. Sannibale,
14
S. Saraf,
25
P. Sarin,
17
B. S. Sathyaprakash,
7
S. Sato,
23
P. R. Saulson,
31
R. Savage,
15
P. Savov,
6
S. Schediwy,
50
R. Schilling,
2
R. Schnabel,
2
R. Schofield,
43
B. F. Schutz,
1, 7
P. Schwinberg,
15
S. M. Scott,
4
A. C. Searle,
4
B. Sears,
14
F. Seifert,
2
D. Sellers,
16
A. S. Sengupta,
7
P. Shawhan,
41
D. H. Shoemaker,
17
A. Sibley,
16
J. A. Sidles,
49
X. Siemens,
14, 6
D. Sigg,
15
S. Sinha,
30
A. M. Sintes,
35, 1
B. J. J. Slagmolen,
4
J. Slutsky,
18
J. R. Smith,
2
M. R. Smith,
14
K. Somiya,
2, 1
K. A. Strain,
40
D. M. Strom,
43
A. Stuver,
32
T. Z. Summerscales,
3
K.-X. Sun,
30
M. Sung,
18
P. J. Sutton,
14
H. Takahashi,
1
D. B. Tanner,
39
M. Tarallo,
14
R. Taylor,
14
R. Taylor,
40
J. Thacker,
16
K. A. Thorne,
32
K. S. Thorne,
6
A. Th ̈uring,
36
K. V. Tokmakov,
40
C. Torres,
33
C. Torrie,
40
G. Traylor,
16
M. Trias,
35
W. Tyler,
14
D. Ugolini,
34
C. Ungarelli,
38
K. Urbanek,
30
H. Vahlbruch,
36
arXiv:0708.3818v2 [gr-qc] 1 Dec 2009
M. Vallisneri,
6
C. Van Den Broeck,
7
M. Varvella,
14
S. Vass,
14
A. Vecchio,
38
J. Veitch,
40
P. Veitch,
37
A. Villar,
14
C. Vorvick,
15
S. P. Vyachanin,
21
S. J. Waldman,
14
L. Wallace,
14
H. Ward,
40
R. Ward,
14
K. Watts,
16
D. Webber,
14
A. Weidner,
2
M. Weinert,
2
A. Weinstein,
14
R. Weiss,
17
S. Wen,
18
K. Wette,
4
J. T. Whelan,
1
D. M. Whitbeck,
32
S. E. Whitcomb,
14
B. F. Whiting,
39
C. Wilkinson,
15
P. A. Willems,
14
L. Williams,
39
B. Willke,
36, 2
I. Wilmut,
26
W. Winkler,
2
C. C. Wipf,
17
S. Wise,
39
A. G. Wiseman,
51
G. Woan,
40
D. Woods,
51
R. Wooley,
16
J. Worden,
15
W. Wu,
39
I. Yakushin,
16
H. Yamamoto,
14
Z. Yan,
50
S. Yoshida,
28
N. Yunes,
32
M. Zanolin,
17
J. Zhang,
42
L. Zhang,
14
C. Zhao,
50
N. Zotov,
19
M. Zucker,
17
H. zur M ̈uhlen,
36
and J. Zweizig
14
(The LIGO Scientific Collaboration, http://www.ligo.org)
1
Albert-Einstein-Institut, Max-Planck-Institut f ̈ur Gravitationsphysik, D-14476 Golm, Germany
2
Albert-Einstein-Institut, Max-Planck-Institut f ̈ur Gravitationsphysik, D-30167 Hannover, Germany
3
Andrews University, Berrien Springs, MI 49104 USA
4
Australian National University, Canberra, 0200, Australia
5
California Institute of Technology, Pasadena, CA 91125, USA
6
Caltech-CaRT, Pasadena, CA 91125, USA
7
Cardiff University, Cardiff, CF24 3AA, United Kingdom
8
Carleton College, Northfield, MN 55057, USA
9
Charles Sturt University, Wagga Wagga, NSW 2678, Australia
10
Columbia University, New York, NY 10027, USA
11
Embry-Riddle Aeronautical University, Prescott, AZ 86301 USA
12
Hobart and William Smith Colleges, Geneva, NY 14456, USA
13
Inter-University Centre for Astronomy and Astrophysics, Pune - 411007, India
14
LIGO - California Institute of Technology, Pasadena, CA 91125, USA
15
LIGO Hanford Observatory, Richland, WA 99352, USA
16
LIGO Livingston Observatory, Livingston, LA 70754, USA
17
LIGO - Massachusetts Institute of Technology, Cambridge, MA 02139, USA
18
Louisiana State University, Baton Rouge, LA 70803, USA
19
Louisiana Tech University, Ruston, LA 71272, USA
20
Loyola University, New Orleans, LA 70118, USA
21
Moscow State University, Moscow, 119992, Russia
22
NASA/Goddard Space Flight Center, Greenbelt, MD 20771, USA
23
National Astronomical Observatory of Japan, Tokyo 181-8588, Japan
24
Northwestern University, Evanston, IL 60208, USA
25
Rochester Institute of Technology, Rochester, NY 14623, USA
26
Rutherford Appleton Laboratory, Chilton, Didcot, Oxon OX11 0QX United Kingdom
27
San Jose State University, San Jose, CA 95192, USA
28
Southeastern Louisiana University, Hammond, LA 70402, USA
29
Southern University and A&M College, Baton Rouge, LA 70813, USA
30
Stanford University, Stanford, CA 94305, USA
31
Syracuse University, Syracuse, NY 13244, USA
32
The Pennsylvania State University, University Park, PA 16802, USA
33
The University of Texas at Brownsville and Texas Southmost College, Brownsville, TX 78520, USA
34
Trinity University, San Antonio, TX 78212, USA
35
Universitat de les Illes Balears, E-07122 Palma de Mallorca, Spain
36
Universit ̈at Hannover, D-30167 Hannover, Germany
37
University of Adelaide, Adelaide, SA 5005, Australia
38
University of Birmingham, Birmingham, B15 2TT, United Kingdom
39
University of Florida, Gainesville, FL 32611, USA
40
University of Glasgow, Glasgow, G12 8QQ, United Kingdom
41
University of Maryland, College Park, MD 20742 USA
42
University of Michigan, Ann Arbor, MI 48109, USA
43
University of Oregon, Eugene, OR 97403, USA
44
University of Rochester, Rochester, NY 14627, USA
45
University of Salerno, 84084 Fisciano (Salerno), Italy
46
University of Sannio at Benevento, I-82100 Benevento, Italy
47
University of Southampton, Southampton, SO17 1BJ, United Kingdom
48
University of Strathclyde, Glasgow, G1 1XQ, United Kingdom
49
University of Washington, Seattle, WA, 98195
50
University of Western Australia, Crawley, WA 6009, Australia
51
University of Wisconsin-Milwaukee, Milwaukee, WI 53201, USA
52
Washington State University, Pullman, WA 99164, USA
(Dated: February 13, 2013)
2
We report on an all-sky search with the LIGO detectors for periodic gravitational waves in the
frequency range 50 – 1000 Hz and with the frequency’s time derivative in the range
1
×
10
8
Hz s
1
to zero. Data from the fourth LIGO science run (S4) have been used in this search. Three different
semi-coherent methods of transforming and summing strain power from Short Fourier Transforms
(SFTs) of the calibrated data have been used. The first, known as “StackSlide”, averages normalized
power from each SFT. A “weighted Hough” scheme is also developed and used, and which also
allows for a multi-interferometer search. The third method, known as “PowerFlux”, is a variant
of the StackSlide method in which the power is weighted before summing. In both the weighted
Hough and PowerFlux methods, the weights are chosen according to the noise and detector antenna-
pattern to maximize the signal-to-noise ratio. The respective advantages and disadvantages of these
methods are discussed. Observing no evidence of periodic gravitational radiation, we report upper
limits; we interpret these as limits on this radiation from isolated rotating neutron stars. The best
population-based upper limit with 95% confidence on the gravitational-wave strain amplitude, found
for simulated sources distributed isotropically across the sky and with isotropically distributed spin-
axes, is 4
.
28
×
10
24
(near 140 Hz). Strict upper limits are also obtained for small patches on the
sky for best-case and worst-case inclinations of the spin axes.
PACS numbers: 04.80.Nn, 95.55.Ym, 97.60.Gb, 07.05.Kf
I. INTRODUCTION
We report on a search with the LIGO (Laser Interfer-
ometer Gravitational-wave Observatory) detectors [1, 2]
for periodic gravitational waves in the frequency range
50 – 1000 Hz and with the frequency’s time derivative in
the range
1
×
10
8
Hz s
1
to zero. The search is carried
out over the entire sky using data from the fourth LIGO
science run (S4). Isolated rotating neutron stars in our
galaxy are the prime target.
Using data from earlier science runs, the LIGO Sci-
entific Collaboration (LSC) has previously reported on
searches for periodic gravitational radiation, using a long-
period coherent method to target known pulsars [3, 4, 5],
using a short-period coherent method to target Scorpius
X-1 in selected bands and search the entire sky in the
160
.
0 – 728
.
8 Hz band [6], and using a long-period semi-
coherent method to search the entire sky in the 200 –
400 Hz band [7]. Einstein@Home, a distributed home
computing effort running under the BOINC architecture
[8], has also been searching the entire sky using a coher-
ent first stage, followed by a simple coincidence stage [9].
In comparison, this paper: 1) examines more sensitive
data; 2) searches over a larger range in frequency and its
derivative; and 3) uses three alternative semi-coherent
methods for summing measured strain powers to detect
excess power from a continuous gravitational-wave signal.
The first purpose of this paper is to present results from
our search for periodic gravitational waves in the S4 data.
Over the LIGO frequency band of sensitivity, the S4 all-
sky upper limits presented here are approximately an or-
der of magnitude better than published previously from
earlier science runs [6, 7]. After following up on outliers
in the data, we find that no candidates survive, and thus
report upper limits. These are interpreted as limits on
radiation from rotating neutron stars, which can be ex-
pressed as functions of the star’s ellipticity and distance,
allowing for an astrophysical interpretation. The best
population-based upper limit with 95% confidence on the
gravitational-wave strain amplitude, found for simulated
sources distributed isotropically across the sky and with
isotropically distributed spin-axes, is 4
.
28
×
10
24
(near
140 Hz). Strict upper limits are also obtained for small
patches on the sky for best-case and worst-case inclina-
tions of the spin axes.
The second purpose of this paper, along with the pre-
vious coherent [6] and semi-coherent [7] papers, is to lay
the foundation for the methods that will be used in fu-
ture searches. It is well known that the search for periodic
gravitational waves is computationally bound; to obtain
optimal results will require a hierarchical approach that
uses coherent and semi-coherent stages [10, 11, 12, 13]. A
fifth science run (S5), which started in November 2005, is
generating data at initial LIGO’s design sensitivity. We
plan to search this data using the best methods possible,
based on what is learned from this and previous analyses.
In the three methods considered here, one searches
for cumulative excess power from a hypothetical periodic
gravitational wave signal by examining successive spec-
tral estimates based on Short Fourier Transforms (SFTs)
of the calibrated detector strain data channel, taking into
account the Doppler modulations of detected frequency
due to the Earth’s rotational and orbital motion with
respect to the Solar System Barycenter (SSB), and the
time derivative of the frequency intrinsic to the source.
The simplest method presented, known as “StackSlide”
[12, 13, 14, 15], averages normalized power from each
SFT. In the Hough method reported previously [7, 10],
referred to here as “standard Hough”, the sum is of bi-
nary zeroes or ones, where an SFT contributes unity if
the power exceeds a normalized power threshold. In this
paper a “weighted Hough” scheme, henceforth also re-
ferred to as “Hough”, has been developed and is simi-
lar to that described in Ref. [16]. This scheme also al-
lows for a multi-interferometer search. The third method,
known as “PowerFlux” [17], is a variant of the StackSlide
method in which the power is weighted before summing.
In both the weighted Hough and PowerFlux methods,
the weights are chosen according to the noise and de-
tector antenna pattern to maximize the signal-to-noise
3
ratio.
The Hough method is computationally faster and more
robust against large transient power artifacts, but is
slightly less sensitive than StackSlide for stationary data
[7, 15]. The PowerFlux method is found in most fre-
quency ranges to have better detection efficiency than
the StackSlide and Hough methods, the exceptions oc-
curring in bands with large non-stationary artifacts, for
which the Hough method proves more robust. However,
the StackSlide and Hough methods can be made more
sensitive by starting with the maximum likelihood statis-
tic (known as the
F
-statistic [6, 10, 18]) rather than SFT
power as the input data, though this improvement comes
with increased computational cost. The trade-offs among
the methods means that each could play a role in our fu-
ture searches.
In brief, this paper makes several important contribu-
tions. It sets the best all-sky upper limits on periodic
gravitational waves to date, and shows that these limits
are becoming astrophysically interesting. It also intro-
duces methods that are crucial to the development of
our future searches.
This paper is organized as follows: Section II briefly
describes the LIGO interferometers, focusing on improve-
ments made for the S4 data run, and discusses the sen-
sitivity and relevant detector artifacts. Section III pre-
cisely defines the waveforms we seek and the associated
assumptions we have made. Section IV gives a detailed
description of the three analysis methods used and sum-
marizes their similarities and differences, while Section
V gives the details of their implementations and the
pipelines used. Section VI discusses the validation of the
software and, as an end-to-end test, shows the detection
of simulated pulsar signals injected into the data stream
at the hardware level. Section VII describes the search
results, and Section VIII compares the results from the
three respective methods. Section IX concludes with a
summary of the results, their astrophysical implications,
and future plans.
II. THE LIGO DETECTOR NETWORK AND
THE S4 SCIENCE RUN
The LIGO detector network consists of a 4-km inter-
ferometer in Livingston Louisiana (called L1) and two
interferometers in Hanford Washington, one 4-km and
another 2-km (H1 and H2, respectively).
The data analyzed in this paper were produced dur-
ing LIGO’s 29.5-day fourth science run (S4) [19]. This
run started at noon Central Standard Time (CST) on
February 22 and ended at midnight CST on March 23,
2005. During the run, all three LIGO detectors had dis-
placement spectral amplitudes near 2
.
5
×
10
19
m Hz
1
/
2
in their most sensitive frequency band near 150 Hz. In
units of gravitational-wave strain amplitude, the sensi-
tivity of H2 is roughly a factor of two worse than that
of H1 and L1 over much of the search band. The typical
strain sensitivities in this run were within a factor of two
of the design goals. Figure 1 shows representative strain
spectral noise densities for the three interferometers dur-
ing the run. As discussed in Section V below, however,
non-stationarity of the noise was significant.
Changes to the interferometers before the S4 run in-
cluded the following improvements [19]:
Installation of active seismic isolation of support
structures at Livingston to cope with high anthro-
pogenic ground motion in the 1-3 Hz band.
Thermal compensation with a CO
2
laser of mirrors
subject to thermal lensing from the primary laser
beam to a greater or lesser degree than expected.
Replacement of a synthesized radio frequency oscil-
lator for phase modulation with a crystal oscillator
before S4 began (H1) and mid-way through the S4
run (L1), reducing noise substantially above 1000
Hz and eliminating a comb of
37 Hz lines. (The
crystal oscillator replacement for H2 occurred after
the S4 run.)
Lower-noise mirror-actuation electronics (H1, H2,
& L1).
Higher-bandwidth laser frequency stabilization
(H1, H2, & L1) and intensity stabilization (H1 &
L1).
Installation of radiation pressure actuation of mir-
rors for calibration validation (H1).
Commissioning of complete alignment control sys-
tem for the L1 interferometer (already implemented
for H1 & H2 in S3 run).
Refurbishment of lasers and installation of photo-
diodes and electronics to permit interferometer op-
eration with increased laser power (H1, H2, & L1).
Mitigation of electromagnetic interference (H1, H2,
& L1) and acoustic interference (L1).
The data were acquired and digitized at a rate of 16384
Hz. Data acquisition was periodically interrupted by dis-
turbances such as seismic transients, reducing the net
running time of the interferometers. The resulting duty
factors for the interferometers were 81% for H1 and H2,
and 74% for L1. While the H1 and H2 duty factors were
somewhat higher than those in previous science runs, the
L1 duty factor was dramatically higher than the
'
40%
typical of the past, thanks to the increased stability from
the installation of the active seismic isolation system at
Livingston.
III. SIGNAL WAVEFORMS
The general form of a gravitational-wave signal is de-
scribed in terms of two orthogonal transverse polariza-
tions defined as “+” with waveform
h
+
(
t
) and “
×
” with
4
FIG. 1: Median amplitude strain noise spectral densities from
the three LIGO interferometers during the S4 run, along with
the Initial LIGO design sensitivity goal.
waveform
h
×
(
t
). The calibrated response seen by an in-
terferometric gravitational-wave detector is then [18]
h
(
t
) =
F
+
(
t,α,δ,ψ
)
h
+
(
t
) +
F
×
(
t,α,δ,ψ
)
h
×
(
t
)
,
(1)
where
t
is time in the detector frame,
α
is the source
right ascension,
δ
is the source declination,
ψ
is the po-
larization angle of the wave, and
F
+
,
×
are the detector
antenna pattern functions for the two orthogonal polar-
izations. For periodic (nearly pure sinusoidal) gravita-
tional waves, which in general are elliptically polarized,
the individual components
h
+
,
×
have the form
h
+
(
t
)
=
A
+
cos Φ(
t
)
,
(2)
h
×
(
t
)
=
A
×
sin Φ(
t
)
,
(3)
where
A
+
and
A
×
are the amplitudes of the two polariza-
tions, and Φ(
t
) is the phase of the signal at the detector.
(One can also define the initial phase of the signal, Φ
0
,
but in this paper it can be taken to be an unknown and
irrelevant constant).
For an isolated quadrupolar gravitational-wave emit-
ter, characterized by a rotating triaxial ellipsoid mass
distribution, the amplitudes
A
+
and
A
×
are related to
the inclination angle of the source,
ι
, and the wave am-
plitude,
h
0
, by:
A
+
=
1
2
h
0
(
1 + cos
2
ι
)
,
(4)
A
×
=
h
0
cos
ι,
(5)
where
ι
is the angle of its spin axis with respect to the
line of sight between source and detector. For such a star,
the gravitational-wave frequency,
f
, is twice the rotation
frequency,
ν
, and the amplitude
h
0
is given by
h
0
=
16
π
2
G
c
4
Iν
2
d
.
(6)
Here
d
is the distance to the star,
I
is the principal mo-
ment of inertia with respect to its spin axis, and

is the
equatorial ellipticity of the star [18]. Assuming that all of
the frequency’s derivative,
̇
f
, is due to emission of grav-
itational radiation and that
I
takes the canonical value
10
38
kgm
2
, we can relate

to
f
and
̇
f
and use Eq. (6) to
obtain
h
sd
= 4
.
54
×
10
24
(
1 kpc
d
)(
250 yr
f/
(4
̇
f
)
)
1
2
,
(7)
by eliminating

, or

sd
= 7
.
63
×
10
5
(
̇
f
10
10
Hz s
1
)
1
2
(
100 Hz
f
)
5
2
,
(8)
by eliminating
d
. These are referred to, respectively,
as the
spin-down limits
on strain and ellipticity. (See
Eqs. (8), (9), and (19) of [6] for more details of the deriva-
tion.)
Note that the methods used in this paper are sen-
sitive to periodic signals from any type of isolated
gravitational-wave source (e.g., freely precessing or os-
cillating neutron stars as well as triaxial ones), though
we present upper limits in terms of
h
0
and

.
Be-
cause we use semi-coherent methods, only the instan-
taneous signal frequency in the detector reference frame,
2
πf
(
t
) =
d
Φ(
t
)
/dt
, needs to be calculated. In the de-
tector reference frame this can, to a very good approx-
imation, be related to the instantaneous SSB-frame fre-
quency
ˆ
f
(
t
) by [7]
f
(
t
)
ˆ
f
(
t
) =
ˆ
f
(
t
)
v
(
t
)
·
ˆn
c
,
(9)
where
v
(
t
) is the detector’s velocity with respect to the
SSB frame, and
ˆn
is the unit-vector corresponding to the
sky-location of the source. In this analysis, we search
for
ˆ
f
(
t
) signals well described by a nominal frequency
ˆ
f
0
at the start of the S4 run
t
0
and a constant first time
derivative
̇
f
, such that
ˆ
f
(
t
) =
ˆ
f
0
+
̇
f
(
t
t
0
)
.
(10)
These equations ignore corrections to the time interval
t
t
0
at the detector compared with that at the SSB and
relativistic corrections. These corrections are negligible
for the one month semi-coherent searches described here,
though the LSC Algorithm Library (LAL) code [20] used
by our searches does provide routines that make all the
corrections needed to provide a timing accuracy of 3
μs
.
(The LAL code also can calculate
f
(
t
) for signals arriving
from periodic sources in binary systems. Including un-
known orbital parameters in the search, however, would
greatly increase the computational cost or require new
methods beyond the scope of this article.)
5
FIG. 2: An illustration of the discrete frequency bins of the
Short Fourier Transform (SFTs) of the data are shown verti-
cally, with the discrete start times of the SFTs shown hori-
zontally. The dark pixels represent a signal in the data. Its
frequency changes with time due to Doppler shifts and in-
trinsic evolution of the source. By sliding the frequency bins,
the power from a source can be lined up and summed after
appropriate weighting or transformation. This is, in essence,
the starting point for all of the semi-coherent search meth-
ods presented here, though the actual implementations differ
significantly.
IV. OVERVIEW OF THE METHODS
A. Similarities and Differences
The three different analysis methods presented here
have many features in common, but also have important
differences, both major and minor. In this Section we
give a brief overview of the methods.
1. The parameter space
All three methods are based on summing measures of
strain power from many SFTs that have been created
from 30-minute intervals of calibrated strain data. Each
method also corrects explicitly for sky-position depen-
dent Doppler modulations of the apparent source fre-
quency due to the Earth’s rotation and its orbital motion
around the SSB, and the frequency’s time derivative, in-
trinsic to the source (see Fig. 2). This requires a search
in a four-dimensional parameter space; a template in the
space refers to a set of values:
λ
=
{
ˆ
f
0
,
̇
f,α,δ
}
. The
third method, PowerFlux, also searches explicitly over
polarization angle, so that
λ
=
{
ˆ
f
0
,
̇
f,α,δ,ψ
}
.
All three methods search for initial frequency
ˆ
f
0
in the
range 50 – 1000 Hz with a uniform grid spacing equal to
the size of an SFT frequency bin,
δf
=
1
T
coh
= 5
.
556
×
10
4
Hz
.
(11)
where
T
coh
is the time-baseline of each SFT. The range
of
ˆ
f
0
is determined by the noise curves of the interferom-
eters, likely detectable source frequencies [21], and limi-
tations due to the increasing computational cost at high
frequencies.
The range of
̇
f
values searched is [
1
×
10
8
, 0] Hz s
1
for the StackSlide and PowerFlux methods and [
2
.
2
×
10
9
, 0] Hz s
1
for the Hough method. The ranges of
̇
f
are determined by the computational cost, as well as by
the low probability of finding an object with
|
̇
f
|
higher
than the values searched—in other words, the ranges of
̇
f
are narrow enough to complete the search in a reasonable
amount of time, yet wide enough to include likely sig-
nals. All known isolated pulsars spin down more slowly
than the two values of
|
̇
f
|
max
used here, and as seen in
the results section, the ellipticity required for higher
|
̇
f
|
is improbably high for a source losing rotational energy
primarily via gravitational radiation at low frequencies.
A small number of isolated pulsars in globular clusters
exhibit slight spin-up, believed to arise from acceleration
in the Earth’s direction; such spin-up values have magni-
tudes small enough to be detectable with the zero-spin-
down templates used in these searches, given a strong
enough signal. The parameter ranges correspond to a
minimum spin-down timescale
f/
|
4
̇
f
|
(the gravitational-
wave spin-down age) of 40 years for a source emitting at
50 Hz and 800 years for a source at 1000 Hz. Since for
known pulsars [22] this characteristic timescale is at least
hundreds of years for frequencies on the low end of our
range and tens of millions of years for frequencies on the
high end, we see again that the ranges of
|
̇
f
|
are wide
enough to include sources from this population.
As discussed in our previous reports [6, 7], the number
of sky points that must be searched grows quadratically
with the frequency
ˆ
f
0
, ranging here from about five thou-
sand at 50 Hz to about two million at 1000 Hz. All three
methods use nearly isotropic grids which cover the en-
tire sky. The PowerFlux search also divides the sky into
regions according to susceptibility to stationary instru-
mental line artifacts. Sky grid and spin-down spacings
and other details are provided below.
2. Upper limits
While the parameter space searched is similar for the
three methods, there are important differences in the way
upper limits are set. StackSlide and Hough both set
population-based frequentist limits on
h
0
by carrying out
Monte Carlo simulations of a random population of pul-
sar sources distributed uniformly over the sky and with
isotropically distributed spin-axes. PowerFlux sets strict
frequentist limits on circular and linear polarization am-
plitudes
h
Circ
limit
0
and
h
Lin
limit
0
, which correspond to
limits on most and least favorable pulsar inclinations,
respectively. The limits are placed separately on tiny
patches of the sky, with the highest strain upper limits
presented here. In this context “strict” means that, re-
gardless of its polarization angle
ψ
or inclination angle
ι
,
regardless of its sky location (within fiducial regions dis-
cussed below), and regardless of its frequency value and
spin-down within the frequency and spin-down step sizes
of the search template, an isolated pulsar of true strain
6
amplitude
h
0
= 2
h
Lin
limit
0
, would have yielded a higher
measured amplitude than what we measure, in at least
95% of independent observations. The circular polariza-
tion limits
h
Circ
limit
0
apply only to the most favorable
inclinations (
ι
0,
π
), regardless of sky location and
regardless of frequency and spin-down, as above.
Due to these different upper limit setting methods,
sharp instrumental lines are also handled differently.
StackSlide and Hough carry out removal of known instru-
mental lines of varying widths in individual SFTs. The
measured powers in those bins are replaced with random
noise generated to mimic the noise observed in neigh-
boring bins. This line cleaning technique can lead to a
true signal being missed because its apparent frequency
may coincide with an instrumental line for a large num-
ber of SFTs. However, population-averaged upper limits
are determined self-consistently to include loss of detec-
tion efficiency due to line removal, by using Monte Carlo
simulations.
Since its limits are intended to be strict, that is, valid
for any source inclination and for any source location
within its fiducial area, PowerFlux must handle instru-
mental lines differently. Single-bin lines are flagged dur-
ing data preparation so that when searching for a partic-
ular source an individual SFT bin power is ignored when
it coincides with the source’s apparent frequency. If more
than 80% of otherwise eligible bins are excluded for this
reason, no attempt is made to set a limit on strain power
from that source. In practice, however, the 80% cutoff is
not used because we have found that all such sources lie
in certain unfavorable regions of the sky, which we call
“skybands” and which we exclude when setting upper
limits. These skybands depend on source frequency and
its derivative, as described in Sec. V D 4.
3. Data Preparation
Other differences among the methods concern the data
windowing and filtering used in computing Fourier trans-
forms and concern the noise estimation. StackSlide and
Hough apply high pass filters to the data above 40Hz, in
addition to the filter used to produce the calibrated data
stream, and use Tukey windowing. PowerFlux applies no
additional filtering and uses Hann windowing with 50%
overlap between adjacent SFT’s. StackSlide and Hough
use median-based noise floor tracking [23, 24, 25]. In
contrast, Powerflux uses a time-frequency decomposition.
Both of these noise estimation methods are described in
Sec. V.
The raw, uncalibrated data channels containing the
strain measurements from the three interferometers are
converted to a calibrated “
h
(
t
)” data stream, following
the procedure described in [26], using calibration refer-
ence functions described in [27]. SFTs are generated di-
rectly from the calibrated data stream, using 30-minute
intervals of data for which the interferometer is operat-
ing in what is known as science-mode. The choice of 30
minutes is a tradeoff between intrinsic sensitivity, which
increases with SFT length, and robustness against fre-
quency drift during the SFT interval due to the Earth’s
motion, source spin-down, and non-stationarity of the
data [7]. The requirement that each SFT contain contigu-
ous data at nominal sensitivity introduces duty factor loss
from edge effects, especially for the Livingston interfer-
ometer (
'
20%) which had typically shorter contiguous-
data stretches. In the end, the StackSlide and Hough
searches used 1004 SFTs from H1 and 899 from L1,
the two interferometers with the best broadband sen-
sitivty. For PowerFlux, the corresponding numbers of
overlapped SFTs were 1925 and 1628. The Hough search
also used 1063 H2 SFTs. In each case, modest require-
ments were placed on data quality to avoid short periods
with known electronic saturations, unmonitored calibra-
tion strengths, and the periods immediately preceding
loss of optical cavity resonance.
B. Definitions And Notation
Let
N
be the number of SFTs,
T
coh
the time-baseline
of each SFT, and
M
the number of uniformly spaced
data points in the time domain from which the SFT is
constructed. If the time series is denoted by
x
j
(
j
=
0
,
1
,
2
...M
1), then our convention for the discrete
Fourier transform is
̃
x
k
= ∆
t
M
1
j
=0
x
j
e
2
π
i
jk/M
,
(12)
where
k
= 0
,
1
,
2
...
(
M
1), and ∆
t
=
T
coh
/M
. For
0
k
M/
2, the frequency index
k
corresponds to a
physical frequency of
f
k
=
k/T
coh
.
In each method, the “power” (in units of spectral den-
sity) associated with frequency bin
k
and SFT
i
is taken
to be
P
i
k
=
2
|
̃
x
i
k
|
2
T
coh
.
(13)
It proves convenient to define a normalized power by
ρ
i
k
=
P
i
k
S
i
k
.
(14)
The quantity
S
i
k
is the single-sided power spectral density
of the detector noise at frequency
f
k
, the estimation of
which is described below. Furthermore, a threshold,
ρ
th
,
can be used to define a
binary count
by [10]:
n
i
k
=
{
1 if
ρ
i
k
ρ
th
0 if
ρ
i
k
< ρ
th
.
(15)
When searching for a signal using template
λ
the de-
tector antenna pattern and frequency of the signal are
found at the midpoint time of the data used to gener-
ate each SFT. Frequency dependent quantities are then
7
Quantity
Description
P
i
Power for SFT
i
& template
λ
ρ
i
Normalized power for SFT
i
& template
λ
n
i
Binary count for SFT
i
& template
λ
S
i
Power spect. noise density for SFT
i
& template
λ
F
i
+
F
+
at midpoint of SFT
i
for template
λ
F
i
×
F
×
at midpoint of SFT
i
for template
λ
TABLE I: Summary of notation used.
evaluated at a frequency index
k
corresponding to the bin
nearest this frequency. To simplify the equations in the
rest of this paper we drop the frequency index
k
and use
the notation given in Table I to define various quantities
for SFT
i
and template
λ
.
C. Basic StackSlide, Hough, and PowerFlux
Formalism
We call the detection statistics used in this search the
“StackSlide Power”,
P
, the “Hough Number Count”,
n
,
and the “PowerFlux Signal Estimator”,
R
. The basic
definitions of these quantities are given below.
Here the simple StackSlide method described in [15]
is used; the “StackSlide Power” for a given template is
defined as
P
=
1
N
N
1
i
=0
ρ
i
,
(16)
This normalization results in values of
P
with a mean
value of unity and, for Gaussian noise, a standard devi-
ation of 1
/
N
. Details about the value and statistics of
P
in the presence and absence of a signal are given in
Appendix B and [15].
In the Hough search, instead of summing the normal-
ized power, the final statistic used in this paper is a
weighted sum of the binary counts, giving the “Hough
Number Count”:
n
=
N
1
i
=0
w
i
n
i
.
(17)
where the Hough weights are defined as
w
i
1
S
i
{
(
F
i
+
)
2
+
(
F
i
×
)
2
}
,
(18)
and the weight normalization is chosen according to
N
1
i
=0
w
i
=
N .
(19)
With this choice of normalization the Hough Number
Count
n
lies within the range [0
,N
]. Thus, we take a
binary count
n
i
to have greater weight if the SFT
i
has a
lower noise floor and if, in the time-interval correspond-
ing to this SFT, the beam pattern functions are larger
for a particular point in the sky. Note that the sensitiv-
ity of the search is governed by the ratios of the different
weights, not by the choice of overall scale. In the next sec-
tion we show that these weights maximize the sensitivity,
averaged over the orientation of the source. This choice
of
w
i
was originally derived in [16] using a different argu-
ment and is similar to that used in the PowerFlux circular
polarization projection described next. More about the
Hough method is given in [7, 10].
The PowerFlux method takes advantage of the fact
that less weight should be given to times of greater noise
variance or smaller detector antenna response to a signal.
Noting that power estimated from the data divided by
the antenna pattern increases the variance of the data at
times of small detector response, the problem reduces to
finding weights that minimize the variance, or in other
words that maximize the signal-to-noise ratio. The re-
sulting PowerFlux detection statistic is [17],
R
=
2
T
coh
N
1
i
=0
W
i
P
i
/
(
F
i
ψ
)
2
N
1
i
=0
W
i
,
(20)
where the PowerFlux weights are defined as
W
i
= [(
F
i
ψ
)
2
]
2
/S
2
i
,
(21)
and where
(
F
i
ψ
)
2
=
{
(
F
i
+
)
2
linear polarization
(
F
i
+
)
2
+ (
F
i
×
)
2
circular polarization
.
(22)
As noted previously, the PowerFlux method searches us-
ing four linear polarization projections and one circular
polarization projection. For the linear polarization pro-
jections, note that (
F
i
+
)
2
is evaluated at the angle
ψ
,
which is the same as (
F
i
×
)
2
evaluated at the angle
ψ
π/
4;
for circular polarization, the value of (
F
i
+
)
2
+ (
F
i
×
)
2
is in-
dependent of
ψ
. Finally note that the factor of 2
/T
coh
in
Eq. (20) makes
R
dimensionless and is chosen to make it
directly related to an estimate of the squared amplitude
of the signal for the given polarization. Thus
R
is also
called in this paper the “PowerFlux Signal Estimator”.
(See [17] and Appendix A for further discussion.)
We have shown in Eqs. (16)-(22) how to compute the
detection statistic (or signal estimator) for a given tem-
plate. The next section gives the details of the imple-
mentation and pipelines used, where these quantities are
calculated for a set of templates
λ
and analyzed.
V. IMPLEMENTATIONS AND PIPELINES
A. Running Median Noise Estimation
The implementations of the StackSlide and Hough
methods described below use a “running median” to esti-
mate the mean power and, from this estimate, the power
8
spectral density of the noise, for every frequency bin of
every SFT. PowerFlux uses a different noise decomposi-
tion method described in its implementation section be-
low.
Note that for Gaussian noise, the single-sided power
spectral density can be estimated using
S
i
k
=
2
〈|
̃
x
i
k
|
2
T
coh
(23)
where the angle brackets represent an ensemble average.
The estimation of
S
i
k
must guard against any biases in-
troduced by the presence of a possible signal and also
against narrow spectral disturbances. For this reason the
mean,
〈|
̃
x
i
k
|
2
, is estimated via the median. We assume
that the noise is stationary within a single SFT, but al-
low for non-stationarities across different SFTs. In every
SFT we calculate the “running median” of
|
̃
x
i
k
|
2
for ev-
ery 101 frequency bins centered on the
k
th
bin, and then
estimate
〈|
̃
x
i
k
|
2
[23, 24, 25] by dividing by the expected
ratio of the median to the mean.
Note, however, that in the StackSlide search, after
the estimated mean power is used to compute
S
i
k
in
the denominator of Eq. (14) these terms are summed in
Eq. (16), while the Hough search applies a cutoff to ob-
tain binary counts in Eq. (15) before summing. This re-
sults in the use of a different correction to get the mean
in the StackSlide search from that used in the Hough
search. For a running median using 101 frequency bins,
the effective ratio of the median to mean used in the
StackSlide search was 0
.
691162 (which was chosen to nor-
malize the data so that the mean value of the StackSlide
Power equals one) compared with the expected ratio for
an exponential distribution of 0
.
698073 used in the Hough
search (which is explained in Appendix A of [7]). It is im-
portant to realize that the results reported here are valid
independent of the factor used, since any overall constant
scaling of the data does not affect the selection of outliers
or the reported upper limits, which are based on Monte
Carlo injections subjected to the same normalization.
B. The StackSlide Implementation
1. Algorithm and parameter space
The StackSlide method uses power averaging to gain
sensitivity by decreasing the variance of the noise [12, 13,
14, 15]. Brady and Creighton [12] first described this ap-
proach in the context of gravitational-wave detection as
a part of a hierarchical search for periodic sources. Their
method consists of averaging the power from a demodu-
lated time series, but as an approximation did not include
the beam pattern response of the detector. In Ref. [15],
a simple implementation is described that averages the
normalized power given in Eq. (14). Its extension to aver-
aging the maximum likelihood statistic (known as the
F
-
statistic) which does include the beam pattern response
FIG. 3: Flow chart for the pipeline used to find the upper
limits presented in this paper using the StackSlide method.
is mentioned in Ref. [15] (see also [6, 10, 18]), and further
extensions of the StackSlide method are given in [13].
As noted above, the simple StackSlide method given
in [15] is used here and the detection statistic, called the
“StackSlide Power”, is defined by Eq. (16). The normal-
ization is chosen so that the mean value of
P
is equal
to 1 and its standard deviation is 1
/
N
for Gaussian
noise alone. For simplicity, the StackSlide Power signal-
to-noise ratio (in general the value of
P
minus its mean
value and then divided by the standard deviation of
P
)
will be defined in this paper as (
P
1)
N
, even for non-
Gaussian noise.
The StackSlide code, which implements the method de-
scribed above, is part of the C-based LSC Algorithms Li-
brary Applications (LALapps) stored in the lscsoft CVS
repository [20]. The code is run in a pipeline with options
set to produce the results from a search and from Monte
Carlo simulations. Parallel jobs are run on computer
clusters within the LSC, in the Condor environment [29],
and the final post processing steps are performed using
Matlab [30]. The specific StackSlide pipeline used to find
the upper limits presented in this paper is shown in Fig. 3.
The first three boxes on the left side of the pipeline can
also be used to output candidates for follow-up searches.
A separate search was run for each successive 0
.
25 Hz
band within 50
1000 Hz. The spacing in frequency used
is given by Eq. (11). The spacing in
̇
f
was chosen as that
which changes the frequency by one SFT frequency bin
during the observation time
T
obs
, i.e., so that
̇
fT
obs
=
δf
.
For simplicity
T
obs
= 2
.
778
×
10
6
seconds
'
32
.
15 days
was chosen, which is greater than or equal to
T
obs
for each
interferometer. Thus, the
̇
f
part of the parameter space
was over-covered by choosing
|
δ
̇
f
|
=
δf
T
obs
=
1
T
coh
T
obs
= 2
×
10
10
Hz s
1
.
(24)
9
Values of
̇
f
in the range [
1
×
10
8
Hz s
1
,
0 Hz s
1
]
were searched. This range corresponds to a search over
51 values of
̇
f
, which is the same as PowerFlux used in
its low-frequency search (discussed in Section. V D).
The sky grid used is similar to that used for the all-
sky search in [6], but with a spacing between sky-grid
points appropriate for the StackSlide search. This grid is
isotropic on the celestial sphere, with an angular spacing
between points chosen for the 50-225 Hz band, such that
the maximum change in Doppler shift from one sky grid
point to the next would shift the frequency by half a bin.
This is given by
δθ
0
=
0
.
5
cδf
ˆ
f
(
v
sin
θ
)
max
= 9
.
3
×
10
3
rad
(
300Hz
ˆ
f
)
,
(25)
where
v
is the magnitude of the velocity
v
of the detector
in the SSB frame, and
θ
is the angle between
v
and the
unit-vector
ˆn
giving the sky-position of the source. Equa-
tions (24) and (25) are the same as Eqs. (19) and (22) in
[7], which represent conservative choices that over-cover
the parameter space. Thus, the parameter space used
here corresponds to that in Ref. [7], adjusted to the S4
observation time, and with the exception that a stere-
ographic projection of the sky is not used. Rather an
isotropic sky grid is used like the one used in [6].
One difficulty is that the computational cost of the
search increases quadratically with frequency, due to the
increasing number of points on the sky grid. To reduce
the computational time, the sky grid spacing given in
Eq. (25) was increased by a factor of 5 above 225 Hz.
This represents a savings of a factor of 25 in computa-
tional cost. It was shown through a series of simulations,
comparing the upper limits in various frequency bands
with and without the factor of 5 increase in grid spacing,
that this changes the upper limits on average by less than
than 0
.
3%, with a standard deviation of 2%. Thus, this
factor of 5 increase was used to allow the searches in the
225
1000 Hz band to complete in a reasonable amount
of time.
It is not surprising that the sky grid spacing can be in-
creased, for at least three reasons. First, the value for
δθ
0
given in Eq. (25) applies to only a small annular region on
the sky, and is smaller than the average change. Second,
only the net change in Doppler shift during the observa-
tion time is important, which is less than the maximum
Doppler shift due to the Earth’s orbital motion during
a one month run. (If the Doppler shift were constant
during the entire observation time, one would not need
to search sky positions even if the Doppler shift varied
across the sky. A source frequency would be shifted by
a constant amount during the observation, and would be
detected, albeit in a frequency bin different from that
at the SSB.) Third, because of correlations on the sky,
one can detect a signal with negligible loss of SNR much
farther from its sky location than the spacing above sug-
gests.
FIG. 4: The StackSlide Power for the 145
155 Hz band with
no sliding. Harmonics of 1 Hz instrumental lines are clearly
seen in H1 (top) and L1 (bottom). These lines are removed
from the data by the StackSlide and Hough searches using the
method described in the text, while PowerFlux search tracks
these lines and avoids them when setting upper limits.
FIG. 5: The L1 amplitude spectral density in a narrow fre-
quency band estimated from 10 SFTs before and after the line
cleaning used by the StackSlide pipeline. In the band shown,
the 150 Hz bin, and one bin either side of this bin have been
replaced with estimates of the noise based on neighboring
bins.
2. Line cleaning
Coherent instrumental lines exist in the data which
can mimic a continuous gravitational-wave signal for pa-
rameter space points that correspond to little Doppler
modulation. Very narrow instrumental lines are removed
(“cleaned”) from the data. In the StackSlide search, a
line is considered “narrow” if its full width is less than
10
IFO
f
start
f
step
Num. ∆
f
left
f
right
Description
Hz
Hz
Hz
Hz
H1
46.7
1
0.0
0.0
Cal. Line
H1 393.1
1
0.0
0.0
Cal. Line
H1 973.3
1
0.0
0.0
Cal. Line
H1 1144.3 —
1
0.0
0.0
Cal. Line
H1
0.0
1.0 1500 0.0006 0.0006 1 Hz Comb
L1
54.7
1
0.0
0.0
Cal. Line
L1 396.7
1
0.0
0.0
Cal. Line
L1 1151.5 —
1
0.0
0.0
Cal. Line
L1
0.0
1.0 1500 0.0006 0.0006 1 Hz Comb
TABLE II: Instrumental lines cleaned during the StackSlide
search. The frequencies cleaned are found by starting with
that given in the first column, and then taking steps in fre-
quency given in the second column, repeating this the num-
ber of times shown in the third column; the fourth and fifth
columns show how many additional Hz are cleaned to the
immediate left and right of each line.
Excluded Bands
Description
Hz
[57
,
63)
Power lines
[
n
60
1
,n
60 + 1)
n
= 2 to 16 Power line harmonics
[340
,
350)
Violin modes
[685
,
690)
Violin mode harmonics
[693
,
696)
Violin mode harmonics
TABLE III: Frequency bands excluded from the StackSlide
search.
5% of the 0
.
25 Hz band, or less than 0
.
0125 Hz. The line
must also have been identified
a priori
as a known instru-
ment artifact. Known lines with less than this width were
cleaned by replacing the contents of bins corresponding
to lines with random values generated by using the run-
ning median to find the mean power using 101 bins from
either side of the lines. This method is also used to esti-
mate the noise, as described in Section V A.
It was found when characterizing the data that a comb
of narrow 1 Hz harmonics existed in the H1 and L1 data,
as shown in Fig. 4. Table II shows the lines cleaned dur-
ing the StackSlide search. As the table shows, only this
comb of narrow 1 Hz harmonics and injected lines used
for calibration were removed. As an example of the clean-
ing process, Fig. 5 shows the amplitude spectral density
estimated from 10 SFTs before and after line cleaning,
for the band with the 1 Hz line at 150 Hz.
The cleaning of very narrow lines has a negligible effect
on the efficiency to detect signals. Very broad lines, on
the other hand, cannot be handled in this way. Bands
with very broad lines were searched without any line
cleaning. There were also a number of highly disturbed
bands, dominated either by the harmonics of 60 Hz power
lines or by the violin modes of the suspended optics, that
were excluded from the StackSlide results. (Violin modes
FIG. 6: Measure confidence vs.
h
0
for an example band (140
140
.
25 Hz in H1). A best-fit straight line is used to find the
value of
h
0
corresponding to 95% confidence and to estimate
the uncertainties in the results (see text).
refer to resonant excitations of the steel wires that sup-
port the interferometer mirrors.) These are shown in
Table III. While these bands can be covered by adjusting
the parameters used to find outliers and set upper limits,
we will wait for future runs to do this.
3. Upper limits method
After the lines are cleaned, the powers in the SFTs
are normalized and the parameter space searched, with
each template producing a value of the StackSlide Power,
defined in Eq. (16). For this paper, only the “loudest”
StackSlide Power is kept, resulting in a value
P
max
for
each 0
.
25 Hz band, and these are used to set upper lim-
its on the gravitational-wave amplitude,
h
0
. (The loud-
est coincident outliers are also identified, but none sur-
vive as candidates after follow-up studies described in
Sec. VII A 1.) The upper limits are found by a series of
Monte Carlo simulations, in which signals are injected
in software with a fixed value for
h
0
, but with otherwise
randomly chosen parameters, and the parameter space
points that surround the injection are searched. The
number of times the loudest StackSlide Power found dur-
ing the Monte Carlo simulations is greater than or equal
to
P
max
is recorded, and this is repeated for a series of
h
0
values. The 95% confidence upper limit is defined to
be the value of
h
0
that results in a detected StackSlide
Power greater than or equal to
P
max
95% of the time.
As shown in Fig. 3, the line cleaning described above
is done after each injection is added to the input data,
which folds any loss of detection efficiency due to line
cleaning into the upper limits self-consistently.
Figure 6 shows the measured confidence versus
h
0
for
11
an example frequency band. The upper limit finding pro-
cess involves first making an initial guess of its value, then
refining this guess using a single set of injections to find
an estimate of the upper limit, and finally using this es-
timate to run several sets of injections to find the final
value of the upper limit. These steps are now described
in detail.
To start the upper limit finding process, first an initial
guess,
h
guess
0
, is used as the gravitational-wave amplitude.
The initial guess need not be near the sought-after upper
limit, just sufficiently large, as explained below. A single
set of
n
injections is done (specifically
n
= 3000 was used)
with random sky positions and isotropically distributed
spin axes, but all with amplitude
h
guess
0
. The output
list of StackSlide Powers from this set of injections is
sorted in ascending order and the 0
.
05
n
’th (specifically
for
n
= 3000 the 150th) smallest value of the StackSlide
Power is found, which we call
P
0
.
05
, Note that the goal is
to find the value of
h
0
that makes
P
0
.
05
=
P
max
, so that
95% of the output powers are greater than the maximum
power found during the search. This is what we call the
95% confidence upper limit. Of course, in general
P
0
.
05
will not equal
P
max
unless our first guess was very lucky.
However, as per the discussion concerning Eq. (B5),
P
1
is proportional to
h
2
0
(i.e, removing the mean value due
to noise leaves on average the power due to the presence
of a signal). Thus, an estimate of the 95%
h
0
confidence
upper limits is given by the following rescaling of
h
guess
0
,
h
est
0
=
P
max
1
P
0
.
05
1
h
guess
0
.
(26)
Thus an estimated upper limit,
h
est
0
, is found from a sin-
gle set of injections with amplitude
h
guess
0
; the only re-
quirement is that
h
guess
0
is chosen loud enough to make
P
0
.
05
>
1.
It is found that using Eq. (26) results in a estimate of
the upper limit that is typically within 10% of the final
value. For example, the estimated upper limit found in
this way is indicated by the circled point in Fig. 6. The
value of
h
est
0
then becomes the first value for
h
0
in a series
of Monte Carlo simulations, each with 3000 injections,
which use this value and 8 neighboring values, measur-
ing the confidence each time. The Matlab [30] polyfit
and polyval functions are then used to find the best-fit
straight line to determine the value of
h
0
corresponding
to 95% confidence and to estimate the uncertainties in
the results. This is the final step of the pipeline shown
in Fig. 3.
C. The Hough Transform Implementation
1. Description of Algorithm
The Hough transform is a general method for pattern
recognition, invented originally to analyze bubble cham-
ber pictures from CERN [31, 32]; it has found many
applications in the analysis of digital images [33]. This
method has already been used to analyze data from the
second science run (S2) of the LIGO detectors [7] and a
detailed description can be found in [10]. Here we present
only a brief description, emphasizing the differences be-
tween the previous S2 search and the S4 search described
here.
The Hough search uses a weighted sum of the binary
counts as its final statistic, as given by Eqs. (15) and (19).
In the standard Hough search as presented in [7, 10], the
weights are all set to unity. The weighted Hough trans-
form was originally discussed in [16]. The software for
performing the Hough transform has been adapted to use
arbitrary weights without any significant loss in compu-
tational efficiency. Furthermore, the robustness of the
Hough transform method in the presence of strong tran-
sient disturbances is not compromised by using weights
because each SFT contributes at most
w
i
(which is of
order unity) to the final number count.
The following statements can be proven using the
methods of [10]. The mean number count in the ab-
sence of a signal is ̄
n
=
Np
, where
N
is the number
of SFTs and
p
is the probability that the normalized
power, of a given frequency bin and SFT defined by
Eq. (14), exceeds a threshold
ρ
th
, i.e.,
p
is the proba-
bility that a frequency bin is selected in the absence of
a signal. For unity weighting, the standard deviation is
simply
σ
=
Np
(1
p
). However, with more general
weighting, it can be shown that
σ
is given by
σ
=
||
w
||
2
p
(1
p
)
,
(27)
where
||
w
||
2
=
N
1
i
=0
w
2
i
. A threshold
n
th
on the number
count corresponding to a false alarm rate
α
H
is given by
n
th
=
Np
+
2
||
w
||
2
p
(1
p
) erfc
1
(2
α
H
)
.
(28)
Therefore
n
th
depends on the weights of the correspond-
ing template
λ
. In this case, the natural detection statis-
tic is not the “Hough Number Count”
n
, but the
signifi-
cance
of a number count, defined by
s
=
n
̄
n
σ
,
(29)
where ̄
n
and
σ
are the expected mean and standard de-
viation for pure noise. Values of
s
can be compared di-
rectly across different templates characterized by differ-
ing weight distributions.
The threshold
ρ
th
(c.f. Eq. 15) is selected to give the
minimum false dismissal probability
β
H
for a given false
alarm rate. In [7] it was shown that the optimal choice for
ρ
th
is 1
.
6 which correspond to a peak selection probability
p
=
e
ρ
th
0
.
2. It can be shown that the optimal choice
is unchanged by the weights and hence
ρ
th
= 1
.
6 is used
once more [34].
Consider a population of sources located at a given
point in the sky, but having uniformly distributed spin
axis directions. For a template that is perfectly matched
12
FIG. 7: The improvement in the significance as a function of
the mismatch in the sky-position. A signal is injected in fake
noise at
α
=
δ
= 0 and the weights are calculated at
α
=
δ
=
δθ
. The curve is the observed significance as a function of
δθ
while the horizontal line is the observed significance when no
weights are used. See main text for more details.
in frequency, spin-down, and sky-position, and given the
optimal peak selection threshold, it can be shown [34]
that the weakest signal that can cross the threshold
n
th
with a false dismissal probability
β
H
has an amplitude
h
0
= 3
.
38
S
1
/
2
(
||
w
||
w
·
X
)
1
/
2
1
T
coh
,
(30)
where
S
= erfc
1
(2
α
H
) + erfc
1
(2
β
H
)
,
(31)
X
i
=
1
S
i
{
(
F
i
+
)
2
+
(
F
i
×
)
2
}
.
(32)
As before,
F
i
+
and
F
i
×
are the values of the beam pat-
tern functions at the mid-point of the
i
th
SFT. To derive
(30) we have assumed that the number of SFTs
N
is
sufficiently large and that the signal is weak [10].
From (30) it is clear that the scaling of the weights
does not matter;
w
i
kw
i
leaves
h
0
unchanged for any
constant
k
. More importantly, it is also clear that the
sensitivity is best, i.e.
h
0
is minimum, when
w
·
X
is
maximum:
w
i
X
i
.
(33)
This result is equivalent to Eq. (18).
In addition to improving sensitivity in single-
interferometer analysis, the weighted Hough method al-
lows automatic optimal combination of Hough counts
from multiple interferometers of differing senstivities.
Ideally, to obtain the maximum increase in sensitiv-
ity, we should calculate the weights for each sky-location
separately. In practice, we break up the sky into smaller
patches and calculate one weight for each sky-patch cen-
ter. The gain from using the weights will be reduced if the
sky patches are too large. From equation (32), it is clear
that the dependence of the weights on the sky-position is
only through the beam pattern functions. Therefore, the
sky patch size is determined by the typical angular scale
over which
F
+
and
F
×
vary; thus for a spherical detector
using the beam pattern weights would not gain us any
sensitivity. For the LIGO interferometers, we have in-
vestigated this issue with Monte-Carlo simulations using
random Gaussian noise. Signals are injected in this noise
corresponding to the H1 interferometer at a sky-location
(
α
0
0
), while the weights are calculated at a mismatched
sky-position (
α
0
+
δθ,δ
0
+
δθ
). The significance values are
compared with the significance when no weights are used.
An example of such a study is shown in Fig. 7. Here, we
have injected a signal at
α
=
δ
= 0, cos
ι
= 0
.
5, zero
spin-down, Φ
0
=
ψ
= 0, and a signal to noise ratio corre-
sponding approximately to a 6-
σ
level without weights.
The figure shows a gain of
10% at
δθ
= 0, decreas-
ing to zero at
δθ
0
.
3 rad. We get qualitatively similar
results for other sky-locations, independent of frequency
and other parameters. There is an additional gain due
to the non-stationarity of the noise itself, which depends,
however, on the quality of the data. In practice, we have
chosen to break the sky up into 92 rectangular patches
in which the average sky patch size is about 0
.
4 rad wide,
corresponding to a maximum sky position mismatch of
δθ
= 0
.
2 rad in Fig. 7.
2. The Hough Pipeline
The Hough analysis pipeline for the search and for
setting upper limits follows roughly the same scheme as
in [7]. In this section we present a short description of
the pipeline, mostly emphasizing the differences from [7]
and from the StackSlide and PowerFlux searches. As
discussed in the previous subsection, the key differences
from the S2 analysis [7] are (i) using the beam-pattern
and noise weights, and (ii) using SFTs from multiple in-
terferometers.
The total frequency range analyzed is 50-1000 Hz, with
a resolution
δf
= 1
/T
coh
as in (11). The resolution in
̇
f
is
2
.
2
×
10
10
Hz s
1
given in (24), and the reference time
for defining the spin-down is the start-time of the obser-
vation. However, unlike StackSlide and PowerFlux, the
Hough search is carried out over only 11 values of
̇
f
, in-
cluding zero, in the range [
2
.
2
×
10
9
Hz s
1
, 0 Hz s
1
].
This choice is driven by the technical design of the cur-
rent implementation, which uses look-up-tables and par-
tial Hough maps as in [7]. This implementation of the
Hough algorithm is efficient when analyzing all resolv-
able points in
̇
f
, as given in (24), but this approach is
incompatible with the larger
̇
f
step sizes used in the other
search methods, which permit those searches to search a
larger
̇
f
range for comparable computational cost.
The sky resolution is similar to that used by the Stack-
Slide method for
f <
225 Hz as given by (25). At fre-
quencies higher than this, the StackSlide sky-resolution
is 5 times coarser, thus the Hough search is analyzing
about 25 more templates at a given frequency and spin-
13
FIG. 8: Two example histograms of the normalized Hough
number count compared to a Gaussian distribution for the H1
detector in the frequency band 150-151 Hz. The upper figure
corresponds to a a patch located at the north pole for the
case in which the weights are used. The number of templates
analyzed in this 1Hz band is of 11
×
10
6
, the number of SFTs
1004, the corresponding mean ̄
n
= 202
.
7 and
σ
= 12
.
94 is
obtained from the weights. The lower figure corresponds to
a patch at the equator using the same data. In this case the
number of templates analyzed in this 1Hz band is of 10
.
5
×
10
6
,
and its corresponding
σ
= 14
.
96.
down value. In each of the 92 sky patches, by means of
the stereographic projection, the sky patch is mapped to
a two dimensional plane with a uniform grid of that res-
olution
δθ
0
. Sky Patches slightly overlap to avoid gaps
among them (see [7] for further details).
Figure 8 shows examples of histograms of the number
counts in two particular sky patches for the H1 detector
in the 150-151 Hz band. In all the bands free of instru-
mental disturbances, the Hough number count distribu-
tions follows the expected theoretical distribution, which
can be approximated by a Gaussian distribution. Since
the number of SFTs for H1 is 1004, the corresponding
IFO
f
start
f
step
n
f
left
f
right
Description
Hz
Hz
Hz
Hz
H1 392.365
1
0.01
0.01
Cal. SideBand
H1 393.835
1
0.01
0.01
Cal. SideBand
H2
54.1
1
0.0
0.0
Cal. Line
H2
407.3
1
0.0
0.0
Cal. Line
H2
1159.7
1
0.0
0.0
Cal. Line
H2 110.934 36.9787 4
0.02
0.02
37 Hz Oscillator
L1 154.6328 8.1386 110 0.01
0.01
8.14 Hz Comb
L1
0.0
36.8725 50 0.02
0.02 37 Hz Oscillator (*)
TABLE IV: Instrumental lines cleaned during the Hough
search that were not listed in Table II (see text). (*) These
lines were removed only in the multi-interferometer search.
mean ̄
n
= 202
.
7 and the standard deviation is given by
Eq. (27). The standard deviation is computed from the
weights
w
and varies among different sky patches because
of varying antenna pattern functions.
The upper limits on
h
0
are derived from the
loud-
est event
, registered over the entire sky and spin-down
range in each 0
.
25 Hz band, not from the highest number
count. As for the StackSlide method, we use a frequentist
method, where upper limits refer to a hypothetical pop-
ulation of isolated spinning neutron stars which are uni-
formly distributed in the sky and have a spin-down rate
̇
f
uniformly distributed in the range [
2
.
2
×
10
9
Hz s
1
,
0 Hz s
1
]. We also assume uniform distributions for the
parameters cos
ι
[
1
,
1],
ψ
[0
,
2
π
], and Φ
0
[0
,
2
π
].
The strategy for calculating the 95% upper limits is
roughly the same scheme as in [7], except for the treat-
ment of narrow instrumental lines.
Known spectral disturbances are removed from the
SFTs in the same way as for the StackSlide search. The
known spectral lines are, of course, also consistently re-
moved after each signal injection when performing the
Monte-Carlo simulations to obtain the upper limits.
The narrow instrumental lines “cleaned” from the SFT
data are the same ones cleaned during the StackSlide
search shown in Table II, together with ones listed in Ta-
ble IV. The additional lines listed in Table IV are cleaned
to prevent large artifacts in one instrument from increas-
ing the false alarm rate of the Hough multi-interferometer
search. Note that the L1 36.8725 Hz comb was eliminated
mid-way through the S4 run by replacing a synthesized
radio frequency oscillator for phase modulation with a
crystal oscillator, and these lines were not removed in
the Hough L1 single-interferometer analysis.
No frequency bands have been excluded from the
Hough search, although the upper limits reported on the
bands shown in Table III, that are dominated by 60 Hz
power line harmonics or violin modes of the suspended
optics, did not always give satisfactory convergence to an
upper limit. In a few of these very noisy bands, upper
limits were set by extrapolation, instead of interpolation,
of the Monte-Carlo simulations. Therefore the results
14