Search for continuous gravitational waves from 20 accreting millisecond
x-ray pulsars in O3 LIGO data
R. Abbott
etal.
*
(LIGO Scientific Collaboration, Virgo Collaboration, and KAGRA Collaboration)
(Received 19 September 2021; accepted 16 December 2021; published 19 January 2022)
Results are presented of searches for continuous gravitational waves from 20 accreting millisecond x-ray
pulsars with accurately measured spin frequencies and orbital parameters, using data from the third
observing run of the Advanced LIGO and Advanced Virgo detectors. The search algorithm uses a hidden
Markov model, where the transition probabilities allow the frequency to wander according to an unbiased
random walk, while the
J
-statistic maximum-likelihood matched filter tracks the binary orbital phase.
Three narrow subbands are searched for each target, centered on harmonics of the measured spin frequency.
The search yields 16 candidates, consistent with a false alarm probability of 30% per subband and target
searched. These candidates, along with one candidate from an additional target-of-opportunity search done
for SAX J
1808
.
4
−
3658
, which was in outburst during one month of the observing run, cannot be
confidently associated with a known noise source. Additional follow-up does not provide convincing
evidence that any are a true astrophysical signal. When all candidates are assumed nonastrophysical, upper
limits are set on the maximum wave strain detectable at 95% confidence,
h
95%
0
. The strictest constraint is
h
95%
0
¼
4
.
7
×
10
−
26
from IGR J
17062
−
6143
. Constraints on the detectable wave strain from each target
lead to constraints on neutron star ellipticity and
r
-mode amplitude, the strictest of which are
ε
95%
¼
3
.
1
×
10
−
7
and
α
95%
¼
1
.
8
×
10
−
5
respectively. This analysis is the most comprehensive and sensitive
search of continuous gravitational waves from accreting millisecond x-ray pulsars to date.
DOI:
10.1103/PhysRevD.105.022002
I. INTRODUCTION
Second generation, ground-based gravitational wave
detectors, specifically the Advanced Laser Interferometer
Gravitational wave Observatory (Advanced LIGO)
[1]
and
Advanced Virgo
[2]
, have detected more than 50 compact
binary coalescence events in recent years
[3
–
5]
.
Continuous gravitational waves from rapidly rotating
neutron stars are also potential sources, e.g., a nonaxisym-
metry due to mountains on the surface, or stellar oscillation
modes in the interior
[6
–
8]
. There are no reported detec-
tions of continuous gravitational waves to date, despite a
number of searches in Advanced LIGO and Advanced
Virgo data
[9
–
30]
.
Low-mass x-ray binaries (LMXBs) are a high-priority
target for continuous gravitational wave searches. LMXBs
are composed of a compact object, such as a neutron star,
1
which accretes matter from a stellar-mass (
≲
1
M
⊙
)
companion
[31]
. The accretion exerts a torque that may
spin up the compact object. Electromagnetic (EM) obser-
vations show that even the pulsar with the highest known
frequency, PSR J
1748
−
2446
ad at 716 Hz
[32]
, rotates
well below the centrifugal break-up frequency, estimated at
∼
1400
Hz
[33]
. Gravitational wave emission may provide
the balancing torque in binary systems such as these,
stopping the neutron star from spinning up to the break-
up frequency
[34,35]
. If so, there should thus be a
correlation between accretion rate (which is inferred via
x-ray flux) and the strength of the continuous gravitational
wave emission
[34
–
37]
. The LMXB Scorpius X-1 is the
brightest extra-Solar x-ray source in the sky, making it a
prime target for searches for continuous gravitational waves
[11,18,38,39]
.
Some LMXBs have EM observations of pulsations
during
“
outburst
”
events lasting days to months, which
allow for measurement of their rotational frequency,
f
⋆
,to
an accuracy of
∼
10
−
8
Hz, and measurement of their binary
ephemerides
[31,40]
. LMXBs that are observed to go into
outburst and have measurable pulsations with millisecond
periods are sometimes called accreting millisecond x-ray
pulsars (AMXPs). If the rotational frequency is known,
computationally cheap narrowband searches are possible.
Six AMXPs were previously searched for continuous
gravitational waves, one in Science Run 6 (S6) using the
TwoSpect algorithm
[41,42]
, and five in Observing Run 2
(O2) using the same Hidden Markov Model (HMM)
*
Full author list given at the end of the article.
1
LMXBs in which the compact object is a stellar-mass black
hole are not expected to function as continuous gravitational
wave sources and are not discussed in this paper.
PHYSICAL REVIEW D
105,
022002 (2022)
2470-0010
=
2022
=
105(2)
=
022002(42)
022002-1
© 2022 American Physical Society
algorithm we use in this work
[12,43]
. No significant
candidates were found in either search. Searches for
continuous gravitational waves from LMXBs are difficult
as the rotation frequency may wander stochastically
on timescales of
≲
1
yr
[44]
, limiting the duration of
coherent integration. A HMM tracks a wandering signal,
and is the search algorithm we use here, following
Refs.
[11,12,43,45]
.
Advanced LIGO and Advanced Virgo began the third
Observing Run (O3) on April 1 2019,
15
∶
00
UTC. There
was a month-long commissioning break between October 1
2019,
15
∶
00
UTC, and November 1 2019,
15
∶
00
UTC,
after which observations resumed until March 27, 2020,
17
∶
00
UTC. This month-long break divides O3 into two
segments: O3a and O3b. In this work we search the full O3
dataset for continuous gravitational wave signals from
AMXPs with known rotational frequencies. The search
is a more sensitive version of an analogous search in O2
data
[12]
, with an expanded target list. We briefly review
the algorithm and O2 search in Sec.
II
. In Secs.
III
and
IV
we describe the targets and the parameter space respec-
tively. We discuss the data used in Sec.
V
. In Sec.
VI
we
describe the vetoes applied to discriminate between terres-
trial and astrophysical candidates. In Sec.
VII
we present
the results of the search. In Sec.
VIII
we describe an
additional target-of-opportunity search performed for one
of the targets that was in outburst during O3a. We provide
upper limits for the detectable wave-strain, and astrophysi-
cal implications thereof, in Sec.
IX
. We conclude in Sec.
X
.
II. SEARCH ALGORITHM
The search in this paper follows the same prescription as
the O2 searches for Scorpius X-1
[11]
and LMXBs with
known rotational frequency
[12]
. It is composed of two parts:
a HMM which uses the Viterbi algorithm to efficiently track
the most likely spin history, and the
J
-statistic, which
calculates the likelihood a gravitational wave is present
given the detector data, and the orbital parameters of both
the Earthandthe LMXB. The HMM formalismisidentical to
that used in Refs.
[11,12,38,43,45]
, and the
J
-statistic was
first introduced in Ref.
[43]
. Below, we provide a brief review
of both the HMM and the
J
-statistic.
A. HMM
In a Markov process, the probability of finding the
system in the current state depends only on the previous
state. In a hidden Markov process the states are not directly
observable and must be inferred from noisy data. In this
paper, the hidden state of interest is the gravitational wave
frequency
f
ð
t
Þ
. Although the rotation frequency
f
⋆
ð
t
Þ
of
every target in this search is measured accurately from EM
pulsations, we allow
f
ð
t
Þ
≠
f
⋆
ð
t
Þ
in general for three
reasons: (i) different emission mechanisms emit at different
multiples of
f
⋆
[46]
; (ii) a small, fluctuating drift may arise
between
f
ð
t
Þ
and
f
⋆
ð
t
Þ
, if the star
’
s core (where the
gravitational-wave-emitting mass or current quadrupole
may reside) decouples partially from the crust (to which
EM pulsations are locked)
[45,47]
; and, (iii) the rotational
frequency of the crust may also drift stochastically due to a
fluctuating accretion torque
[31,44]
. The gravitational-
wave frequency is therefore hidden even though the EM
measurement of
f
⋆
helps restrict the searched frequency
space, as described in Sec.
IV
.
Following the notation of Refs.
[11,12]
we label the
hidden state variable as
q
ð
t
Þ
. In our model, it transitions
between a discrete set of allowed values
f
q
1
;
...
;q
N
Q
g
at
discrete times
f
t
0
;
...
;t
N
T
g
. The probability of the state
transitioning from
q
i
at time
t
n
to
q
j
at time
t
n
þ
1
is
determined by the transition matrix
A
q
j
q
i
. In this search, as
in previous searches of LMXBs
[11,12,38]
, the transition
matrix is
A
q
j
q
i
¼
1
3
ð
δ
q
j
q
i
þ
1
þ
δ
q
j
q
i
þ
δ
q
j
q
i
−
1
Þ
;
ð
1
Þ
where
δ
ij
is the Kronecker delta. Equation
(1)
corresponds
to allowing
f
ð
t
Þ
to move 0, or
1
frequency bins, with
equal probability, at each discrete transition. It implicitly
defines the signal model for
f
ð
t
Þ
to be a piece-wise constant
function, with jumps in frequency allowed at the discrete
times
f
t
0
;
...
;t
N
T
g
. This is a well-tested approximation for
an unbiased random walk
[43,45]
.
The total duration of the search is
T
obs
, which we split
into
N
T
coherent equal chunks of length
T
drift
, where
N
T
¼b
T
obs
=T
drift
c
, and
b
...
c
indicates rounding down to
the nearest integer. We justify our choice of
T
drift
in Sec.
IV
.
In essence, it needs to be short enough to ensure that
f
⋆
ð
t
Þ
does not wander by more than one frequency bin during
each time segment, but ideally no shorter in order to
maximize the signal-to-noise ratio in each segment. For
each time segment the likelihood that the observation
o
j
is
related to the hidden state
q
i
is given by the emission matrix
L
o
j
q
i
. We calculate
L
o
j
q
i
from the data via a frequency
domain estimator, e.g., the
J
-statistic, as discussed in
Sec.
II B
.
The probability that the hidden path is
Q
¼
f
q
ð
t
0
Þ
;
...
;q
ð
t
N
T
Þg
given a set of observations
O
¼
f
o
ð
t
0
Þ
;
...
;o
ð
t
N
T
Þg
is
P
ð
Q
j
O
Þ¼
Π
q
ð
t
0
Þ
A
q
ð
t
1
Þ
q
ð
t
0
Þ
L
o
ð
t
1
Þ
q
ð
t
1
Þ
...
×
A
q
ð
t
N
T
Þ
q
ð
t
N
T
−
1
Þ
L
o
ð
t
N
T
Þ
q
ð
t
N
T
Þ
;
ð
2
Þ
where
Π
q
ð
t
0
Þ
is the prior probability of starting in the state
q
ð
t
0
Þ
, and is taken to be uniform within a certain range
guided by EM measurements of
f
⋆
. The Viterbi algorithm
is a computationally efficient way to find the path
Q
that
maximizes Eq.
(2) [48]
.
R. ABBOTT
et al.
PHYS. REV. D
105,
022002 (2022)
022002-2
The detection statistic we use in this work is
L
¼
ln
P
ð
Q
j
O
Þ
, i.e., the log-likelihood of the most likely path
given the data. The search outputs one
P
ð
Q
j
O
Þ
value per
frequency bin, corresponding to the optimal path
Q
terminating in that frequency bin.
B.
J
-statistic
Any long-lived gravitational wave signal from an LMXB
observed by the detectors is Doppler modulated by the
orbital motion of the detectors around the Solar System
barycenter, and by the orbital motion of the compact object
in its binary. The
F
-statistic is a frequency domain
estimator originally designed for isolated neutron stars,
and accounts for the Earth
’
s annual orbital motion (as well
as the amplitude modulation caused by the Earth
’
s diurnal
rotation)
[49]
. Algorithms that implement the
F
-statistic,
such as
lalapps_ComputeFstatistic_v2
[50]
,
have subsequently added functionality to account for
modulation of the signal due to binary motion.
The
J
-statistic accounts for the binary modulation via a
Jacobi-Anger expansion of the orbit
[43]
. It ingests
F
-statistic
“
atoms
”
as calculated for an isolated source
as an input, assumes the binary is in a circular orbit,
2
and
requires three binary orbital parameters: the period
P
, the
projected semimajor axis
a
0
, and the time of passage
of the ascending node
T
asc
. We use the
J
-statistic as the
frequency domain estimator
L
o
j
q
i
in this paper, as in
Refs.
[11,12]
. The
J
-statistic is a computationally efficient
algorithm, as it reuses
F
-statistic atoms when searching
over a template bank of binary orbital parameters.
III. TARGETS
The AMXPs chosen as targets for this search, along with
their positions, orbital elements, and pulsation frequencies
are listed in Table
I
. These 20 targets constitute all known
AMXPs with observed coherent pulsations and precisely
measured orbital elements as of April 2021.
3
For details on
the relevant EM observations, principally in the x-ray band,
see Refs.
[31,40,53,54]
.
Most AMXPs are transient, with
“
active
”
(outburst) and
“
quiescent
”
phases. Pulsations, and therefore
f
⋆
, are only
observed during the active phase. Active phases are
typically associated with accretion onto the neutron star,
however accretion can also happen during quiescence
[92]
.
The frequency derivatives,
_
f
⋆
, in the active phase and in the
quiescent phase are set by the accretion torque and
magnetic dipole braking respectively
[92,93]
. The value
of
_
f
⋆
has implications for the continuous gravitational
wave signal strength (see Sec.
IX C
), as well as the choice
of
T
drift
(see Sec.
IVA
).
One target, SAX J
1808
.
4
−
3658
, went into outburst
during O3a
[84,94,95]
. It may be the case that continuous
gravitational waves are only emitted when an AMXP is in
outburst
[96]
. If so, we increase our signal-to-noise ratio by
searching only data from the times that it was in outburst,
compared to searching the entirety of O3 data. To inves-
tigate this possibility, we perform in Sec.
VIII
an additional
target-of-opportunity search for continuous gravitational
waves from SAX J
1808
.
4
−
3658
while it is in outburst.
IV. SEARCH PARAMETERS
The
J
-statistic matched filter requires specification of
the source sky position [right ascension (RA) and decli-
nation (Dec)], the orbital period
P
, the projected semimajor
axis
a
0
, and the orbital phase
φ
a
at the start of the search.
The orbital phase can be equivalently specified via a time of
passage through the ascending node,
T
asc
. EM observations
constrain all of these parameters, as well as the spin
frequency
f
⋆
. These measurements, along with their
associated uncertainties, are listed in Table
I
.
There are several mechanisms that could lead to con-
tinuous gravitational wave emission from an AMXP, in its
active or quiescent phase.
“
Mountains
”
on the neutron star
surface, be they magnetically or elastically supported, emit
at
2
f
⋆
and potentially
f
⋆
[97]
. The dominant continuous
gravitational wave emission from
r
-mode oscillations
(Rossby waves excited by radiation-reaction instabilities)
is predicted to be at
∼
4
f
⋆
=
3
[98
–
101]
. Thus, we search
frequency subbands centered on
f
1
;
4
=
3
;
2
g
f
⋆
for each
target. As in Refs.
[11,12]
we choose a subband width
of
∼
0
.
61
Hz.
4
Recent work indicates that the continuous gravitational
wave signal from
r
-modes could emit at a frequency far
from
4
f
⋆
=
3
due to equation-of-state-dependent relativistic
corrections, and so comprehensive searches for
r
-modes
may need to cover hundreds of Hz for the targets listed in
Table
I [103,104]
. The exact range of frequencies to search
is a nonlinear function of
f
⋆
, and does not necessarily
include
4
f
⋆
=
3
(see Eq. (17) of Ref.
[104]
). However, these
estimates are still uncertain. We deliberately search
∼
0
.
61
Hz subbands centered on
4
f
⋆
=
3
, as an exhaustive
2
This assumption is justified as none of the targets described in
Sec.
III
have measurable eccentricity with sufficient precision
[31,40]
.
3
We do not include the AMXPAquila X-1
[51,52]
in our target
list as there is a large uncertainty on all three binary orbital
elements, compared to the other 20 AMXPs. One would need to
search
>
10
10
binary orbital templates, an order of magnitude
more than the rest of the targets combined. The number of binary
orbital templates is calculated as a function of the uncertainty in
orbital elements in Sec.
IV B
.
4
Other narrowband searches, such as Refs.
[10,102]
, search
subbands whose width,
∼
10
−
3
f
, scales with frequency. We note
that 0.61 Hz is comparable to
10
−
3
f
for the harmonics of
f
⋆
that
we search in this paper, but is
2
20
Δ
f
, where
Δ
f
is the frequency
bin size defined in Sec.
IVA
. Having the number of frequency
bins in the subband equal a power of two speeds up the Fourier
transform
[11]
.
SEARCH FOR CONTINUOUS GRAVITATIONAL WAVES FROM 20
...
PHYS. REV. D
105,
022002 (2022)
022002-3
TABLE I. Target list: position (RA and Dec), orbital period (
P
), projected semimajor axis in light-seconds (
a
0
), time of passage through the ascending node as measured near the
time of the most recent outburst (
T
asc
), the time of passage through the ascending node as propagated to the start of O3 (
T
asc
;
O
3
), as described in Sec.
IV B
, and frequency of
observed pulsations (
f
⋆
). Numbers in parentheses indicate reported
1
σ
errors (68% confidence level), unless otherwise noted. All objects have positional uncertainty
≤
1
sinRA
and
≤
0
.
5
00
in Dec.
Target
RA
Dec
P/s
a
0
=
lt-s
T
asc
=
GPS time
T
asc
;
O
3
=
GPS time
f
⋆
=
Hz
Refs.
IGR J
00291
þ
5934
00h29m03.05s
þ
59
°
34
0
18
.
93
00
8844.07673(9) 0.064993(2)
1122149932.93(5) 1238157687(1)
598.89213099(6)
[55,56]
MAXI J
0911
−
655
09h12m02.46s
−
64
°
52
0
06
.
37
00
2659.93312(47) 0.017595(9)
1145507148.0(9)
1238165918(16)
339.9750123(3)
[57,58]
XTE J
0929
−
314
09h29m20.19s
−
31
°
23
0
03
:
2
00
2614.746(3)
0.006290(9)
705152406.1(9)
1238165763(612) 185.105254297(9)
[59,60]
IGR J
16597
−
3704
16h59m32.902s
−
37
°
07
0
14
:
3
00
2758.2(3)
0.00480(3)
1193053416(9)
1238163777(4907) 105.1758271(3)
[61,62]
IGR J
17062
−
6143
17h06m16.29s
−
61
°
42
0
40
:
6
00
2278.21124(2) 0.003963(6)
1239389342(4)
1238165942(4)
163.656110049(9)
[63]
IGR J
17379
−
3747
17h37m58.836s
−
37
°
46
0
18
.
35
00
6765.8388(17) 0.076979(14) 1206573046.6(3)
1238162748(8)
468.083266605(7)
[64,65]
SAX J
1748
.
9
−
2021
17h48m52.161s
−
20
°
21
0
32
.
406
00
31555.300(3)
0.38757(2)
1109500772.5(8)
1238151731(12)
442.3610957(2)
[40,66]
NGC 6440 X-2
17h48m52.76s
−
20
°
21
0
24
:
0
00
3457.8929(7)
0.00614(1)
956797704(2)
1238166449(57)
205.89221(2)
[67,68]
IGR J
17494
−
3030
17h49m23.62s
−
30
°
29
0
58
.
999
00
4496.67(3)
0.015186(12)
1287797911(1)
1238163668(331) 376.05017022(4)
[69]
Swift J
1749
.
4
−
2807
17h49m31.728s
−
28
°
08
0
05
.
064
00
31740.8417(27) 1.899568(11) 1298634645.85(12) 1238136602(5)
517.92001385(6)
[70
–
72]
IGR J
17498
−
2921
17h49m56.02s
−
29
°
19
0
20
:
7
00
13835.619(1)
b
0.365165(5)
b
997147537.43(7)
b
1238164020(6)
400.99018734(9)
b
[73,74]
IGR J
17511
−
3057
17h51m08.66s
−
30
°
57
0
41
:
0
00
12487.5121(4)
0.2751952(18) 936924316.03(3) 1238160570(10)
244.83395145(9)
[75,76]
XTE J
1751
−
305
17h51m13.49s
−
30
°
37
0
23
:
4
00
2545.3414(38)
a
0.010125(5)
a
701914663.57(3)
a
1238164644(487) 435.31799357(3)
a
[77,78]
Swift J
1756
.
9
−
2508
17h56m57.43s
−
25
°
06
0
27
:
4
00
3282.40(4)
0.00596(2)
1207196675(9)
1238166119(378) 182.06580377(11)
[79]
IGR J
17591
−
2342
17h59m02.86s
−
23
°
43
0
08
:
3
00
31684.7503(5)
1.227714(4)
1218341207.72(8) 1238144176.7(3) 527.425700578(9)
[80,81]
XTE J
1807
−
294
18h06m59.8s
−
29
°
24
0
30
00
2404.4163(3)
0.004830(3)
732384720.7(3)
1238165711(63)
190.62350702(4)
[82,83]
SAX J
1808
.
4
−
3658
18h08m27.647s
−
36
°
58
0
43
.
90
00
7249.155(3)
0.062809(7)
1250296258.5(2)
1238161173(5)
400.97521037(1)
[84]
XTE J
1814
−
338
18h13m39.02s
−
33
°
46
0
22
:
3
00
15388.7229(2)
a
0.390633(9)
a
739049147.41(8)
a
1238151597(4)
314.35610879(1)
a
[85,86]
IGR J
18245
−
2452
18h24m32.51s
−
24
°
52
0
07
:
9
00
39692.812(7)
0.76591(1)
1049865088.37(9) 1238128096(33)
254.3330310(1)
[87,88]
HETE J
1900
.
1
−
2455
19h00m08.65s
−
24
°
55
0
13
:
7
00
4995.2630(5)
0.01844(2)
803963262.3(8)
1238161513(43)
377.296171971(5)
[89
–
91]
a
90% confidence level.
b
3
σ
error.
R. ABBOTT
et al.
PHYS. REV. D
105,
022002 (2022)
022002-4
broadband search lies outside the scope of this paper, which
aims to conduct fast, narrowband searches at astrophysi-
cally motivated harmonics of
f
⋆
while accommodating
frequency wandering within those subbands, a challenge in
its own right.
A.
T
drift
and frequency binning
Another key parameter for the search algorithm
described in Sec.
II
is the coherence time
T
drift
.Asin
Refs.
[11,12]
we fix
T
drift
¼
10
d for each target.
5
This
choice of
T
drift
is guided by observations of Scorpius X-1
[44]
. Quantitative studies of how x-ray flux variability in
AMXPs impacts searches for continuous gravitational
waves are absent from the literature. The choice to use
T
drift
¼
10
d balances the increased sensitivity achieved via
longer coherence times with the knowledge that the
gravitational wave frequency may wander stochastically,
e.g., due to fluctuations in the mass accretion rate. The
particular value
T
drift
¼
10
d has been adopted in all
previous Viterbi LMXB searches
[11,12,38]
and is justified
approximately with reference to a simple random-walk
interpretation of fluctuations in the x-ray flux of Scorpius
X-1
[44,105,106]
, but other values are reasonable too.
We remind the reader that the choice of
T
drift
implicitly
fixes the proposed signal model as one in which the
frequency may wander step-wise zero, plus or minus one
frequency bin every
T
drift
¼
10
d. The size of the frequency
bins,
Δ
f
, is fixed by the resolution implied by the coherence
time, i.e.,
Δ
f
¼
1
=
ð
2
T
drift
Þ¼
5
.
787037
×
10
−
7
Hz, for
T
drift
¼
10
d. As
Δ
f
depends on
T
drift
, changing the coher-
ence time explicitly changes the signal model, e.g., if
T
drift
is
halved and
T
obs
is kept constant, then both
N
T
and
Δ
f
double; thus thesignal can move up to a factoroffourmore in
frequency in the same
T
obs
. The connection between the
coherence time and signal model features in all semicoherent
search methods. However, for a HMM-based search such as
this, the choice of coherence time is not limited by computa-
tional cost, as it is in all-sky searches or searches based on the
F
-statistic
[17,107]
.
This analysis does not search over any frequency
derivatives. The maximum absolute frequency derivative,
j
_
f
max
j
, that does not change the frequency more than one
frequency bin over the course of one coherent chunk is
j
_
f
max
j¼
Δ
f
T
drift
≈
6
.
7
×
10
−
13
Hz s
−
1
:
ð
3
Þ
When measured, the long-term secular frequency derivative
is well below this value for all of our targets, see Sec.
IX C
for details.
B. Number of orbital templates
The orbital elements are known to high precision, with the
uncertainty in
P
satisfying
σ
P
≲
10
−
3
s, the uncertainty in
a
0
satisfying
σ
a
0
≲
10
−
4
light-seconds(lt-s),andtheuncertainty
in
T
asc
satisfying
σ
T
asc
≲
1
s. However,
T
asc
is measured
relative to the target
’
s most recent outburst, which is often
years before the start of O3 (
T
O
3
;
start
¼
1238166483
GPS
time). We need to propagate it forward in time. This
propagation compounds the uncertainty in
T
asc
, viz.
[11,12,39]
σ
T
asc
;
O
3
¼½
σ
2
T
asc
þð
N
orb
σ
P
Þ
2
1
=
2
;
ð
4
Þ
where
N
orb
is the number of orbits between the observed
T
asc
and
T
asc
;
O
3
. Henceforth
T
asc
and
σ
T
asc
symbolize their values
when propagated to
T
O
3
;
start
.
To conduct the search over the orbital elements for each
target and subband we construct a rectangular grid in
the parameter space defined by
ð
P
3
σ
P
;a
0
3
σ
a
0
;
T
asc
3
σ
T
asc
Þ
. For three targets, XTE J
0929
−
314
,IGR
J
16597
−
3704
, and IGR J
17494
−
3030
, the range
ð
T
asc
P=
2
Þ
is smaller than
ð
T
asc
3
σ
T
asc
Þ
and we use the former.
We assume that
P
and
a
0
remain within the same bin for the
entire search. While some targets have a nonzero measure-
ment of
_
PT
obs
(
_
a
0
T
obs
), in all cases it ismuch smaller than the
template spacing in
P
(
a
0
)
[56,63,108]
.
It is unlikely that the true source parameters lie exactly
on a grid point in the parameter space. Thus the grid is
spaced such that the maximum mismatch,
μ
max
,isnever
more than an acceptable level. The mismatch is defined as
the fractional loss in signal-to-noise ratio between the
search executed at the true parameters and at the nearest
grid point
[109]
. We calculate the number of grid points
required for
P
,
a
0
and
T
asc
using Eq. (71) of Ref.
[109]
, i.e.,
N
P
¼
π
2
ffiffiffi
6
p
μ
−
1
=
2
max
fa
0
γ
T
drift
P
2
σ
P
;
ð
5
Þ
N
a
0
¼
3
π
ffiffiffi
2
p
μ
−
1
=
2
max
f
σ
a
0
;
ð
6
Þ
N
T
asc
¼
6
π
2
ffiffiffi
2
p
μ
−
1
=
2
max
fa
0
1
P
σ
T
asc
;
ð
7
Þ
where
γ
is a refinement factor defined in general in Eq. (67)
of Ref.
[109]
. In the case of O3, the semicoherent segments
are contiguous so we have
γ
¼
N
T
¼
36
.Wefix
μ
max
¼
0
.
1
. A set of software injections into O3 data verifies that a
template grid constructed with
μ
max
¼
0
.
1
results in a
maximum fractional loss in signal-to-noise ratio of 10%.
We make the conservative choice of rounding
N
P
,
N
a
0
, and
N
T
asc
up to the nearest integer, after setting
f
to the highest
frequency in each 0.61 Hz subband. As in Ref.
[12]
we find
N
a
0
¼
1
for each target and subband, and so hold
a
0
constant at its central value while searching over
P
and
T
asc
.
5
We consider additional
T
drift
durations for the target-of-
opportunity search for continuous gravitational waves from
SAX J
1808
.
4
−
3658
during its O3a outburst in Sec.
VIII
.
SEARCH FOR CONTINUOUS GRAVITATIONAL WAVES FROM 20
...
PHYS. REV. D
105,
022002 (2022)
022002-5
Table
II
shows
N
P
,
N
T
asc
, and
N
tot
¼
N
P
N
T
asc
for each
target and subband. When Eq.
(5)
or
(7)
predicts only two
templates for a given subband we round up to three,
ensuring that the central value of
P
or
T
asc
from EM
observations is included in the template bank. Note that the
EM observations are sufficiently precise that
<
5
×
10
4
templates are required across all targets and subbands. This
is in contrast to the O2 search for continuous gravitational
waves from Scorpius X-1, for which
∼
10
9
templates were
needed, mainly due to the large uncertainty in
a
0
, and the
unknown rotation frequency
[11]
.
C. Thresholds
The output of the search algorithm outlined in Sec.
II
is a
L
value corresponding to the most likely path through each
subband for each orbital template
ð
P; a
0
;T
asc
Þ
. We flag a
template for further follow-up if
L
exceeds a threshold,
L
th
,
given an acceptable probability of false alarm. To determine
L
th
we need to know how often pure noise yields
L
>
L
th
.
The distribution of
L
in noise-only data is unknown
analytically, but depends on
P
,
a
0
, and the frequency, so
Monte Carlo simulations are used to determine
L
th
in each
subband for each target.
We estimate the distribution of
L
in noise via two
methods: (i) using realizations of synthetic Gaussian noise
generated using the
lalapps_Makefakedata_v5
pro-
gram in the LIGO Scientific Collaboration Algorithm
Library (LALSuite)
[50]
, and (ii) searching O3 data in
off-target locations to simulate different realizations of true
detector noise. As in Refs.
[11,12]
we generate realizations
for each target and subband, and apply the search algorithm
described in Sec.
II
to each realization to recover samples
from the noise-only distribution of
L
. Details on how we
use these samples to find
L
th
for each subband are given in
the Appendix
A
. Unless otherwise noted,
L
th
refers to the
lower of the two thresholds derived from the methods listed
above to minimize false dismissals.
To define
L
th
we must also account for a
“
trials factor
”
due to the number of templates searched in each subband.
We assume that in noise-only data the spacing between
templates is sufficiently large such that each template
TABLE II. Starting frequencies,
f
s
, for each
∼
0
.
61
Hz-wide subband, number of templates needed to cover the
P
and
T
asc
domains in
that subband,
N
P
and
N
T
asc
respectively, and the total number of templates for each subband,
N
tot
¼
N
P
N
T
asc
. The projected semimajor
axis
a
0
is known precisely enough that we have
N
a
0
¼
1
for each subband.
Target
f
s
(Hz)
N
P
N
T
asc
N
tot
Target
f
s
(Hz)
N
P
N
T
asc
N
tot
IGR J
00291
þ
5934
598.6
1
3
3
IGR J
17498
−
2921
400.7
1
17
17
798.5
1
3
3
534.7
1
22
22
1197.8
1
3
3
802.0
3
33
99
MAXI J
0911
−
655
339.7
1
10
10
IGR J
17511
−
3057
244.5
1
14
14
453.3
3
14
42
326.4
1
19
19
679.9
3
20
60
489.7
1
28
28
XTE J
0929
−
314
184.8
3
52
156
XTE J
1751
−
305
435.0
4
195
780
246.8
3
69
207
580.4
5
260
1300
370.2
3
104
312
870.6
8
390
3120
IGR J
16597
−
3704
104.9
49
23
1127
Swift J
1756
.
9
−
2508
181.8
10
34
340
140.2
65
31
2015
242.8
13
45
585
210.4
97
46
4462
364.1
20
67
1340
IGR J
17062
−
6143
163.4
1
1
1
IGR J
17591
−
2342
527.1
1
3
3
218.2
1
1
1
703.2
3
3
9
327.3
1
1
1
1054.9
3
4
12
IGR J
17379
−
3747
467.8
4
12
48
XTE J
1807
−
294
190.3
1
7
7
624.1
5
15
75
254.2
1
9
9
936.2
7
23
161
381.2
1
13
13
SAX J
1748
.
9
−
2021
442.1
3
18
54
SAX J
1808
.
4
−
3658
400.7
4
5
20
589.8
3
24
72
534.6
5
7
35
884.7
3
36
108
802.0
7
10
70
NGC 6440 X-2
205.6
1
6
6
XTE J
1814
−
338
314.1
1
9
9
274.5
1
8
8
419.1
1
12
12
411.8
1
12
12
628.7
1
17
17
IGR J
17494
−
3030
375.7
21
112
2352
IGR J
18245
−
2452
254.0
3
44
132
501.4
27
150
4050
339.1
3
58
174
752.1
41
224
9184
508.7
5
87
435
Swift J
1749
.
4
−
2807
517.6
7
43
301
HETE J
1900
.
1
−
2455
377.0
1
17
17
690.6
9
57
513
503.1
1
22
22
1035.8
13
85
1105
754.6
1
33
33
R. ABBOTT
et al.
PHYS. REV. D
105,
022002 (2022)
022002-6
returns a statistically independent
L
. We can therefore
relate the false alarm probability for a search of a subband
with
N
tot
templates,
α
N
tot
, to the probability of a false alarm
for a single template,
α
, viz.
α
N
tot
¼
1
−
ð
1
−
α
Þ
N
tot
:
ð
8
Þ
Previous comparable searches have set
α
N
tot
between 0.01
and 0.3
[11,12,38,39]
. In this search, we fix
α
N
tot
¼
0
.
3
, i.e.,
set the acceptable probability of false alarm at 30% per
subband. As we search a total of
20
×
3
¼
60
subbands, we
expect
∼
18
candidates above
L
th
due to noise alone (i.e.,
false alarms), a reasonable number on which to perform
more exhaustive follow-up. Looking ahead to the results in
Sec.
VII
we recover 4611 candidates above
L
th
. While this
number is much higher than the
∼
18
false alarms expected,
almost all of these candidates are non-Gaussian noise
artifacts in one (or both) of the detectors. All but 16 of
the 4611 candidates are eliminated by the vetoes outlined in
Sec.
VI
. We reiterate that
L
th
in each subband is the lower
of the two thresholds described in Appendix
A
, lowering
conservatively the probability of false dismissal.
D. Computing resources
A mix of central processing unit (CPU) and graphical
processing unit (GPU) resources are used. The GPU
implementation of the
J
-statistic is identical to that used
in Refs.
[11,12]
. The entire search across all targets and
subbands takes
∼
30
CPU-hours and
∼
40
GPU-hours when
using compute nodes equipped with Xeon Gold 6140 CPUs
and NVIDIA P100 12 GB PCIe GPUs. Producing
L
th
for
each subband, as described in Sec.
IV C
, takes an additional
∼
5
×
10
2
CPU-hours and
∼
4
×
10
3
GPU-hours to perform
the search on different noise realizations. The additional
follow-up in Appendix
B1
requires an additional
∼
10
3
CPU-hours and
∼
10
2
GPU-hours.
V. O3 DATA
We use the full dataset from O3, spanning from April 1,
2019,
15
∶
00
UTC to March 27, 2020,
17
∶
00
UTC, from
the LIGO Livingston and Hanford observatories. We do not
use any data from the Virgo interferometer in this analysis,
due to its lower sensitivity compared to the two LIGO
observatories in the frequency subbands over which we
search
[110]
. The data products ingested by the search
algorithm described in Sec.
II
are short Fourier transforms
(SFTs) lasting 1800s. Times when the detectors were
offline, poorly calibrated, or were impacted by egregious
noise, are excluded from analysis by using
“
Category 1
”
vetoes as defined in section V. 2 of Ref.
[110]
. The SFTs are
generated from the
“
C01 calibrated self-gated
”
dataset,
which is the calibrated strain data with loud transient
glitches removed
[111]
. Transient glitches otherwise
impact the noise floor, as described in section
VI
.1 of
Ref.
[110]
. The median systematic error of the strain
magnitude across O3 is
<
2%
[112,113]
.
The coherence time
T
drift
¼
10
d splits the data into
N
T
¼
36
segments. However, due to the month-long
commissioning break between O3a and O3b there are
two segments without any SFTs. These two segments,
starting at October 8, 2019,
15
∶
00
UTC and October 15,
2019,
15
∶
00
UTC, are replaced with a uniform log-like-
lihood for all frequency bins, which allows the HMM to
effectively skip over them while still allowing spin wander-
ing. When generating synthetic data in Secs.
IV C
and
IX
the same two data segments are also replaced with uniform
log-likelihoods to emulate the real search.
VI. VETOES
When a candidate is returned with
L
>
L
th
we must
decide whether there are reasonable grounds to veto the
candidate as nonastrophysical. We use three of the vetoes
from Ref.
[12]
: the known line veto, detailed in Sec.
VI A
,
the single interferometer veto, detailed in Sec.
VI B
, and the
off-target veto, detailed in Sec.
VI C
. The false dismissal
rate of these vetoes is less than 5% (see detailed safety
investigations in Sec. IV B of Ref.
[38]
and Sec. IV B
of Ref.
[11]
).
A. Known line veto
As part of the detector characterization process many
harmonic features are identified as instrumental
“
known
lines
”
[110,114]
. However, the exact source of these
harmonic features is sometimes unidentified, and their
impact cannot always be mitigated through isolating
hardware components or post-processing the data
[110,114]
. We use the vetted known lines list in Ref.
[115]
.
Any candidate close to a known line at frequency
f
line
is
vetoed. Precisely, if for any time
0
≤
t
≤
T
obs
the candi-
date
’
s frequency path
f
ð
t
Þ
satisfies
j
f
ð
t
Þ
−
f
line
j
<
2
π
a
0
f
line
=P;
ð
9
Þ
then the candidate is vetoed.
6
B. Single interferometer veto
Aninstrumentalartifactisunlikely tobecoincidentinboth
detectors, so the candidate
’
s
L
should be dominated by only
one of the detectors if the signal is nonastrophysical. On the
6
One might consider an additional Doppler broadening factor
of
2
π
a
⊕
=
1
yr, where
a
⊕
is the mean Earth-Sun distance, as
stationary lines in the detector frame get Doppler shifted when
transforming the data to the frame of reference of the source. We
opt not to apply this factor for simplicity in this search, as the
exact pattern of Doppler modulation depends strongly on the sky
location of the target. Looking ahead to the results in Sec.
VII
,we
note that none of the 16 surviving candidates is within
2
π
fa
⊕
=
1
yr of any known line.
SEARCH FOR CONTINUOUS GRAVITATIONAL WAVES FROM 20
...
PHYS. REV. D
105,
022002 (2022)
022002-7
other hand, an astrophysical signal may need data from both
detectors to be detected, or if it is particularly strong may be
seen in both detectors individually.
We label the original log-likelihood as
L
∪
, and we also
calculate the two single interferometer log-likelihoods
L
a
and
L
b
(where the higher
L
is labeled with
b
for definite-
ness). There are four possible outcomes for this veto:
(1) If the
L
value in one detector is subthreshold, while
the other is above the two-detector
L
value, i.e., one
has
L
a
<
L
th
and
L
b
>
L
∪
and
f
b
ð
t
Þ
, the frequency
path associated with
L
b
, is close to the frequency
path of the candidate when using data from both
detectors,
f
∪
ð
t
Þ
, i.e.,
j
f
∪
ð
t
Þ
−
f
b
ð
t
Þj
<
2
π
a
0
f
∪
=P;
ð
10
Þ
then the candidate is likely to be a noise artifact in
detector
b
, and is vetoed.
(2) If one has
L
a
<
L
th
and
L
b
>
L
∪
, but Eq.
(10)
does
not hold then the candidate signal cannot be vetoed,
as the single-interferometer searches did not find the
same candidate. This could indicate that the candi-
date is a weak astrophysical signal that needs data
from both detectors to be detectable.
(3) If one has
L
a
>
L
th
and
L
b
>
L
th
, the candidate
could represent a strong astrophysical signal that is
visible in data from both detectors independently, or
it could represent a common noise source. Candi-
dates in this category cannot be vetoed.
(4) If one has
L
a
<
L
th
and
L
b
<
L
∪
, data from both
detectors is needed for the candidate to be above
threshold, possible indicating a weak astrophysical
signal. Candidates in this category cannot be vetoed.
C. Off-target veto
The third veto we apply to a candidate is to search an off-
target sky position with the same orbital template. If the off-
target search returns
L
>
L
th
then the candidate is likely
instrumental rather than astrophysical. For this veto, off-
target corresponds to shifting the target sky position
þ
40
m
in RA and
þ
10
° in Dec.
VII. O3 SEARCH RESULTS
The results of the search of all 20 targets are summarized
in Fig.
1
, with
α
N
tot
¼
0
.
3
, i.e., a nominal probability of
false alarm per subband of 30%. Each symbol indicates,
for all templates with
L
>
L
th
, the terminating frequency
bin and
p
noise
, the probability that a search of that
candidate
’
s subband in pure noise would return at least
one candidate at least as loud as the one seen.
Equation
(A6)
in Appendix
A5
defines
p
noise
explicitly.
Each candidate is colored according to
L
. We note that high
L
does not always correspond to low
p
noise
due to the
FIG. 1. Summary of search results across all targets and subbands with
L
>
L
th
. The different symbols correspond to candidates from
different targets. The ordinate shows
p
noise
for each candidate, the probability that a search of that candidate
’
s subband in pure noise
would return at least one candidate at least as loud as the one seen. The color of each candidate indicates
L
(see color bar at right).
Candidates that are eliminated by the vetoes outlined in Sec.
VI
are not shown for clarity. Details on the search results are in Sec.
VII
and
Appendix
B
.
R. ABBOTT
et al.
PHYS. REV. D
105,
022002 (2022)
022002-8