of 21
Search for gravitational waves from Scorpius X-1 in the second Advanced
LIGO observing run with an improved hidden Markov model
B. P. Abbott
etal.
*
(LIGO Scientific Collaboration and Virgo Collaboration)
(Received 28 June 2019; published 4 December 2019)
We present results from a semicoherent search for continuous gravitational waves from the low-mass x-ray
binary Scorpius X-1, using a hidden Markov model (HMM) to track spin wandering. This search improves
on previous HMM-based searches of LIGO data by using an improved frequency domain matched
filter, the
J
-statistic, and by analyzing data from Advanced LIGO
s second observing run. In the frequency
range searched, from 60 to 650 Hz, we find no evidence of gravitational radiation. At 194.6 Hz, the most
sensitive search frequency, we report an upper limit on gravitational wave strain (at 95% confidence) of
h
95%
0
¼
3
.
47
×
10
25
when marginalizing over source inclination angle. This is the most sensitive search for
Scorpius X-1, to date, that is specifically designed to be robust in the presence of spin wandering.
DOI:
10.1103/PhysRevD.100.122002
I. INTRODUCTION
Rotating neutron stars with nonaxisymmetric deforma-
tions are predicted to emit persistent, periodic gravitational
radiation. They are a key target for continuous-wave searches
performed with gravitational wave (GW) detectors such as
the second-generation Advanced Laser Interferometer
Gravitational-Wave Observatory (Advanced LIGO)
[1
5]
and Virgo
[4]
. The time-varying quadrupole moment neces-
sary for GW emission may result from thermal
[6,7]
,or
magnetic
[8
10]
gradients,
r
-modes
[11
13]
, or nonaxisym-
metric circulation of the superfluid interior
[14
17]
.These
mechanisms produce signals at certain multiples of the spin
frequency
f
[1]
. Of particular interest are accreting
low-mass x-ray binaries (LMXB), such as Scorpius X-1
(Sco X-1), where a neutron star is spun up by accretion from
its stellar companion. Electromagnetic observations of
LMXBs to date imply
f
620
Hz
[18]
, well short of the
theoretical centrifugal break-up limit
f
1
.
5
kHz
[19]
.
Regardless of the exact GW mechanism, the latter observa-
tion suggests an equilibrium between the spin-up accretion
torque, and GW spin-down torque
[20
22]
. Torque balance
also implies a relation between x-ray luminosity and the GW
strain, making Sco X-1, the brightest LMXB x-ray source,
the most promising known target.
Initial LIGO, a first-generation detector, started taking
science data in 2002. It reached its design sensitivity in
Science Run 5 (S5) starting 2005
[23]
, and exceeded it in
Science Run 6 (S6)
[24]
. Following detector upgrades,
the second-generation Advanced LIGO interferometer
[2]
began taking science data during Observing Run 1 (O1),
which ran from September 2015 to January 2016. The
strain noise in O1 is 3 to 4 times lower than S6 between
100 and 300 Hz
[25]
. During this period, LIGO observed
three binary black hole mergers, GW150914
[26]
,
GW151012, and GW151226
[27]
. Observing Run 2
(O2) began in November 2016, and ran until August 26,
2017. From August 1, 2017, the two LIGO detectors were
joined by Virgo, resulting in a three-detector network. As
well as further binary black hole mergers
[28]
, LIGO and
Virgo made the first gravitational wave observation of a
binary neutron-star merger during O2
[29]
.
No search has yet reported a detection of a continuous-
wave source. To date, four searches for Sco X-1 have been
conducted on Initial LIGO data, and three on Advanced
LIGO data. The first search coherently analyzed the most-
sensitive six hour segment from Science Run 2 (S2) using
the
F
-statistic
[30]
, a maximum-likelihood detection sta-
tistic
[31]
. The second was a directed, semicoherent
analysis using the
C
-statistic
[32]
. The third, also a directed
analysis, used the TwoSpect algorithm on doubly Fourier
transformed S5 data
[33
35]
. The fourth applied the
radiometer algorithm
[36]
to conduct a directed search
on S4
[37]
,S5
[38]
, and later O1
[39]
data. Three LMXB
searches have been performed with Advanced LIGO data,
comprising the radiometer search
[39]
, an analysis based
on a hidden Markov model (HMM)
[40]
, and a cross-
correlation analysis
[41
43]
. The upper limits established
by these searches are summarized in Table
I
.
Astrophysical modeling and x-ray observations suggest
that the spin frequency of a LMXB wanders stochastically
in response to fluctuations in the hydromagnetic accretion
torque
[46
49]
. As no electromagnetic measurements of
f
are available to guide a gravitational wave search for
Sco X-1, such searches must either account for spin
wandering or limit their observing times and/or coherence
*
Full author list given at the end of the article.
PHYSICAL REVIEW D
100,
122002 (2019)
2470-0010
=
2019
=
100(12)
=
122002(21)
122002-1
© 2019 American Physical Society
times in accordance with the anticipated timescale and
amplitude of the spin wandering
[50]
. For example, the
sideband search described in Ref.
[32]
is restricted to data
segments no longer than ten days. The HMM tracker, first
applied to the search for Sco X-1 in Ref.
[40]
, is an effective
technique for detecting the most probable underlying spin
frequency,
f
ð
t
Þ
and thus accounting for spin wandering.
The signal from a binary source is Doppler shifted, as the
neutron star revolves around the barycenter of the binary,
dispersing power into orbital sidebands near the source-
frame emission frequency. The separation of these side-
bands and the source-frame frequency depends on the
binary orbital parameters and
f
, but is typically within
0.05% of the gravitational wave frequency for a source
such as Sco X-1. Four maximum-likelihood matched
filters have been developed to detect these sidebands:
the
C
-statistic, which weights sidebands equally
[32]
, the
binary modulated
F
-statistic
[51]
, the Bessel-weighted
F
-
statistic
[52]
, and the
J
-statistic, which extends the Bessel-
weighted
F
-statistic to account for the phase of the binary
orbit
[53]
. Any of these matched filters can be combined
with the HMM to conduct a search for signals from a binary
source that accounts for spin wandering.
In this paper, we combine the
J
-statistic described in
Ref.
[53]
with the HMM described in that paper and
Refs. [40, 50], and perform a directed search of Advanced
LIGO O2 data for evidence of a gravitational wave signal
from Sco X-1. In the search band 60
650 Hz, we find no
evidence of a gravitational wave signal. The paper is
organized as follows. In Sec.
II
, we briefly review the
HMM and the
J
-statistic. In Sec.
III
, we discuss the search
strategy and parameter space. In Sec.
IV
, we report on the
results from the search and veto candidates corresponding
to instrumental artifacts. In Sec.
V
, we discuss the search
sensitivity and consequent upper limits on the gravitational
wave strain.
II. SEARCH ALGORITHM
In this section, we outline the two key components of
the search algorithm: the HMM used to recover the most
probable spin history
f
0
ð
t
Þ
, and the
J
-statistic, the matched
filter that accounts for the Doppler shifts introduced by the
orbital motions of Earth and the LMXB. The HMM
formalism is the same as used in Refs.
[40,52,53]
,sowe
review it only briefly. The
J
-statistic is described fully in
Ref.
[53]
; again, we review it briefly.
A. HMM formalism
A Markov model describes a stochastic process in terms
of a state variable
q
ð
t
Þ
, which transitions between allowable
states
f
q
1
;
...
;q
N
Q
g
at discrete times
f
t
0
;
...
;t
N
T
g
. The
transition matrix
A
q
j
q
i
represents the probability of jumping
from state
q
i
at the time
t
¼
t
n
to
q
j
at
t
¼
t
n
þ
1
depending
only on
q
ð
t
n
Þ
. A HMM extends the Markov model to
situations where direct observation of
q
ð
t
Þ
is impossible
[
q
ð
t
Þ
is called the hidden state]. Instead one measures an
observable state
o
ð
t
Þ
selected from
f
o
1
;
...
;o
N
o
g
, which is
related to the hidden state by the emission matrix
L
o
j
q
i
,
which gives the likelihood that the system is in state
q
i
given the observation
o
j
. In gravitational wave searches
for LMXBs like Sco X-1, where the spin frequency cannot
be measured electromagnetically, it is natural to map
q
ð
t
Þ
to
f
0
ð
t
Þ
and
o
ð
t
Þ
to the raw interferometer data, some
equivalent intermediate data product (e.g., short Fourier
transforms), or a detection statistic (e.g.,
F
-statistic,
J
-statistic).
In a LMXB search, we divide the total observation
(duration
T
obs
) into
N
T
equal segments of length
T
drift
¼
T
obs
=N
T
. In practice,
T
drift
is chosen on astrophysi-
cal grounds to give
N
T
¼d
T
obs
=T
drift
e
based on an esti-
mates of plausible spin-wandering timescales
[50]
; in this
paper we follow Ref.
[40]
in choosing
T
drift
¼
10
d. The
tracker is able to track the signal even if the spin frequency
occasionally jumps by two bins as it can catch up to the
signal path, although with an attendant loss of sensitivity as
the recovered must include a step that contains only noise.
In each segment, the emission probability
L
o
j
q
i
is
computed from some frequency domain estimator
G
ð
f
Þ
such as the maximum likelihood
F
-or
J
-statistic (dis-
cussed in Sec.
II B
). The frequency resolution of the
TABLE I. Summary of indicative upper limits achieved in previous searches for Sco X-1. VSR2 and VSR3 are Virgo Science Runs 2
and 3, respectively. Where applicable, the upper limits refer to signals of unknown polarization.
Search
Data
Upper limit
Reference
F
-statistic
S2
h
95%
0
2
×
10
22
at 464
484 Hz, 604
626 Hz
[31]
C
-statistic
S5
h
95%
0
8
×
10
25
at 150 Hz
[32]
TwoSpect
S6, VSR2, VSR3
h
95%
0
2
×
10
23
at 20
57.25 Hz
[34]
Radiometer
S4, S5
h
90%
0
2
×
10
24
at 150 Hz
[38, 44]
TwoSpect
S6
h
95%
0
1
.
8
×
10
24
at 165 Hz
[45]
Radiometer
O1
h
90%
0
6
.
7
×
10
25
at 130
175 Hz
[39]
Viterbi 1.0
O1
h
95%
0
8
.
3
×
10
25
at 106 Hz
[40]
Cross-correlation
O1
h
95%
0
2
.
3
×
10
25
at 175 Hz
[43]
B. P. ABBOTT
et al.
PHYS. REV. D
100,
122002 (2019)
122002-2
estimator is
Δ
f
drift
¼
1
=
ð
2
T
drift
Þ
. The probability that an
observation
O
¼f
o
ð
t
1
Þ
;
...
;o
ð
t
N
T
Þg
is associated with a
particular hidden path
Q
¼f
q
ð
t
0
Þ
;
...
;q
ð
t
N
T
Þg
is then
given by
P
ð
Q
j
O
Þ
L
o
ð
t
N
T
Þ
q
ð
t
N
T
Þ
A
q
ð
t
N
T
Þ
q
ð
t
N
T
1
Þ

