of 38
All-sky search for periodic gravitational waves in LIGO S4 data
B. Abbott,
15
R. Abbott,
15
R. Adhikari,
15
J. Agresti,
15
P. Ajith,
2
B. Allen,
2,50
R. Amin,
19
S. B. Anderson,
15
W. G. Anderson,
50
M. Arain,
39
M. Araya,
15
H. Armandula,
15
M. Ashley,
4
S. Aston,
38
P. Aufmuth,
14
C. Aulbert,
1
S. Babak,
1
S. Ballmer,
15
H. Bantilan,
8
B. C. Barish,
15
C. Barker,
16
D. Barker,
16
B. Barr,
40
P. Barriga,
49
M. A. Barton,
40
K. Bayer,
18
K. Belczynski,
25
J. Betzwieser,
18
P. T. Beyersdorf,
28
B. Bhawal,
15
I. A. Bilenko,
22
G. Billingsley,
15
R. Biswas,
50
E. Black,
15
K. Blackburn,
15
L. Blackburn,
18
D. Blair,
49
B. Bland,
16
J. Bogenstahl,
40
L. Bogue,
17
R. Bork,
15
V. Boschi,
15
S. Bose,
51
P. R. Brady,
50
V. B. Braginsky,
22
J. E. Brau,
43
M. Brinkmann,
2
A. Brooks,
37
D. A. Brown,
15,6
A. Bullington,
31
A. Bunkowski,
2
A. Buonanno,
41
O. Burmeister,
2
D. Busby,
15
R. L. Byer,
31
L. Cadonati,
18
G. Cagnoli,
40
J. B. Camp,
23
J. Cannizzo,
23
K. Cannon,
50
C. A. Cantley,
40
J. Cao,
18
L. Cardenas,
15
M. M. Casey,
40
G. Castaldi,
46
C. Cepeda,
15
E. Chalkley,
40
P. Charlton,
9
S. Chatterji,
15
S. Chelkowski,
2
Y. Chen,
1
F. Chiadini,
45
D. Chin,
42
E. Chin,
49
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,
16
T. Corbitt,
18
D. Coward,
49
D. Coyne,
15
J. D. E. Creighton,
50
T. D. Creighton,
15
R. P. Croce,
46
D. R. M. Crooks,
40
A. M. Cruise,
38
A. Cumming,
40
J. Dalrymple,
32
E. D’Ambrosio,
15
K. Danzmann,
14,2
G. Davies,
7
D. DeBra,
31
J. Degallaix,
49
M. Degree,
31
T. Demma,
46
V. Dergachev,
42
S. Desai,
33
R. DeSalvo,
15
S. Dhurandhar,
13
M. Dı
́
az,
34
J. Dickson,
4
A. Di Credico,
32
G. Diederichs,
14
A. Dietz,
7
E. E. Doomes,
30
R. W. P. Drever,
5
J.-C. Dumas,
49
R. J. Dupuis,
15
J. G. Dwyer,
10
P. Ehrens,
15
E. Espinoza,
15
T. Etzel,
15
M. Evans,
15
T. Evans,
17
S. Fairhurst,
7,15
Y. Fan,
49
D. Fazi,
15
M. M. Fejer,
31
L. S. Finn,
33
V. Fiumara,
45
N. Fotopoulos,
50
A. Franzen,
14
K. Y. Franzen,
39
A. Freise,
38
R. Frey,
43
T. Fricke,
44
P. Fritschel,
18
V. V. Frolov,
17
M. Fyffe,
17
V. Galdi,
46
J. Garofoli,
16
I. Gholami,
1
J. A. Giaime,
17,19
S. Giampanis,
44
K. D. Giardina,
17
K. Goda,
18
E. Goetz,
42
L. M. Goggin,
15
G. Gonza
́
lez,
19
S. Gossler,
4
A. Grant,
40
S. Gras,
49
C. Gray,
16
M. Gray,
4
J. Greenhalgh,
27
A. M. Gretarsson,
11
R. Grosso,
34
H. Grote,
2
S. Grunewald,
1
M. Guenther,
16
R. Gustafson,
42
B. Hage,
14
D. Hammer,
50
C. Hanna,
19
J. Hanson,
17
J. Harms,
2
G. Harry,
18
E. Harstad,
43
T. Hayler,
27
J. Heefner,
15
I. S. Heng,
40
A. Heptonstall,
40
M. Heurs,
2
M. Hewitson,
2
S. Hild,
14
E. Hirose,
32
D. Hoak,
17
D. Hosken,
37
J. Hough,
40
E. Howell,
49
D. Hoyland,
38
S. H. Huttner,
40
D. Ingram,
16
E. Innerhofer,
18
M. Ito,
43
Y. Itoh,
50
A. Ivanov,
15
D. Jackrel,
31
B. Johnson,
16
W. W. Johnson,
19
D. I. Jones,
47
G. Jones,
7
R. Jones,
40
L. Ju,
49
P. Kalmus,
10
V. Kalogera,
25
D. Kasprzyk,
38
E. Katsavounidis,
18
K. Kawabe,
16
S. Kawamura,
24
F. Kawazoe,
24
W. Kells,
15
D. G. Keppel,
15
F. Ya. Khalili,
22
C. Kim,
25
P. King,
15
J. S. Kissel,
19
S. Klimenko,
39
K. Kokeyama,
24
V. Kondrashov,
15
R. K. Kopparapu,
19
D. Kozak,
15
B. Krishnan,
1
P. Kwee,
14
P. K. Lam,
4
M. Landry,
16
B. Lantz,
31
A. Lazzarini,
15
B. Lee,
49
M. Lei,
15
J. Leiner,
51
V. Leonhardt,
24
I. Leonor,
43
K. Libbrecht,
15
P. Lindquist,
15
N. A. Lockerbie,
48
M. Longo,
45
M. Lormand,
17
M. Lubinski,
16
H. Lu
̈
ck,
14,2
B. Machenschalk,
1
M. MacInnis,
18
M. Mageswaran,
15
K. Mailand,
15
M. Malec,
14
V. Mandic,
15
S. Marano,
45
S. Ma
́
rka,
10
J. Markowitz,
18
E. Maros,
15
I. Martin,
40
J. N. Marx,
15
K. Mason,
18
L. Matone,
10
V. Matta,
45
N. Mavalvala,
18
R. McCarthy,
16
D. E. McClelland,
4
S. C. McGuire,
30
M. McHugh,
21
K. McKenzie,
4
J. W. C. McNabb,
33
S. McWilliams,
23
T. Meier,
14
A. Melissinos,
44
G. Mendell,
16
R. A. Mercer,
39
S. Meshkov,
15
C. J. Messenger,
40
D. Meyers,
15
E. Mikhailov,
18
S. Mitra,
13
V. P. Mitrofanov,
22
G. Mitselmakher,
39
R. Mittleman,
18
O. Miyakawa,
15
S. Mohanty,
34
G. Moreno,
16
K. Mossavi,
2
C. MowLowry,
4
A. Moylan,
4
D. Mudge,
37
G. Mueller,
39
S. Mukherjee,
34
H. Mu
̈
ller-Ebhardt,
2
J. Munch,
37
P. Murray,
40
E. Myers,
16
J. Myers,
16
T. Nash,
15
G. Newton,
40
A. Nishizawa,
24
K. Numata,
23
B. O’Reilly,
17
R. O’Shaughnessy,
25
D. J. Ottaway,
18
H. Overmier,
17
B. J. Owen,
33
Y. Pan,
41
M. A. Papa,
1,50
V. Parameshwaraiah,
16
P. Patel,
15
M. Pedraza,
15
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,
16
D. Rabeling,
4
H. Radkins,
16
R. Rahkola,
43
N. Rainer,
2
M. Rakhmanov,
33
M. Ramsunder,
33
K. Rawlins,
18
S. Ray-Majumder,
50
V. Re,
38
H. Rehbein,
2
S. Reid,
40
D. H. Reitze,
39
L. Ribichini,
2
R. Riesen,
17
K. Riles,
42
B. Rivera,
16
N. A. Robertson,
15,40
C. Robinson,
7
E. L. Robinson,
38
S. Roddy,
17
A. Rodriguez,
19
A. M. Rogan,
51
J. Rollins,
10
J. D. Romano,
7
J. Romie,
17
R. Route,
31
S. Rowan,
40
A. Ru
̈
diger,
2
L. Ruet,
18
P. Russell,
15
K. Ryan,
16
S. Sakata,
24
M. Samidi,
15
L. Sancho de la Jordana,
36
V. Sandberg,
16
V. Sannibale,
15
S. Saraf,
26
P. Sarin,
18
B. S. Sathyaprakash,
7
S. Sato,
24
P. R. Saulson,
32
R. Savage,
16
P. Savov,
6
S. Schediwy,
49
R. Schilling,
2
R. Schnabel,
2
R. Schofield,
43
B. F. Schutz,
1,7
P. Schwinberg,
16
S. M. Scott,
4
A. C. Searle,
4
B. Sears,
15
F. Seifert,
2
D. Sellers,
17
A. S. Sengupta,
7
P. Shawhan,
41
D. H. Shoemaker,
18
A. Sibley,
17
X. Siemens,
15,6
D. Sigg,
16
S. Sinha,
31
A. M. Sintes,
36,1
B. J. J. Slagmolen,
4
J. Slutsky,
19
J. R. Smith,
2
M. R. Smith,
15
K. Somiya,
2,1
K. A. Strain,
40
D. M. Strom,
43
A. Stuver,
33
T. Z. Summerscales,
3
K.-X. Sun,
31
M. Sung,
19
P. J. Sutton,
15
H. Takahashi,
1
D. B. Tanner,
39
M. Tarallo,
15
R. Taylor,
15
R. Taylor,
40
J. Thacker,
17
K. A. Thorne,
33
K. S. Thorne,
6
A. Thu
̈
ring,
14
K. V. Tokmakov,
40
C. Torres,
34
C. Torrie,
40
G. Traylor,
17
M. Trias,
36
W. Tyler,
15
D. Ugolini,
35
C. Ungarelli,
38
K. Urbanek,
31
H. Vahlbruch,
14
M. Vallisneri,
6
C. Van Den Broeck,
7
M. Varvella,
15
S. Vass,
15
A. Vecchio,
38
J. Veitch,
40
PHYSICAL REVIEW D
77,
022001 (2008)
1550-7998
=
2008
=
77(2)
=
022001(38)
022001-1
©
2008 The American Physical Society
P. Veitch,
37
A. Villar,
15
C. Vorvick,
16
S. P. Vyachanin,
22
S. J. Waldman,
15
L. Wallace,
15
H. Ward,
40
R. Ward,
15
K. Watts,
17
D. Webber,
15
A. Weidner,
2
M. Weinert,
2
A. Weinstein,
15
R. Weiss,
18
S. Wen,
19
K. Wette,
4
J. T. Whelan,
1
D. M. Whitbeck,
33
S. E. Whitcomb,
15
B. F. Whiting,
39
C. Wilkinson,
16
P. A. Willems,
15
L. Williams,
39
B. Willke,
14,2
I. Wilmut,
27
W. Winkler,
2
C. C. Wipf,
18
S. Wise,
39
A. G. Wiseman,
50
G. Woan,
40
D. Woods,
50
R. Wooley,
17
J. Worden,
16
W. Wu ,
39
I. Yakushin,
17
H. Yamamoto,
15
Z. Yan,
49
S. Yoshida,
29
N. Yunes,
33
M. Zanolin,
18
J. Zhang,
42
L. Zhang,
15
C. Zhao,
49
N. Zotov,
20
M. Zucker,
18
H. zur Mu
̈
hlen,
14
and J. Zweizig
15
(LIGO Scientific Collaboration)
1
Albert-Einstein-Institut, Max-Planck-Institut fu
̈
r Gravitationsphysik, D-14476 Golm, Germany
2
Albert-Einstein-Institut, Max-Planck-Institut fu
̈
r Gravitationsphysik, D-30167 Hannover, Germany
3
Andrews University, Berrien Springs, Michigan 49104 USA
4
Australian National University, Canberra, 0200, Australia
5
California Institute of Technology, Pasadena, California 91125, USA
6
Caltech-CaRT, Pasadena, California 91125, USA
7
Cardiff University, Cardiff, CF24 3AA, United Kingdom
8
Carleton College, Northfield, Minnesota 55057, USA
9
Charles Sturt University, Wagga Wagga, NSW 2678, Australia
10
Columbia University, New York, New York 10027, USA
11
Embry-Riddle Aeronautical University, Prescott, Arizona 86301 USA
12
Hobart and William Smith Colleges, Geneva, New York 14456, USA
13
Inter-University Centre for Astronomy and Astrophysics, Pune - 411007, India
14
Leibnitz Universita
̈
t Hannover, D-30167 Hannover, Germany
15
LIGO - California Institute of Technology, Pasadena, California 91125, USA
16
LIGO Hanford Observatory, Richland, Washington 99352, USA
17
LIGO Livingston Observatory, Livingston, Louisiana 70754, USA
18
LIGO - Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, USA
19
Louisiana State University, Baton Rouge, Louisiana 70803, USA
20
Louisiana Tech University, Ruston, Louisiana 71272, USA
21
Loyola University, New Orleans, Louisiana 70118, USA
22
Moscow State University, Moscow, 119992, Russia
23
NASA/Goddard Space Flight Center, Greenbelt, Maryland 20771, USA
24
National Astronomical Observatory of Japan, Tokyo 181-8588, Japan
25
Northwestern University, Evanston, Illinois 60208, USA
26
Rochester Institute of Technology, Rochester, New York 14623, USA
27
Rutherford Appleton Laboratory, Chilton, Didcot, Oxon OX11 0QX United Kingdom
28
San Jose State University, San Jose, California 95192, USA
29
Southeastern Louisiana University, Hammond, Louisiana 70402, USA
30
Southern University and A&M College, Baton Rouge, Louisiana 70813, USA
31
Stanford University, Stanford, California 94305, USA
32
Syracuse University, Syracuse, New York 13244, USA
33
The Pennsylvania State University, University Park, Pennsylvania 16802, USA
34
The University of Texas at Brownsville and Texas Southmost College, Brownsville, Texas 78520, USA
35
Trinity University, San Antonio, Texas 78212, USA
36
Universitat de les Illes Balears, E-07122 Palma de Mallorca, Spain
37
University of Adelaide, Adelaide, SA 5005, Australia
38
University of Birmingham, Birmingham, B15 2TT, United Kingdom
39
University of Florida, Gainesville, Florida 32611, USA
40
University of Glasgow, Glasgow, G12 8QQ, United Kingdom
41
University of Maryland, College Park, Maryland 20742 USA
42
University of Michigan, Ann Arbor, Michigan 48109, USA
43
University of Oregon, Eugene, Oregon 97403, USA
44
University of Rochester, Rochester, New York 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 Western Australia, Crawley, WA 6009, Australia
50
University of Wisconsin-Milwaukee, Milwaukee, Wisconsin 53201, USA
B. ABBOTT
et al.
PHYSICAL REVIEW D
77,
022001 (2008)
022001-2
51
Washington State University, Pullman, Washington 99164, USA
(Received 11 September 2007; published 10 January 2008; publisher error corrected 29 February 2008)
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
semicoherent 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, 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.
DOI:
10.1103/PhysRevD.77.022001
PACS numbers: 04.80.Nn, 07.05.Kf, 95.55.Ym, 97.60.Gb
I. INTRODUCTION
We report on a search with the LIGO (Laser
Interferometer 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
Scientific Collaboration (LSC) has previously reported on
searches for periodic gravitational radiation, using a long-
period coherent method to target known pulsars [
3
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 coherent
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 semicoherent
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 order 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 expressed 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 inclinations
of the spin axes.
The second purpose of this paper, along with the pre-
vious coherent [
6
] and semicoherent [
7
] papers, is to lay
the foundation for the methods that will be used in future
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 semicoherent stages [
10
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 spectral
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
15
], averages nor-
malized power from each SFT. In the Hough method
reported previously [
7
,
10
], referred to here as ‘‘standard
Hough,’’ the sum is of binary zeroes or ones, where a SFT
contributes unity if the power exceeds a normalized power
threshold. In this paper a ‘‘weighted Hough’’ scheme,
ALL-SKY SEARCH FOR PERIODIC GRAVITATIONAL
...
PHYSICAL REVIEW D
77,
022001 (2008)
022001-3
henceforth also referred to as ‘‘Hough,’’ has been devel-
oped and is similar to that described in Ref. [
16
]. This
scheme also allows 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 detector antenna pattern to maximize the
signal-to-noise 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 frequency ranges
to have better detection efficiency than the StackSlide and
Hough methods, the exceptions occurring in bands with
large nonstationary 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 statistic (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 mean that each
could play a role in our future 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 introduces
methods that are crucial to the development of our future
searches.
This paper is organized as follows: Sec. II briefly de-
scribes the LIGO interferometers, focusing on improve-
ments made for the S4 data run, and discusses the
sensitivity 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 Sec. 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
Sec. 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 interfer-
ometer in Livingston, Louisiana (called L1) and two inter-
ferometers in Hanford, Washington, one 4-km and another
2-km (H1 and H2, respectively).
The data analyzed in this paper were produced during
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 displacement
spectral amplitudes near
2
:
5

10

19
mHz

1
=
2
in their
most sensitive frequency band near 150 Hz. In units of
gravitational-wave strain amplitude, the sensitivity of H2 is
roughly a factor of 2 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 2 of the design goals.
Figure
1
shows representative strain spectral noise den-
sities for the three interferometers during the run. As
discussed in Sec. V below, however, nonstationarity of
the noise was significant.
Changes to the interferometers before the S4 run in-
cluded the following improvements [
19
]:
(i) Installation of active seismic isolation of support
structures at Livingston to cope with high anthropo-
genic ground motion in the 1– 3 Hz band.
(ii) 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.
(iii) Replacement of a synthesized radio frequency os-
cillator for phase modulation with a crystal oscil-
lator before S4 began (H1) and midway 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.)
(iv) Lower-noise mirror-actuation electronics (H1, H2,
and L1).
(v) Higher-bandwidth laser frequency stabilization (H1,
H2, and L1) and intensity stabilization (H1 and L1).
(vi) Installation of radiation pressure actuation of mir-
rors for calibration validation (H1).
(vii) Commissioning of complete alignment control
system for the L1 interferometer (already imple-
mented for H1 and H2 in S3 run).
100
1000
10
−23
10
−22
10
−21
10
−20
10
−19
10
−18
LIGO Detector Sensitivities During S4 Science Run
Frequency (Hz)
Strain noise amplitude spectral density (Hz
−1/2
)
H2
L1
H1
LIGO SRD goal (4 km)
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.
B. ABBOTT
et al.
PHYSICAL REVIEW D
77,
022001 (2008)
022001-4
(viii) Refurbishment of lasers and installation of photo-
diodes and electronics to permit interferometer
operation with increased laser power (H1, H2,
and L1).
(ix) Mitigation of electromagnetic interference (H1,
H2, and L1) and acoustic interference (L1).
The data were acquired and digitized at a rate of
16 384 Hz. Data acquisition was periodically interrupted
by disturbances 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 polarizations
defined as ‘‘

’’ with waveform
h


t

and ‘‘

’’ with
waveform
h


t

. The calibrated response seen by an inter-
ferometric 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 polarization
angle of the wave, and
F

;

are the detector antenna
pattern functions for the two orthogonal polarizations.
For periodic (nearly pure sinusoidal) gravitational 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 polar-
izations, 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 irrele-
vant constant).
For an isolated quadrupolar gravitational-wave emitter,
characterized by a rotating triaxial-ellipsoid mass distribu-
tion, the amplitudes
A

and
A

are related to the inclina-
tion angle of the source,

, and the wave amplitude,
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 fre-
quency,

, 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 moment
of inertia with respect to its spin axis, and

is the equato-
rial ellipticity of the star [
18
]. Assuming that all of the
frequency’s derivative,
_
f
, is due to emission of gravita-
tional radiation and that
I
takes the canonical value
10
38
kg m
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 derivation.)
Note that the methods used in this paper are sensitive to
periodic signals from any type of isolated gravitational-
wave source (e.g., freely precessing or oscillating neutron
stars as well as triaxial ones), though we present upper
limits in terms of
h
0
and

