of 36
All-sky search for continuous gravitational waves from isolated neutron stars using
Advanced LIGO and Advanced Virgo O3 data
The LIGO Scientific Collaboration, The Virgo Collaboration, and The KAGRA Collaboration
(Dated: January 4, 2022)
We present results of an all-sky search for continuous gravitational waves which can be produced
by spinning neutron stars with an asymmetry around their rotation axis, using data from the third
observing run of the Advanced LIGO and Advanced Virgo detectors. Four different analysis methods
are used to search in a gravitational-wave frequency band from 10 to 2048 Hz and a first frequency
derivative from
10
8
to 10
9
Hz/s. No statistically-significant periodic gravitational-wave signal
is observed by any of the four searches. As a result, upper limits on the gravitational-wave strain
amplitude
h
0
are calculated. The best upper limits are obtained in the frequency range of 100 to 200
Hz and they are
1
.
1
×
10
25
at 95% confidence-level. The minimum upper limit of 1
.
10
×
10
25
is
achieved at a frequency 111.5 Hz. We also place constraints on the rates and abundances of nearby
planetary- and asteroid-mass primordial black holes that could give rise to continuous gravitational-
wave signals.
I. INTRODUCTION
The Advanced LIGO [1] and Advanced Virgo [2] de-
tectors have made numerous detections of gravitational
waves (GW), to date consisting of short-duration (tran-
sient) GW emitted during the inspirals and mergers of
compact binary systems of black holes (BH), neutron
stars (NS), [3, 4], as well as mixed NS-BH binaries [5].
Among still undiscovered types of GW radiation are long-
lasting, almost-monochromatic continuous waves (CW),
whose amplitudes and frequencies change much more
slowly compared to those of transient sources (on the
timescale of years rather than seconds). Astrophysically,
promising sources of CW are rotating, non-axisymmetric
NS, emitting GW at a frequency close to, or related to,
their spin frequency. Deviations from the symmetry (a
NS ‘deformation’) may be caused by fluid instabilities,
such as in the case of r-modes, or by elastic, thermal or
magnetic stresses in the crust and/or core of NS, and
may be acquired at various stages of stars’ isolated evo-
lution, or during an interaction with a companion in a
binary system (for recent reviews on sources of CW, see
e.g., [6–8]). Discovery of CW emitted by NS would allow
to probe their still mysterious interiors, study properties
of dense matter in conditions distinct from those occur-
ring in inspirals and mergers of binary NS systems, as
well as carry out additional tests of the theory of gravity
[9]. Due to intrinsically smaller GW amplitude of CW
in comparison to the already-detected transient sources,
searches for CW from rotating non-axisymmetric NS are
essentially limited to the Galaxy.
The search presented here is not limited to
gravitational-wave signals from deformed rotating neu-
tron stars. Another source of quasi-monochromatic, per-
sistent GWs are very light, planetary- and asteroid-mass,
inspiraling primordial black holes (PBHs), which could
comprise a fraction or the totality of dark matter [10].
Full author list given at the end of the article.
Such signals would arise from inspiraling PBHs whose
chirp masses are less than
O
(10
5
)
M
and whose GW
frequencies are less than
250 Hz, and would be indis-
tinguishable from those arising from non-axisymmetric
rotating NSs spinning up.
Recent detections of black holes made by the LIGO-
Virgo-KAGRA Collaboration have revived interest in
PBHs: low spin measurements and the rate inferences
are consistent with those expected for BHs that formed
in the early universe [11]. Existence of light PBHs is
well-motivated theoretically and experimentally: recent
detections of star and quasar microlensing events [12–
14] suggest compact objects or PBHs with masses be-
tween 10
6
and 10
5
M
could constitute a fraction of
dark matter of order
f
PBH
0
.
01, which is consistent
within the unified scenario for PBH formation presented
in [15], but greater than expected for free-floating (i.e.
not bound to an orbit) planets [16] (e.g. the hypothet-
ical Planet 9 could be a PBH with a mass of 10
6
M
that was captured by the solar system [17]). PBHs may
also collide with NS and be responsible for the origin of
NS-mass BHs, potentially detectable in the LIGO-Virgo-
KAGRA searches [18]. However, constraints arising from
such observations [10], even those that come from the
LIGO-Virgo merging rate inferences [19, 20] and stochas-
tic background searches [21, 22], rely on modelling as-
sumptions, and can be evaded if, for example, PBHs
formed in clusters [23–28]. It is therefore important to
develop complementary probes of these mass regimes to
test different PBH formation models [29, 30], which is
possible by searching for continuous GWs.
Searches for continuous waves are usually split in three
different domains:
targeted searches
look for signals from
known pulsars;
directed searches
look for signals from
known sky locations;
all-sky searches
look for signals from
unknown sources. All-sky searches for
a priori
unknown
CW sources have been carried out in the Advanced LIGO
and Advanced Virgo data previously [31–43]. A recent
review on pipelines for wide parameter-space searches can
be found in [44].
Here we report on results from an all-sky, broad fre-
arXiv:2201.00697v1 [gr-qc] 3 Jan 2022
2
quency range search using the most-sensitive data to
date, the LIGO-Virgo O3 observing run, employing four
different search pipelines: the
FrequencyHough
[45],
Sky-
Hough
[46],
Time-Domain
F
-statistic
[47, 48], and
SOAP
[49]. Each pipeline uses different data analysis meth-
ods and covers different regions of the frequency and fre-
quency time derivative parameter space, although there
exist overlaps between them (see Table I and Fig. 1 for
details). The search is performed for frequencies between
10 Hz and 2048 Hz and for a range of frequency time
derivative between -10
8
Hz/s and 10
9
Hz/s, covering
the whole sky. We note here that the search is generally-
agnostic to the type of the GW source, so the results are
not actually limited to signals from non-axisymmetric
rotating NS in our Galaxy.
A comprehensive multi-
stage analysis of the signal outliers obtained by the four
pipelines has not revealed any viable candidate for a con-
tinuous GW signal. However we improve the broad-range
frequency upper limits with respect to previous O1 and
O2 observing run and also with respect to the recent
analysis of the first half of the O3 run [39]. This is also
the first all-sky search for CW sources that uses the Ad-
vanced Virgo detector’s data.
The article is organized as follows: in Section II we
describe the O3 observing run and provide details about
the data used. Section III we present an overview of
the pipelines used in the search. Section IV, details of
the data-analysis pipelines are described. Section V, we
describe the results obtained by each pipeline, namely
the signal candidates and the sensitivity of the search
whereas Section VI contains a discussion of the astro-
physical implications of our results.
II. DATA SETS USED
The data set used in this analysis was the third ob-
serving run (O3) of the Advanced LIGO and Advanced
Virgo GW detectors [1, 2]. LIGO is made up of two
laser interferometers, both with 4 km long arms. One is
at the LIGO Livingston Observatory (L1) in Louisiana,
USA and the other is at the LIGO Hanford Observa-
tory (H1) in Washington, USA. Virgo (V1) consists of
one interferometer with 3 km arms located at European
Gravitational Observatory (EGO) in Cascina, Italy. The
O3 run took place between the 2019 April 1 and the
2020 March 27. The run was divided into two parts,
O3a and O3b, separated by one month commissioning
break that took place in October 2019. The duty fac-
tors for this run were
76%,
71%,
76% for L1, H1,
V1 respectively. The maximum uncertainties (68% confi-
dence interval) on the calibration of the LIGO data were
of 7%/11% in magnitude and 4 deg/9 deg in phase for
O3a/O3b data ([50, 51]). For Virgo, it amounted to 5%
in amplitude and 2 deg in phase, with the exception of
the band 46 - 51 Hz, for which the maximum uncertainty
was estimated as 40% in amplitude and 34 deg in phase
during O3b. For the smaller range 49.5 - 50.5 Hz, the
calibration was unreliable during the whole run [52].
III. OVERVIEW OF SEARCH PIPELINES
In this section we provide a broad overview of the four
pipelines used in the search. The three pipelines:
Fre-
quencyHough
,
SkyHough
, and
Time-Domain
F
-statistic
have been used before in several all-sky searches of the
LIGO data. The
SOAP
pipeline is a new pipeline ap-
plied for the first time to an all-sky search. It uses novel
algorithms.
SOAP
aims at a fast, preliminary search
of the data before more sensitive but much more time
consuming methods are applied (see [44] for a review on
pipelines for wide parameter-space searches). The indi-
vidual pipelines are described in more detail in the fol-
lowing section.
A. Signal model
The GW signal in the detector frame from an isolated,
asymmetric NS spinning around one of its principal axis
of inertia is given by [47]:
h
(
t
) =
h
0
[
F
+
(
t,α,δ,ψ
)
1 + cos
2
ι
2
cos
φ
(
t
)
+
F
×
(
t,α,δ,ψ
) cos
ι
sin
φ
(
t
)]
,
(1)
where
F
+
and
F
×
are the antenna patterns of the de-
tectors dependent on right ascension
α
, declination
δ
of
the source and polarization angle
ψ
,
h
0
is the amplitude
of the signal,
ι
is the angle between the total angular
momentum vector of the star and the direction from the
star to the Earth, and
φ
(
t
) is the phase of the signal. The
amplitude of the signal is given by:
h
0
=
4
π
2
G
c
4
I
zz
f
2
d
1
.
06
×
10
26
(

10
6
)
×
(
I
zz
10
38
kg m
2
)(
f
100 Hz
)
2
(
1 kpc
d
)
,
(2)
where
d
is the distance from the detector to the source,
f
is the GW frequency (assumed to be twice the rotation
frequency of the NS),

