Search for continuous gravitational waves from ten H.E.S.S. sources
using a hidden Markov model

The catalogue of gravitational wave (GW) transients [1] observed during the first two observing runs of the Advanced Laser Interferometer Gravitational$-$Wave Observatory (aLIGO) and advanced Virgo detectors is a testament to the recent advances in the field of GW astronomy [2, 3]. The signals observed to date primarily come from the coalescence of binary black holes [4, 5, 6] and binary neutron stars [7, 8]. An isolated neutron star whose mass is distributed asymmetrically about its rotation axis is also a promising source of persistent, quasi$-$monochromatic continuous gravitational waves (CW) [9, 10]. Emission from such sources is predicted to be at multiples of the star’s spin frequency $f_{*}$ [11]. For example, a mass quadrupole caused by a thermoelastic or magnetic mountain emits a CW signal at $f_{*}$ and $2f_{*}$ [12, 13], r$-$modes emit at roughly $4f_{*}$/3 [14] and a current quadrupole produced by non$-$axisymmetric circulation of neutron super$-$fluid pinned to the crust emits a signal at $f_{*}$ [15].

Extracting a CW signal from noisy detector data often involves matched filtering against a bank of templates. Managing the size of this template bank and thus the computational cost of the search is an important part of designing a practical search [16, 17]. For a coherent search over a total observation time $T_{\rm{obs}}$, the required number of templates grow as $T_{\rm{obs}}^{n}$, whereas the sensitivity increases as $T_{\rm{obs}}^{1/2}$ [18]. The value of the exponent $n$ is determined by the number of unknown parameters. For example, one has $n\gtrsim 5$ for an all$-$sky search, which typically involves searching over sky position, frequency, and frequency derivatives [18].

In a search directed at a particular object, such as the one described here, $n$ typically has a smaller value as the sky position, and in many cases, some frequency derivatives are known. Nevertheless the computational challenge is acute, and it is common to use a semi$-$coherent approach, which involves dividing the search into $N_{s}$ short blocks of length $T_{\rm{s}}$, coherently searching for signals in each block, and then incoherently combining these blocks to cover the entire search space. The sensitivity of a semi$-$coherent search than scales as $\sim N_{\rm{s}}^{1/4}T_{\rm{s}}^{1/2}$ [18]. Although this makes it less sensitive than coherent matched filtering, it allows more data to be used, which in turn improves the sensitivity of the search while allowing it to remain computationally tractable [18].

Phase templates for CW signals from isolated pulsars are often expressed as Taylor expansions about some reference time $t_{0}$ of form (e.g., see Eq. (18) in Ref. [19])

where $\phi(t_{0})$ is the phase at the reference time, ${f_{0}^{k}}(t_{0})$ is the $k^{\rm{th}}$ time derivative of $f_{0}$ evaluated at $t_{0}$ and $k_{\rm max}$ is the highest$-$order time derivative included in the search. One therefore has $n=k_{\rm max}(k_{\rm max}-1)/2$, which grows rapidly with $k_{\rm max}$ [20].

The challenge increases dramatically when the phase of the signal wanders stochastically on timescales of days to weeks. Radio and X$-$ray timing of isolated pulsars show “spin$-$wandering” to be a widespread phenomenon, often going by the name timing noise [21, 22, 23]. It is attributed to various mechanisms including changes in the star’s magnetosphere [24], spin micro$-$jumps [25], super$-$fluid dynamics in the stellar interior [26, 27] and fluctuations in the spin$-$down torque [28, 29]. It is especially pronounced in young pulsars with characteristic ages $\lesssim 10$kyr [30, 21], among which include three targets in this paper. While spin$-$wandering could, in principal, be modelled with higher$-$order derivatives, the number of derivatives required would make $k_{\rm{max}}$ impractically large and thus make the search computationally infeasible.

