Draft version November 29, 2021
Typeset using L
A
T
E
X
twocolumn
style in AASTeX631
Searches for Gravitational Waves from Known Pulsars at Two Harmonics in the Second and Third
LIGO-Virgo Observing Runs
The LIGO Scientific Collaboration, the Virgo Collaboration, and the KAGRA Collaboration
(Dated: November 29, 2021)
ABSTRACT
We present a targeted search for continuous gravitational waves (GWs) from 236 pulsars using data
from the third observing run of LIGO and Virgo (O3) combined with data from the second observing
run (O2). Searches were for emission from the
l
=
m
= 2 mass quadrupole mode with a frequency
at only twice the pulsar rotation frequency (single harmonic) and the
l
= 2
,m
= 1
,
2 modes with a
frequency of both once and twice the rotation frequency (dual harmonic). No evidence of GWs was
found so we present 95% credible upper limits on the strain amplitudes
h
0
for the single harmonic
search along with limits on the pulsars’ mass quadrupole moments
Q
22
and ellipticities
ε
. Of the
pulsars studied, 23 have strain amplitudes that are lower than the limits calculated from their elec-
tromagnetically measured spin-down rates. These pulsars include the millisecond pulsars J0437
−
4715
and J0711
−
6830 which have spin-down ratios of 0.87 and 0.57 respectively. For nine pulsars, their
spin-down limits have been surpassed for the first time. For the Crab and Vela pulsars our limits are
factors of
∼
100 and
∼
20 more constraining than their spin-down limits, respectively. For the dual
harmonic searches, new limits are placed on the strain amplitudes
C
21
and
C
22
. For 23 pulsars we
also present limits on the emission amplitude assuming dipole radiation as predicted by Brans-Dicke
theory.
1.
INTRODUCTION
To date, the LIGO and Virgo observatories have made
detections of numerous sources of gravitational radia-
tion. These detections have been of transient grav-
itational waves (GWs) from the inspiral and subse-
quent mergers of compact binary objects including bi-
nary black holes and binary neutron stars (Abbott
et al. 2021a). Recently, the list of observed events ex-
panded to include neutron star-black hole binaries (Ab-
bott et al. 2021b). There remain other types of GW
sources that are yet to be observed such as continuous
GW (CW) sources. Unlike transients, CW signals are
almost monochromatic, with their amplitude and fre-
quency changing very slowly over year-long timescales.
The mass quadrupoles of these sources, such as deformed
neutron stars, are expected to be far smaller than those
involved in compact binaries and therefore only local
galactic sources are likely to produce detectable signals.
Likely candidates for producing such signals are neu-
tron stars spinning with some non-axisymmetric defor-
mation (Zimmermann & Szedenits 1979). This may be
a solid deformation such as mountains on the crust pro-
duced during cooling (Ushomirsky et al. 2000), from bi-
nary accretion (Gittins & Andersson 2021) or due to
strong magnetic fields (Bonazzola & Gourgoulhon 1996;
Cutler 2002). GW radiation can also be caused by
fluid modes of oscillation beneath the crust such as r-
modes (Andersson 1998; Friedman & Morsink 1998). By
detecting CWs, light can be shed on the structure of
the star. Additionally, detections of such GWs can be
used to test general relativity via the constraint of non-
standard GW polarization (Isi et al. 2017; Abbott et al.
2019a). A more thorough discussion of various methods
of GW emission from neutron stars can be found in Riles
(2017); Glampedakis & Gualtieri (2018).
The structure of this paper is as follows. Section 1.1
outlines the types of CW searches. Section 1.2 describes
the types of signal models used in this analysis. Sec-
tion 2 describes the search methods used. Section 3
covers both the GW and EM data used. We present our
results in Section 4 with conclusions in Section 5.
1.1.
Continuous-wave searches
There are broadly three types of CW searches. Tar-
geted searches look for signals from known pulsars for
which their rotational phases can be accurately deter-
arXiv:2111.13106v1 [astro-ph.HE] 25 Nov 2021
2
mined from electromagnetic (EM) observations (e.g.,
Abbott et al. 2017a, 2019b; Nieder et al. 2019, 2020;
Abbott et al. 2020, 2021c). This simplifies the search
as the EM observations can be used to derive a timing
solution and it is assumed that the GW phase evolu-
tion is locked to the EM evolution. This means the
search is over a small parameter space, generally lim-
ited to the unknown signal amplitude and orientation
of the source, which allows a more sensitive search than
other methods. In some targeted searches, the assump-
tion that the GW evolution follows the EM evolution is
relaxed and the search is performed in a narrow band
around the expected frequency and spin-down rate (Ab-
bott et al. 2017b, 2019c). In this case, the search is more
computationally expensive due to the larger parameter
range being searched and slightly less sensitive because
of a higher trials factor. To overcome this, narrow-band
searches often look at fewer targets. Directed searches
look for signals from small sky regions that are believed
to have a high probability of containing a neutron star,
such as supernova remnants (e.g., Abbott et al. 2019d;
Ming et al. 2019; Lindblom & Owen 2020; Papa et al.
2020; Abbott et al. 2021d). As the timing solution can-
not be derived from EM observations, a wide range of ro-
tational parameters must be searched. All-sky searches
look for signals in all sky directions and over a wide
range of rotational parameters (e.g., Abbott et al. 2018,
2019e; Dergachev & Papa 2020; Abbott et al. 2021e;
Steltner et al. 2021). Both these methods suffer increas-
ing computational costs and decreasing sensitivity of the
searches as parameter space increases.
Searches of all three types have been performed and so
far no convincing evidence for CWs has been observed.
However, searches have probed new regimes, such as pro-
viding upper limits on emission that are more stringent
(i.e., smaller) than those based on energetics arguments.
For example, for several pulsars including the Crab pul-
sar, Vela pulsar (Abbott et al. 2019b), J0537
−
6910 (Ab-
bott et al. 2021c,f) and two millisecond pulsars (Abbott
et al. 2020) the direct upper limits set on the GW am-
plitude are more constraining than limits based on the
assumption that all the pulsars’ spin-down luminosity is
radiated through GWs, known as the spin-down limit.
In this paper we report the new results of a tar-
geted search for CW signals from 236 pulsars using
the most recent LIGO and Virgo datasets including
the second and third observing runs (O2 and O3 re-
spectively). LIGO (Aasi et al. 2015) consists of two
gravitational-wave detectors situated in Hanford, Wash-
ington (H1) and Livingston, Louisiana (L1) while Virgo
(Acernese et al. 2015) is located in Cascina, Pisa (V1).
The ephemerides for the pulsars have been derived from
observations using the CHIME, Hobart, Jodrell Bank,
MeerKAT, Nancay, NICER and UTMOST observato-
ries. More details on these observations can be found in
Section 3.2.
1.2.
Signal models
We assume that the gravitational-wave (GW) emis-
sion is locked to the rotational phase of the pulsar. For
the ideal case of a triaxial star rotating steadily about a
principal moment of inertia axis, the GW emission is at
twice the star’s spin frequency,
f
rot
. However, there are
mechanisms that can produce variations to this 2
f
rot
frequency. For example, a superfluid component with
a misaligned spin axis within the star could produce
a dual-harmonic emission at both once and twice the
rotation frequency without leaving an imprint on the
EM emission (Jones 2010). Therefore we perform two
searches: one at just twice the pulsar rotation frequency
and one at both one and two times the frequency, which
is referred to as a dual harmonic search.
The waveform used in the dual harmonic search is
detailed in Jones (2010) and used in Pitkin et al. (2015);
Abbott et al. (2017a, 2019b, 2020). The signals
h
21
and
h
22
at once and twice the pulsar rotation frequency can
be defined as
h
21
=
−
C
21
2
[
F
D
+
(
α,δ,ψ
;
t
) sin
ι
cos
ι
cos
(
Φ(
t
) + Φ
C
21
)
+
F
D
×
(
α,δ,ψ
;
t
) sin
ι
sin
(
Φ(
t
) + Φ
C
21
)
]
,
(1)
h
22
=
−
C
22
[
F
D
+
(
α,δ,ψ
;
t
)(1 + cos
2
ι
) cos
(
2Φ(
t
) + Φ
C
22
)
+
2
F
D
×
(
α,δ,ψ
;
t
) cos
ι
sin
(
2Φ(
t
) + Φ
C
22
)
]
,
(2)
where
C
21
and
C
22
are the dimensionless constants that
give the component amplitudes, the angles (
α,δ
) are the
right ascension and declination of the source, while the
angles (
ι,ψ
) describe the orientation of the source’s spin
axis with respect to the observer in terms of inclination
and polarization, Φ
C
21
and Φ
C
22
are phase angles at a
defined epoch and Φ(
t
) is the rotational phase of the
source. The antenna functions
F
D
+
and
F
D
×
describe
how the two polarization components (plus and cross)
are projected onto the detector.
For the ideal case of a steadily spinning triaxial star
emitting GWs only at twice the rotation frequency, the
equatorial ellipticity can be defined as
ε
≡
|
I
xx
−
I
yy
|
I
zz
,
(3)
where (
I
xx
,I
yy
,I
zz
) are the source’s principal moments
of inertia, with the star rotating about the
z
-axis. The
mass quadrupole of the source
Q
22
is often quoted and
3
is related to the ellipticity as
Q
22
=
I
zz
ε
√
15
8
π
.
(4)
For single harmonic emission,
C
21
from equation (1) can
be set as 0, leaving only
C
22
in equation (2). The am-
plitude can then be parameterised as the dimensionless
h
0
: the amplitude of the circularly polarised signal that
would be observed if the source lay directly above or
below the plane of the detector and had its spin axis
pointed directly towards or away from the detector. The
following equations are defined in Aasi et al. (2014).
h
0
= 2
C
22
=
16
π
2
G
c
4
I
zz
εf
2
rot
d
,
(5)
where
d
is the distance of the source. The spin-down
limit
h
sd
0
of a source is given by:
h
sd
0
=
1
d
(
5
GI
zz
2
c
3
|
̇
f
rot
|
f
rot
)
1
/
2
,
(6)
where
̇
f
rot
is the first derivative of the rotational fre-
quency, i.e., the spin-down rate, and provides an am-
plitude limit assuming that all of the rotational energy
lost by the pulsar is converted to gravitational-wave en-
ergy (Owen 2005). When the directly observed
h
0
is
smaller than
h
sd
0
, the spin-down limit can be said to
have been surpassed. This information is most often
represented by quoting the “spin-down ratio”, i.e., the
ratio between
h
0
and
h
sd
0
. If assuming that there is no
mechanism (e.g., accretion) providing some additional
spin-up torque, the direct amplitude constraints probe
a new physical regime only when the spin-down limit is
surpassed. There are two types of spin-down rate: in-
trinsic and observed. The observed spin-down rate can
be affected by the transverse velocity of the source, e.g.,
the Shlovskii effect (Shklovskii 1970), so where possi-
ble the intrinsic spin-down rate is used to calculate the
spin-down limit.
1.2.1.
Non-GR polarization signal
In this paper we also perform a search for GWs with
polarizations as predicted in a modification of the stan-
dard general relativity (GR) proposed by Brans and
Dicke. The Brans-Dicke theory (Brans & Dicke 1961)
predicts three independent polarization states: two ten-
sor polarizations, as in GR, and an additional scalar
polarization. The dominant scalar radiation in Brans-
Dicke theory originates from the time-dependent dipole
moment (see Verma 2021, for details). The signal
h
11
due to dipole radiation is given by (see Eqs. (63) - (67)
and (70) - (71) of Verma 2021)
h
11
=
−
h
d
0
c
(
α,δ
;
t
) sin
ι
sin(Φ(
t
) + Φ
0
)
,
(7)
where
c
(
α,δ
;
t
) is the amplitude modulation function
and Φ
0
is the phase angle at time
t
= 0. The explicit for-
mula for
c
(
α,δ
;
t
) is given by Eq. (64) of Verma (2021).
We see that the dipole radiation comes at the rotational
frequency of the pulsar. We assume that the only non-
vanishing component
D
of the dipole moment in the
pulsar’s frame is in the
x
-direction. Then the amplitude
h
d
0
of the signal is given by
h
d
0
=
4
πG
c
3
ζ
Df
rot
d
,
(8)
where
ζ
is the parameter of the Brans-Dicke theory (see
Section 3 of Verma 2021, for details).
2.
ANALYSIS
As in Abbott et al. (2017a, 2019b), three largely in-
dependent analysis methods were used for the searches
in this paper: the
Time-domain Bayesian method
(Sec-
tion 2.1), the
F
/
G
/
D
-statistic method (Section 2.2),
and the
5n-vector method
(Section 2.3). The
F
/
G
/
D
-
statistic and the 5n-vector methods are only used in
searches for pulsars deemed to be high-value: those
which surpass their spin-down limits in the Bayesian
analysis. The Bayesian and
F
/
G
/
D
-statistic methods
search for two signal models: a single harmonic signal
emitted by the
l
=
m
= 2 mass quadrupole mode at
twice the pulsar rotation frequency and a dual harmonic
signal emitted by the
l
=
m
= 2
and
l
= 2
,m
= 1 modes
at twice and once the frequency. The 5n-vector method
restricts the latter search to the
l
= 2
,m
= 1 mode
only. Only one method, the
D
-statistic, is used for the
Brans-Dicke polarization search. The GW emission is
assumed to precisely follow the phase evolution deter-
mined through EM observations, although uncertainties
in values from the EM observations are not accounted
for here.
2.1.
Time-domain Bayesian method
The raw GW strain data are heterodyned using their
expected phase evolution, which includes corrections for
the relative motion of the source with respect to the
detector and relativistic effects (Dupuis & Woan 2005).
They are then low-pass filtered using a cut-off frequency
of 0
.
25 Hz and then down-sampled to one sample per
minute (1/60 Hz bandwidth) centred about the expected
signal frequency now at 0 Hz. For the dual harmonic
search this is repeated so that a time series is obtained
centred at both
f
rot
and 2
f
rot
. Bayesian inference is
used to estimate the remaining unknown signal parame-
ters and the evidence for the signal model (Pitkin et al.
2017). For the parameter inference, the prior distribu-
tions used were those given in Appendix 2 of Abbott
4
et al. (2017a). They were uninformative uniform pri-
ors for the orientation angles, unless restricted ranges
were appropriate as discussed in Section 2.1.2. For the
amplitude parameters, Fermi-Dirac distribution priors
were used (see Section 2.3.5 of Pitkin et al. 2017). The
Fermi-Dirac distributions for each pulsar were set such
that they were close to flat over the bulk of the likeli-
hood while penalizing very large values. To avoid basing
the priors on current detector data the priors were con-
structed by choosing Fermi-Dirac parameters that gave
distributions for which the 95% probability upper bound
was equivalent to the estimated upper limit sensitivity
of the combined LIGO S6 and Virgo VSR4 science runs
at the particular pulsar GW signal frequency.
1
This method also considers the effect of glitches on the
pulsars (Section 2.1.1) and can perform searches with
restrictions on the pulsar orientation (Section 2.1.2).
2.1.1.
Glitches
Although their frequency is usually very stable, pul-
sars occasionally experience a transient increase in ro-
tation frequency and frequency derivative. Such events
are called glitches and are most common in young, non-
recycled pulsars, although they do rarely occur in mil-
lisecond pulsars (Cognard & Backer 2004; McKee et al.
2016). Some of our sample of pulsars glitched during
the course of O2 and O3. We assume that glitches af-
fect the GW phase identically to the EM phase, but with
the addition of an unknown phase offset at the time of
the glitch. This phase offset is included in the parame-
ter inference. For glitches that occur before or after the
range of the data, no phase offset is needed. The pulsars
which experienced glitching during the course of this
analysis are J0534+2200, also known as the Crab pulsar
(Shaw et al. 2021), J0908
−
4913 (Lower et al. 2019) and
J1105
−
6107. They are shown in Table 1 along with the
time (MJD) of the glitch.
2.1.2.
Restricted orientations
Occasionally, the orientation of a pulsar can be deter-
mined from modeling of X-ray observations of its pul-
sar wind nebula (Ng & Romani 2004, 2008). In such
cases, these values can be included as narrow priors on
inclination and polarization angle rather than using an
uninformative uniform prior. Results still assuming uni-
form priors are also included. Such pulsars are shown
in Table 2 below along with their restricted prior ranges
1
For the dual harmonic search for pulsars with signal components
below 20 Hz, the
C
21
priors were constructed without using the
estimated VSR4 Virgo sensitivity. This was to prevent prior from
dominating over the likelihood in this frequency region.
(Abbott et al. 2017a), which are assumed to be Gaussian
about the given mean and standard deviation.
In the case of the Crab pulsar, which both experienced
a glitch and has sufficient observations for restricted
priors, four individual analyses are performed. Each
analysis accounts for the glitch, with combinations of
dual/single harmonic search and restricted/unrestricted
priors.
2.2.
F
/
G
/
D
-statistic method
The time-domain
F
/
G
/
D
-statistic method uses the
F
-
statistic derived in Jaranowski et al. (1998) and the
G
statistic derived in Jaranowski & Kr ́olak (2010). The in-
put data for the analysis using this method are the het-
erodyned data used in the time-domain Bayesian anal-
ysis. The
F
-statistic is used when the amplitude, phase
and polarization of the signal are unknown, whereas the
G
-statistic is applied when only amplitude and phase are
unknown, and the polarization of the signal is known
(as described in Section 2.1.2). The methods have been
used in several analyses of LIGO and Virgo data (Abadie
et al. 2011; Aasi et al. 2014; Abbott et al. 2017a, 2020).
The method also uses the
D
-statistic developed in Verma
(2021) to search for dipole radiation in Brans-Dicke the-
ory. The
D
-statistic search is the first search of LIGO
and Virgo data for dipole radiation as predicted by
Brans-Dicke theory.
In this method a signal is detected in the data if
the value of the
F
-,
G
- or
D
-statistic exceeds a certain
threshold corresponding to an acceptable false alarm
probability. We consider the false alarm probability of
1% for the signal to be significant. The
F
-,
G
-, and
D
-statistics are computed for each detector and each
inter-glitch period separately. The results from different
detectors or different inter-glitch periods are then com-
bined incoherently by adding the respective statistics.
When the values of the statistics are not statistically
significant, we set upper limits on the amplitude of the
gravitational-wave signal.
2.3.
The 5
n
-vector method
The 5
n
-vector method, derived in Astone et al. (2012),
is a multi-detector matched filter in the frequency do-
main, based on the sidereal modulation of the expected
signal amplitude and phase. The method has been used
in several analysis of LIGO and Virgo detector data
(Abadie et al. 2011; Aasi et al. 2014; Abbott et al. 2017a,
2020), and recently Abbott et al. (2021a) combined with
the Band-Sample Data (BSD) framework (Piccinni et al.
2018). This is based on the construction of BSD files,
i.e., complex time series that cover 10 Hz and 1 month
of the original dataset. Using the BSD files, the compu-
5
Table 1.
Pulsars with glitches occurring during the course of the runs used in this analysis.
PSR
Epoch (MJD)
J0534+2200 (Crab)
58687.6448
±
0.0033
J0908
−
4913
58767.34
±
4.5
J1105
−
6107
58582.24
Table 2.
Pulsars with observations sufficient to restrict their orientation priors in terms of inclination and polarization and the
values used as the constraints.
PSR
Ψ (rad)
ι
1
(rad)
ι
2
(rad)
J0534+2200 (Crab)
2.1844
±
0.0016
1.0850
±
0.0149
2.0566
±
0.0149
J0835-4510 (Vela)
2.2799
±
0.0015
1.1048
±
0.0105
2.0368
±
0.0105
J1952+3252
0.2007
±
0.1501
J2229+6114
1.7977
±
0.0454
0.8029
±
0.1100
2.3387
±
0.1100
tational cost of the analysis is reduced to a few CPU-
minutes per source per detector.
The 5
n
-vector method uses a frequentist approach:
the significance of a certain candidate, characterised by
a value of the detection statistic, is established through
the p-value, that is the probability to obtain a larger
value for the statistics in the hypothesis of noise only.
The statistic distribution can be inferred considering a
range of off-source frequencies, as in Astone et al. (2014).
In case of no detection, upper limits on the amplitude
for the single harmonic search are set following Abbott
et al. (2019b). For the dual harmonic search, for simplic-
ity, we only consider the emission at once the rotation
frequency, so we set upper limit on the
C
21
parameter
alone.
For the pulsars which experienced glitching, each
inter-glitch period is analysed independently and then
the resulting statistics are combined incoherently. For
pulsars for which an estimation of the polarization pa-
rameters can be derived from EM observations, see Ta-
ble 2, two searches are carried on: one assuming un-
informative uniform priors on
ψ
and
ι
and one using
restricted Gaussian priors. Only O3 data from LIGO
and Virgo detectors have been used in this search.
3.
DATA SETS USED
3.1.
Gravitational-wave data
The data set used O2 and O3 runs. The O2 run
took place between 2016 October 30 (MJD: 57722.667)
and 2017 August 25 (MJD: 57990.917). Virgo joined
O2 on 2017 August 1. The duty factors for L1, H1,
and V1, were 57%, 59%, and 80%, respectively. O3
took place between 2019 April 1 (MJD: 58574.625) and
2020 March 27 (MJD: 58935.708). Virgo was opera-
tional for the whole of the O3 run. The duty factors for
this run were 76%, 71%, and 76% for L1, H1, and V1,
respectively. The uncertainties on the amplitude and
phase calibration of the detectors for O2 can be found
in Cahillane et al. (2017); Acernese et al. (2018) and
those for O3 can be found in Sun et al. (2020, 2021);
Acernese et al. (2021). For O2, the maximum 1
σ
am-
plitude uncertainties over the range 10-2000 Hz were be-
tween about [
−
2
.
5
,
+7
.
5]% and [
−
8
,
+4]% for H1 and
L1, respectively, and for Virgo the the maximum uncer-
tainty was 5.1%. For O3, the maximum 1
σ
amplitude
uncertainties over the range 10-2000 Hz were between
about [
−
5
,
+7]% and [
−
5
.
5
,
+5
.
5]% for H1 and L1, re-
spectively, and for Virgo the maximum uncertainty was
5%. These ranges are the maximum upper and lower
bound over the full frequency range and over different
measurement epochs over the run, so at specific frequen-
cies/times the uncertainty can be far smaller.
The data used underwent cleaning processes (Davis
et al. 2019; Viets & Wade 2021; Acernese et al. 2021),
specifically the removal of narrowband spectral artifacts
at the calibration line frequencies and power line fre-
quencies. A discussion on the consequences of perform-
ing a search using LIGO data with the narrowband
cleaning of Viets & Wade (2021) applied compared to
that without it applied can be found in Appendix A.
3.2.
Electromagnetic data
EM observations of pulsars produce the timing so-
lutions used as input to the GW searches. These ob-
servations have been made in radio and X-ray wave-
lengths. The observatories which have contributed to
the data set are: the Canadian Hydrogen Intensity Map-
ping Experiment (CHIME) (as part of the CHIME Pul-
sar Project; Amiri et al. 2021), the Mount Pleasant Ob-
servatory 26 m telescope, the 42 ft telescope and Lovell
telescope at Jodrell Bank, the MeerKAT project (as
part of the MeerTime project; Bailes et al. 2020), the
Nan ̧cay Decimetric Radio Telescope, the Neutron Star
Interior Composition Explorer (NICER) and the Molon-
6
glo Observatory Synthesis Telescope (as part of the UT-
MOST pulsar timing programme; Jankowski et al. 2019;
Lower et al. 2020). Pulsar timing solutions were deter-
mined from this data using
Tempo
(Nice et al. 2015) or
Tempo2
(Edwards et al. 2006; Hobbs et al. 2006, 2009)
to fit the model parameters.
We select pulsars whose rotation frequency is greater
than 10 Hz so they are within the sensitivity band of
the GW detectors. This leads to primarily targeting
millisecond pulsars and fast spinning young pulsars. Of
the 236 pulsars in this analysis, 74 are different from the
221 used in the O2 analysis (Abbott et al. 2019b). There
are 168 pulsars in binary systems and 161 are millisec-
ond pulsars with frequencies above about 100 Hz. The
pulsar J0537
−
6910 is not included due to the recently
published searches for it in Abbott et al. (2021f,c).
For some pulsars, ephemerides were only available for
the course of O3. In such cases, only GW data from
O3 was used. This was the case for 102 out of the 236
pulsars in this analysis.
4.
ANALYSIS RESULTS
No evidence for gravitational-wave signals from any
of the included pulsars was found. The results for all
except the high-value targets are shown in Table 3. The
high-value pulsars are shown in Table 4 and discussed in
Section 4.1. As no CWs were observed, we present the
95% credible upper limits on the gravitational-wave am-
plitudes
C
22
and
C
21
for the dual harmonic run (search-
ing for the mass quadrupole modes
l
= 2,
m
= 1
,
2) and
the gravitational-wave amplitude
h
0
for the single har-
monic (
l
= 2,
m
= 2) search. These were all calculated
using coherently combined data from all three detectors
over the O2 and O3 observing runs or just the O3 run,
as appropriate. Due to the calibration amplitude sys-
tematic uncertainties for the detectors, the amplitude
estimates can have uncertainties of up to
∼
8%.
Figure 1 shows the 95% credible upper limits on the
gravitational-wave amplitudes
C
22
and
C
21
for all pul-
sars using the time-domain Bayesian method. In addi-
tion, it shows joint detector sensitivity estimates for the
two amplitudes based on the representative power spec-
tral densities for the detectors over the course of O3. For
an explanation of how these estimates were calculated,
see Appendix C in Abbott et al. (2019b).
The 95% credible upper limits for the GW amplitude
h
0
from the single harmonic analysis for all pulsars,
again using the results from the time-domain Bayesian
method, are shown as stars in Figure 2. The spin-down
limit for each pulsar is represented as a grey triangle. If
the observed upper limit for a pulsar is below the spin-
down limit, this is shown via a dotted green line from the
spin-down limit to the
h
0
limit which is plotted within
a shaded circle. The solid line gives the joint detector
sensitivity estimate over the course of O3.
2
Figure 3 shows a histogram of the spin-down ratio
h
95%
0
/h
sd
0
from the Bayesian analysis for every pulsar for
which calculating a spin-down rate was possible.
3
These
values rely on the pulsar distance and period derivative
which can have associated errors. These are not taken
into account in this study but their presence should be
kept in mind. The mass quadrupole
Q
22
(see equation 4)
also relies on these values and the ellipticity
ε
(equa-
tion (5)) additionally requires the fiducial moment of
inertia
I
fid
zz
which is taken to be 10
38
kg m
2
. The values
for distance and period derivative used can be found in
Table 3 and Table 4.
The single harmonic search was used to place lim-
its on the mass quadrupole
Q
95%
22
which can be used
to find the pulsar’s ellipticity
ε
95%
using equations (5)
and (4). However, for pulsars that did not surpass their
spin-down limits these
Q
22
and
ε
values would lead to
spin-down rates
̇
P
rot
that are greater than (and thus are
in violation of) their measured values. The results for
the Bayesian analysis are shown in terms of the mass
quadrupole
Q
22
and ellipticity
ε
in Figure 4. Also in-
cluded are histograms of the upper limits and spin-down
limits as well as contour lines of equal characteristic age
τ
calculated under the assumption that all spin-down is
due to energy loss through GW emission, i.e., the brak-
ing index is
n
= 5.
From the Bayesian search twenty-three pulsars have
direct upper limits that are below their spin-down limit,
with 89 pulsars within a factor of 10 of their spin-
down limit. There were 90 millisecond pulsars with a
spin-down ratio less than 10. For the dual harmonic
search, the most constraining upper limit for
C
21
was
J2302+4442 with 7
.
05
×
10
−
27
. The smallest
C
22
upper
limit was 2
.
05
×
10
−
27
for J1537
−
5312. As physically
meaningful constraints for the single harmonic search
are only achieved once the spin-down limit has been
surpassed, the following best limits are taken from the
23 pulsars that had
h
95%
0
/h
sd
0
<
1. The smallest spin-
down limit was 0.009 for J0534+2200 (the Crab pul-
sar). The pulsar with the smallest upper limit on
h
0
was J1745
−
0952 with 4
.
72
×
10
−
27
. The best
Q
22
upper
limit was achieved by J0711
−
6830 with 4
.
07
×
10
29
kg m
2
which led to the best limit on ellipticity of 5
.
26
×
10
−
9
.
2
The sensitivity estimate for O3 alone is used as it dominates
compared to the estimate for the O2 run.
3
Spin-down rates cannot be calculated for pulsars with insuffi-
cient distance, frequency or frequency derivative data (see equa-
tion (6)).
7
10
1
10
2
10
3
Gravitational-wave Frequency (Hz)
10
−
27
10
−
26
10
−
25
10
−
24
10
−
23
C
21
/
22
Strain Sensitivity
C
21
Sensitivity estimate
C
21
Results
C
22
Sensitivity estimate
C
22
Results
Figure 1.
The 95% credible upper limits on the gravitational-wave amplitudes for all 236 pulsars using the time-domain
Bayesian method. The pink stars and green crosses show the 95% credible upper limits for the GW amplitudes (
C
22
and
C
21
)
for the dual harmonic search. The solid lines show the estimated sensitivity of all three detectors combined over the course of
O3.
For each pulsar, we performed a model comparison be-
tween the assumption of the data being consistent with
a coherent signal compared to the assumption of an in-
coherent signal
or
noise. This was calculated for both
the dual harmonic (
l
= 2,
m
= 1
,
2) and single harmonic
(
l
= 2,
m
= 2) searches. Specifically, the base-10 loga-
rithm of the Bayesian odds between models is calculated,
which will be referred to as
O
. Of all the pulsars in this
search, none had
O
>
0, meaning in all cases incoherent
noise was more likely than a coherent signal. The pul-
sar with the highest odds overall was J2010
−
1323 with
−
0
.
77.
4.1.
High-value targets
Table 4 shows the results for the analyses on the
high-value targets for the
Bayesian
, the
F
/
G
-statistic
and the
5
n
-vectors
analyses. In this case, Statistic
a
and Statistic
b
change depending on which analysis
method was used. For the Bayesian analysis Statistic
a
and Statistic
b
represent the base-10 logarithm of the
Bayesian odds,
O
, for the dual- and single-harmonic
searches, respectively. This is the same as for the results
in Table 3. For the
F
/
G
-statistic
and
5
n
-vector
analysis
methods they represent the false-alarm probabilities at
the
l
= 2
,m
= 1 and
l
=
m
= 2 modes respectively.
By definition, all high-value pulsars surpassed their
spin-down limits in the Bayesian analysis. Several pul-
sars glitched during the course of the runs: J0534+2200
(Crab pulsar), J0908
−
4913 and J1105
−
6107. The times
of the glitches are shown in Table 1 and the process
for dealing with them is outlined in Section 2.1.1. Ad-
ditionally, some have sufficient information from EM
observations on their orientation to restrict their pri-
ors: J0534+2200 (Crab pulsar), J0835
−
4510 (Vela),
J1952+3252 and J2229+6114. This is discussed in 2.1.2
and the pulsars’ restricted ranges are quoted in Table 2.
For each pulsar with either a glitch or restricted pri-
ors, individual analyses were performed assuming GW
emission at both 2
f
rot
and
f
rot
and just 2
f
rot
. In the