L
o
ð
t
1
Þ
q
ð
t
1
Þ
×
A
q
ð
t
1
Þ
q
ð
t
0
Þ
Π
q
ð
t
0
Þ
;
ð
1
Þ
where
Π
q
i
is the prior, i.e., the probability the system starts
in state
q
i
at
t
¼
t
0
. For this search, we take a flat prior.
[Note that there is no initial observation
o
ð
t
0
Þ
as the initial
state of the system is captured by the prior.] The task, then,
is to find the optimal hidden path
Q
, that is, the path
Q
that maximizes
P
ð
Q
j
O
Þ
given
O
. We find
Q
efficiently
with the recursive Viterbi algorithm
[54]
, which is dis-
cussed in detail in Appendix A of Ref.
[40]
.
In this paper, we follow the convention in Ref.
[40]
of
defining the Viterbi detection score
S
for a path as the
number of standard deviations by which that path
s log-
likelihood exceeds the mean log-likelihood of all paths.
Mathematically we have
S
¼
ln
δ
q
ð
t
N
T
Þ
μ
ln
δ
ð
t
N
T
Þ
σ
ln
δ
ð
t
N
T
Þ
;
ð
2
Þ
where
μ
ln
δ
ð
t
N
T
Þ
¼
N
1
Q
X
N
Q
i
¼
1
ln
δ
q
i
ð
t
N
T
Þ
;
ð
3
Þ
σ
2
ln
δ
ð
t
N
T
Þ
¼
N
1
Q
X
N
Q
i
¼
1
½
ln
δ
q
i
ð
t
N
T
Þ
μ
ln
δ
ð
t
N
T
Þ