In this paper, we use a hidden Markov model (HMM) solved by the classic Viterbi algorithm to search for CW signals from isolated neutron stars. This scheme combines an existing, efficient and thoroughly tested maximum likelihood detection statistic called the $\mathcal{F}-$statistic [19] with a HMM. Loosely speaking, the $\mathcal{F}-$statistic coherently searches for a constant$-$frequency signal within a block of data, and the HMM tracks the stochastic wandering of the signal frequency from one block to the next [31]. The method shares some common elements with other semi$-$coherent methods such as StackSlide [20, 32] and facilitates the use of data from the entire second observing (O2) run while remaining computationally manageable. Finally, it can accommodate the situation where the gravitational$-$wave$-$emitting quadrupole is not phase locked to the stellar crust, even when the ephemeris is accurately measured via electromagnetic observation [33].

To sharpen the astrophysical focus of the search, we select 10 pulsars whose pulsations and/or pulsar wind nebulae have been observed by very high$-$energy (VHE) $\gamma-$ray surveys such as the High Energy Stereoscopic System (H.E.S.S.) and the Very Energetic Radiation Imaging Telescope Array System (VERITAS). Previous CW searches have searched for signals from eight of these targets [34, 35, 33, 36, 37, 38, 39, 40, 41, 42, 43]. Leptonic models are often used to completely describe the TeV $\gamma-$rays originating from these sources [44, 45, 46]. The $\gamma-$ray emission is unlikely to be tied directly to a gravitational$-$wave$-$emitting mass or current quadrupole, as it comes from low$-$density plasma accelerated in magnetosphere vacuum gaps [47, 48]. However, it is entirely possible that $\gamma-$ray emission is a reliable co$-$signature of gravitational radiation, because it is tied to young pulsars, whose mass and current quadrupoles at birth have had relatively little time to relax via slow dissipative processes [49, 50, 51]. Moreover such objects are surrounded by active magnetospheric and particle production processes, which may react back on the star. For example, the heating of polar caps by downward$-$moving particles bombarding the stellar surface can induce non$-$axisymmetric variations in the stellar mass density and thus yield a gravitational$-$wave$-$emitting mass quadrupole [52, 53, 54, 55].

In anticipation of a future detection, it is crucial to note that the TeV spectrum provides a comprehensive picture of the electron population around the pulsar [56]. We can therefore obtain better estimates for the age and magnetic field strength around the progenitors, while multi$-$wavelength observations can yield accurate estimates of the distance to the source. Therefore, the sources in this paper offer unusually rich opportunities to do high$-$impact astrophysics in the event of a detection. Likewise, in the absence of a detection, the accurate age, distance, and magnetisation measurements can be converted into accurate GW spin$-$down upper limits [9, 57, 41] to be interpreted in the light of the upper limits from the GW search.

This paper is organised as follows. We outline the targets selected for this study and briefly review their GeV and TeV properties in Section II. In Section III, we review the results from previous CW searches targeting the selected pulsars. Section IV explains the search algorithm which includes a review of the signal model, the detection statistic known as the $\mathcal{F}-$statistic and the HMM formulation. In Section V, we define the search parameter space for each pulsar as well as the detection threshold and data quality vetoes. The results are presented in Section VI. We conclude with a brief statement on the astrophysical implications in Section VII.

The targets for this search are selected from the publicly available TeV catalogue ¹¹1http://tevcat2.uchicago.edu/. We restrict the analysis to TeV$-$emitting pulsars, sources firmly associated with isolated neutron stars and located within the galactic disk (i.e., up to a distance of 7kpc). With these restrictions in place, the catalogue includes two $\gamma-$ray$-$emitting pulsars and eight H.E.S.S. sources firmly identified as pulsar wind nebulae (PWNe) powered by known pulsars. Table 1 gives the astrophysical parameters for the ten targets. Note that a small discrepancy exists between the position of the H.E.S.S. sources and their associated pulsar due to the pulsar being offset from the centre of the TeV emission region, possibly due to its proper motion, asymmetric crushing of the PWN by the surrounding supernova remnant or by asymmetric pulsar outflow [46, 59]. Below, we briefly discuss the GeV and TeV properties of the H.E.S.S. targets.