is the ellipticity or asymmetry of
the star, given by (
I
xx
I
yy
)
/I
zz
, and
I
zz
is the moment
of inertia of the star with respect to the principal axis
aligned with the rotation axis.
We assume that the phase evolution of the GW signal
can be approximated with a second order Taylor expan-
sion around a fiducial reference time
τ
r
:
φ
(
τ
) =
φ
o
+ 2
π
[
f
(
τ
τ
r
) +
̇
f
2!
(
τ
τ
r
)
2
]
,
(3)
where
φ
o
is an initial phase and
f
and
̇
f
are the frequency
and first frequency derivative at the reference time. The
3
10
16
5×10
12
2×10
11
10
10
10
9
10
8
FrequencyHough
SkyHough
TD Fstat
SOAP
10
20
65
100
200
350
750
1000
2048
f
GW
[
Hz
]
10
8
10
9
2×10
10
10
11
10
16
0
f
GW
[
Hz
/
s
]
FIG. 1. Frequency and frequency derivative search ranges
of the four pipelines: the
FrequencyHough
pipeline ranges
marked in grey,
SkyHough
in red,
Time-Domain
F
-statistic
in blue, and
SOAP
in magenta. See Table I for details.
relation between the time at the source
τ
and the time
at the detector
t
is given by:
τ
(
t
) =
t
+
~r
(
t
)
·
~n
c
+ ∆
E
S
,
(4)
where
~r
(
t
) is the position vector of the detector in the
Solar System Barycenter (SSB) frame, and
~n
is the unit
vector pointing to the NS; ∆
E
and ∆
S
are respec-
tively the relativistic Einstein and Shapiro time delays.
In standard equatorial coordinates with right ascension
α
and declination
δ
, the components of the unit vector
~n
are given by (cos
α
cos
δ,
sin
α
cos
δ,
sin
δ
).
B. Parameter space analyzed
All the four pipelines perform an all-sky search, how-
ever the frequency and frequency derivative ranges ana-
lyzed are different for each pipeline. The detailed ranges
analyzed by the four pipelines are summarized in Table I
and presented in Fig. 1. The
FrequencyHough
pipeline
analyzes a broad frequency range between 10 Hz and
2048 Hz and a broad frequency time derivative range be-
tween -10
8
Hz/s and 10
9
Hz/s. A very similar range
of
f
and
̇
f
is analyzed by
SOAP
pipeline. The
SkyHough
pipeline analyzes a narrower frequency range where the
detectors are most sensitive whereas
Time-Domain
F
-
statistic
pipeline analyzes
f
and
̇
f
ranges of the bulk of
the observed pulsar population (see Fig. 2 in Sect. IV C).
C. Detection statistics
As all-sky searches cover a large parameter space they
are computationally very expensive and it is computa-
tionally prohibitive to analyze coherently the data from
the full observing run using optimal matched-filtering.
As a result each of the pipelines developed for the anal-
ysis uses a semi-coherent method. Moreover to reduce
the computer memory and to parallelize the searches the
data are divided into narrow bands. Each analysis be-
gins with sets of
short Fourier transforms
(SFTs) that
span the observation period, with coherence times rang-
ing from 1024s to 8192s. The
FrequencyHough
,
SkyHough
and
SOAP
pipelines compute measures of strain power
directly from the SFTs and create detection statistics
by stacking those powers with corrections for frequency
evolution applied. The
FrequencyHough
and
SkyHough
pipelines use
Hough
transform to do the stacking whereas
SOAP
pipeline uses the
Viterbi
algorithm. The
Time-
Domain
F
-statistic
pipeline extracts band-limited 6-day
long time-domain data segments from the SFT sets and
applies frequency evolution corrections coherently to ob-
tain the
F
-statistic ([47]). Coincidences are then required
among multiple data segments with no stacking.
D. Outlier follow-up
All four pipelines perform a follow-up analysis of
the statistically significant candidates (outliers) obtained
during the search. All pipelines perform vetoing of the
outliers corresponding to narrow, instrumental artifacts
(lines) in the advanced LIGO detectors ([53]). Several
other consistency vetoes are also applied to eliminate
outliers. The
FrequencyHough
,
SkyHough
, and
Time-
Domain
F
-statistic
pipelines perform follow-up of the
candidates by processing the data with increasing long
coherence times whereas
SOAP
pipeline use
convolu-
tional neural networks
to do the post processing.
E. Upper limits
No periodic gravitational wave signals were observed
by any of the four pipelines and and all the pipelines ob-
tain upper limits on their strength. The three pipelines
SkyHough
,
Time-Domain
F
-statistic
and
SOAP
obtain
the upper limits by injections of the signals according to
the model given in Section III A above for an array of sig-
nal amplitudes
h
0
and randomly choosing the remaining
parameters. The
FrequencyHough
pipeline obtains upper
limits using an analytic formula (see Eq. 6) that depends
on the spectral density of the noise of the detector. The
formula was validated by a number of tests consisting of
injecting signals to the data.
4
Pipeline
Frequency [Hz] Frequency derivative [Hz/s]
FrequencyHough 10
2048
-10
8
10
9
SkyHough
65
350
-10
9
5
×
10
12
SOAP
40
1000
-10
9
10
9
1000
2000
-10
8
10
8
TD Fstat
20
200
-3
.
2
×
10
9
f/
100
0
200
750
-2
×
10
10
2
×
10
11
TABLE I. Frequency and frequency derivative search ranges of the four pipelines.
Band [Hz]
T
FFT
[s]
δf
[Hz]
δ
̇
f
[Hz/s]
10–128
8192
1
.
22
×
10
4
3
.
92
×
10
12
128–512
4096
2
.
44
×
10
4
7
.
83
×
10
12
512–1024
2048
4
.
88
×
10
4
1
.
57
×
10
11
1024–2048
1024
9
.
76
×
10
4
3
.
13
×
10
11
TABLE II. Properties of the FFTs used in the
Frequency-
Hough
pipeline. The time duration
T
FFT
refers to the length
in seconds of the data chunks on which the FFT is computed.
The frequency bin width is the inverse of the time duration,
while the spin-down bin width is computed as
δ
̇
f
=
δf/T
obs
,
where
T
obs
is the total run duration.
IV. DETAILS OF SEARCH METHODS
A.
FrequencyHough
The
FrequencyHough
pipeline is a semi-coherent proce-
dure in which interesting points (i.e. outliers) are selected
in the signal parameter space, and then are followed-up
in order to confirm or reject them. This method has
been used in several past all-sky searches of Virgo and
LIGO data [31, 34, 35, 54]. A detailed description of
the methodology can be found in [45]. In the following,
we briefly describe the main analysis steps and specific
choices used in the search.
Calibrated detector data are used to build “short du-
ration” and cleaned [55] Fast Fourier Transform (FFTs),
with duration
T
FFT
which depends on the frequency band
being considered, see Table II.
Next, local maxima are selected based on the square
root of the equalized power of the data
1
passing a dimen-
sionless threshold of Θ = 1
.
58. The collection of these
time-frequency peaks forms the so-called
peakmap
.
The peakmap is cleaned of the strongest disturbances
using a
line persistency
veto [45].
The time-frequency peaks of the peakmap are prop-
erly shifted, for each sky position
2
, to compensate
1
Computed as the ratio of the squared modulus of each FFT of
the data and an auto-regressive estimation of the average power
spectrum, see [55] for more details.
2
Over a suitable grid, which bin size depends on the frequency
and sky location.
the Doppler effect due to the detector motion [45].
The shifted peaks are then fed to the
FrequencyHough
algorithm [45], which transforms each peak to the
frequency/spin-down plane of the source. The frequency
and spin-down bins (which we will refer to as
coarse
bins
in the following) depend on the frequency band, as shown
in Table II, and are defined, respectively, as
δf
= 1
/T
FFT
and
δ
̇
f
=
δf/T
obs
, where
T
obs
is the total run duration.
In practice, the nominal frequency resolution has been
increased by a factor of 10 [45], as the
FrequencyHough
is not computationally bounded by the width of the fre-
quency bin. The algorithm, moreover, adaptively weights
any noise non-stationarity and the time-varying detector
response [56].
The whole analysis is split into tens of thousands of
independent jobs, each of which covers a small portion
of the parameter space. Moreover, for frequencies above
512 Hz a GPU-optimized implementation of the
Frequen-
cyHough
transform has been used [57].
The output of a
FrequencyHough
transform is a 2-
D histogram in the frequency/spin-down plane of the
source.
Outliers, that is significant points in this plane, are
selected by dividing each 1 Hz band of the corresponding
histogram into 20 intervals and taking, for each interval,
and for each sky location, the one or (in most cases) two
candidates with the highest histogram number count [45].
All the steps described so far are applied separately to
the data of each detector involved in the analysis.
As in past analyses [31, 34], candidates from each
detector are clustered and then coincident candidates
among the clusters of a pair of detectors are found us-
ing a distance metric
3
d
FH
built in the four-dimensional
parameter space of sky position (
λ, β
) (in ecliptic coordi-
nates), frequency
f
and spin-down
̇
f
. Pairs of candidates
with distance
d
FH
3 are considered coincident. In the
current O3 analysis, coincidences have been done only
3
The metric is defined as
d
FH
=
(
f
δf
)
2
+
(
̇
f
δ
̇
f
)
2
+
(
λ
δλ
)
2
+
(
β
δβ
)
2
,
f
, ∆
̇
f
, ∆
λ
, and ∆
β
are the differences, for each parameter,
among pairs of candidates of the two detectors, and
δf
,
δ
̇
f
,
δλ
,
and
δβ
are the corresponding bin widths.
5
among the two LIGO detectors for frequencies above 128
Hz, while also coincidences H1 - Virgo and L1 - Virgo
have been considered for frequencies below 128 Hz, where
the difference in sensitivity (especially in the very low fre-
quency band) is less pronounced.
Coincident candidates are ranked according to the
value of a statistic built using the distance and the
Fre-
quencyHough
histogram weighted number count of the
coincident candidates [45]. After the ranking, the eight
outliers in each 0.1 Hz band with the highest values of
the statistic are selected and subject to the follow-up.
1. Follow-up
The
FrequencyHough
follow-up runs on each outlier of
each coincident pair. It is based on the construction of a
new peakmap, over
±
3
coarse
bins around the frequency
of the outlier, with a longer
T
FFT
. This new peakmap is
built after the removal of the signal frequency variation
due to the Doppler effect for a source located at the out-
lier sky position.
A new refined grid on the sky is built around this point,
covering
±
3
coarse
bins, in order to take into account the
uncertainty on the outlier parameters. For each point of
this grid we remove the residual Doppler shift from the
peakmap by properly shifting the frequency peaks. Each
new corrected peakmap is the input for the
Frequency-
Hough
transform to explore the frequency and the spin-
down range of interest (
±
3
coarse
bins for the frequency
and the spin-down). The most significant peak among all
the
FrequencyHough
histograms, characterized by a set
of
refined
parameters, is selected and subject to further
post-processing steps.
First, the
significance
veto (V1) is applied. It con-
sists in building a new peakmap over 0
.
2 Hz around the
outlier refined frequency, after correcting the data with
its refined parameters. The corrected peakmap is then
projected on the frequency axis. Its frequency range is
divided in sub-bands, each covering
±
2 coarse frequency
bins. The maximum of the projection in the sub-band
containing the outlier is compared with the maxima se-
lected in the remaining
off-source
intervals. The outlier
is kept if it ranks as first or second for both detectors.
Second, a
noise line
veto (V2) is used, which discards out-
liers whose frequency, after the removal of the Doppler
and spin-down corrections, overlaps a band polluted by
known instrumental disturbances.
The
consistency
test (V3) discards pairs of coincident
outliers if their Critical Ratios (CRs), properly weighted
by the detector noise level, differ by more than a factor
of 5. The CR is defined as
CR
=
x
μ
σ
,
(5)
where
x
is the value of the peakmap projection in a given
frequency bin,
μ
is the average value and
σ
the standard
deviation of the peakmap projection.
The
distance
veto (V4) consists in removing pairs of coin-
cident outliers with distance
d
FH
>
6 after the follow-up.
Finally, outliers with distance
d
FH
<
3 from
hardware
injections
are also vetoed (V5).
Outliers which survive all these vetoes are scrutinized
more deeply, by applying a further follow-up step, based
on the same procedures just described, but further in-
creasing the segment duration
T
FFT
.
2. Parameter space
The
FrequencyHough
search covers the frequency range
[10, 2048] Hz, a spin-down range between -10
8
Hz/s to
10
9
Hz/s and the whole sky. The frequency and spin-
down resolutions are given in Tab. II. The sky resolution,
on the other hand, is a function of the frequency and of
the sky position and is defined in such a way that for two
nearby sky cells the maximum frequency variation, due
to the Doppler effect, is within one frequency bin, see [45]
for more details.
3. Upper limits
Upper limits are computed for every 1 Hz sub-band
in the range of 20–2048 Hz
4
, considering only the LIGO
detectors, as Virgo sensitivity is worse for most of the
analyzed frequency band. First, for each detector we use
the analytical relation [45]
h
UL,
95%
4
.
97
N
1
/
4
S
n
(
f
)
T
FFT
CR
max
+ 1
.
6449
,
(6)
where
N
=
T
obs
/T
FFT
,
S
n
(
f
) is the detector average
noise power spectrum and
CR
max
is the maximum out-
lier CR
5
, in the given 1 Hz band. For each 1 Hz band,
the final upper limit is the worse among those computed
separately for Hanford and Livingston.
As verified through a detailed comparison based on
LIGO and Virgo O2 and O3 data, this procedure pro-
duces conservative upper limits with respect to those ob-
tained through the injection of simulated signals, which
is computationally much heavier [58].
Moreover, it has been shown that the upper limits ob-
tained through injections are always above those based
on Eq. 6 when the minimum CR in each 1 Hz sub-band is
used. The two curves based, respectively, on the highest
and the smallest CR delimit a region containing both a
4
Although the search starts at 10 Hz, we decided to compute
upper limit starting from 20 Hz, due to the unreliable calibration
at lower frequency.
5
Defined by Eq. 5 and where in this case the various quantities
are computed over the Frequency-Hough map
6
Parameter Resolution
δf
1
.
4
×
10
4
Hz
δ
̇
f
5
×
10
12
Hz/s
δθ
0
.
69 Hz
/f
TABLE III. Parameter-space resolutions employed by the
SkyHough
pipeline.
more stringent upper limit estimate and the search sen-
sitivity estimate, that is the minimum strain of a de-
tectable signal. Any astrophysical implication of our re-
sults, discussed in Sec. V will be always based on the
most conservative estimate.
B.
SkyHough
SkyHough
[46, 59] is a semicoherent pipeline based on
the Hough transform to look for CW signals from isolated
neutron stars. Several versions of this pipeline have been
used throughout the initial [60, 61] and advanced [31,
32] detector era, as well as to look for different kinds of
signals such as CW from neutron stars in binary systems
[40, 41, 62] or long-duration GW transients [63]. The
current implementation of
SkyHough
closely follows that
of [32] and includes an improved suite of post-processing
and follow-up stages [64–66].
1. Parameter space
The
SkyHough
pipeline searches over the standard four
parameters describing a CW signal from isolated NS: fre-
quency
f
, spin-down
̇
f
and sky position, parametrized
using equatorial coordinates
α,δ
.
Parameter-space resolutions are given in [46]
δf
=
1
T
SFT
, δ
̇
f
=
δf
T
obs
, δθ
=
c/v
T
SFT
P
f
f
,
(7)
where
θ
represents either of the sky angles,
v/c
'
10
4
represents the average detector velocity as a fraction of
the speed of light, and the pixel factor P
f
= 2 is a tun-
able overresolution parameter. Table III summarizes the
numerical values employed in this search.
The
SkyHough
all-sky search covers the most sensitive
frequency band of the advanced LIGO detectors, between
65 Hz and 350 Hz. This band is further sub-divided
into ∆
f
= 0
.
025 Hz sub-bands, resulting in a total of
11400 frequency bands. Spin-down values are covered
from
1
×
10
9
Hz
/
s to 5
×
10
12
Hz
/
s, which include
typical spin-up values associated to CW emission from
the evaporation of boson clouds around black holes [67].
2. Description of the search
The first stage of the
SkyHough
pipeline performs
a multi-detector search using H1 and L1 SFTs with
T
SFT
= 7200
s
. Each 0
.
025 Hz sub-band is analyzed sep-
arately using the same two step strategy as in [32, 41]:
parameter-space is efficiently analyzed using
SkyHough
’s
look-up table approach; the top 0.1% most significant
candidates are further analyzed using a more sensitive
statistic. The result for each frequency sub-band is a
toplist containing the 10
5
most significant candidates
across the sky and spin-down parameter-space.
Each toplist is then clustered using a novel approach
presented in [64] and firstly applied in [41]. A parameter-
space distance is defined using the average mismatch
in frequency evolution between two different parameter-
space templates
d
(
~
λ,
~
λ
) =
T
SFT
N
SFT
N
SFT
α
=0
f
(
t
α
;
~
λ
)
f
(
t
α
;
~
λ
)
,
(8)
where
f
(
t
;
~
λ
) is defined as
f
(
t
;
~
λ
) =
[
f
+ (
t
t
ref
)
·
̇
f
]
·
[
1 +
~v
(
t
)
·
~n
c
]
(9)
and
~
λ
=
{
f,
̇
f,α,δ
}
refers to the phase-evolution param-
eters of the template.
Clusters are constructed by pairing together templates
in consecutive frequency bins such that
d
(
~
λ,
~
λ
)
1.
Each cluster is characterized by its most significant el-
ement (the
loudest
element). From each 0
.
025 Hz sub-
band, we retrieve the forty most significant clusters for
further analysis. This results in a total of 456000 candi-
dates to follow-up.
The loudest cluster elements are first sieved through
the
line veto
, a standard tool to discard clear instrumen-
tal artifacts using the list of known, narrow, instrumen-
tal artifacts (lines) in the advanced LIGO detectors [53]:
If the instantaneous frequency of a candidate overlaps
with a frequency band containing an instrumental line of
known origin, the candidate is ascribed an instrumental
origin and consequently ruled out.
Surviving candidates are then followed-up using
PyFstat
, a Python package implementing a Markov-
chain Monte Carlo (MCMC) search for CW signals
[65, 68]. The follow-up uses the
F
-statistic as a (log)
Bayes factor to sample the posterior probability distribu-
tion of the phase-evolution parameters around a certain
parameter-space region
P(
~
λ
|
x
)
e
F
(
~
λ
;
x
)
·
P(
~
λ
)
,
(10)
where P(
~
λ
) represents the prior probability distribution
of the phase-evolution parameters. The
F
-statistic, as
opposed to the
SkyHough
number count, allows us to use
longer coherence times, increasing the sensitivity of the
follow-up with respect to the main search stage.
7
Stage
0
1
2
3
4
5
N
seg
660
330
92
24
4
1
T
coh
0.5 day 1 day 4 days 15 days 90 days 360 days
TABLE IV. Coherence-time configuration of the multi-stage
follow-up employed by the
SkyHough
pipeline. The data
stream is divided into a fix number of segments of the same
length; the reported coherence time is an approximate value
obtained by dividing the observation time by the number of
segments at each stage.
As initially described in [68], the effectiveness of an
MCMC follow-up is tied to the number of templates cov-
ered by the initial prior volume, suggesting a hierarchi-
cal approach: coherence time should be increased follow-
ing a ladder so that the follow-up is able to converge
to the true signal parameters at each stage. We follow
the proposal in [66] and compute a coherence-time lad-
der using
N
= 10
3
(see Eq. (31) of [68]) starting from
T
coh
= 1 day including an initial stage of
T
coh
= 0
.
5 days.
The resulting configuration is collected in Table IV.
The first follow-up stage is similar to that employed
in [40, 41]: an MCMC search around the loudest candi-
date of the selected clusters is performed using a coher-
ence time of
T
coh
= 0
.
5 days. Uniform priors containing
4 parameter-space bins in each dimension are centered
around the loudest candidate. A threshold is calibrated
using an injection campaign: any candidate whose loud-
est 2
F
value over the MCMC run is lower than 2
F
= 3450
is deemed inconsistent with CW signal.
The second follow-up stage is a variation of the method
described in [66], previously applied to [69, 70]. For each
outlier surviving the initial follow-up stage (stage 0 in
Table IV), we construct a Gaussian prior using the me-
dian and inter-quartile range of the posterior samples
and run the next-stage MCMC follow-up. The resulting
maximum 2
F
is then compared to the expected 2
F
in-
ferred from the previous MCMC follow-up stage. Highly-
discrepant candidates are deemed inconsistent with a
CW signal and hence discarded.
Given an MCMC stage using
ˆ
N
segments from which
a value of 2
ˆ
F
is recovered, the distribution of 2
F
values
using
N
segments is well approximated by
P(2
F|
N,
2
ˆ
F
,
ˆ
N
) = Gauss(2
F
;
μ,σ
)
,
(11)
where
μ
=
ρ
2
0
+ 4
N ,
(12)
σ
2
= 8
·
(
N
+
ˆ
N
+
ρ
2
0
)
,
(13)
and
ρ
2
0
= 2
ˆ
F−
4
ˆ
N
is a proxy for the (squared) SNR [71].
Equation (11) is exact in the limit of
N,
ˆ
N