2
;
ð
4
Þ
δ
q
i
ð
t
N
T
Þ
denotes the likelihood of the most likely path
ending in state
q
i
at step
N
T
, and
δ
q
ð
t
N
T
Þ¼
max
i
δ
q
i
ð
t
N
T
Þ
is the likelihood of the optimal path overall.
B.
J
-statistic
The frequency domain estimator
G
ð
f
Þ
converts the
interferometer data into the likelihood that a signal is
present at frequency
f
. For a continuous-wave search for an
isolated neutron star, the maximum-likelihood
F
-statistic
[30]
is a typical choice for
G
ð
f
Þ
. The
F
-statistic accounts
for the diurnal rotation of Earth, and its orbit around the
Solar System barycenter. It is an almost optimal matched
filter for a biaxial rotor
[55]
.
For a neutron star in a binary system, such as a LMXB, the
signal isfrequency (Doppler)modulated bythebinary orbital
motion as well. Reference
[40]
used the Bessel-weighted
F
-statistic to account for this modulation, without using
information about the orbital phase. Reference
[53]
intro-
duced the
J
-statistic, which is a matched filter that extends
the
F
-statistic to include orbital phase in the signal model.
The orbital Doppler effect distributes the
F
-statistic power
into approximately
2
m
þ
1
orbital sidebands separated by
P
1
,with
m
¼d
2
π
f
a
0
e
, where
d
·
e
denotes rounding
up to the nearest integer,
P
is the orbital period, and
a
0
¼ð
a
sin
i
Þ
=c
is the light travel time across the projected
semimajor axis (where
a
is semimajor axis and
i
is the
inclination angle of the binary). For a zero-eccentricity
Keplerian orbit, the Jacobi-Anger identity may be used to
expand the signal
h
ð
t
Þ
in terms of Bessel functions, sug-
gesting a matched filter of the form
[52,53]
G
ð
f
Þ¼
F
ð
f
Þ
B
ð
f
Þ
;
ð
5
Þ
with
B
ð
f
Þ¼
X
m
s
¼
m
J
s
ð
2
π
f
0
a
0
Þ
e
is
φ
a
δ
ð
f
s=P
Þ
;
ð
6
Þ
where
J
s
ð
z
Þ
is the Bessel function of the first kind of order
s
,
φ
a
is the orbital phase at a reference time, and
δ
is the Dirac
delta function.
All else being equal, using the
J
-statistic instead of the
Bessel-weighted
F
-statistic improves sensitivity by a factor
of approximately 4. Reference
[53]
, particularly Sec.
IV
of that paper, examines the difference between the two
estimators in depth.
The Bessel-weighted
F
-statistic requires a search over
a
0
but does not depend on
φ
a
. By contrast, the more-
sensitive
J
-statistic involves searching over
φ
a
too. In this
paper we apply the
J
-statistic to search for Sco X-1.
Details of the search and priors derived from electromag-
netic measurements are discussed in Sec.
III
.
III. LIGO O2 SEARCH
A. Sco X-1 parameters
The matched filter described in Sec.
II B
depends on
three binary orbital parameters: the period
P
, the projected
semimajor axis
a
0
, and the phase
φ
a
. The
F
-statistic
depends on the sky location
α
(right ascension) and
δ
(declination), and optionally the source frequency deriva-
tives. For this search, we assume there is no secular
evolution in frequency. The other parameters have been
measured electromagnetically for Sco X-1 and are pre-
sented in Table
II
.
For
α
,
δ
, and
P
, the uncertainties in the electromagnetic
measurements are small enough that they have no appreci-
able effect on the sensitivity of the search
[51,60,61]
, and a
single, central value can be assumed. However, the uncer-
tainties in
a
0
and
φ
a
cannot be neglected. The time spent
searching orbital parameters scales as the number of
(
a
0
,
φ
a
) pairs. Careful selection of the ranges of
a
0
and
φ
a
is essential to keep computational costs low.
SEARCH FOR GRAVITATIONAL WAVES FROM SCORPIUS X-1
...
PHYS. REV. D
100,
122002 (2019)
122002-3
The previous analysis described in Ref.
[40]
used the
Bessel-weighted
F
-statistic in place of the
J
-statistic, and
searched over a uniformly gridded range of
a
0
, where the
grid resolution did not depend on frequency. However, the
J
-statistic is more sensitive to mismatch in the binary
orbital parameters, so a finer grid is required. We must also
choose an appropriate grid for
φ
a
. (The Bessel-weighted
F
-statistic is independent of
φ
a
.)
As the
J
-statistic has a similar overall response to
parameter mismatches as the binary
F
-statistic, we follow
the formalism in Ref.
[51]
to select an appropriate param-
eter space gridding. We choose a grid which limits the
maximum loss in signal-to-noise ratio (mismatch)
μ
max
to
μ
max
¼
0
.
1
. Equation (71) in Ref.
[51]
gives a general
equation for the number of grid points needed for each
search parameter. For the particular search considered in
this paper, the number of choices for
a
0
and
φ
a
are
N
a
0
¼