PSR J0534$+$2200 (Crab pulsar): Fermi Large Area Telescope (Fermi$-$LAT) observations of the $\gamma-$rays emitted by the Crab pulsar reveal a light curve with two main peaks, labelled P1 and P2, both of which are extremely stable in position across the $\gamma-$ray energy band (See Fig. 1 of Ref. [64]). The first $\gamma-$ray pulse leads the radio pulses by 281$\mu$s [64]. This combined with the measurement of photons up to $\sim$20GeV constraints the production site of non$-$thermal emission in the pulsar magnetosphere [64]. In addition, observations taken by VERITAS between 2007$-$2011 reveal $\gamma-$ray pulses at energies $>$100GeV [65] which are narrower than those recorded by Fermi$-$LAT at 100MeV. Combined analysis of Fermi$-$LAT and VERITAS data yields important information about the shape and location of the particle acceleration site and TeV emission mechanism [65].

PSR J0835$-$4510 (Vela pulsar): Fermi$-$LAT observations of this pulsar reveal two narrow and widely separated pulses (labelled P1 and P2) typically observed in high energy pulsars [66]. A complex bridge emission (labelled P3) is also observed between the two main peaks (See Fig. 2 of Ref. [67]). H.E.S.S. II$-$CT5 observations from 2013 to 2015 confirm the characteristics previously known from Fermi$-$LAT studies and reveal a change in the P2 pulse morphology as well as the onset of a new component near P2 at higher energies [68]. Fermi$-$LAT and H.E.S.S. measurements of the phase$-$averaged spectrum also contains important information about the location and origin of TeV emission [67].

HESS J1018$-$589$-$B: This is one of two regions of TeV $\gamma-$ray emission around the supernova remnant G284.3$-$1.8. The TeV emission has a radius of $0.15\pm 0.03^{\circ}$ and is likely to be associated with a pulsar wind nebula powered by PSR J1016$-$5857 [69].

HESS J1356-645: The morphology of this source in radio and X$-$rays, combined with its spectral energy distribution in TeV, as measured by Fermi$-$LAT in 2013, leads to its classification as a clearly identified PWN [70]. The emission is most likely associated with the young and energetic pulsar PSR J1357$-$6429, which is located at a distance of $\sim$5pc from the center of the TeV emission [71].

HESS J1420$-$607 (Kookaburra PWN): This is one of the two TeV sources in a complex of compact and extended radio and X$-$ray sources, called Kookaburra. It is centered north of a young and energetic pulsar PSR J1420$-$6048, believed to be one of the sources powering the complex [72]. Pulsed $\gamma-$ray and extended X$-$ray emissions from the source provide evidence for its PWN nature [73].

HESS J1514$-$591 (MSH 15$-$52): This composite supernova remnant contains a bright X$-$ray PWN powered by PSR J1513$-$5908 which is surrounded by a shell that dominates the emission in the radio band. Correlation between the H.E.S.S. and X$-$ray data and pulsar jets in the X$-$ray domain lead to its classification as a clearly identified PWN [74]. The TeV morphology also suggests asymmetric expansion of the supernova remnant, displacing the PWN in the direction away from the high photon density region [75].

HESS J1718$-$385: Observation of this H.E.S.S. source by the space based X$-$ray Multi$-$Mirror Mission (XMM$-$Newton) reveals a hard X$-$ray source at the position of PSR J1718$-$3825. Diffuse emission in the vicinity of the pulsar suggests the existence of a synchrotron nebula and strengthens its association with the pulsar. However, the relationship between the X$-$ray and $\gamma-$ray emission is not straight forward. The overall asymmetry of the nebula with respect to the pulsar is consistent with the idea of supernova remnant expanding into a non$-$uniform environment [76].

HESS J1825-137: TeV emissions from the PWN extend out to 1.5^∘ from the pulsar’s location, which makes it one the most luminous and potentially the largest of all firmly identified PWNs in the Milky Way [77]. The XMM$-$Newton’s observation of the region, including the pulsar, reveals an elongated non$-$thermal$-$emitting nebula with the pulsar located on the northern border and a tail towards the peak of the H.E.S.S. source. This is also consistent with TeV $\gamma-$rays originating from the inverse Compton mechanism within the PWN powered by PSR J1826$-$1256 [78].