1 or
ρ
2
0

1. In this search, however, we calibrate a bracket
on (2
F−
μ
)
for each follow-up stage using an injection
campaign, shown in table V. Candidates outside of the
bracket are deemed inconsistent with a CW signal.
Comparing stages (2
F −
μ
)
bracket
Stage 0 v.s. Stage 1
(-1.79, 1.69)
Stage 1 v.s. Stage 2
(-1.47, 1.35)
Stage 2 v.s. Stage 3
(-0.94, 0.80)
Stage 3 v.s. Stage 4
(-0.63, 0.42)
Stage 4 v.s. Stage 5
(-0.34, 0.11)
TABLE V. 2
F
consistency brackets employed in the multi-
stage follow-up of the
SkyHough
pipeline. Brackets were com-
puted using a campaign of 500 software-injected signals rep-
resenting an isotropic population of uniformly sky-distributed
NS at 150 representative frequency bands with an amplitude
corresponding to the
h
95%
0
sensitivity estimation. The implied
false dismissal probability is
.
1
/
(150
×
500)
'
1
.
3
×
10
5
.
Stages correspond to those described in Table IV.
Any surviving candidates are subject to manual in-
spection in search for obvious instrumental causes such
as hardware-injected artificial signals or narrow instru-
mental artifacts.
C.
Time-Domain
F
-statistic
The
Time-Domain
F
-statistic
search method has been
applied to an all-sky search of VSR1 data [48] and all-sky
searches of the LIGO O1 and O2 data [31, 32, 34]. The
main tool of the pipeline is the
F
-statistic [47] with which
one can coherently search the data over a reduced param-
eter space consisting of signal frequency, its derivatives,
and the sky position of the source. However, a coherent
all-sky search over the long data set like the whole data of
O3 run is computationally prohibitive. Thus the data are
divided into shorter time domain segments. Moreover, to
reduce the computer memory required to do the search,
the data are divided into narrow-band segments that are
analyzed separately. As a result the
Time-Domain
F
-
statistic
pipeline consists of two parts. The first part
is the coherent search of narrowband, time-domain seg-
ments. The second part is the search for coincidences
among the parameters of the candidates obtained from
the coherent search of all the time domain segments.
The algorithms to calculate the
F
-statistic in the co-
herent search are described in Sec. 6.2 of [48]. The time
series is divided into segments, called frames, of six side-
real days long each. Moreover the data are divided into
sub-bands of 0.25 Hz overlapped by 0.025 Hz. The O3
data has a number of non-science data segments. The
values of these bad data are set to zero. For our analysis,
we choose only segments that have a fraction of bad data
less than 60% both in H1 and L1 data and there is an
overlap of more than 50% between the data in the two de-
tectors. This requirement results in forty-one 6-day-long
data segments for each sub-band. For the search we use a
four-dimensional grid of templates (parameterized by fre-
quency, spin down rate, and two more parameters related
to the position of the source in the sky) constructed in
Sec. 4 of [72] with grid’s minimal match parameter MM
8
chosen to be
3
/
2. This choice of the grid spacing led
to the following resolution for the four parameters of the
space that we search
f
'
1
.
9
×
10
6
Hz
,
(14a)
̇
f
'
1
.
1
×
10
11
Hz/s
,
(14b)
α
'
7
.
4
×
10
2
(
100 Hz
f
)
rad
,
(14c)
δ
'
1
.
5
×
10
2
(
100 Hz
f
)
rad
.
(14d)
We set a fixed threshold of 15.5 for the
F
-statistic and
record the parameters of all threshold crossings, together
with the corresponding values of the
F
-statistic. In the
second stage of the analysis we use exactly the same co-
incidence search algorithm as in the analysis of VSR1
data and described in detail in Sec. 8 of [48] with only
one change. We use a different coincidence cell from that
described in [48]. In [48] the coincidence cell was con-
structed from Taylor expansion of the autocorrelation
function of the
F
-statistic. In the search performed here
the chosen coincidence cell is a suitably scaled grid cell
used in the coherent part of the pipeline. We scale the
four dimensions of the grid cell by different factors given
by [16 8 2 2] corresponding to frequency, spin down rate
(frequency derivative), and two more parameters related
to the position of the source in the sky respectively. This
choice of scaling gives optimal sensitivity of the search.
We search for coincidences in each of the bands analyzed.
Before identifying coincidences we veto candidate signals
overlapping with the instrumental lines identified by in-
dependent analysis of the detector data. To estimate
the significance of a given coincidence, we use the for-
mula for the false alarm probability derived in the ap-
pendix of [48]. Sufficiently significant coincidences are
called outliers and are subject to a further investigation.
1. Parameter space
Our
Time-Domain
F
-statistic
analysis is a search over
a 4-dimensional space consisting of four parameters: fre-
quency, spin-down rate and sky position. As we search
over the whole sky the search is very computationally in-
tensive. Given that our computing resources are limited,
to achieve a satisfactory sensitivity we have restricted
the range of frequency and spin-down rates analyzed to
cover the frequency and spin-down ranges of the bulk of
the observed pulsars. Thus we have searched the gravita-
tional frequency band from 20 Hz to 750 Hz. The lower
frequency of 20 Hz is chosen due to the low sensitivity of
the interferometers below 20 Hz. In the frequency 20 Hz
to 130 Hz range, assuming that the GW frequency is
twice the spin frequency, we cover young and energetic
pulsars, such as Crab and Vela. In the frequency range
from 80 Hz to 160 Hz we can expect GW signal due to
r-mode instabilities [73, 74]. In the frequency range from
160 Hz to 750 Hz we can expect signals from most of the
recycled millisecond pulsars, see Fig. 3 of [75].
For the GW frequency derivative
̇
f
we have chosen
a frequency dependent range. Namely, for frequencies
less than 200 Hz we have chosen
̇
f
to be in the range
[
f/τ
min
,
0], where
τ
min
is a limit on pulsar’s charac-
teristic age, and we have taken
τ
min
= 1000 yr. For
frequencies greater than 200 Hz we have chosen a fixed
range for the spin-down rate. As a result, the following
ranges of
̇
f
were searched in our analysis:
0
>
̇
f >
3
.
2
×
10
9
f
100 Hz
Hz/s
,
for
f <
200 Hz
,
(15a)
2
×
10
11
Hz/s
>
̇
f >
2
×
10
10
Hz/s
,
for
f >
200 Hz
.
(15b)
In Fig. 2 we plot GW frequency derivatives against
GW frequencies (assuming the GW frequency is twice
the spin frequency of the pulsar) for the observed pulsars
from the ATNF catalogue [76]. We show the range of the
GW frequency derivative selected in our search, and one
can see that the expected frequency derivatives of the
observed pulsars are well within this range. Note, finally,
that we have made the conservative choice of including
positive values of the frequency derivative (‘spin-up’), in
order to search as wide a range as possible. In most
cases, however, the pulsars that appear to spin-up are
in globular clusters, for which the local forces make the
measurement unreliable [77].
2. Sensitivity of the search
In order to assess the sensitivity of the
F
-statistic
search, we set upper limits on the intrinsic GW amplitude
h
0
in each 0.25 Hz bands. To do so, we generate signals
for an array of 8 amplitudes
h
0
and for randomly selected
sky positions (samples drawn uniformly from the sphere).
For each amplitude, we generate 100 signals with
f
,
̇
f
,
the polarization angle
ψ
and cosine of the inclination an-
gle
ι
are chosen from uniform random distributions in
their respective ranges. The signals are added to the
real data segments, and searches are performed with the
same grids and search set-up as for the real data search,
in the neighbourhood of injected signal parameters. We
search
±
6 grid points for
̇
f
and
±
1 grid points for the
sky positions away from the true values of the signal’s
parameters. We consider a signal
detected
if
coinci-
dence multiplicity for the injected signal is higher than
the highest signal multiplicity in a given sub-band and in
a given hemisphere in the real data search
. The detec-
tion efficiency is the fraction of recovered signals. We
estimate the
h
95%
0
, i.e., 95% confidence upper limit on
9
20
100
750
f
GW
[
Hz
]
10
18
10
17
10
16
10
15
10
14
10
13
10
12
10
11
10
10
10
9
10
8
10
7
|
f
GW
| [
Hz
/
s
]
J0534+2200 (Crab)
J0835
4510 (Vela)
J0537
6910
J0058
7218
FIG. 2. Frequency time derivative for tentative emission of
GWs (
̇
f
2
̇
f
rot
) as a function of the frequency of emitted
GWs (
f
2
f
rot
), where
f
rot
and
̇
f
rot
are rotational fre-
quency and frequency time derivative for known pulsars, ob-
tained from the Australia Telescope National Facility (ATNF)
database [76]. The vertical axis shows the absolute value for
both negative values of the frequency time derivative (“spin-
down”, blue dots) and positive values (“spinup”, red plus
symbols). Blue dashed lines represent spin-down limits used
in the
Time-Domain
F
-statistic
search: for
f <
200 Hz, 0
>
̇
f >
f/τ
min
, where
τ
min
= 1000 yr denotes a limit on pul-
sar’s characteristic age; for
f >
200 Hz,
̇
f >
2
×
10
10
Hz/s.
For
f >
200 Hz in the case of spinning-up objects, in the
F
-statistic search we admit a positive range of values to
̇
f <
2
×
10
11
Hz/s. The boundary of this range is marked
by a red continuous line.
the GW amplitude
h
0
, by fitting
6
a sigmoid function to
a range of detection efficiencies
E
as a function of in-
jected amplitudes
h
0
,
E
(
h
0
) =
(
1 +
e
k
(
x
0
h
0
)
)
1
, with
k
and
x
0
being the parameters of the fit. Figure 3 presents
an example fit to the simulated data with 1
σ
errors on
the
h
95%
0
estimate marked in red.
D.
SOAP
SOAP
[49] is a fast, model-agnostic search for long
duration signals based on the Viterbi algorithm [81]. It
is intended as both a rapid initial search for isolated NSs,
6
For the
h
95%
0
fitting procedure, we use the
python 3
[78]
scipy-optimize
[79]
curve
fit
package, implementing the
Levenberg-Marquardt least squares algorithm, to obtain the best
fitted parameters,
x
0
and
k
, to the sigmoid function. Errors of
parameters
δx
0
and
δk
are obtained from the covariance matrix
and used to calculate the standard deviation
σ
e
of the detection
efficiency as a function of
h
0
i.e., the confidence bands around the
central values of the fit. In practice, we use the
uncertainties
package [80] to obtain the
±
1
σ
standard deviation on the
h
0
value.
0.5
1.0
1.5
2.0
2.5
3.0
h
o
1e
25
0.0
0.2
0.4
0.6
0.8
1.0
detection efficiency E [%]
data
fit
1 fit error at 95% UL
E
efficiency errors
FIG. 3. Example sigmoid function fit (green solid line) to the
injected data efficiencies (blue dots), representing the detec-
tion efficiency
E
as a function of injected GW amplitude
h
0
used in
Time-Domain
F
-statistic
search. Pale red and blue
curves mark the 1
σ
confidence band obtained from the un-
certainty of the fit. Red error bar marks the
±
1
σ
standard
deviation on the
h
95%
0
value, corresponding to the efficiency
of 0.95 (indicated by the horizontal dashed gray line). Verti-
cal errors for each efficiency represent 1
σ
standard binomial
errors related to detection rate,
σ
E
=
E
(1
E
)
/N
i
, where
E
is the efficiency and
N
i
= 100 is the number of injections
for each GW amplitude. The data shown relates to the sub-
band with frequency of the lower edge of the band equal to
725.95 Hz.
quickly providing candidates for other search methods
to investigate further, as well as a method to identify
long duration signals which may not follow the standard
Continuous Wave
(CW) frequency evolution. In its most
simple form
SOAP
analyzes a spectrogram to find the
continuous time-frequency track which gives the highest
sum of fast Fourier transform power. If there is a signal
present within the data then this track is the most likely
to correspond to that signal. The search pipeline consists
of three main stages, the initial
SOAP
search [49], the
post processing step using convolutional neural networks
[82] and a parameter estimation stage.
1. Data preparation
The data used for this search starts as calibrated de-
tector data which is used to create a set of fast Fourier
transforms (FFTs) with a coherence time of 1800 s. The
power spectrum of these FFTs are then summed over
one day, i.e. every 48 FFTs. Assuming that the signal
remains within a single bin over the day, this averages
out the antenna pattern modulation and increases the
SNR in a given frequency bin. As the frequency of a
CW signal increases, the magnitude of the daily Doppler
modulation also increases, therefore the assumption that
a signal remains in a single frequency bin within one day
no longer holds. Therefore, the analysis is split into 4
separate bands (40-500 Hz, 500-1000 Hz, 1000-1500 Hz,
1500-2000 Hz) where for each band the Doppler modula-
tions are accounted for by taking the sum of the power
in adjacent frequency bins. For the bands starting at