π
ffiffiffi
2
p
2
μ
1
=
2
max
f
0
Δ
a
0

;
ð
7
Þ
N
φ
a
¼

1
2
μ
1
=
2
max
f
0
a
0

2
π
P

Δ
φ
a

;
ð
8
Þ
where
Δ
a
0
and
Δ
φ
a
are the widths of the search ranges
for
a
0
and
φ
a
respectively. The number of orbital parameters
to be searched depends on the search frequency. Accordingly
for each search subband, we adopt a different grid resolution,
with the grid refined at higher frequencies. In the subband
beginning at 60 Hz, we have
N
a
0
¼
768
and
N
φ
a
¼
78
;in
the subband beginning at 650 Hz, we have
N
a
0
¼
8227
and
N
φ
a
¼
824
. In principle we could achieve further
computational savings by noting that
N
φ
a
also depends
on
a
0
, but for safety we use the largest
a
0
.
The search range for
a
0
is
1
.
45
a
0
=
ð
1
s
Þ
3
.
25
, which
matches the most recent electromagnetic measurement
[57]
and widens the error bars on the widely cited and previous
best published measurement,
a
0
¼
1
.
44

0
.
18
s
[62]
.
The orbital phase
φ
a
can be related to the electromag-
netically measured time of ascension,
T
asc
, given in
Table
II
,by
φ
a
¼
2
π
T
asc
=P
ð
mod
2
π
Þ
:
ð
9
Þ
The one-sigma uncertainty in the published value for
T
asc
is

50
s
[57,59]
for a time of ascension at GPS time
974416624 s (in November 2010). As O2 took place
significantly after this time, to make a conservative estimate
on appropriate error bars for
T
asc
, we advance
T
asc
by
adding 3135 orbital periods to the time of ascension taken
from Ref.
[57]
. As there is uncertainty associated with the
measured orbital period, this widens the one-sigma uncer-
tainty of
T
asc
to

144
s, which we round up to

150
s.
To cover a significant portion of the measured
T
asc
range
while keeping the search computationally feasible, we
search a two-sigma range around the central
T
asc
, namely,
1164 543 014
T
asc
=
ð
1
s
Þ
1 164 543 614
(expressed for
presentation purposes as the time of the last ascension
before the start of O2).
As there is no electromagnetic measurement of
f
for
Sco X-1, we search the band
60
f
=
ð
1
Hz
Þ
650
,
where LIGO is most sensitive, again adopting a uniform
prior (see Sec.
II A
for a discussion of the HMM prior). The
same band is analyzed in Ref.
[40]
. For computational
convenience, we split the band into blocks of approxi-
mately 0.61 Hz (discussed further in Sec.
III B
).
The final electromagnetically measured parameter is the
polarization angle,
ψ
. Because the
F
-statistic components
of the
J
-statistic are maximized over the polarization
angle, the
J
-statistic is insensitive to
ψ
.
TABLE II. Electromagnetic measurements of the sky position and binary orbital parameters of Sco X-1. The uncertainties represent
one-sigma (68%) confidence intervals, except for
a
0
, for which hard limits are given. The search resolution for
a
0
and
T
asc
is different in
each frequency subband, as discussed in Sec.
III A
. The search range for the time of ascension is the observed time of ascension
propagated forward to the start of O2.
Observed parameter
Symbol
Value
Reference
Right ascension
α
16 h 19 m 55.0850 s
[56]
Declination
δ
15
°
38
0
24
:
9
00
[56]
Orbital period
P
68023
.
86048

0
.
0432
s
[57]
Projected semimajor axis
a
0
[1.45, 3.25] s
[57]
Polarization angle
ψ
234

3
°
[58]
Orbital inclination angle
i
44

6
°
[58]
Time of ascension
T
asc
974416624