HESS J1831$-$098: Detected by the H.E.S.S. survey in 2011, this source is proposed to be a member of a new class of TeV emitters known as extended TeV halos [79]. A 67ms pulsar, PSR J1831$-$0952, lies at a small angular offset of $\sim$0.05^∘ from the best$-$fit position of the H.E.S.S. source. A conversion efficiency from spin$-$down power to 1$-$20 TeV $\gamma-$rays of $\epsilon=1\%$ would be required to power this H.E.S.S. source, a value typically inferred for other TeV/PWN systems. This combined with the lack of other plausible counterparts in multi$-$wavelength searches provides evidence for a PWN as the most favorable scenario [80].

HESS J1849-000: This extended TeV source is coincident with a hard X$-$ray source IGR J18490$-$0000 [81]. Follow$-$up observations of the source by XMM$-$Newton and Rossi X$-$ray Timing Explorer confirm HESS J1849$-$000 as a PWN powered by a young and energetic pulsar, PSR J1849$-$0001 [63]. This study also provides clear evidence of a PWN that is asymmetric in X$-$rays [82].

Numerous CW searches have previously looked for signals from eight of the targets chosen for this study. In particular, PSR J0534$+$2200 (Crab) has been extensively studied using data from the LIGO science runs (S) 2$-$5 [34, 35, 33, 36] and Virgo science run (VSR) 4 [33]. Most of these searches exploited a Bayesian method, except Refs. [36] and [33] which used Markov Chain Monte Carlo (MCMC) and 5n$-$vector techniques, respectively. Although none of these searches found strong evidence of a CW signal from the Crab pulsar, they improved the upper limits (ULs) on the GW strain amplitude from $4.1\times 10^{-23}$ to $2.4\times 10^{-25}$ over a period of five years.

Data from VSR2 [37] and VSR4 [39] was used to search for a CW signal from PSR J0835$-$4510 (Vela). These searches were carried out using three largely independent methods; Bayesian and $\mathcal{F}/\mathcal{G}$-statistic methods in Ref. [37] and a 5n-vector method in Ref. [39]. No CW signal was detected and ULs were set on the GW strain amplitude. In 2014, Aasi et al. combined observations from S6, VSR2 and VSR4 to look for CW signals from 179 pulsars [38], four of which are included in this HMM search. For seven high$-$value targets, the search was performed using three largely independent methods; Bayesian, $\mathcal{F}/\mathcal{G}$-statistic, and 5n-vector methods. Only the Bayesian method was applied to the remaining 172 pulsars. The search found no credible evidence for GW emission from any pulsar and produced ULs on the emission amplitude.

In early 2017, Abbott et al. [41] repeated a similar analysis for 200 pulsars, five of which are considered here, using the data from aLIGO’s first observing run (O1). For 11 high$-$value targets, the search was performed using the Bayesian, $\mathcal{F}/\mathcal{G}$-statistic, and 5n-vector methods. Only the Bayesian method was used for the remaining 189 pulsars. No detections were reported, and ULs were set on the GW strain amplitude and ellipticity of each target.

Following the second observing run (O2), Abbott et al. [42] performed a CW search for 33 pulsars using the 5n-vector pipeline. These included three of the targets chosen for this study. The narrow$-$band search found no evidence for a GW signal from any of the targets, but significantly improved on ULs for pulsars with rotation frequency $>$30Hz. Later that year, Abbott et al. [43] combined data from O1 and O2 to search for CW signals from 222 pulsars at $f_{0}$ and $2f_{0}$. As in Ref. [40], they used Bayesian, $\mathcal{F}/\mathcal{G}$-statistic and 5n-vector methods for 34 high$-$value targets and only the Bayesian method for the remaining 188 pulsars. Although the search found no strong evidence for GW emission from any pulsar, it did allow for updated ULs on the GW amplitude, mass quadrupole moment and ellipticity for 167 pulsars, and first such limits for 55 others.

We summarise the ULs set by the aforementioned studies for targets relevant to this search in Table 2. Although these targets have been studied in the past, the HMM enables the examination of a much more inclusive set of spin$-$wandering templates than the Taylor$-$expanded phase model assumed in the above references. Young pulsars with TeV $\gamma-$ray emission, such as the ones considered here, exhibit rapid spin$-$wandering. Hence they derive maximum benefit from being searched again with a HMM.

