Draft version June 29, 2022
Typeset using L
A
T
E
X
twocolumn
style in AASTeX631
Narrowband searches for continuous and long-duration transient gravitational waves from known
pulsars in the LIGO-Virgo third observing run
The LIGO Scientific Collaboration, the Virgo Collaboration, and the KAGRA Collaboration
(Dated: June 29, 2022)
ABSTRACT
Isolated neutron stars that are asymmetric with respect to their spin axis are possible sources of
detectable continuous gravitational waves. This paper presents a fully-coherent search for such signals
from eighteen pulsars in data from LIGO and Virgo’s third observing run (O3). For known pulsars,
efficient and sensitive matched-filter searches can be carried out if one assumes the gravitational ra-
diation is phase-locked to the electromagnetic emission. In the search presented here, we relax this
assumption and allow the frequency and frequency time-derivative of the gravitational waves to vary
in a small range around those inferred from electromagnetic observations. We find no evidence for
continuous gravitational waves, and set upper limits on the strain amplitude for each target. These
limits are more constraining for seven of the targets than the spin-down limit defined by ascribing
all rotational energy loss to gravitational radiation. In an additional search we look in O3 data for
long-duration (hours–months) transient gravitational waves in the aftermath of pulsar glitches for six
targets with a total of nine glitches. We report two marginal outliers from this search, but find no
clear evidence for such emission either. The resulting duration-dependent strain upper limits do not
surpass indirect energy constraints for any of these targets.
1.
INTRODUCTION
Continuous gravitational waves (CWs) are quasi-
monochromatic signals expected to be ever-present in
the data of gravitational-wave (GW) detectors such
as Advanced LIGO (Aasi et al. 2015a) and Advanced
Virgo (Acernese et al. 2015). While the observation of
transient GWs from compact binary coalescences has
become nearly commonplace (Abbott et al. 2021b), CWs
have yet to be detected as of the third observing run
(O3). One of the most enticing and commonly sought
after sources of CWs is a rapidly spinning, asymmetric
neutron star (NS); see Sieniawska & Bejger (2019) and
Haskell & Schwenzer (2020) for recent reviews. In the
case of a triaxial NS, CW emission occurs at twice the
rotation frequency of the star.
Many NSs are observed as pulsars by radio, X-ray or
γ
-ray telescopes (Lorimer & Kramer 2012). Pulsars can
often be timed extremely precisely—in the best cases
the arrival of new pulses can be predicted to within tens
of nano-seconds. This precision can enable exciting sci-
ence, including sensitive tests of general relativity (Wex
& Kramer 2020), placing constraints on the equation
of state of the dense matter inside NSs (Lattimer &
Prakash 2004; Kramer & Wex 2009; Ho et al. 2015), and
using deviations in timing residuals of pulsars to search
for a stochastic gravitational-wave background (Verbiest
et al. 2016; Arzoumanian et al. 2020; Goncharov et al.
2021). Detecting GWs from spinning NSs would add
a completely new messenger to the study of these ex-
treme objects (Glampedakis & Gualtieri 2018; Haskell
& Schwenzer 2020).
In this paper, we search in LIGO–Virgo O3 data
(taken in 2019–2020) for CWs from NSs that have been
observed as pulsars and precisely timed in either the ra-
dio or X-ray bands. We identify eighteen promising can-
didates for which the observed spin-down (negative fre-
quency derivative) implies indirect limits on CW emis-
sion that fall within a factor of 3 of the expected sen-
sitivity of the search. Some other analyses assume the
phase of CWs to be locked to the rotational phase of the
crust of the star as observed by electromagnetic (EM)
telescopes (Aasi et al. 2014; Abbott et al. 2017a, 2018,
2019a, 2020b, 2021c,a; Ashok et al. 2021; Nieder et al.
2020). Here we relax that assumption, and allow for the
frequency of rotation, and its derivative, to differ from
the EM-observed values by a small factor: the so-called
“narrowband” search approach (Abbott et al. 2008; Aasi
et al. 2015b; Abbott et al. 2017b, 2019b; Ashok et al.
2021; Nieder et al. 2020). We use two separate anal-
ysis pipelines, the 5
n
-vector (Astone et al. 2014; Mas-
trogiovanni et al. 2017) and the frequency-domain
F
-
statistic (Jaranowski et al. 1998; Wette et al. 2018b)
arXiv:2112.10990v2 [gr-qc] 27 Jun 2022
2
pipelines, to perform phase-coherent searches for CWs
on O3 data over our widened parameter space. Using
two separate pipelines allows us to cross-check results
between pipelines, compare limits set by the two meth-
ods, and increase confidence in any potential detection
by requiring both pipelines to see the same signal.
A scenario where GW and EM emission have simi-
lar, but slightly different phase evolution is plausible,
e.g.
, when there is differential rotation between the rigid
crust and superfluid parts of the star. Possible observa-
tional evidence for this comes, for example, from pulsar
glitches (Link et al. 1992, 2000; Lyne et al. 2000; Fuentes
et al. 2017; Haskell 2017): sudden spin-up events that
are often followed by an exponential relaxation back to
the simple spin-down scenario. In many models (Haskell
& Melatos 2015), the superfluid and non-superfluid com-
ponents build up a lag with respect to one another, and
a glitch represents the sudden re-coupling of the two
components (Anderson & Itoh 1975; Lyne et al. 2000).
Glitches are also directly relevant for GW searches:
first, some of our CW search targets glitched during O3.
For these, we perform separate phase-coherent searches
covering the parts of O3 before and after the glitches.
Second, it is also possible that glitches trigger increased
GW emission in their aftermath (van Eysden & Melatos
2008; Bennett et al. 2010; Prix et al. 2011; Melatos et al.
2015; Singh 2017; Yim & Jones 2020). Hence, we also
perform additional searches for long-duration transient
CW-like signals with durations from a few hours to four
months after observed glitches during (or shortly before)
O3, covering nine glitches from six pulsars.
The rest of this paper is structured as follows. We
discuss our selection of target pulsars, and the EM ob-
servations that we use to guide our search in Sec. 2 and
the O3 GW data set in Sec. 3. In Sec. 4 we describe
our analysis methods. We present the results of the two
CW searches in Sec. 5. We then cover the detailed
setup and the results of the search for GW emission in
the aftermath of pulsar glitches in Sec. 6. In Sec. 7
we conclude the paper with a discussion of the results
obtained by all three analyses and their astrophysical
implications. Throughout the rest of the paper we will
often refer to the analyses using the
F
-statistic and 5
n
-
vector pipelines search over all data as the two “CW”
searches, and the search for long-duration transients af-
ter glitches as the “transient search.”
2.
ELECTROMAGNETIC DATA AND TARGET
SELECTION
We start with timing solutions from EM observations
of pulsars (also referred to as
ephemerides
), and search
in a narrow band in frequency and spin-down, as further
described below in Sec. 4. The widths of the frequency
and spin-down search bands are usually larger than the
uncertainty in the timing solutions. The observations we
use were made in radio and X-ray bands, and were pro-
vided by the following observatories: the Canadian Hy-
drogen Intensity Mapping Experiment (CHIME) (Ban-
dura et al. 2014; Amiri et al. 2021), University of Tas-
mania’s Mount Pleasant Observatory 26 m telescope,
the 42 ft telescope and Lovell telescope at Jodrell Bank,
the MeerKAT telescope (observations made as part of
the MeerTime project, Bailes et al. 2020), the Nan ̧cay
Decimetric Radio Telescope, the Neutron Star Interior
Composition Explorer (NICER) (Gendreau 2016) and
the UTMOST timing program with the Molonglo Ob-
servatory Synthesis Telescope (MOST) (Bailes et al.
2017; Jankowski et al. 2018; Lower et al. 2020). The
Tempo
(Nice et al. 2015) and
Tempo2
(Edwards et al.
2006) timing packages were used to fit the model pa-
rameters and provide timing solutions.
We select our targets for the CW searches in a similar
manner to Abbott et al. (2019b), based on the sensitive
frequency band of the Advanced LIGO and Virgo de-
tectors and availability of precise ephemerides over the
duration of O3 from EM observations. We have ana-
lyzed all isolated pulsars, except for one, with a rota-
tion frequency between 10 and 350 Hz and for which the
spin-down limit falls within a factor of 3 of the expected
sensitivity of the full network over the course of O3.
This frequency range includes all the high-value targets
identified in Abbott et al. (2021a) for which it could be
possible to go below the spin-down limit. Pulsars in this
frequency range would produce CWs between 20 Hz and
700 Hz. The only pulsar we do not analyze that satisfies
these criteria is PSR J0537–6910, which was analyzed in
detail using a narrowband approach to search for
r
-mode
emission in Abbott et al. (2021d) and using a targeted
approach in Abbott et al. (2021c). Searches lasting up
to 120 days during inter-glitch periods for this pulsar are
performed by the post-glitch transient search presented
in Sections 4.4 and 6.1.
We estimate the expected sensitivity as
h
sens
= Θ
√
S
(
f
)
T
obs
,
(1)
where
S
(
f
) is the power spectral density as a function of
frequency for the LIGO Hanford, LIGO Livingston and
Virgo detectors (H1, L1 and V1 respectively) and
T
obs
is the observing time assuming 11 months of data with
a duty cycles of 75%, 77% and 76% respectively. The
factor Θ
∼
30 encodes the scaling of typical narrowband
searches but is a function of the number of templates
employed in the analysis (Astone et al. 2014).
3
On the other hand, the spin-down limit is an indirect
upper limit (UL) on the GW amplitude assuming all of
a pulsar’s rotational energy loss comes in the form of
GWs (Jaranowski et al. 1998; Prix 2009). It is given by
h
sd
≈
2
.
55
×
10
−
25
1 kpc
d
√
I
38
|
̇
f
spin
|
10
−
13
Hz
/
s
1 Hz
f
spin
,
(2)
where
I
38
is the star’s moment of inertia in units of
10
38
kg m
2
,
f
spin
is its rotation frequency, and
|
̇
f
spin
|
is
the absolute value of its spin-down rate (Abbott et al.
2019b). We have computed spin-down limits accord-
ing to the most recent distance estimates given in the
ATNF catalog (Manchester et al. 2005), version 1.64,
and extrapolating the rotational frequencies and spin-
down rates to O3. For several pulsars, we have used
more recent distance estimates from the literature.
In Table 1 we present a list of targets along with the
ranges of GW frequency and spin-down parameters and
corresponding number of templates covered for the two
pipelines used to conduct the CW searches. We discuss
our parameter range choices for each pipeline in Sec. 4.
The observatory yielding the timing solution we use for
each source is noted in the far right column.
Two pulsars on our list of CW search targets,
PSR J0534+2200 (the Crab) and PSR J1105–6107,
glitched during O3 (Shaw et al. 2021), on 2019 July 23
and 2019 April 9 respectively. One other, PSR J1813–
1749, shows marginal evidence for a glitch on or close to
2019 Aug 03 (Ho et al. 2020a).
For this reason, we perform two separate searches, be-
fore and after the estimated glitch epoch, which are iden-
tified in Table 1 with suffixes “bg” and “ag” for “before
glitch” and “after glitch” respectively.
For PSR J0534+2200, we use O3 data until the last
observation before the glitch. The analysis after the
glitch uses data from ten days after the glitch epoch,
accounting for the estimated relaxation time, until the
end of O3.
For PSR J1105–6107, we only search after the glitch
since the glitch occurs nearly at the start of the run. The
after-glitch search uses data from two days after the esti-
mated glitch until the end of O3. For PSR J1813–1749,
we perform two separate searches: one fully-coherent
search across the full O3 duration, assuming the pul-
sar did not glitch, and a search before the glitch. No
search after the glitch is performed since the glitch oc-
curs near the end of the EM timing model (although a
post-glitch transient search is conducted assuming the
timing model remains valid within uncertainties; see be-
low). Further details about these glitches will also be
discussed in Sec. 6.
For the post-glitch transient search, we have used
information from the glitch catalogs maintained at
ATNF (Hobbs et al. 2021) and Jodrell Bank (Espinoza
et al. 2011; Basu et al. 2022; Shaw et al. 2021). We have
selected six pulsars with GW frequencies
f >
15 Hz and
with glitches observed during (or shortly before) O3:
PSR J0534+2200, PSR J0537–6910, PSR J0908–4913,
PSR J1105–6107, PSR J1813–1749, and PSR J1826–
1334. Another pulsar within our frequency band of in-
terest, PSR J2021+3651, was observed to glitch during
O3 by Jodrell Bank, but the glitch time uncertainty is
too large to make our transient search setup feasible. As
mentioned before, for PSR J1813–1749 it is not certain
if a glitch actually occurred (Ho et al. 2020a), but we
perform an opportunistic search here with its assumed
parameters.
The
ephemerides
for
PSR
J0534+2200
and
PSR J1826–1334 were provided by Jodrell Bank; a de-
tailed discussion of the PSR J0534+2200 glitch is given
in Shaw et al. (2021). Ephemerides for PSR J0908–4913
and PSR J1105–6107 are derived from UTMOST radio
observations; these are discussed further in Lower et al.
(2019, 2020). A potential ambiguity in timing solutions
from periodically-scheduled surveys like UTMOST, par-
ticularly affecting glitch size estimates, was discussed
by Dunn et al. (2021), but for these two targets it has
been confirmed that the provided values are correct. For
PSR J1813–1749 and PSR J0537–6910 we use NICER
X-ray observations, as reported in Ho et al. (2020a,b)
and Abbott et al. (2021c). Details for all targets will be
listed in Sec. 6.1.
3.
GW DATA
We analyze GW strain data taken at the LIGO Han-
ford Observatory (H1), LIGO Livingston Observatory
(L1) and Virgo (V1), during their O3 observing run.
O3 consisted of two separate sections, separated by a
month-long commissioning break. O3a ran from 2019
April 1, 15:00 UTC until 2019 October 1, 15:00 UTC.
O3b ran from 2019 November 1, 15:00 UTC, to 2020
March 27, 17:00 UTC. See Buikema et al. (2020) and
Acernese et al. (2019) for general descriptions of the per-
formances of the LIGO and Virgo detectors respectively
during O3. The calibration of this data set and its un-
certainty budget are described in Sun et al. (2020, 2021)
and Acernese et al. (2022a).
Data preparation for all searches starts with removing
times of large transient noise whose cause is known, re-
ferred to as “CAT1” vetoes (Davis et al. 2021; Acernese
et al. 2022b) from calibrated strain data.
4
Table 1
. Set-up parameters for fully coherent CW search pipelines.
Name
R.A.
DEC.
f
∆
f
(5v) ∆
f
(
F
)
̇
f
∆
̇
f
(5v)
∆
̇
f
(
F
)
̈
f
∆
̈
f n
5v
total
n
F
total
Ref.
×
10
−
13
×
10
−
15
×
10
−
15
×
10
−
23
×
10
−
23
×
10
7
×
10
7
Hz
Hz
Hz
Hz s
−
1
Hz s
−
1
Hz s
−
1
Hz s
−
2
Hz s
−
2
J0534+2200 bg 05
h
34
m
31
.
97
s
+22
◦
00
′
52
.
07
′′
59
.
241
...
0
.
24
−
7370
.
0
...
2900
.
0
2360
.
0
5
.
2
...
983
.
7
a
J0534+2200 ag
05
h
34
m
31
.
97
s
+22
◦
00
′
52
.
07
′′
59
.
241
0
.
24
0
.
24
−
7370
.
0
2300
.
0
2900
.
0
2360
.
0
5
.
2
498
.
0 9362
.
0
a
J0711–6830
07
h
11
m
54
.
18
s
−
68
◦
30
′
47
.
37
′′
364
.
234
0
.
72
1
.
5
−
0
.
00989
2
.
1
0
.
004
0
.
0
0
.
0
4
.
391
45
.
85
b
J0835–4510
08
h
35
m
20
.
52
s
−
45
◦
10
′
34
.
28
′′
22
.
371
0
.
089
0
.
089
−
313
.
0
130
.
0
120
.
0
504
.
0
0
.
0
32
.
9
281
.
6
c
J1101–6101
11
h
01
m
44
.
96
s
−
61
◦
01
′
39
.
6
′′
31
.
846
0
.
13
0
.
13
−
45
.
3
19
.
0
18
.
0
0
.
0
0
.
0
7
.
026
61
.
52
d
J1105–6107
11
h
05
m
25
.
71
s
−
61
◦
07
′
55
.
63
′′
31
.
644
0
.
13
0
.
13
−
79
.
7
32
.
0
32
.
0
554
.
0
35
.
0
11
.
64
266
.
0
e
J1809–1917
18
h
09
m
43
.
13
s
−
19
◦
17
′
38
.
2
′′
24
.
166
0
.
097
0
.
097
−
74
.
4
32
.
0
30
.
0
3
.
7
0
.
14
8
.
886
152
.
3
a
J1813-1749 bg
18
h
13
m
35
.
11
s
−
17
◦
49
′
57
.
57
′′
44
.
703
...
0
.
18
−
1290
.
0
...
510
.
0
0
.
0
0
.
0
...
92
.
08
d
J1813-1749 full
18
h
13
m
35
.
11
s
−
17
◦
49
′
57
.
57
′′
44
.
703
0
.
18
0
.
18
−
1290
.
0
520
.
0
510
.
0
0
.
0
0
.
0
266
.
4 2270
.
0
d
J1828–1101
18
h
28
m
18
.
85
s
−
11
◦
01
′
51
.
72
′′
27
.
754
0
.
11
0
.
11
−
57
.
0
23
.
0
23
.
0
4
.
23
9
.
4
7
.
481 1162
.
0
a
J1833–0827
18
h
33
m
40
.
26
s
−
08
◦
27
′
31
.
53
′′
23
.
449
0
.
094
0
.
094
−
25
.
2
11
.
0
10
.
0
−
0
.
463
0
.
082
2
.
873
242
.
2
e
J1838–0655
18
h
38
m
3
.
13
s
−
06
◦
55
′
33
.
4
′′
28
.
363
0
.
11
0
.
11
−
199
.
0
83
.
0
80
.
0
67
.
1
7
.
8
27
.
13
555
.
4
d
J1856+0245
18
h
56
m
50
.
91
s
+02
◦
45
′
53
.
17
′′
24
.
714
0
.
099
...
−
189
.
0
79
.
0
...
0
.
0
0
.
0
22
.
42
...
a
J1913+1011
19
h
13
m
20
.
34
s
+10
◦
11
′
22
.
97
′′
55
.
694
0
.
22
0
.
22
−
52
.
5
23
.
0
21
.
0
−
0
.
321
2
.
6
15
.
02 1951
.
0
a
J1925+1720
19
h
25
m
27
.
06
s
+17
◦
20
′
27
.
42
′′
26
.
434
0
.
11
0
.
11
−
36
.
6
15
.
0
15
.
0
3
.
93
6
.
3
4
.
538
469
.
2
a
J1928+1746
19
h
28
m
42
.
55
s
+17
◦
46
′
29
.
67
′′
29
.
097
0
.
12
0
.
12
−
55
.
7
23
.
0
22
.
0
0
.
0
0
.
0
7
.
845
68
.
41
a
J1935+2025
19
h
35
m
41
.
94
s
+20
◦
25
′
40
.
1
′′
24
.
955
0
.
092
0
.
1
−
189
.
0
79
.
0
76
.
0
95
.
0
3
.
4
20
.
99
581
.
3
a
J1952+3252
19
h
52
m
58
.
21
s
+32
◦
52
′
40
.
51
′′
50
.
587
0
.
2
0
.
2
−
74
.
9
32
.
0
30
.
0
2
.
92
0
.
012
18
.
61 3908
.
0
f
J2124–3358
21
h
24
m
43
.
84
s
−
33
◦
58
′
45
.
06
′′
405
.
588
0
.
91
1
.
6
−
0
.
0169
2
.
1
0
.
0068
0
.
0
0
.
0
5
.
595
48
.
06
g
J2229+6114
22
h
29
m
6
.
57
s
+61
◦
14
′
10
.
9
′′
38
.
709
0
.
15
0
.
16
−
590
.
0
240
.
0
260
.
0
1170
.
0
0
.
0
107
.
3 1021
.
0
f
Note
— This table lists frequency
f
and spin-down (
̇
f
,
̈
f
) parameters corresponding to twice the rotational parameters in the
pulsar ephemerides at the start of the O3 run (2019 April 1 15:00 UTC, MJD 58574.000), as discussed in Sec. 4.1. Right
ascension (RA) and Declination (DEC) are given in J2000 coordinates. In the column headings,
F
stands for the
F
-statistic
search discussed in Sec. 4.3 and 5v stands for the 5
n
-vector method discussed in Sec. 4.2. The suffixes “bg” and “ag” refer to
“before” and “after” glitch respectively, while “full” corresponds to the full run in the case of PSR J1813–1749.
(a) Jodrell Bank 42ft telescope and Lovell telescope (b) MeerKAT telescope (part of the MeerTime project) (Bailes et al.
2020) (c) University of Tasmania’s Mt. Pleasant Observatory 26m telescope (d) NICER (Gendreau 2016) (e) MOST (Bailes
et al. 2017; Jankowski et al. 2018; Lower et al. 2020) (f) CHIME (part of the CHIME pulsar project, Bandura et al. 2014;
Amiri et al. 2021) (g) Nan ̧cay Decimetric Radio Telescope. This table is available online in a machine-readable format (LIGO
Scientific Collaboration et al. 2022).
5
This procedure removed 0
.
7 days out of 246 days of H1
data, 2
.
2 days out of 256 days of L1 data, and 0
.
46 days
out of 251 days of V1 data.
The data are then cleaned to remove some known
monochromatic or quasi-monochromatic detector arti-
facts that can be effectively monitored by external sen-
sors, like signals injected for calibration, or contam-
ination from power mains (Davis et al. 2019; Viets
& Wade 2021; Sun et al. 2020, 2021; Acernese et al.
2022a,b). The H1 and L1 data near the 60 Hz power
mains frequency are also cleaned using a non-linear fil-
tering method outlined in Vajente et al. (2020).
The data from L1 and H1 are then “gated” to remove
further large transient artifacts, which removes an ad-
ditional 1
.
0 days of data from L1 and 0
.
62 days of data
from H1 (Zweizig & Riles 2021; Davis et al. 2021).
There are some differences in how each search pipeline
uses these data, as will also be discussed in detail in
Sec. 4.
The 5
n
-vector search creates Short Fourier
Databases (SFDBs) from the time domain data. SFDBs
consist of a collection of short-duration (1024 s long)
fast Fourier transforms (FFTs), overlapped by half and
Tukey-windowed with a window shape parameter of
α
= 0
.
001. During the construction of the SFDBs, an
extra time-domain cleaning (Astone et al. 2005 , in ad-
dition to those discussed in the previous paragraph) is
applied to identify delta-like time-domain disturbances.
The
F
-statistic and transient pipelines create 1800 s
Tukey-windowed short Fourier transforms (SFTs) with
a window shape parameter of
α
= 0
.
001, and then ex-
tract a narrow frequency range around the frequency
band for each pulsar.
Additional narrow lines with identified instrumental
origin are listed in Goetz et al. (2021) and Piccinni et al.
(2021); we use these for pipeline-level cleaning or post-
processing vetoes, as discussed in later sections.
4.
SIGNAL MODEL AND SEARCH METHODS
4.1.
Signal model
With all three pipelines, we search for quasi-
monochromatic GW signals with fixed intrinsic ampli-
tude at approximately twice the pulsar rotation fre-
quency,
f
≈
2
f
spin
, allowing for some deviation in the
narrowband search approach. The general frequency
evolution model for such GWs, following the standard
model for a pulsar’s
f
spin
as a function of time
t
′
at the
pulsar, is
f
(
t
′
) =
f
Taylor
(
t
′
) +
f
glitch
(
t
′
)
.
(3)
For most pulsars, no glitches have been observed during
O3, and hence only the first term for the long-time spin-
down evolution is relevant:
f
Taylor
(
t
′
) =
N
∑
k
=0
f
(
k
)
(
t
′
−
T
ref
)
k
k
!
,
(4)
with a reference time
T
ref
and up to
N
frequency deriva-
tives
f
(
k
)
. We will also use
̇
f
,
̈
f
and
...
f
for the first three
derivatives in this paper. The three search methods in
this paper have different maximum
N
they can handle,
as discussed in Appendix A.
For glitching pulsars, the additional term is
f
glitch
(
t
′
) = Θ(
t
′
−
T
gl
)
[
M
∑
k
=0
∆
f
(
k
)
gl
(
t
′
−
T
gl
)
k
k
!
+
δf
R
e
−
(
t
′
−
T
gl
)
/τ
R
]
(5)
(or multiple such terms for pulsars with multiple
glitches), where Θ(
x
) is the Heaviside step function,
T
gl
is the glitch time, ∆
f
(
k
)
gl
with
M
≤
N
are the perma-
nent jumps in the frequency and its derivatives,
δf
R
is
the part of the frequency jump that relaxes back, and
τ
R
is the relaxation time. Relaxation is not necessarily
observed for all glitches. The glitch term is typically a
small correction on top of
f
Taylor
(
t
′
).
Building on this frequency evolution model, the full
GW signal received by a detector is
h
(
t
) =
F
+
(
t
;
α,δ,ψ
)
h
+
(
t
) +
F
×
(
t
;
α,δ,ψ
)
h
×
(
t
)
,
(6)
where
F
+
,
×
(
t
;
α,δ,ψ
) is the detector response to + and
×
polarized GWs received at the detector at time
t
,
coming from right ascension
α
and declination
δ
with
polarization angle
ψ
(Jaranowski et al. 1998). Timing
corrections for translating from
t
′
to
t
are discussed be-
low.
The amplitude of each GW polarization can be written
as
h
+
(
t
) =
1
2
h
0
(1 + cos
2
ι
) cos
[
φ
0
+ Φ(
t
;
{
f
(
k
)
}
,α,δ
)
]
,
(7a)
h
×
(
t
) =
h
0
cos
ι
sin
[
φ
0
+ Φ(
t
;
{
f
(
k
)
}
,α,δ
)
]
,
(7b)
where
ι
is the inclination angle of the pulsar,
φ
0
is the
initial phase,
h
0
is an overall dimensionless strain am-
plitude as discussed below, and Φ(
t
;
{
f
(
k
)
}
,α,δ
) is the
integrated phase from the timing model of Eq. (3).
To convert from
t
′
at the source to
t
at the detec-
tor, several timing corrections typically need to be ad-
dressed: proper motion of the source, binary orbital mo-
tion, and the motion of the detectors around the solar
6
system barycenter. Proper motions between the pul-
sar and the solar system can usually be ignored (Co-
vas 2020) for these searches, and in this analysis we do
not cover sources in binary orbits as in,
e.g.
, Abbott
et al. (2021a,e, 2022); Ashok et al. (2021). However, the
sky position dependent correction between solar system
barycenter and detector frame is a crucial part of the
analysis methods described in the following, and is im-
plemented in all cases (Jaranowski et al. 1998).
For GWs from a non-axisymmetric star, the amplitude
h
0
can be written as
h
0
=
4
π
2
G
c
4
I
zz
f
2
d
,
(8)
where
d
is the distance to the NS,
f
is the GW fre-
quency, and
I
zz
is the principal moment of inertia. The
ellipticity,
, is given by
=
|
I
xx
−
I
yy
|
I
zz
.
(9)
In summary, our CW signal model depends on two
sets of parameters, the frequency-evolution parame-
ters
λ
=
{
f
(
k
)
,α,δ
}
, and the amplitude parameters
A
=
{
h
0
,
cos
ι,ψ,φ
0
}
. The signal model for the long-
duration transient search is also based on the one dis-
cussed here, but with an additional window function in
time, as will be discussed below in Sec. 4.4.
For both CW pipelines, we search over a narrow range
of GW frequencies and first spin-down terms centered
on their respective central estimates from the source
ephemeris:
f
= 2
f
spin
(1
±
κ
) and
̇
f
= 2
̇
f
spin
(1
±
κ
),
where we take
κ
= 2
×
10
−
3
. This value is intended to
allow for physical offsets between the frequency of EM-
observed pulsar rotation and the GW emission process,
as discussed in the introduction and in more detail in
Abbott et al. (2008, 2019b).
The frequency and spin-down changes due to glitches
offer evidence of how large of a physical offset we might
expect. While glitches with ∆
f
(1)
gl
/
̇
f
spin
&
2
×
10
−
3
are
observed occasionally, this value of
κ
encompasses most
glitches that have been observed (see,
e.g.
, Fuentes et al.
2017; Lower et al. 2021) and increasing this value signif-
icantly to encompass all observed glitch sizes would not
be computationally feasible. The full widths we allow for
both parameters, ∆
f
= 4
κf
spin
and ∆
̇
f
= 4
κ
̇
f
spin
, are
shown in Table 1 for each target we consider. Higher-
order frequency derivatives are handled differently by
each pipeline, which we discuss in the next two sub-
sections and in Appendix A.
4.2. 5
n
-vector narrowband pipeline
The 5
n
-vector narrowband pipeline
uses the 5
n
-vector
method of Astone et al. (2010, 2012) to search for CWs
in a narrow frequency and spin-down range around the
source ephemerides. This method searches for the char-
acteristic modulations imprinted by the Earth’s rotation
on a GW signal, splitting it into 5 harmonics. The search
is applied for the
n
detectors present in the network. For
more details on the implementation of the method we
use here, see Mastrogiovanni et al. (2017).
The starting point of the analysis is a collection of
1024 s long FFTs built as explained in the previous sec-
tion. For each target, we extract a frequency band large
enough to contain the Doppler modulation of every tem-
plate that we consider. When higher-order frequency
derivatives are used to fit the EM data, the correspond-
ing GW spin-down terms (including
̈
f
) are fixed to twice
the values provided in the source ephemerides.
For each target, we then correct for the Doppler mod-
ulation in the time domain using a non-uniform down-
sampling to 1 Hz. Using this time-series, we then apply
two matched filters (for the two CW polarizations) and
construct a detection statistic
S
combining coherently
the matched filter results from all detectors (Mastrogio-
vanni et al. 2017) for each template in the
f
and
̇
f
space.
The detection statistic is defined as
S
=
|
H
+
|
2
|
A
+
|
4
+
|
H
×
|
2
|
A
×
|
4
,
(10)
where
A
+
/
×
are
the
Fourier
components
of
F
+
/
×
(
t
; ˆ
n,ψ
= 0) and
H
+
/
×
are GW amplitude es-
timators built from the Fourier components of the data
and
A
+
/
×
. The template bank bin size in frequency is
given by the inverse of the time-series duration and the
spin-down bin size is given by its inverse squared.
The final step is to select the local maximum of the
detection statistic in every 10
−
4
Hz band over all spin-
down values considered. Within this set of points in
the parameter space, we select as outliers those with a
p
-value (defined as the tail probability of the detection
statistic noise-only distribution) below a 1% threshold,
taking into account the number of templates applied in
the
f
–
̇
f
space. The noise-only distribution of the detec-
tion statistic used to estimate the threshold for candi-
date selection is built using all the templates excluded
from the analysis in each 10
−
4
Hz band. The detec-
tion statistic value corresponding to the threshold for
candidate selection is extrapolated from the noise-only
distribution using an exponential fit of the tail.
Finally, to compute the 95% confidence level ULs
h
95%
0
on the CW amplitude in the case of no detection, we
perform software injections with simulated GW signals
in each 10
−
4
Hz sub-band to estimate the value of
h
0
at which we achieve 95% detection efficiency at a false-
alarm probability of 1%.
7
4.3.
Frequency-domain
F
-statistic pipeline
We also use the frequency-domain
F
-statistic pipeline
to search for CWs in a narrow frequency and spin-down
range around the source ephemerides. The
F
-statistic is
a matched filter which is maximized over the amplitude
parameters
A
of the quasi-monochromatic signals (Jara-
nowski et al. 1998; Cutler & Schutz 2005); see Tenorio
et al. (2021a) for a review covering its applications.
For each target we search over the range of GW fre-
quency
f
and spin-down
̇
f
, with matched-filter tem-
plates chosen using the multi-detector
F
-statistic met-
ric (Prix 2007; Wette et al. 2008), which we call through
the
LALSuite
(LIGO Scientific Collaboration 2021) pro-
gram
lalapps
Weave
(Wette et al. 2018b). If the second
spin derivative
̈
f
spin
is non-zero, we also search over the
range
̈
f
= 2(
̈
f
spin
±
3
σ
̈
f
spin
) where
σ
̈
f
spin
is the 1
σ
un-
certainty given in the source ephemeris. Higher order
frequency derivatives, up to
f
(4)
, are fixed to the EM-
measured values. Because of constraints on the pipeline
infrastructure, if the source ephemeris for a target in-
cludes derivatives above
f
(4)
, we refit the arrival times
with fewer frequency derivatives. This results in slightly
different ephemerides being used for this search, and for
the 5
n
-vector search. We discuss targets for which this
is the case in Appendix A.
We place templates with a maximum matched-filter
mismatch from a putative source of 0.02, which leads to
a very dense grid compared to the one used by the 5
n
-
vector search. This can be seen by comparing the final
two (non-reference) columns of Table 1. The analysis
uses a set of 1800 s Tukey-windowed SFTs over the full
data span for all three detectors, and produces a detec-
tion statistic, which we denote 2
F
, for each matched-
filter template. See Sec. 3 for a full discussion of data
processing and data quality cuts made before generation
of the SFTs.
Following Tenorio et al. (2022), we construct 10
4
randomly-chosen batches of templates from our search
results for each pulsar, and fit a Gumbel distribution
to the distribution of the maxima of these batches. We
can then propagate this Gumbel distribution based on
the total number of templates to estimate the tail prob-
ability (
p
-value) of the largest template across the full
search for that target. We outline our implementation
of this method more fully in Appendix B.1.
If any templates have a
p
-value of less than 1% and
corresponding large 2
F
values, we perform follow-up
studies of these templates. A list of frequencies of known
narrow spectral artifacts,
f
l
, and their nominal widths,
w
l
, are collated for all three detectors (Goetz et al. 2021;
Piccinni et al. 2021). If
|
f
l
−
f
large
|
< w
l
, where
f
large
is the frequency of the large-2
F
template, then we veto
that template. We also look at the single-detector 2
F
value calculated for each detector. If the individual 2
F
value for such a template in a single detector is larger
than the 2
F
value calculated using the whole network,
then we also veto that large-2
F
template (Keitel et al.
2014; Leaci 2015). Any remaining large-2
F
templates
with
p <
1% are then flagged for further follow-up, and
could be considered candidates for detection (subject to
further studies).
In the absence of a detection, we set ULs at 95% confi-
dence,
h
95%
0
, for each target using the software-injection
scheme set out in,
e.g.
, Abbott et al. (2007).
4.4.
Post-glitch transient search
The search for long-duration transient CW-like sig-
nals from glitching pulsars is motivated by the idea that
some of the excess energy liberated in a glitch could cor-
respond to transient changes in the quadrupole moment
of the star and be radiated away in GWs (van Eysden
& Melatos 2008; Bennett et al. 2010; Prix et al. 2011;
Melatos et al. 2015; Singh 2017; Yim & Jones 2020) dur-
ing the post-glitch relaxation phase.
To search for such signals, we use the transient
F
-
statistic method introduced by Prix et al. (2011) and
previously applied to LIGO O2 data in searches for GWs
from glitches in the Crab and Vela pulsars (Keitel et al.
2019). The idea is to search for long-duration transients
that are “CW-like” in the sense of following the standard
quasi-monochromatic CW signal model from Sec. 4.1
with only an additional window function
$
(
t
;
t
0
,τ
) in
time applied:
h
(
t
;
λ,
A
,
T
) =
$
(
t
;
t
0
,τ
)
h
(
t
;
λ,
A
)
,
(11)
where the transient parameters
T
consist of the window
shape, signal start time
t
0
, and a duration parameter
τ
. As in Keitel et al. (2019), we limit ourselves here to
the simplest case of rectangular windows, meaning the
signal is exactly a standard CW that starts at time
t
0
and is cut off at time
t
0
+
τ
, with no additional ampli-
tude evolution. Some form of amplitude decay would be
expected for a realistic signal linked to post-glitch re-
laxation (Yim & Jones 2020). As demonstrated in Prix
et al. (2011) and Keitel et al. (2019), the signal-to-noise
ratio (S/N) loss from using rectangular windows in a
search for exponentially decaying signals is mild, while
using exponentially windowed search templates would
mean a much higher computational cost (Keitel & Ash-
ton 2018).
The search uses the same underlying
F
-statistic
library
functions
in
LALSuite
(LIGO
Scien-
tific Collaboration 2021) as the CW search de-
scribed in Sec. 4.3, analyzing SFT data with the
8
10
2
10
3
Frequency [Hz]
10
−
26
10
−
25
10
−
24
h
0
H1 sensitivity
L1 sensitivity
V1 sensitivity
5
n
-vector 95% CL UL
F
-stat 95% CL UL
Spin-down limit
Figure 1.
The red solid, blue dashed, and purple dotted curves show the expected sensitivities for H1, L1, and V1 respectively,
using Eq. 1. The blue pentagons indicate the median 95% CL ULs from the 5
n
-vector search across all 10
−
4
Hz sub-bands for
each source. The black crosses indicate 95% CL ULs from the
F
-statistic search, which are set across the full search range for
each target. The orange triangles indicate the spin-down limit,
h
sd
, with error bars that reflect uncertainty in the distance to
each source. In a few cases the error bars are smaller than the size of the markers. We discuss and compare these limits in more
detail in Sec. 7. We do not show uncertainties on ULs here, although we discuss uncertainties due to the UL-setting method
as well as calibration uncertainties in Sec. 7. The data for this figure are available online (LIGO Scientific Collaboration et al.
2022).
lalapps
ComputeFstatistic
v2
program which was
used before,
e.g.
, in searches for CWs from supernova
remnants (Aasi et al. 2015c; Abbott et al. 2019c; Lind-
blom & Owen 2020). For each target, we perform a
search at fixed sky location over a grid in
f
(
k
)
param-
eters, with the number of spin-down terms depending
on the pulsar ephemerides, going up to at most a
...
f
term. At each grid point
λ
, the standard multi-detector
F
-statistic is calculated over the maximum duration
of interest. Intermediate per-SFT quantities from this
calculation, the so-called
F
-statistic atoms, are kept.
As described in Prix et al. (2011), partial sums of these
atoms are evaluated in a loop over a range of
{
t
0
,τ
}
.
The resulting matrix of
F
mn
=
F
(
t
0
m
,τ
n
) statistics
gives the likelihoods of transient CW-like signals in
each time range. We marginalize
F
mn
over uniform
priors on
t
0
and
τ
, obtaining the Bayes factor
B
tS
/
G
for transient CW-like signals against Gaussian noise as
our detection statistic. As demonstrated by Prix et al.
(2011),
B
tS
/
G
improves detection efficiency compared
to taking the maximum of
F
mn
, and as demonstrated
by Tenorio et al. (2022), it also produces cleaner back-
ground distributions on real data. The computational
cost of these searches scales linearly with the number of
λ
templates and with the product of the numbers of
t
0
and
τ
values, and is dominated by the partial summing
steps over the latter two parameters. See Prix et al.
(2011) and Keitel & Ashton (2018) for details.
As this is only the second search thus far for long-
transient GWs of this type (after Keitel et al. 2019), we
describe its practical setup in more detail in Sec. 6.1.
5.
CW SEARCH RESULTS
9
Table 2.
UL results on CW strain from the 5
n
-vector (5v) and
F
-statistic (
F
) pipelines.
Name
f h
95
0
(5v)
h
95
0
(
F
)
h
sd
h
95
0
/h
sd
h
95
0
/h
sd
(5v)
(
F
)
sd
d
σ
d
Hz
×
10
−
26
×
10
−
26
×
10
−
26
(5v)
(
F
)
×
10
−
5
×
10
−
5
×
10
−
5
kpc
kpc
J0534+2200 bg
59.24
...
8.79
143.0
...
0.061
...
4.6
76.0
2.00
a
0.5
J0534+2200 ag
59.24
5.2
6.09
143.0
0.036
0.043
2.7
3.2
76.0
2.00
a
0.5
J0711–6830
364.23
2.29
2.69
1.25
1.8
2.2
0.0017
0.002
0.00095
0.11
0.041
J0835–4510
22.37
36.2
39.5
341.0
0.11
0.12
19.0
21.0
180.0
0.28
b
0.02
J1101–6101
31.85
8.96
10.5
4.25
2.1
2.5
59.0
68.0
28.0
7.00
2.7
J1105–6107
31.64
8.94
10.1
17.1
0.52
0.59
20.3
23.0
38.0
2.36
0.92
J1809–1917
24.17
21.7
26.1
13.7
1.6
1.9
110.0
140.0
73.0
3.27
1.3
J1813–1749 full
44.70
4.61
5.85
21.9
0.21
0.27
10.0
13.0
64.0
6.20
c
2.4
J1813–1749 bg
44.70
...
9.61
21.9
...
0.44
...
21.0
64.0
6.20
c
2.4
J1828–1101
27.75
14.8
16.3
7.67
1.9
2.1
87.0
96.0
45.0
4.77
1.9
J1833–0827
23.45
22.9
29.3
5.88
3.9
5.0
180.0
230.0
46.0
4.50
1.8
J1838–0655
28.36
13.2
15.3
10.2
1.3
1.5
100.0
120.0
79.0
6.60
d
0.9
J1856+0245
19.35
20.5
...
11.2
1.8
...
200.0
...
110.0
6.32
2.46
J1913+1011
55.69
3.47
4.37
5.36
0.65
0.81
4.9
6.1
7.5
4.61
1.8
J1925+1720
26.43
14.3
17.6
5.93
2.4
3.0
98.0
120.0
41.0
5.06
2.0
J1928+1746
29.10
10.5
13.5
8.15
1.3
1.7
51.0
66.0
39.0
4.34
1.7
J1935+2025
24.95
19.2
23.7
15.3
1.3
1.5
130.0
170.0
110.0
4.60
1.8
J1952+3252
50.59
3.73
4.93
10.3
0.36
0.48
4.1
5.5
11.0
3.00
e
2.0
J2124–3358
405.59
2.4
2.9
0.374
6.4
7.7
0.0057
0.0068
0.00095 0.44
f
0.05
J2229+6114
38.71
5.78
7.23
33.1
0.17
0.22
11.0
14.0
63.0
3.00
g
2.0
Note
— We show the ULs on strain amplitude at 95% confidence for both pipelines,
h
95
0
, the spin-down limit for each target
h
sd
along with the implied limit on ellipticity
for each pipeline, and the limit on ellipticity implied by the spin-down
limit
sd
. Finally, we also show the ratio of the 95% confidence UL and the spin-down limit. We surpass the spin-
down limit for seven targets: PSR J0534+2200, PSR J0835-4510, PSR J1105–6107, PSR J1813–1749, PSR J1913+1011,
PSR J1952+3252, and PSR J2229+6114. We do not show uncertainty on ULs here, although we discuss uncertainty due
to the UL-setting method as well as calibration uncertainty in Sec. 7. Distance estimates
d
and their uncertainties
σ
d
are
from dispersion measures (with a fiducial uncertainty of 40% from Yao et al. 2017) unless noted with one of the following:
(a) Kaplan et al. (2008) (b) Dodson et al. (2003) (c) Camilo et al. (2021) (d) Gotthelf & Halpern (2008) (e) Verbiest
et al. (2012) (f) Reardon et al. (2021) (g) Halpern et al. (2001). This table and a machine-readable file are available
online (LIGO Scientific Collaboration et al. 2022).
As described below, we find no significant candi-
dates, although both the
F
-statistic and 5
n
-vector CW
searches find outliers requiring further follow up after
the vetoes in Sec. 4 were applied. Hence, we set obser-
vational ULs on the GW strain from each target, con-
straining GW emission below the spin-down limit on
seven of the eighteen target pulsars. All UL results are
shown in Fig. 1, and listed in Table 2. We give details of
the results from both pipelines, including outliers that
were followed up, in the rest of this section.
When discussing ULs in the following section, it
should be noted that physically meaningful constraints
on the GW energy emission are set when the ULs are
lower than the spin-down limit (i.e., when the limit
has been
surpassed
). Ellipticity constraints for pulsars
whose spin-down limit is not surpassed would imply a
|
̇
f
spin
|
higher than the one observed.
Finally, distance estimates used for the spin-down lim-
its are given in Table 2. These are either from the ATNF
pulsar catalogue (based on dispersion measures), or from
the literature. For the pulsar PSR J1813–1749, different
models of the electron density in the Galaxy yield dif-
ferent distance estimates (Camilo et al. 2021). We have
used the more optimistic estimate
d
= 6
.
2
±
2
.
4 kpc,
although a different model gives a more pessimistic es-
timate of
d
= 12 kpc.
5.1. 5
n
-vector narrowband pipeline