Searches for Gravitational Waves from Known Pulsars at Two Harmonics in the Second and Third LIGO-Virgo Observing Runs

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 electromagnetically 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 ﬁ rst 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.


INTRODUCTION
To date, the LIGO and Virgo observatories have made detections of numerous sources of gravitational radiation. These detections have been of transient gravitational waves (GWs) from the inspiral and subsequent mergers of compact binary objects including binary black holes and binary neutron stars (Abbott et al. 2021a). Recently, the list of observed events expanded to include neutron star-black hole binaries (Abbott 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 frequency 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 neutron stars spinning with some non-axisymmetric deformation (Zimmermann & Szedenits 1979). This may be a solid deformation such as mountains on the crust produced 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 rmodes (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 nonstandard GW polarization 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. Section 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.

Continuous-wave searches
There are broadly three types of CW searches. Targeted searches look for signals from known pulsars for which their rotational phases can be accurately deter-mined from electromagnetic (EM) observations (e.g., Abbott et al. 2017aAbbott et al. , 2019bNieder et al. 2019Nieder et al. , 2020Abbott et al. 2020Abbott et al. , 2021cAshok et al. 2021). This simplifies the search as the EM observations can be used to derive a timing solution and it is assumed that the GW phase evolution is locked to the EM evolution. This means the search is over a small parameter space, generally limited to the unknown signal amplitude and orientation of the source, which allows a more sensitive search than other methods. In some targeted searches, the assumption that the GW evolution follows the EM evolution is relaxed and the search is performed in a narrow band around the expected frequency and spindown rate (Abbott et al. 2017b(Abbott et al. , 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. As the timing solution cannot be derived from EM observations, a wide range of rotational parameters must be searched. All-sky searches look for signals in all sky directions and over a wide range of rotational parameters. A review of directed and all-sky methods and previous searches can be found in Tenorio et al. (2021). Both these methods suffer increasing 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 providing 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 pulsar, Vela pulsar (Abbott et al. 2019b), J0537−6910 (Abbott et al. 2021c,d) and two millisecond pulsars (Abbott et al. 2020) the direct upper limits set on the GW amplitude 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 targeted 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 respectively). LIGO (Aasi et al. 2015) consists of two gravitational-wave detectors situated in Hanford, Washington (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.

Signal models
We assume that the gravitational-wave (GW) emission 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 2f 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. (2017aAbbott et al. ( , 2019bAbbott et al. ( , 2020. The signals h 21 and h 22 at once and twice the pulsar rotation frequency can be defined as h 22 = −C 22 F D + (α, δ, ψ; t)(1 + cos 2 ι) cos 2Φ(t) + Φ C 22 + 2F D × (α, δ, ψ; t) cos ι sin 2Φ(t) + Φ C 22 , 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 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 is related to the ellipticity as For single harmonic emission, C 21 from equation (1) can be set as 0, leaving only C 22 in equation (2). The amplitude 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).
where d is the distance of the source. The spin-down limit h sd 0 of a source is given by: whereḟ rot is the first derivative of the rotational frequency, i.e., the spin-down rate, and provides an amplitude limit assuming that all of the rotational energy lost by the pulsar is converted to gravitational-wave energy (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: intrinsic 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 possible the intrinsic spin-down rate is used to calculate the spin-down limit.

Non-GR polarization signal
In this paper we also perform a search for GWs with polarizations as predicted in a modification of the standard general relativity (GR) proposed by Brans and Dicke. The Brans-Dicke theory (Brans & Dicke 1961) predicts three independent polarization states: two tensor 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) where c(α, δ; t) is the amplitude modulation function and Φ 0 is the phase angle at time t = 0. The explicit formula 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 nonvanishing 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 where ζ is the parameter of the Brans-Dicke theory (see Section 3 of Verma 2021, for details).

ANALYSIS
As in Abbott et al. (2017aAbbott et al. ( , 2019b, three largely independent analysis methods were used for the searches in this paper: the Time-domain Bayesian method (Section 2.1), the F/G/D-statistic method (Section 2.2), and the 5n-vector method (Section 2.3). The F/G/Dstatistic 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 determined through EM observations, although uncertainties in values from the EM observations are not accounted for here.

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 2f rot . Bayesian inference is used to estimate the remaining unknown signal parameters and the evidence for the signal model . For the parameter inference, the prior distributions used were those given in Appendix 2 of Abbott et al. (2017a). They were uninformative uniform priors 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 likelihood while penalizing very large values. This choice of prior means that the amplitude posteriors will be dominated by the likelihood even when no signal is observed. Any upper limits derived from the posteriors will be more conservative than those that would be found from using an uninformative Jeffreys prior that was uniform in the logarithm of the amplitude, i.e., p(h 0 ) ∝ 1/h 0 . To avoid basing the priors on current detector data the priors were constructed 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).

Glitches
Although their frequency is usually very stable, pulsars occasionally experience a transient increase in rotation frequency and frequency derivative. Such events are called glitches and are most common in young, nonrecycled pulsars, although they do rarely occur in millisecond 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 affect 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 parameter 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.

Restricted orientations
Occasionally, the orientation of a pulsar can be determined from modeling of X-ray observations of its pulsar 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 uniform priors are also included. Such pulsars are shown in Table 2 below along with their restricted prior ranges (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.

F/G/D-statistic method
The time-domain F/G/D-statistic method uses the Fstatistic derived in Jaranowski et al. (1998) and the G statistic derived in Jaranowski & Królak (2010). The input data for the analysis using this method are the heterodyned data used in the time-domain Bayesian analysis. 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. 2017aAbbott et al. , 2020. The method also uses the D-statistic developed in Verma (2021) to search for dipole radiation in Brans-Dicke theory. 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 combined 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.

The 5n-vector method
The 5n-vector method, derived in Astone et al. (2012), is a multi-detector matched filter in the frequency domain, 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. 2017aAbbott et al. , 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 computational cost of the analysis is reduced to a few CPUminutes per source per detector. The 5n-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 simplicity, 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 parameters can be derived from EM observations, see Table 2, two searches are carried on: one assuming uninformative uniform priors on ψ and ι and one using restricted Gaussian priors. Only O3 data from LIGO and Virgo detectors have been used in this search. 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 operational 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. (2020Sun et al. ( , 2021; Acernese et al. (2022). For O2, the maximum 1σ amplitude uncertainties over the range 10-2000 Hz were between about [−2.5, +7.5]% and [−8, +4]% for H1 and L1, respectively, and for Virgo the the maximum uncertainty 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, respectively, 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 frequencies/times the uncertainty can be far smaller.
The data used underwent cleaning processes (Davis et al. 2019;Viets & Wade 2021;Acernese et al. 2022), specifically the removal of narrowband spectral artifacts at the calibration line frequencies and power line frequencies. A discussion on the consequences of performing 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.

Electromagnetic data
EM observations of pulsars produce the timing solutions used as input to the GW searches. These observations have been made in radio and X-ray wavelengths. The observatories which have contributed to the data set are: the Canadian Hydrogen Intensity Mapping Experiment (CHIME) (as part of the CHIME Pulsar Project; Amiri et al. 2021), the Mount Pleasant Observatory 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çay Decimetric Radio Telescope, the Neutron Star Interior Composition Explorer (NICER) and the Molonglo Observatory Synthesis Telescope (as part of the UT-MOST pulsar timing programme; Jankowski et al. 2019;Lower et al. 2020). Pulsar timing solutions were determined from this data using Tempo (Nice et al. 2015) or Tempo2 (Edwards et al. 2006;Hobbs et al. 2006aHobbs et al. , 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 millisecond 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. (2021d,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.

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 amplitudes C 22 and C 21 for the dual harmonic run (searching for the mass quadrupole modes l = 2, m = 1, 2) and the gravitational-wave amplitude h 0 for the single harmonic (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 systematic 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 pulsars using the time-domain Bayesian method. In addition, it shows joint detector sensitivity estimates for the two amplitudes based on the representative power spectral 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 spindown 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, frequency derivative and principal moment of inertia, which all have associated uncertainties. These are not taken into account in this study, for which we use the best-fit values listed in Table 3 and 4 and a fiducial moment of inertia I fid zz of 10 38 kg m 2 , but their presence should be kept in mind. Distance errors are primarily based on uncertainties in the galactic free electron distribution (Yao et al. 2017), which can lead to distance errors on the order of a factor of two. Nearby pulsars, for which parallax measurements are possible, will generally have smaller distance uncertainties. Table 3 provides a reference for the distance to each pulsar, which can be used to find an estimate of the associated error as required. The relative uncertainties in frequency derivative are generally much smaller than the distance uncertainties. The principal moment of inertia is equation of state dependent and could range between ∼ (1 − 3) × 10 38 kg m 2 (see, e.g., Abbott et al. 2007). The mass quadrupole Q 22 and the ellipticity ε limits also rely on these values; for example, our given ellipticity upper limits are inversely proportional to I zz /10 38 kg m 2 .
The single harmonic search was used to place limits 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Ṗ 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 included 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 braking 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 spindown 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 spindown limit was 0.009 for J0534+2200 (the Crab pulsar). 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 . This pulsar has a dispersion measure distance of 0.11 kpc, which makes it relatively close by. However, its high ecliptic latitude makes it very insensitive to parallax measurement (Reardon et al. 2021).
For each pulsar, we performed a model comparison between the assumption of the data being consistent with a coherent signal compared to the assumption of an incoherent 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 logarithm 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 pulsar with the highest odds overall was J2010−1323 with −0.77. Table 4 shows the results for the analyses on the highvalue targets for the Bayesian, the F/G-statistic and the 5n-vectors analyses. In this case, the two columns named "Statistic" have different values depending on  Table 3 and assuming the canonical moment of inertia). For those pulsars which surpass their spin-down limits, their results are plotted within shaded circles. The pink curve gives an estimate of the expected strain sensitivity of all three detectors combined during the course of O3. The highlighted pulsars are those with the best h0, Q22 and spin-down ratio out of the pulsars which surpassed their spin-down limit, as well as the best h0 limit out of the whole sample. The Vela pulsar is highlighted and the pulsar J0537−6910 upper and spin-down limits calculated in Abbott et al. (2021c) are also included for completeness. which analysis method was used. For the Bayesian analysis they give 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 5n-vector analysis methods they represent the false-alarm probabilities at the l = 2, m = 1 and l = m = 2 modes respectively.

High-value targets
By definition, all high-value pulsars surpassed their spin-down limits in the Bayesian analysis. Several pulsars 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. Additionally, some have sufficient information from EM observations on their orientation to restrict their priors: 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 priors, individual analyses were performed assuming GW emission at both 2f rot and f rot and just 2f rot . In the case of the Crab pulsar, which both experienced a glitch and has sufficient observations for restricted priors, four individual analyses were performed. Each analysis accounted for the glitch, with combinations of dual/single harmonic search and restricted/unrestricted priors. The values shown in Table 4 are from the searches with glitches accounted for via an unknown phase offset. When a pulsar had a restricted prior search, the results are shown in parentheses next to the unrestricted results. The Crab pulsar is of interest due to its high spindown luminosity. For the single harmonic Bayesian analysis and with the glitch accounted for by a phase offset, its upper limit as a fraction of the spin-down limit is only 0.0094(0.0085) meaning that GWs contribute to less than 0.009% of the available spin-down luminosity. This is consistent with previous studies that also surpassed the spin-down limit (Abbott et al. 2019b(Abbott et al. , 2017b. Its h 95% 0 upper limit was found to be 1.3(1.2) × 10 −26 . With a distance of 2 kpc and period derivative of 4.2 × 10 −13 s s −1 , the upper limits on the mass quadrupole and ellipticity were calculated to be Q 95% 22 = 5.6(5.0)×10 32 kg m 2 and ε 95% = 7.2(6.5)×10 −6 . The base-10 logarithm of the Bayesian odds for this analysis favouring a coherent signal over incoherent noise is -2.6(-2.7).
The Vela pulsar also has a very high spin-down luminosity and is considered another source of interest. Unlike the Crab pulsar, the Vela pulsar did not experience any glitches over the course of this analysis. In the single harmonic Bayesian analysis, the spin-down limit was surpassed with a ratio of 0.052(0.051), with h 95% 0 = 1.8(1.7) × 10 −25 . This ratio corresponds to a maximum of 0.27% of the spin-down luminosity being emitted by GWs. The previous known pulsar search by Abbott et al. (2019b) found the spin-down ratio to be 0.042 with h 95% 0 = 1.4 × 10 −25 which is lower than in this analysis. This is due to significant noise in the LIGO Hanford detector at twice Vela's rotational frequency, with an angular sensing control dither line be-ing the most likely culprit. 4 However, this analysis is an improvement on the more recent measurement of the spin-down ratio of 0.067 and h 95% 0 = 2.2 × 10 −25 produced in Abbott et al. (2020). The upper limits on the mass quadrupole Q 95% 22 and ellipticity ε were calculated to be 7.2(7.1)×10 33 kg m 2 and 9.3(9.2)×10 −5 , respectively. These values were calculated assuming a distance of 0.28 kpc and a period derivative of 1.2×10 −13 s s −1 . The base-10 logarithm of the Bayesian odds for this pulsar in the single harmonic analysis was -1.1(-1.0).
The pulsar J0537−6910 has the highest spin-down luminosity but has not been included in this search due to recently published searches for it in Abbott et al. (2021c,d). The limits, which can be found in Table 3 of Abbott et al. (2021c), are shown for comparison. They found h 95% Table 5 shows the results for the analyses on the highvalue targets using the D-statistic analyses to search for dipole radiation predicted by Brans-Dicke theory. No significant signal was detected and consequently upper limits are obtained. The 95% confidence upper limits are given in the second last column of Table 5 and the last column gives the false alarm probability, i.e., the probability that the obtained values of the D-statistic result only from the noise in the data. The most constraining upper limit for dipole radiation is obtained for the millisecond pulsar J0437−4715.

CONCLUSIONS
In this paper, we have searched for evidence of GWs from 236 pulsars over the course of the LIGO and Virgo O2 and O3 runs and across all three detectors (LIGO Hanford, LIGO Livingston and Virgo). This is an improvement on the 221 pulsars from the O1 and O2 analysis in Abbott et al. (2019b). Searches were carried out for two different emission models. One assumed GW emission from the l = m = 2 mass quadrupole mode and the other assumed emission from the l = 2, m = 1, 2 modes. For the single harmonic search, new upper limits on h 0 were produced and a total of 23 pulsars surpassed their spin-down limits (24 if one includes J0537-6910 from Abbott et al. (2021c)). This is an improvement from the 20 pulsars in Abbott et al. (2019b) and includes 9 pulsars for which their spin-down limit had not previously been surpassed. For the dual harmonic Pulsars for which our direct upper limits have surpassed their spin-down limits are highlighted within a shaded circle with a dotted green line linking the limit to its spin-down limit. Also included are pink contour lines of equal characteristic age τ = P/4Ṗ assuming that GW emission alone is causing spin-down. To the right of the plot, histograms of both these direct limits and spin-down limits are shown by filled and empty bars respectively.
search, new limits on C 21 and C 22 are found. For pulsars which were deemed high-value, two more analysis methods were included for robustness: the F/G-statistic method and the 5n-vector method. The millisecond pulsars that surpassed their spindown limits, J0437−4715 and J0711−6830, have ellipticity upper limits of 8.5×10 −9 and 5.3×10 −9 , respectively. Comparing these values to the left-hand panel of Figure 3 in Gittins & Andersson (2021) finds that our results are lower than the maximum values predicted for a variety of neutron star equations of state. Our results are therefore providing new constraints in physically realistic parts of the ellipticity parameter space.
No search found strong evidence of GW emission from any of the pulsars. However, with so many pulsars now surpassing their spin-down limit, including the millisecond pulsars J0437−4715 and J0711−6830 (Abbott et al. 2020), the next observing run O4 could add more pulsars to this count and bring us closer to observing CWs from pulsars for the first time. In addition to the search for CW signals consisting of the tensorial polarizations predicted by GR, this paper provides the first search ex-plicitly targeting emission of scalar polarization modes predicted by Brans-Dicke theory.

ACKNOWLEDGEMENTS
This material is based upon work supported by NSF's LIGO Laboratory which is a major facility fully funded by the National Science Foundation. The authors also gratefully acknowledge the support of the Science and Technology Facilities Council (STFC) of the United Kingdom, the Max-Planck-Society (MPS), and the State of Niedersachsen/Germany for support of the construction of Advanced LIGO and construction and operation of the GEO600 detector. Additional support for Advanced LIGO was provided by the Australian Research Council. The authors gratefully acknowledge the Italian Istituto Nazionale di Fisica Nucleare (INFN), the French Centre National de la Recherche Scientifique (CNRS) and the Netherlands Organization for Scientific Research (NWO), for the construction and operation of the Virgo detector and the creation and support of the EGO consortium. The authors also gratefully acknowledge research support from these agencies  The Nançay Radio Observatory is operated by the Paris Observatory, associated with the French Centre National de la Recherche Scientifique (CNRS). We acknowledge financial support from "Programme National de Cosmologie et Galaxies" (PNCG) and "Programme National Hautes Energies" (PNHE) of CNRS/INSU, France.
We are grateful to the staff of the Dominion Radio Astrophysical Observatory, which is operated by the National Research Council of Canada. CHIME is funded by a grant from the Canada Foundation for Innovation (CFI) 2012 Leading Edge Fund (Project 31170) and by contributions from the provinces of British Columbia, Québec and Ontario. The CHIME/FRB Project, which enabled development in common with the CHIME/Pulsar instrument, is funded by a grant from the CFI 2015 Innovation Fund (Project 33213) and by contributions from the provinces of British Columbia and Québec, and by the Dunlap Institute for Astronomy and Astrophysics at the University of Toronto. Additional support was provided by the Canadian D.An. acknowledges support from an EPSRC/STFC fellowship (EP/T017325/1).
W.C.G.H. acknowledges support through grants 80NSSC19K1444 and 80NSSC21K0091 from NASA. This work is supported by NASA through the NICER mission and the Astrophysics Explorers Program and uses data and software provided by the High Energy Astrophysics Science Archive Research Center (HEASARC), which is a service of the Astrophysics Science Division at NASA/GSFC and High Energy Astrophysics Division of the Smithsonian Astrophysical Observatory.
We would like to thank all of the essential workers who put their health at risk during the COVID-19 pandemic, without whom we would not have been able to complete this work.

Software:
Astropy ( The data used in this analysis was subject to a cleaning process described in Viets & Wade (2021) which focused on the removal of various narrow-band spectral artifacts at calibration line frequencies. For any pulsars very close to these lines, this cleaning would be expected to provide an improvement in sensitivity. In this appendix we present the comparison of results using this cleaned data against results using data without this cleaning process (which we will refer to as "uncleaned") for a sample of pulsars.
Uncleaned O3 LIGO data is used for a dual harmonic Bayesian analysis of 95 pulsars which had ephemeris data only overlapping with O3. This is compared to the O3-only analysis performed in this paper using the cleaned data. The Virgo data used was the same in both cases. For comparison, the ratio of h 0 upper limits for each pulsar using uncleaned h 0,uncleaned versus cleaned h 0,cleaned data are shown in Figure 5.
The mean ratio of upper limit for uncleaned data versus cleaned data was 0.9966 (with a standard deviation of 0.0486) which suggests no major effect of the line cleaning on the majority results. It should be noted that for this analysis there is a statistical uncertainty on the upper limits of around 1% due to the use of a finite number of posterior samples when calculating them (Abbott et al. 2020). If performing parameter estimation on h 0 using multiple data sets consisting of independent noise realizations drawn from the same distribution, empirically it is found that the resultant upper limits will vary by on order of 30%. Due to the cleaning, the cleaned and uncleaned datasets will contain different, albeit highly correlated, noise. So, a spread of upper limit ratios that are larger than expected from the pure statistical uncertainty on each limit, but smaller that one would get from independent data is to be expected.
As the upper limit ratio spread can be explained as being consistent with expectations from statistical fluctuations, it suggests that very few pulsars are close enough to the cleaned lines for the cleaning to have a significant effect overall. However, it makes sense to apply consistency in using the same dataset for all pulsars being analyzed. In this analysis we chose to use the narrow-band cleaned data.   Table 3 continued  Table 3 continued  Table 3 continued   Table 3 continued    Table 3 continued  Table 3 continued  References-The following is a list of references for pulsar ephemeris data used in this analysis: CHIME: α, Hobart: β, Jodrell Bank: γ, MeerKAT: δ, Nancay: , NICER: ζ, UTMOST: η, .
a The observedṖ has been corrected to account for the relative motion between the pulsar and observer. b The corrected pulsarṖ value is negative, so no value is given and no spin-down limit has been calculated.
c This is a globular cluster pulsar for which a proxy period derivative has been derived assuming a characteristic age of 10 9 years and a braking index of n = 5. · · · · · · · · · · · · · · · · · · · · · Table 4 continued  Table 4 continued  Note-For references and other notes see Table 3. Values in parentheses are those produced using the restricted orientation priors described in Section 2.1.2 a For the Bayesian method this column shows the base-10 logarithm of the Bayesian odds, O, comparing a coherent signal model at both the l = 2, m = 1, 2 modes to incoherent signal models. For the F -/G-statistic method this column shows the false-alarm probability for a signal just at the l = 2, m = 1 mode, assuming that the 2F value has a χ 2 distribution with 4 degrees-of-freedom and the 2G value has a χ 2 distribution with 2 degrees of freedom. For the 5n-vector method this column shows the p-value for a search for a signal at just the l = 2, m = 1 mode, where the null hypothesis being tested is that the data are consistent with pure Gaussian noise. b This is the same as in footnote a , but for all the methods the assumed signal model is from the l = m = 2 mode.
c The observedṖ has been corrected to account for the relative motion between the pulsar and observer. * The pulsar J0537−6910 was not analysed in this study but is included for completion. The values listed here are taken from Abbott et al.
(2021c).  Note-For references and other notes see Table 3. Values in parentheses are those produced using the restricted orientation priors described in Section 2.1.2. The last column shows the false-alarm probability (FAP) for a signal, assuming that the 2D value has a χ 2 distribution with 2 degrees-of-freedom.