Search for continuous gravitational waves from neutron stars in five globular clusters with a phase-tracking hidden Markov model in the third LIGO observing run

HMMs model the evolution and observation of systems fulfilling two important criteria: a) the evolution of the internal state of the system is Markovian, and b) the internal state is hidden from the observer but may be probed probabilistically through measurement. This section briefly describes the formalism involved; for a more detailed review of HMMs see Rabiner [72].

The hidden state of the system is assumed here to be discretised and bounded, so that each state is in the set $\mathcal{Q}=\{q_{1},q_{2},\ldots,q_{N_{Q}}\}$, where $N_{Q}$ is the total number of hidden states. Time is likewise discretised, with the HMM occupying the hidden state $q(t_{i})\in\mathcal{Q}$ at each timestep $t_{i}\in\{t_{0},t_{1},t_{2},\ldots,t_{N_{T}}\}$ ,where $t_{0}$ is an “initial timestep” before any observations have been made. A sequence of possible hidden states occupied during a particular realisation of the HMM is denoted by $Q=[q(t_{0}),q(t_{1}),q(t_{2}),\ldots,q(t_{N_{T}})]$. The dynamics of the system are captured by the transition matrix $A(q_{j},q_{i})$, whose entries are defined as

At each timestep we assume that the system is observed. For a given realisation we write the sequence of observations as $O=[o(t_{1}),o(t_{2}),\ldots,o(t_{N_{T}})]$. The set of possible observations need not be finite, discretised, or bounded. All we require is that there exists an emission probability function $L(o,q_{i})$ defined as

i.e. the probability of observing $o$ at the timestep $t_{n}$, if the hidden state $q_{i}$ is occupied at $t_{n}$.

Finally we specify a belief about the state of the system at the beginning of the analysis, i.e.

In this paper we have no reason to impose anything other than a uniform prior, viz. $\Pi(q_{i})=N_{Q}^{-1}$.

Given the above ingredients, it is straightforward to compute the probability of a hidden state sequence $Q$, given a sequence of observations $O$:

The goal is to compute $\operatorname*{arg\,max}_{Q}\Pr(Q\mid O)$.

In this paper, the full observation is divided into segments of length $T_{\text{coh}}=T_{\mathrm{obs}}/N_{T}$, which are analysed coherently. Each coherent segment is a timestep in the HMM. The hidden states are pairs of frequencies and phases at the start of a segment, $q(t_{n})=[f(t_{n}),\Phi(t_{n})]$. The phase-dependent $\mathcal{B}$-statistic, defined in Section II.2, plays the role of $\ln L(o,q_{i})$. It then remains to specify the transition matrix. We follow Melatos et al. [42] and adopt the following model for the evolution of the signal frequency $f$ and phase $\Phi$:

where $\xi(t)$ is a zero-mean white noise term, i.e.

and $\sigma$ parametrises the noise strength. The forward transition matrix $A[(f_{j},\Phi_{j}),(f_{i},\Phi_{i})]$ is readily obtained by solving the Fokker-Planck equation corresponding to Eqs. (15)–(18). For practical reasons it is useful to also know the backward transition matrix

which is similarly obtained by solving the backward Fokker-Planck equation for this system.

Details of the transition matrices are given in Appendix A, and an illustrative figure is shown in Fig. 1.

The top panel of Fig. 1 shows a heatmap of $A(q_{j},q_{i})$ as a function of $\Delta\Phi=\Phi(t_{n+1})-\Phi(t_{n})$ and $\Delta f=f(t_{n+1})-f(t_{n})$. As shown in Melatos et al. [42] and Appendix A, it is approximately a truncated Gaussian centred on $(\Delta\Phi,\Delta f)=(0,0)$, with positive covariance between $\Delta\Phi$ and $\Delta f$. The bottom panel shows slices through the distribution at $\Delta f=-1\,,0\,,+1\,\mathrm{bins}$. These slices are the distributions used to compute the transition probabilities in the HMM. The $\Delta f=\pm 1$ bin distributions are identical up to the location of the peak and each contain approximately 20% of the total probability, while the $\Delta f=0$ bin contains the other 60%. The transitions probabilities are broad in phase. All three distributions have similar widths. From this figure we see that the probability is concentrated in the same bin as the previous timestep, but some leakage into neighboring bins does occur.