In this paper, we subdivide the total observation time $T_{\rm{obs}}$ into blocks of duration $T_{\rm{drift}}$, during which the signal frequency is assumed to remain constant to a good approximation. Within each block, we search for a signal predicted to be emitted by a biaxial rotor by coherently combining the detector data. This data is available in the form of short Fourier transforms (SFTs), each with length $T_{\rm{SFT}}$ = 30min. In Sections IV.1 and IV.2, we review the signal model and the maximum likelihood detection statistic called the $\mathcal{F}-$statistic, as defined in Ref. [19]. Finally, we incoherently combine the output of each block and recover the optimal frequency path through the data using a HMM. We briefly review the HMM framework in Section IV.3.

We use publicly available data from the O2 run of the aLIGO detectors, specifically the GWOSC$-$4KHz$\_$R1$\_$STRAIN channel, to search for CW signals from the pulsars listed in Table 1 [93]. This observing run took place between 30$-$11$-$2016 and 25$-$08$-$2017. The advanced Virgo detector was also operating for the last month of O2. However, we do not use the data from this detector due to its relatively short observation period and lower sensitivity [93]. During O2, there was a break in the operation of both detectors from 22$-$12$-$2016 23:00:00 UTC to 04$-$01$-$2017 16:00:00 UTC. In addition, a commissioning break took place at the Livingston detector between 08$-$05$-$2017 and 26$-$05$-$2017, while at the Hanford detector it lasted from 08$-$05$-$2017 to 08$-$06$-$2017 [93]. In order to perform a joint search between the two aLIGO detectors, we choose a common period from 04$-$01$-$2017 to 25$-$08$-$2017, with a total $T_{\rm{obs}}$ $\approx$ 234 days.

In Sections V.1 and V.2, we outline the procedure for setting the search parameters, namely the frequency range and coherence timescale. In Section V.3, we briefly discuss whether pulsar glitches have any affect on this search. Section V.4 outlines the procedure for setting the detection threshold. Finally, we describe the vetoes used to distinguish a real signal from a non$-$Gaussian noise feature in Section V.5.

All of the targets listed in Table 1 return candidates with $\mathcal{L}\geq\mathcal{L}_{\mathrm{th}}$. This is not surprising as the detector noise is not Gaussian. We apply the three data quality vetoes outlined in Section V.5 to each candidate to separate a real signal from a non$-$Gaussian noise artifact. The results are summarised in Table 4 while Figs. 13$-$22 in Appendix D show $\mathcal{L}$ against the terminating frequency of the Viterbi paths for each sub$-$band of each target. The search returns 5,256 candidates across 24 sub$-$bands, only 12 of which survive the three data quality vetoes.

In this paper, we present a search for CW signals from ten pulsars selected on the basis of their high$-$energy emissions and association with TeV sources. The search uses publicly available data from the O2 run of the aLIGO detectors. We use a method which combines the maximum$-$likelihood $\mathcal{F}-$statistic with a HMM to efficiently track the secular spin$-$down and stochastic spin wandering of a CW signal. In addition to being fast, HMM has the ability to search a much wider set of $f_{0}(t)$ histories than a coherent search with a Taylor series phase model, noting that the electromagnetic pulsations and CW signal are not necessarily locked together.

For each target, we scan 1Hz sub$-$bands around $f_{*},\ 4f_{*}/3$ and $2f_{*}$ and omit bands below 10Hz. Additionally, we choose a coherence timescale $T_{\rm{drift}}={\rm{min}}(T_{\rm{drift}}^{\prime},T_{\rm{drift}}^{\prime\prime})$ to ensure that the signal frequency falls within one frequency bin over a single coherent step. For all targets chosen for this search, the duration of the coherent segments is limited by the secular spin$-$down ($T_{\rm{drift}}^{\prime}$) and not the stochastic spin wandering ($T_{\rm{drift}}^{\prime\prime}$). Despite this, it is vital to include stochastic spin wandering in the GW phase model because the CW signal may not be locked to the electromagnetic pulsations; the spin of the gravitational$-$wave$-$emitting quadrupole may wander invisibly a lot more than the electromagnetically timed solid crust.