50
s
[57,59]
Search parameter
Symbol
Search range
Resolution
Frequency
f
0
60
650 Hz
5
.
787037
×
10
7
Hz
Projected semimajor axis
a
0
1.450
3.250 s
Variable
Time of ascension
T
asc
1164543014
1164543614 s
Variable
B. P. ABBOTT
et al.
PHYS. REV. D
100,
122002 (2019)
122002-4
A summary of the search ranges flowing from the
electromagnetically measured parameters of Sco X-1 is
presented in Table
II
.
B. Workflow
The workflow for the search is displayed as a flowchart
in Fig.
1
.
The data from the detector are provided as short Fourier
transforms (SFTs), each covering
T
SFT
¼
1800
s. We
divide the search into subbands, both to facilitate managing
the volume of data, and to ensure that replacing the search
frequency
f
with the midpoint of the subband,
̄
f
, is a good
approximation in Eq.
(6)
. To achieve best performance
from the fast Fourier transforms used to compute the
convolution in
(6)
, it is desirable to have a power of 2
number of frequency bins in the band, so we set the
subband width to be
Δ
f
band
¼
2
20
Δ
f
drift
¼
0
.
6068148
Hz.
This in turn sets the number of hidden states per subband
per binary orbital parameter to be
N
Q
¼
2
20
.
For each subband, we divide the data into
N
T
blocks,
each with duration
T
drift
¼
10
d. We then compute, from
the SFTs, the
F
-statistic
atoms
[63]
(
F
a
,
F
b
) for each
block using the fixed parameters (
α
,
δ
,
P
) in Table
II
.
The next step is to compute the
J
-statistic for the
(
a
0
,
φ
a
) search grid described in Sec.
III A
. The
F
-statistic
atoms do not depend on the binary orbital parameters so
they are not recomputed when calculating the
J
-statistic.
The code to compute the
J
-statistic is based on the
F
-statistic subroutines contained in the LIGO Scientific
Collaboration Algorithm Library
[64]
.
After computing the
J
-statistic, we use the Viterbi
algorithm to compute the optimal paths through the
HMM trellis, i.e., the set of vectors
Q
. In principle, the
tracking problem is three dimensional (over
f
0
,
a
0
, and
φ
a
),
but
a
0
does not vary significantly over
T
obs
1
yr and
φ
a
varies deterministically, with the phase at time step
n
given
by
φ
a
ð
t
n
Þ¼
φ
a
ð
t
n
1
Þþ
2
π
T
drift
=P
. Thus, it is convenient
to search independently over
f
0
and pairs
ð
a
0
;
φ
a
Þ
. This
allows searches over
ð
a
0
;
φ
a
Þ
pairs to be performed in
parallel.
The result of this procedure is one log-likelihood for the
optimal path through the trellis terminating at every 3-tuple
ð
f
0
;a
0
;
φ
a
Þ
. Equation
(2)
converts these log-likelihoods to
Viterbi scores. As the noise power spectral density (PSD) of
the detector is a function of
f
0
, we compute
μ
and
σ
separately for each band. By contrast, the PSD is not a
function of
a
0
and
φ
a
. Therefore, we can recalculate
μ
and
σ
for every
ð
a
0
;
φ
a
Þ
pair (rather than calculating
μ
and
σ
using every log-likelihood across the entire search), thereby
considerably reducing memory use. This has no significant
impact on the Viterbi scores.
For each subband that produces a best Viterbi score
lower than the detection threshold (chosen in Sec.
III C
), we
compute an upper limit on the gravitational wave strain for
a source in that subband. For Viterbi scores that exceed the
threshold, we apply the veto tests described in Sec.
IVA
.
We claim a detection, if a candidate survives all vetoes.
For performance reasons, the most computationally
intensive parts of the search (computing the
J
-statistic,
FIG. 1. Flowchart of the
J
-statistic search pipeline for each
subband. Note that the
F
-statistic atoms are computed once per
block and per subband, then the
J
-statistic is recalculated for
each
ð
a
0
;
φ
a
Þ
pair. The gray ovals are the start and end of the
algorithm, the green rectangles are procedures, the blue (red)
parallelograms are intermediate (input) data, the yellow diamonds
are decision points, and the gray dashed line represents a loop
repeated once for each choice of parameter. The rectangles with a
dashed boundary were run on graphical processing units, while
those with a solid boundary were run on central processing units.
SEARCH FOR GRAVITATIONAL WAVES FROM SCORPIUS X-1
...
PHYS. REV. D
100,
122002 (2019)
122002-5
and the Viterbi tracking) were run using NVIDIA P100
graphical processing units (GPUs). Other steps were run
using CPU codes on Intel Xeon Gold 6140 CPUs.
C. Threshold and false alarm probability
It remains to determine a detection score threshold
S
th
corresponding to the desired false alarm probability.
Consider the probability density function (PDF)
p
n
ð
S
Þ
of
the Viterbi score in noise. For a given threshold
S
th
and a
fixed search frequency and set of binary orbital parameters,
the probability that the score will exceed this threshold (i.e.,
produce a false alarm) is
α
¼
Z
S
th
d
Sp
n
ð
S
Þ
:
ð
10
Þ
In general, the search covers many frequency bins and
choices of binary parameters. The probability
α
N
of a false
alarm over a search covering
N
parameter choices (number
of frequency bins multiplied by number of binary param-
eter choices) is
α
N
¼
1
ð
1
α
Þ
N
:
ð
11
Þ
This equation assumes that the Viterbi score in noise is an
independent random variable at each point in the parameter
space, which is not necessarily true, as the
J
-statistic
calculated for two points nearby in parameter space is
correlated to some degree. However, for
μ
max
¼
0
.
1
as used
in this search, these correlations do not have a significant
impact
[65]
. In practice, we fix
α
N
and
N
and solve
(10)
and
(11)
for
α
and hence
S
th
.
As the noise-only PDF
p
n
ð
S
Þ
of the Viterbi score is
unknown analytically
[40]
, we resort to Monte Carlo sim-
ulations. We generate
10
2
Gaussian noise realizations in
seven subbands of width
Δ
f
band
, namely those starting at 55,
155, 255, 355, 455, 555, and 650 Hz. The noise is generated
using the standard LIGO tool
lalapps_
M
akefakedata_v4
. These are
the same subbands used in Sec. III C of Ref.
[40]
, and the
one-sided noise PSD
S
h
ð
f
Þ
is set to match the O2 data. We
then perform the search described in Sec.
III B
(including
scanning over
a
0
and
φ
a
).
The results of this search produce an empirical version of
p
n
ð
S
Þ
. Plotting the tail of this distribution on a logarithmic
plot suggests that a fit to a function of the form
e
λ
S
is an
appropriate choice to allow the PDF to be extrapolated in
order to solve
(11)
.
We first analyze each band independently to ensure that
there is no frequency dependence in
p
n
ð
S
Þ
. Table
III
gives
the best-fit
λ
, and the threshold
S
th
obtained, for each band
analyzed in isolation. We find that there is no significant
dependence on the subband searched, nor any identifiable
trend in
λ
or
S
th
. Combining the realizations for all bands
produces
λ
¼
3
.
28
and hence
S
th
¼
13
.
66
for
α
¼
0
.
01
.
The empirical PDF and fitted exponential are shown
in Fig.
2
.
D. Sensitivity
After selecting
S
th
, it remains to determine the lowest (as
a function of frequency) characteristic wave strain,
h
95%
0
,
that can be detected with 95% efficiency (i.e., a 5% false
dismissal rate). To do this, we generate Monte Carlo
realizations of Gaussian noise with Sco X-1
like signals
injected. We determine the proportion of signals recovered
as a function of
h
0
and double-check the false alarm
probability quoted above.
For O2, the most sensitive subband of width
Δ
f
band
¼
0
.
6068148
Hz is the one beginning at 194.6 Hz. Following
a typical procedure used to find upper limits for continuous
gravitational wave searches
[66]
, we generate
10
2
noise
realizations and inject signals using the source parameters
in Table
II
, with
T
obs
¼
230
d (the duration of O2),
T
drift
¼
10
d,
N
T
¼
23
,
ffiffiffiffiffi
S
h
p
¼
7
.
058
×
10
24
Hz
1
=
2
,
and cos
ι
¼
1
. The remaining range-bound parameters,
namely
f
0
inj
,
a
0
inj
,
T
asc
inj
, and
ψ
inj
are chosen from a uniform
FIG. 2. Tail of the PDF of the Viterbi score
S
in noise. The
purple histogram shows the empirical PDF derived from
10
2
realizations of the noise analyzed in the seven 0.61 Hz subbands
starting at 55, 155, 255, 355, 455, 555, and 650 Hz. The green
curve is an exponential fitted to the histogram.
TABLE III. Results of investigating the empirical PDF of the
Viterbi score in seven subbands in Gaussian noise. The second
column is
λ
obtained by fitting the PDF to
e
λ
S
. The third column
is the threshold
S
th
obtained by solving Eq.
(11)
.
Start of band (Hz)
λ
S
th
55
3
.
02
14.12
155
3
.
24
13.63
255
3
.
26
13.58
355
3
.
27
13.61
455
3
.
30
13.62
555
3
.
29
13.66
650
3
.
29
13.63
B. P. ABBOTT
et al.
PHYS. REV. D
100,
122002 (2019)
122002-6
distribution within the range given by their one
σ
error bars.
The source frequency
f
0
inj
is chosen from a uniform
distribution on the interval [194.6 Hz, 194.7 Hz]. For each
realization, the signal is injected with progressively lower
h
0
until it can no longer be detected. We denote by
h
0
;
min
;
i
the lowest
h
0
that can be detected in realization
i
. To obtain
h
95%
0
, we take the 95th highest
h
0
;
min
;
i
. The simulations
return the threshold
h
95%
0
¼
1
.
46
×
10
25
at 194.6 Hz.
In general, the signal-to-noise ratio is strongly
affected by the inclination angle
ι
, not just
h
0
. We follow
Ref.
[59]
and define an effective
h
0
that absorbs the
dependence on
ι
:
h
eff
0
¼
h
0
2
1
=
2
f½ð
1
þ
cos
2
ι
Þ
=
2