Melatos et al. [42] included a damping term $-\gamma f$ in Eq. (15). This approximates the tendency of the spin frequency of the neutron star to approach zero in the long term, i.e. it represents the long-term spin down process. In this search we explicitly search over the spin-down rate $f^{(1)}$ as a parameter in the deterministic part of the signal model (see Sections III.5.1 and III.6.1), and assume that once the deterministic spin down is taken into account, the mean drift of the signal frequency is negligible. As such, we impose the condition $\gamma\ll(2fT_{\text{coh}}^{2})^{-1}$ which implies $|f(t_{n+1})-f(t_{n})|\ll(2T_{\mathrm{coh}})^{-1}$ [42]. This renders the effect of $-\gamma f$ in Eq. (15) negligible, and in this work we do not include it. The approximation is justified in more detail in Appendix A.

With the HMM fully specified, we turn to the question of to use it in a CW search. To find signal candidates we use the Viterbi algorithm [72, 40] to find the most likely sequence of hidden states,

More precisely, for each possible terminating frequency bin $f_{i}$ we find the most likely sequence of hidden states which terminates in that bin,

These paths are sorted according to their probability $\Pr[Q^{*}(f_{i})\mid O]$. Paths exceeding a predetermined threshold are marked for follow-up analysis. The determination of the threshold as a function of the desired false alarm probability is discussed in Section IV.1. To maximize numerical accuracy, we work with log probabilities so that the products in Eq. (14) become sums. Hence we threshold candidates according to their log likelihood $\mathcal{L}=\ln\Pr[Q^{*}(f_{i})\mid O]$.

To maximize computational performance the HMM tracking step is performed on GPUs, as with the computation of the $\mathcal{B}$-statistic described in Section II.2.