. Because we use semicoherent
methods, only the instantaneous signal frequency in the
detector reference frame,
2
f

t

d


t

=dt
, needs to be
calculated. In the detector reference frame this can, to a
very good approximation, be related to the instantaneous
SSB-frame frequency
^
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 semicoherent 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.)
ALL-SKY SEARCH FOR PERIODIC GRAVITATIONAL
...
PHYSICAL REVIEW D
77,
022001 (2008)
022001-5
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 differ-
ences, 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 dependent Doppler
modulations of the apparent source frequency due to the
Earth’s rotation and its orbital motion around the SSB, and
the frequency’s time derivative, intrinsic to the source (see
Fig.
2
). This requires a search in a four-dimensional pa-
rameter space; a template in the space refers to a set of
values:

f
^
f
0
;
_
f;;
g
. The third method, PowerFlux,
also searches explicitly over polarization angle, so that

f
^
f
0
;
_
f;;;
g
.
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 a 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 interferome-
ters, likely detectable source frequencies [
21
], and limita-
tions 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
j
_
f
j
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 signals. All known isolated pulsars spin down more
slowly than the two values of
j
_
f
j
max
used here, and as seen
in the results section, the ellipticity required for higher
j
_
f
j
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 ex-
hibit slight spin-up, believed to arise from acceleration in
the Earth’s direction; such spin-up values have magnitudes
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 time scale
f=
j
4
_
f
j
(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 time scale 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
j
_
f
j
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 thousand at
50 Hz to about
2

10
6
at 1000 Hz. All three methods use
nearly isotropic grids which cover the entire sky. The
PowerFlux search also divides the sky into regions accord-
ing to susceptibility to stationary instrumental line arti-
facts. 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 pulsar
sources distributed uniformly over the sky and with iso-
tropically distributed spin axes. PowerFlux sets strict fre-
quentist
limits
on
circular
and
linear-polarization
amplitudes
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, regard-
less 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 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-
FIG. 2 (color online).
An illustration of the discrete frequency
bins of the short Fourier transform (SFTs) of the data are shown
vertically, with the discrete start times of the SFTs shown
horizontally. The dark pixels represent a signal in the data. Its
frequency changes with time due to Doppler shifts and intrinsic
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 semicoherent search methods presented here,
though the actual implementations differ significantly.
B. ABBOTT
et al.
PHYSICAL REVIEW D
77,
022001 (2008)
022001-6
polarization limits
h
Circ
-
limit
0
apply only to the most favor-
able inclinations (

0
,

), regardless of sky location and
regardless of frequency and spin-down, as above.
Because of 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 neighbor-
ing 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 number of
SFTs. However, population-averaged upper limits are de-
termined self-consistently to include loss of detection ef-
ficiency 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 instrumental lines
differently. Single-bin lines are flagged during data prepa-
ration so that when searching for a particular 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 40 Hz, 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
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 reference
functions described in [
27
]. SFTs are generated directly
from the calibrated data stream, using 30-minute intervals
of data for which the interferometer is operating in what is
known as science-mode. The choice of 30 minutes is a
trade-off between intrinsic sensitivity, which increases
with SFT length, and robustness against frequency drift
during the SFT interval due to the Earth’s motion, source
spin-down, and nonstationarity of the data [
7
]. The require-
ment that each SFT contain contiguous data at nominal
sensitivity introduces duty factor loss from edge effects,
especially for the Livingston interferometer (
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 interfer-
ometers with the best broadband sensitivity. For
PowerFlux, the corresponding numbers of overlapped
SFTs were 1925 and 1628. The Hough search also used
1063 H2 SFTs. In each case, modest requirements were
placed on data quality to avoid short periods with known
electronic saturations, unmonitored calibration 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 con-
structed. 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
X
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
j
~
x
i
k
j
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 detec-
tor antenna pattern and frequency of the signal are found at
the midpoint time of the data used to generate each SFT.
Frequency dependent quantities are then evaluated at a
frequency index
k
corresponding to the bin nearest this
ALL-SKY SEARCH FOR PERIODIC GRAVITATIONAL
...
PHYSICAL REVIEW D
77,
022001 (2008)
022001-7
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 defi-
nitions 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
X
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 deviation of
1
=

N
p
. 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 normalized
power, the final statistic used in this paper is a weighted
sum of the binary counts, giving the ‘‘Hough number
count’’:
n

X
N

1
i

0
w
i
n
i
;
(17)
where the Hough weights are defined as
w
i
/
1
S
i
f
F
i


2

F
i


2
g
;
(18)
and the weight normalization is chosen according to
X
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 corresponding to this SFT,
the beam-pattern functions are larger for a particular point
in the sky. Note that the sensitivity of the search is gov-
erned by the ratios of the different weights, not by the
choice of overall scale. In the next section 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 argument 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
resulting
PowerFlux detection statistic is [
17
],
R

2
T
coh
P
N

1
i

0
W
i
P
i
=

F
i

2
P
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 using
four linear-polarization projections and one circular-
polarization projection. For the linear-polarization projec-
tions, 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 inde-
pendent 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 implemen-
tation 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
TABLE I. Summary of notation used.
Quantity
Description
P
i
Power for SFT
i
and template

i
Normalized power for SFT
i
and template

n
i
Binary count for SFT
i
and template

S
i
Power spect. noise density for SFT
i
and template

F
i

F

at midpoint of SFT
i
for template

F
i

F

at midpoint of SFT
i
for template

B. ABBOTT
et al.
PHYSICAL REVIEW D
77,
022001 (2008)
022001-8
spectral density of the noise, for every frequency bin of
every SFT. PowerFlux uses a different noise decomposition
method described in its implementation section below.
Note that for Gaussian noise, the single-sided power
spectral density can be estimated using
S
i
k

2
hj
~
x
i
k
j
2
i
T
coh
;
(23)
where the angle brackets represent an ensemble average.
The estimation of
S
i
k
must guard against any biases intro-
duced by the presence of a possible signal and also against
narrow spectral disturbances. For this reason the mean,
hj
~
x
i
k
j
2
i
, is estimated via the median. We assume that the
noise is stationary within a single SFT, but allow for non-
stationarities across different SFTs. In every SFT we cal-
culate the ‘‘running median’’ of
j
~
x
i
k
j
2
for every 101
frequency bins centered on the
k
th
bin, and then estimate
hj
~
x
i
k
j
2
i
[
23
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 de-
nominator of Eq. (
14
) these terms are summed in Eq. (
16
),
while the Hough search applies a cutoff to obtain binary
counts in Eq. (
15
) before summing. This results 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.691 162 (which was chosen to normalize the data so
that the mean value of the StackSlide Power equals one)
compared with the expected ratio for an exponential dis-
tribution of 0.698 073 used in the Hough search (which is
explained in Appendix A of [
7
]). It is important 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
15
].
Brady and Creighton [
12
] first described this approach 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 demodulated 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 averaging the
maximum likelihood statistic (known as the
F
-statistic)
which does include the beam-pattern response is men-
tioned in Ref. [
15
] (see also [
6
,
10
,
18
]), and further exten-
sions 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
p
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
p
, even for non-Gaussian noise.
The StackSlide code, which implements the method
described above, is part of the C-based LSC Algorithms
Library 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 [
28
],
and the final post processing steps are performed using
Matlab [
29
]. 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
FIG. 3.
Flow chart for the pipeline used to find the upper limits
presented in this paper using the StackSlide method.
ALL-SKY SEARCH FOR PERIODIC GRAVITATIONAL
...
PHYSICAL REVIEW D
77,
022001 (2008)
022001-9