2
þ
cos
2
ι
g
1
=
2
;
ð
12
Þ
allowing us to generalize results from the simulations
above, where all injections were done with cos
ι
¼
1
.
Thus, the result obtained above corresponds to circular
polarization. The electromagnetically measured inclination
of Sco X-1
s orbit is
i
44
°

6
°
[58]
. Although it is not
necessarily the case, if we assume that the orbital
inclination equals the inclination angle
ι
of the putative
neutron star
s spin axis, we obtain
h
ι
44
°
;
95%
0
¼
1
.
35
h
eff
;
95%
0
.
The search in Ref.
[40]
found a scaling relation of the
form
h
95%
0
S
1
=
2
h
f
1
=
4
0
to hold for fixed
T
obs
. The
f
1
=
4
0
dependence arises because the latter search added side-
bands incoherently. In the case of the
J
-statistic, which
adds sidebands coherently, we expect the scaling to depend
just on
h
0
, with
h
95%
0
S
1
=
2
h
:
ð
13
Þ
We verify this scaling in Gaussian noise by repeating the
injection procedure described above in frequency bands
beginning at 55, 355, and 650 Hz. The scaling is the final
ingredient needed to produce the blue dashed curve in
Fig.
3
, which shows the expected sensitivity of a search
over the full search band, assuming Gaussian noise, a 100%
duty cycle, and a circularly polarized signal.
There is no simple scaling similar to
(13)
that can be used
to account for the effect of non-Gaussian noise and the
detector duty cycle. Hence we introduce a multiplicative
correction factor
κ
j
for a selection of subbands indexed
FIG. 3. Sensitivity of a search for Sco X-1 across the frequency band searched in this work. The horizontal axis shows the search
frequency. The vertical axis shows the wave strain
h
95%
0
needed for a 95% detection efficiency, on the assumption cos
ι
¼
1
, namely, a
circularly polarized signal. The blue dashed curve is based on simulations of Gaussian noise, while the red solid curve is corrected for
the non-Gaussian statistics of the noise and the interferometer duty cycle, through multiplying by
κ
freq
ð
f
Þ
(see Sec.
III D
). The diamonds
show
h
95%
0
derived through injections into subbands, again assuming a circularly polarized signal (see Sec.
III D
).
SEARCH FOR GRAVITATIONAL WAVES FROM SCORPIUS X-1
...
PHYS. REV. D
100,
122002 (2019)
122002-7
by
j
, following Ref.
[40]
. We determine
κ
j
by doing
10
2
injections (drawing parameters as described above) into the
detector data for the
j
th subband, again using progressively
lower
h
0
until we determine the minimum
h
0
detected.
Then,
κ
j
equals
h
eff
;
95%
0
for injections into real noise,
divided by
h
eff
;
95%
0
for injections into Gaussian noise.
Producing
κ
j
in this way for a random selection of
subbands in the search band suggests that
κ
depends
weakly on frequency, most likely due to the
J
-statistic
not perfectly summing sidebands
[40]
. A linear fit to the
computed
κ
j
values suggests a frequency-dependent cor-
rection factor
κ
freq
ð
f
Þ¼
1
.
944
þ
4
.
60
×
10
4
f=
ð
1
Hz
Þ
:
ð
14
Þ
We use
κ
freq
ð
f
Þ
to adjust the blue dashed curve in Fig.
3
,
producing the red solid curve in that figure, which represents
the expected sensitivity across the full search band, where the
noise is realistic (i.e., not Gaussian). The 50 subbands
sampled are shown on the plot as gray diamonds.
IV. O2 ANALYSIS
We now analyze the data from LIGO
s O2, using the full
dataset from November 30, 2016 to August 26, 2017,
including data from the LIGO Livingston (L1) and Hanford
(H1) observatories. The Virgo interferometer also partici-
pated in the last two months of O2, but we do not use any
Virgo data in this analysis.
There are two notable pauses in data gathering: an end-
of-year break starting on December 22, 2016 lasting for
13 days, and a commissioning break starting on May 7,
2017 lasting for 19 (L1) or 32 (H1) days.
Data stretches shorter than
T
SFT
are discarded, as is a
period of approximately one month where much of the
band was contaminated due to a blinking light in the
power system and a digital camera (used for detector
diagnostics) that was inadvertently left on. A detailed
discussion of Advanced LIGO detector noise can be
found in Ref.
[67]
. Taking all these factors into account,
the overall duty cycle (i.e., proportion of time spent
gathering science-quality data) for O2 was 51.9% (L1)
and 46.2% (H1).
Because of the commissioning break, one ten-day block
has no data. We fill this block with a uniform log-
likelihood, so that the HMM has no preference for
remaining in the same frequency bin, or moving by one
bin, during the break, while still allowing a maximum drift
of
Δ
f
drift
every ten days. An alternative, but equivalent,
approach would be to remove the break entirely, and alter
the transition matrix
A
q
i
q
j
for that step to allow the HMM to
wander up to two frequency bins. The end-of-year break is
also longer than ten days, but it is covered by two blocks.
Both of the blocks that overlap with the end-of-year break
contain data.
We search the same frequency band as Ref.
[40]
, namely
60
650 Hz. The lower limit is set by LIGO
s poor
sensitivity for signals
25
Hz and the significant contami-
nation from instrumental noise in the band 25
60 Hz. The
sensitivity of the search falls as frequency increases, while
compute time rises dramatically. We terminate the search at
650 Hz, as in Ref.
[40]
.
The results of the search are presented in Fig.
4
, which
shows the frequency and recovered orbital parameters
a
0
and
φ
a
for every path with
S>S
th
. The color of the points
shows the Viterbi score associated with that path. As the
FIG. 4. Candidates identified by the search. The left-hand panel plots the detection score
S
(indicated by color; see color bar) as a
function of final frequency
f
0
ð
t
N
T
Þ
(horizontal axis) and orbital semimajor axis
a
0
(vertical axis) recovered by the HMM. The right-hand
panel plots the candidates with
T
asc
on the vertical axis. Undecorated candidates are eliminated by the known line veto, candidates
marked by blue circles are eliminated by the single-interferometer veto, candidates marked by orange squares are eliminated by the
T
obs
=
2
veto, and the candidates marked by green triangle survive the veto process.
B. P. ABBOTT
et al.
PHYS. REV. D
100,
122002 (2019)
122002-8
most a signal can wander during the observation is
N
T
Δ
f
drift
1
.
3
×
10
5
Hz, which is small compared to
Δ
f
band
(and what can be visually discerned on Fig.
4
), we
define
f
0
for a given path to be equal to
f
0
ð
t
¼
N
T
Þ
for
convenience.
To rule out false alarms, we apply the hierarchy of vetoes
first described in Ref.
[40]
. The vetoes are (1) the known
instrumental lines veto (described in Sec.
IVA 1
below),
(2) the single-interferometer veto (Sec.
IVA 2
), (3) the
T
obs
=
2
veto (Sec.
IVA 3
), and (4) the
T
drift
veto (ultimately
not used, but discussed in Sec. IVA 4 of Ref.
[40]
). To
ensure that the vetoes are unlikely to falsely dismiss a
true signal, we perform the search on a dataset with
synthetic signals injected into it, and ensure that those
injections are not vetoed. These veto safety tests are
described in Sec.
IV B
.
The number of candidates found in the initial search, and
then vetoed at each step, are listed in Table
IV
.
A. Vetoes
1. Known lines veto
There are a large number of persistent instrumental noise
lines identified as part of LIGO
s detector characterization
process
[67,68]
. These lines can arise from a number of
sources, including interference from equipment around the
detector, resonant modes in the suspension system, and
external environmental causes (e.g., the electricity grid).
A noise line generally produces high
j
F
a
j
and
j
F
b
j
values. The convolution in
(6)
reduces the impact of this
somewhat by summing bins near and far from the line, but
in practice the noise lines are strong enough that they
contaminate any candidate nearby. Accordingly, we veto
any candidate whose Viterbi path
f
0
ð
t
Þ
satisfies
j
f
0
ð
t
Þ
f
line
j
<
2
π
a
0
f
0
=P
, for any time
t
along the path
and for any line frequency
f
line
. This veto is efficient,
excluding 14 of the 20 candidates.
2. Single interferometer veto
During O2, L1 was slightly more sensitive than H1, but
overall the sensitivities of the two interferometers were
similar. Accordingly, any astrophysical signal that can be
detected in the combined dataset should either be detected by
the individual detector datasets when analyzed separately
(for stronger signals) or in neither (for weaker signals).
A signal that is detectable in one interferometer only is likely
to be a noise artifact, so we veto it.
Following Ref.
[40]
, we compare the Viterbi scores
obtained from individual detectors to the original combined
score
S
to classify survivors of the known line veto into
four categories discussed below, one of which is vetoed.
Category A
.
One detector returns
S<S
th
, while the
other detector returns
S>S
, and the frequency estimated
by the latter detector is close to that of the original
candidate
f
0
, that is,
j
f
0
f
0
j
<
2
π
a
0
f
0
=P
, where
the subscript
denotes a quantity estimated by the search
in both detectors. This category and the next represent
signals where the score is dominated by one detector. We
veto candidates in category A.
Category B
.
As with category A, one detector returns
S<S
th
, while the other detector returns
S>S
. Unlike
category A, the frequency estimated by the latter detector
is far from the original candidate, i.e.,
j
f
0
f
0
j
>
2
π
a
0
f
0
=P
. In this case, it is possible that there is signal
at
f
0
which is detectable when combining the data from
both detectors but not from one detector, because an artifact
masks its presence. Hence we keep the candidate for
follow-up.
Category C
.
The candidate is seen with
S>S
th
in both
detectors. This could either be a relatively strong signal, or
an artifact from a noise source common to both detectors.
The single-interferometer veto cannot distinguish these
possibilities. Again, we keep the candidate for follow-up.
Category D
.
The candidate is not seen by either
detector, with
S<S
th
in both detectors. This could be a
signal that is too weak to see in either detector individually.
We keep the candidate for follow-up.
Category A of the single-interferometer veto eliminates
two of the remaining six candidates. The two eliminated
candidates were stronger in H1 compared to L1.
3. T
obs
=
2
veto
We divide the observing run into two segments, the first
covering 140 days from November 30, 2016 (GPS time-
stamp 1164562334) to April 19, 2017 (GPS timestamp
1176658334), and the second covering 90 days from
January 19, 2017 (GPS timestamp 1168882334) to
August 25, 2017 (GPS timestamp 1187731792). This
division is chosen to get approximately equal effective
observing time in the two segments. There is no forceful
evidence to suggest that the gravitational wave strength of a
LMXB varies significantly with time (and a signal with
time-varying strength is likely to have a considerably more
complicated form than assumed here); thus we do not
expect a signal to appear preferentially in either segment.
We search the segments separately for the candidates which
survived both preceding vetoes. To determine whether to
veto candidates at this stage, we apply the same set of
categories as in veto 2.
TABLE IV. Number of candidates found in the first pass,
and number remaining after applying the vetoes described in
Sec.
IVA
.
After veto
Survivors
First pass
20
Line
6
Single interferometer
4
T
obs
=
2
3
SEARCH FOR GRAVITATIONAL WAVES FROM SCORPIUS X-1
...
PHYS. REV. D
100,
122002 (2019)
122002-9