The phase-tracking implementation of the HMM assumes that the evolution of the rotational phase is driven by a random walk in the frequency, with a strength parametrised by a control parameter $\sigma$ [see Equations (15)–(18)]. Our approach here is to fix $T_{\mathrm{coh}}$ based on the particular details of the search, whether that involves prior belief about the nature of the spin wandering, or a desire to keep the number of sky position/spindown templates under control (see Section (III.3). With $T_{\mathrm{coh}}$ fixed, we then set $\sigma$ to be consistent with the signal model implied by our choice of $T_{\mathrm{coh}}$. Setting $\sigma$ appropriately is a subtle task. If $\sigma$ is too large, the transition probabilities become uniform in both frequency and phase, the dynamical constraints imposed by the noise model are lost, and the method reduces approximately to the phase-maximized $\mathcal{F}$-statistic [56]. If $\sigma$ is too small then the HMM struggles to track stochastic spin wandering or other unexpected phase evolution. Melatos et al. [42] recommended

but did not validate Eq. (24) empirically.

There is no unique prescription for choosing $\sigma$. The uncertainty about the physical origin of the time-varying mass-quadrupole moment makes it difficult to argue in general that one choice of $\sigma$ (reflecting the assumed nature of the underlying noise process) is optimal for a given search. However, $\sigma$ endows the signal model with flexibility and allows the HMM to absorb small errors in parameters determining the secular spin evolution, e.g. sky position and spin-down rate. This in turn reduces the number of templates and hence the computational cost. As the HMM is restricted by construction to moving by at most one frequency bin per timestep, there is a limit to the parameter error which can be absorbed.³³3For sufficiently large mismatches in sky position or spindown a significant fraction of the total power is also lost in each coherent step, whether or not the HMM is allowed to track the frequency evolution of the signal across coherent steps. The aim, then, is to check whether Eq. (24) strikes an acceptable balance, allowing for both secular and stochastic deviations away from a Taylor-expanded signal model, as reflected in our choice of $T_{\mathrm{coh}}$, without discarding the advantage gained by phase tracking.

We test the suitability of Eq. (24) via synthetic data injections into both Gaussian noise and the LIGO O3 data. For a range of $h_{0}$ values we compute the detection probability $P_{\text{d}}$ for $\sigma=0.2\sigma^{*}$, $\sigma^{*}$, and $5\sigma^{*}$ as well as for no phase tracking, at a fixed false alarm probability of $10^{-13}$ per frequency bin. We inject moderate frequency wandering via stepwise changes in $\ddot{f}$; in each of 52 7-day chunks of injected data $\ddot{f}$ is randomly chosen from a uniform distribution spanning $-1.13\times 10^{-18}\leq\ddot{f}/(1\,\mathrm{Hz}\,\mathrm{s}^{-2})\leq 1.13\times 10^{-18}$. The values of $\phi$, $f$, and $\dot{f}$ are updated self-consistently as the signal evolves. With the chosen $\ddot{f}$ distribution, $f(t)$ wanders by approximately 28 frequency bins on average, after correcting for the mean $\dot{f}=T_{\mathrm{obs}}^{-1}[f(T_{\mathrm{obs}})-f(0)]$. Note that the timescale of $\ddot{f}$ evolution is deliberately chosen not to coincide with $T_{\text{coh}}=5\,\mathrm{d}$ to challenge the algorithm. The sky position (right ascension $\alpha$, declination $\delta$) and initial value of $\dot{f}$ is chosen at random. The value of $f$ is chosen uniformly over the full search range when injecting into synthetic Gaussian noise, but fixed at $645.1442\,\mathrm{Hz}$ when injecting into LIGO O3 data to avoid intersecting with non-Gaussian artefacts in the data. Signal recovery is performed at the injected sky position, with a $10^{-4}\,\mathrm{Hz}$ frequency band centred on $f=f(0)$ and $\dot{f}=\left[f(T_{\text{obs}})-f(0)\right]/T_{\text{obs}}$. The injection and recovery setup is specified in Table 2.

In order to find the detection probability $P_{\text{d}}$ at a given probability of false alarm $P_{\text{fa}}$ we require an estimate of the noise-only distribution of the log likelihoods returned by the HMM. To obtain this we fit an exponential tail to the distribution of log likelihoods returned from a large number of simulations, using the procedure outlined in Appendix B and adopted in previous searches [53, 88]. The injected amplitude of the signal is characterised by $h_{0,\mathrm{eff}}$, a combination of $h_{0}$ and $\cos\iota$ which determines approximately the signal-to-noise ratio for many CW search methods and is defined as [89]

For $\cos\iota=1$ we have $h_{0,\text{eff}}=h_{0}$, for $\cos\iota=0$ we have $h_{0,\text{eff}}=h_{0}/\sqrt{8}$, and averaging over a uniform distribution in $\cos\iota$ gives $h_{0,\text{eff}}=(2/5)^{1/2}h_{0}$. Figure 2 graphs $P_{\text{d}}$ as a function of $h_{0,\mathrm{eff}}$ for the choices of $\sigma$ outlined above, showing both the case where the signals are injected into simulated Gaussian noise and the case where the signals are injected into the LIGO O3 data.

In both cases the choice $\sigma=\sigma^{*}$ performs as well as or better than the alternatives as well as no phase tracking across the full $h_{0,\mathrm{eff}}$ range. We adopt $\sigma=\sigma^{*}$ throughout the rest of this paper.

Another important tunable parameter in the HMM is the number of bins used to track the signal phase, $N_{\Phi}$. Melatos et al. [42] fixed $N_{\Phi}=32$. It is worth checking how reducing $N_{\Phi}$ affects performance. We seek to reduce $N_{\Phi}$ without sacrificing sensitivity, because the computational cost is proportional to $N_{\Phi}$.

To this end we test $N_{\Phi}=4$, $8$, and $32$ (c.f. Section III.4.1). Sensitivity curves showing $P_{\text{d}}$ as a function of $h_{0,\mathrm{eff}}$ are displayed in Figure 3.

The differences are relatively minor, so we proceed with $N_{\Phi}=8$ in what follows.

In principle the possible frequency of CW emission from neutron stars spans a wide range. Observed pulsar spin frequencies range from less than $1\,\mathrm{Hz}$ [90, 91] up to $716\,\mathrm{Hz}$ [92]. Although the relationship between the rotation frequency $f_{\text{rot}}$ and the CW signal frequency $f$ depends on the precise nature of the emission mechanism, one normally assumes $f_{\text{rot}}\leq f\leq 2f_{\text{rot}}$ [1]. We remind the reader that the notional target for this search is a recycled neutron star, whose rotation frequency is likely to fall between $0.1\,\mathrm{kHz}$ and $0.7\,\mathrm{kHz}$, as for observed millisecond pulsars. Note that ${\sim}80\%$ of pulsars above $0.1\,\mathrm{kHz}$ fall in the range $0.1\leq f_{\text{rot}}/1\,\mathrm{kHz}\leq 0.4$ [93]. There are also some practical considerations when it comes to selecting the frequency range. The LIGO interferometers are polluted by Newtonian noise below ${\sim}20\,\mathrm{Hz}$, and narrowband spectral features (instrumental lines) below $\sim 0.1\,\mathrm{kHz}$ [94, 95].

The spin-down rate of a CW signal is weakly constrained a priori. Under the optimistic assumption that CW emission is responsible for all of the spindown of the emitting body, the strain amplitude is equal to

where $I_{zz}$ is the component of the moment of inertia tensor about the rotation axis. The expected sensitivity of this search is roughly [56, 78, 96]

which implies $\dot{f}\sim-10^{-10}\,\mathrm{Hz}\,\mathrm{s}^{-1}$ for $h_{0}^{\mathrm{sens}}=h_{0}^{\mathrm{sd}}$. To beat the spin-down limit we require $h_{0}^{\mathrm{sens}}<h_{0}^{\mathrm{sd}}$. We therefore wish to pick a bound on $\dot{f}$ such that the cross-over point $h_{0}^{\mathrm{sens}}=h_{0}^{\mathrm{sd}}$ occurs at a high enough frequency that we beat the spin-down limit for a significant fraction of plausible signal frequencies. Figure 4 shows the $h_{0}^{\mathrm{sens}}$ curve with three $h_{0}^{\mathrm{sd}}$ curves overplotted with $-\dot{f}=10^{-11}\,\mathrm{Hz}\,\mathrm{s}^{-1}$, $-\dot{f}=1\times 10^{-10}\,\mathrm{Hz}\,\mathrm{s}^{-1}$, and $-\dot{f}=5\times 10^{-10}\,\mathrm{Hz}\,\mathrm{s}$. We set $D=5\,\mathrm{kpc}$ in all cases.

The $\dot{f}=-5\times 10^{-10}\,\mathrm{Hz}\,\mathrm{s}^{-1}$ curve crosses over $h_{0}^{\mathrm{sens}}$ at approximately $f=0.8\,\mathrm{kHz}$. Hence for $0.1\leq f/1\,\mathrm{kHz}\leq 0.8$ and $-5\leq\dot{f}/10^{-10}\,\mathrm{Hz}\,\mathrm{s}^{-1}\leq 0$ we expect to approach or beat the spin-down limit across the full frequency range for each cluster. The range includes the bulk of the astrophysical prior on the signal frequency, recalling that the majority of MSPs have $f_{\mathrm{rot}}<0.4\,\mathrm{kHz}$ and thus $f<0.8\,\mathrm{kHz}$. We therefore adopt $0.1\leq f/1\,\mathrm{kHz}\leq 0.8$ and $-5\times 10^{-10}\leq\dot{f}/1\,\mathrm{Hz}\,\mathrm{s}^{-1}\leq 0$ in the present paper. To facilitate data management, the total frequency range is divided into $0.5\,\mathrm{Hz}$ sub-bands which are searched independently.

We note that electromagnetically timed MSPs have spin-down rates orders of magnitude lower than our sensitivity limits according to Fig. 4 [93]. Their spin down is typically attributed to magnetic dipole braking, with the recycling process suppressing the magnetic dipole moment [97, 98, 9, 10, 99, 100]. Hence the signals to which this search is sensitive are most likely to come from a new class of objects with significantly higher spin-down torques, possibly dominated by gravitational radiation reaction [10, 101].

Four of the clusters targeted in this search (NGC 6325, NGC 6397, NGC 6544, and Terzan 6) are core-collapsed; their density profiles peak sharply at the core. At the same time, our cluster selection criterion favors clusters which produce disrupted binaries where the neutron star may be ejected from the core. Indeed clusters similar to those selected are observed to contain millisecond pulsars many core radii away from the centre of the cluster [76]. Therefore, as a precaution, we extend the search to cover sky positions outside the core. Specifically, we extend the search out to the tidal radius as recorded in Harris [75], ensuring that every point within this radius is covered by at least one sky position template.

We follow previous HMM-based searches and take the frequency resolution to be $(2T_{\mathrm{coh}})^{-1}$. Other choices are possible. In principle any desired frequency resolution may be achieved by padding or decimating the data, e.g. $\delta f=(10T_{\mathrm{coh}})^{-1}$ [48, 34], but $\delta f=(2T_{\mathrm{coh}})^{-1}$ has been successfully used many times previously (see e.g. Refs. [40, 43, 48, 53, 52]) and is assumed in Eq. (24).

The HMM accommodates some displacement between the $\dot{f}$ value used to compute the $\mathcal{B}$-statistic and the true $\dot{f}$ by letting the signal wander into adjacent frequency bins at the boundaries between coherent segments. To determine the $\dot{f}$ resolution we empirically estimate the mismatch, i.e. the fractional loss of signal power which is incurred as the value of $\dot{f}$ used in the search is displaced away from the true value, $\dot{f}_{\mathrm{inj}}$, with

where we define $\mathcal{L}^{\prime}=\mathcal{L}-\mathcal{L}_{\text{noise}}$ and $\mathcal{L}_{\text{noise}}$ is the mean likelihood returned in noise over a $0.5\,\mathrm{Hz}$ sub-band. Note that $\mathcal{L}^{\prime}(\dot{f})$ is the maximum $\mathcal{L}$ value over a $0.5\,\mathrm{Hz}$ sub-band centred on the injected frequency, and so $m(\dot{f})$ does not approach $1$ but plateaus at approximately $0.8$, due to the maximisation over noise-only frequency paths (whereas $\mathcal{L}_{\mathrm{noise}}$ is the mean noise-only $\mathcal{L}$ value in a $0.5\,\mathrm{Hz}$ sub-band, not the maximum). We compute $m(\dot{f})$ by injecting moderately loud signals into Gaussian noise, with $h_{0}=6\times 10^{-26}$ and $S_{h}^{1/2}=4\times 10^{-24}\,\mathrm{Hz}^{-1/2}$. An example of $m(\dot{f})$ is displayed in Fig. 5, generated from a single realisation of noise and signal. All mismatch curves are smoothed using a Savitzky-Golay filter.

The shape of $m(\dot{f})$ is consistent across noise realisations. Hence we take the $\dot{f}$ resolution to be the median difference between the two $\dot{f}$ values satisfying $m(\dot{f})=0.2$, drawn from 100 noise realisations. This gives an $\dot{f}$ resolution of $\delta\dot{f}=3.7\times 10^{-12}\,\mathrm{Hz}\,\mathrm{s}^{-1}$. We verify that the resolution is independent of the injected $f$ and $\dot{f}$ values and adopt it throughout the search.

A mismatch in sky location induces a small apparent variation in the signal frequency due to a residual Doppler shift. Although the locations of globular clusters are well-known, they are extended objects on the sky. In this search we take a conservative approach and search out to the tidal radius of each cluster, as discussed in Section III.5.2. In general a single sky pointing does not suffice to cover the full sky area of each cluster. As with a mismatch in $\dot{f}$, we expect the transition matrix in Equation (15) to track some of the residual Doppler variation. Nevertheless it is important to check how much variation can be accommodated before a significant amount of signal power is lost.

We specify the sky location by right ascension ($\alpha$) and declination ($\delta$), referenced to the J2000 epoch. We empirically estimate the mismatch profiles $m(\alpha)$ and $m(\delta)$ and thence derive resolutions in $\alpha$ and $\delta$. As with $\dot{f}$, we smooth the mismatch profiles using a Savitzky-Golay filter. The treatment is somewhat more involved than in Section III.6.1, as $m(\alpha)$ and $m(\delta)$ depend on the injected values $\alpha_{\mathrm{inj}}$ and $\delta_{\mathrm{inj}}$ as well as $f$, i.e. their shape is not a universal function of sky position offset, because the Doppler modulation varies across the sky. In particular, sky positions at declinations near the ecliptic poles experience relatively little Doppler modulation, and hence proportionally smaller residual Doppler modulations [24]. In this paper, we target five relatively small regions of the sky so we do not produce a full skymap of $m(\alpha,\delta)$. Instead we measure the $\alpha$ and $\delta$ resolutions independently for each cluster. The clusters themselves are fairly small (diameter $\lesssim 10\,\mathrm{arcmin}$), so we approximate the $\alpha$ and $\delta$ resolutions to be uniform within each cluster.

Each cluster is covered by ellipses arranged in a hexagonal grid, with their semimajor and semiminor axes given by the resolutions in $\delta$ and $\alpha$ corresponding to $m(\delta)\leq 0.2$ and $m(\alpha)\leq 0.2$. The centre of one ellipse coincides with the cluster core, and the whole cluster out to the tidal radius is covered by at least one ellipse. The centres of the ellipses are taken as the sky position templates in the search. For simplicity we do not consider degeneracy between sky position and $\dot{f}$, which may reduce modestly the number of templates, but adds to the complexity of implementing the search.

Examples of $m(\alpha)$ and $m(\delta)$ are shown in the top two panels of Fig. 6, calculated for Terzan 6 at a frequency of $600\,\mathrm{Hz}$. They are qualitatively similar, with well-defined troughs, but the resolutions, again defined as the distances between the two abcissae where $m=0.2$, differ significantly — the resolution inferred from the $m(\alpha)$ profile is $1.9\,\mathrm{arcmin}$, compared to $15\,\mathrm{arcmin}$ for $m(\delta)$. The bottom panel of 6 shows the measured resolution in $\alpha$ and $\delta$ as a function of signal frequency for a series of injections at the sky location of Terzan 6.

Clearly there is a strong frequency dependence, and it would be wasteful to simply apply the finest resolution across a wide band. Consequently, we produce curves like those in the bottom panel of Fig. 6 for each target. We interpolate by fitting a power-law function $Af^{k}+B$ for $A$, $k$, and $B$ using the lmfit package [106] to find the resolution in each $0.5\,\mathrm{Hz}$ sub-band.

Some of the $0.5\,\mathrm{Hz}$ sub-bands are so disturbed by non-Gaussian artifacts that they are unusable in the absence of noise cancellation [107, 108, 109, 110]. In view of this, we impose a blanket veto on sub-bands which produce output files larger than $100\,\mathrm{MB}$ (corresponding to more than $1.67\times 10^{5}$ candidates per file). Between $3.9$% and $5.4$% of the full search band is removed thus, with the smallest number of removed sub-bands in NGC 6544 (55 out of 1400), and the highest in NGC 6397 (75 out of 1400). These percentages are roughly in line with the all-sky O3a search reported by Abbott et al. [111], for which 13% of the observing band was removed due to an unmanageable number of candidates (noting that the band in that search extended to $2000\,\mathrm{Hz}$). After this veto, approximately $4.9\times 10^{5}$ candidates remain.

The known-lines veto is simple and efficient, rejecting many candidates with relatively little computational cost. The LIGO Scientific Collaboration maintains a list of known instrumental lines [112, 95]. When assessing coincidence we must take into account that the line frequencies and bandwidths in Ref. [112] are quoted in the detector frame, while the HMM frequency path of a candidate is recorded in the Solar System barycentre frame. The resulting Doppler shift amounts to approximately $\pm 10^{-4}f$ [113]. For each candidate we check whether any portion of the HMM frequency path overlaps with a known line in either detector, after correcting for the Doppler shift. In the event of any overlap the candidate is vetoed. After this veto, $18252$ candidates remain.

The cross-cluster veto excludes candidates which have frequencies overlapping (up to Doppler modulation) with candidates in another cluster. Candidates which appear at the same frequency in multiple clusters are consistent with non-Gaussian noise artifacts.

Given the targeted false alarm rate of one per cluster, we expect five candidate which are due to Gaussian noise fluctuations, and a candidate which overlaps with one of these will be falsely dismissed under this veto. Here we briefly show that the probability of a collision between a Gaussian false alarm and an astrophysical candidate is low. Our criterion for coincidence between two candidates is $|f_{1}-f_{2}|<10^{-4}f_{1}$, where $f_{1}$ and $f_{2}$ are the terminal frequencies of the two candidate frequency paths, and $10^{-4}$ is the Doppler broadening factor. To evaluate the probability of a collision, we simulate a search by picking a candidate frequency uniformly over the search range $100\,\mathrm{Hz}\leq f\leq 800\,\mathrm{Hz}$, which is designated as the frequency of the astrophysical candidate, and then pick five more frequencies, again uniformly, which represent the frequencies of the Gaussian false alarm candidates. We then check the coincidence criterion for each Gaussian false alarm against the astrophysical candidate. Simulating $10^{5}$ searches in this way, we find $67$ collisions, i.e. the probability of falsely dismissing an astrophysical candidate due to coincidence with a Gaussian false alarm is approximately $7\times 10^{-4}$, which we find to be acceptably low. We therefore veto any candidate for which there is a candidate in another cluster with respective frequencies $f_{1}$ and $f_{2}$ satisfying $|f_{1}-f_{2}|<10^{-4}f_{1}$.

$17977$ candidates are vetoed this way. We find that all of the vetoed candidates have at least 4 cross-cluster counterparts, i.e. distinct combinations of $\alpha$, $\delta$, and $\dot{f}$ with at least one frequency path with $\mathcal{L}>\mathcal{L}_{\mathrm{th}}$ overlapping with the vetoed candidate — multiple cross-cluster counterparts may be associated with the same cluster. After this veto, $275$ candidates remain.

If a candidate is due to a non-Gaussian feature which is present in the data from only one interferometer, we expect the significance of the candidate to increase when data from only that interferometer is used. We thus veto candidates satisfying $\mathcal{L}_{X}>\mathcal{L}_{\cup}$, where $\mathcal{L}_{X}$ is the log likelihood of the candidate template using only data from detector $X$ and $\mathcal{L}_{\cup}$ is the log likelihood computed using data from both detectors. Before carrying out this veto we group candidates that have frequencies within $10^{-2}\,\mathrm{Hz}$ of each other, in order to reduce the computation required and enable manual verification of the results of the veto procedure. There are 10 candidate groups in total, appearing in two clusters, NGC 6397 and NGC 6544. To compute $\mathcal{L}_{X}$ over the relevant frequency band we use the loudest sky position and $\dot{f}$ template in each candidate group, and $\mathcal{L}_{\cup}$ is taken to be the loudest log likelihood value in the candidate group. Plots showing $\mathcal{L}_{X}$ in comparison to $\mathcal{L}_{\cup}$ for each candidate are collected in Appendix C. In all but one case, one $\mathcal{L}_{X}$ value peaks above $\mathcal{L}_{\cup}$, indicating that the candidate group is likely to be due to a disturbance in the corresponding interferometer and is not astrophysical. These candidate groups are vetoed.

After the single-interferometer veto, two candidates remain.

The line list used for the known-lines veto [112] includes only lines for which there is good evidence of a non-astrophysical origin. Other line-like features are visible by eye in plots of the detector ASDs but have not yet been confidently associated with a terrestrial source.

As a further check on the final remaining candidate group, we manually inspect the ASD in each detector to check whether the candidate group overlaps with a loud non-Gaussian feature in the data, taking into account the Doppler modulation. Candidates which overlap with a clear feature in an ASD are vetoed, on the grounds that we expect any astrophysical signal to be at most marginally visible in the detector ASDs at distances $\gtrsim 2\,\mathrm{kpc}$ [114]. A plot showing the detector ASDs and corresponding $\mathcal{L}$ values for the loudest candidate in the remaining candidate group, in the cluster NGC $6397$ at $f\approx 652.756\,\mathrm{Hz}$, is shown in Figure 7.

There is a clear, broad peak in the H1 ASD within the Doppler modulation range of the candidate, with no corresponding feature visible in the L1 ASD. We therefore veto this candidate group.

After the unknown lines veto, no candidates remain.

For completeness, in Appendix D we present plots analogous to Fig. 7 for all ten candidate groups which survive the cross-cluster veto.