High-frequency continuous gravitational waves searched in LIGO O3 public data with Einstein@Home

TL;DR

This paper reports the most sensitive search to date for high-frequency continuous gravitational waves in LIGO O3 data using Einstein@Home, setting upper limits on gravitational wave amplitude, neutron star ellipticity, and r-mode amplitude.

Contribution

It introduces a highly sensitive search method for high-frequency gravitational waves in LIGO O3 data utilizing volunteer computing, and provides new upper limits on neutron star properties.

Findings

No detection of gravitational waves was made.

Set upper limits on gravitational wave amplitude at 800 Hz: 1.32×10⁻²⁵.

Excluded neutron stars with ellipticity ≥ 1.96×10⁻⁸ within 100 pc.

Abstract

We search for nearly-monochromatic gravitational wave signals with frequencies $800.0~\textrm{Hz} \leq f \leq 1686.0~\textrm{Hz}$ and spin-down $-2.7\times10^{-9}~\textrm{Hz}\,\textrm{s}^{-1} \leq \dot f \leq 0.2\times 10^{-9}~\textrm{Hz}\,\textrm{s}^{-1}$. We use LIGO O3 public data from the Hanford and Livingston detectors and deploy this search on the Einstein@Home volunteer-computing project. This is the most sensitive search carried out to date in this parameter space. Our results are consistent with a non-detection. We set upper limits on the gravitational wave amplitude $h_{0}$ and translate these to upper limits on neutron star ellipticity and on r-mode amplitude. The most stringent upper limits are at $800~\textrm{Hz}$ with $h_{0} = 1.32\times10^{-25}$, at the $90\%$ confidence level. Searching in the high frequency bands allows us to probe astrophysically interesting ellipticities with our results excluding isolated neutron stars rotating faster than $2.5~\textrm{ms}$ with ellipticities $\epsilon \geq 1.96 \times 10^{-8}\left[\frac{d}{100~\textrm{pc}}\right]$ within a distance $d$ from Earth. Our results also exclude r-mode amplitudes $\alpha \geq 7 \times 10^{-7}\left[\frac{d}{100~\textrm{pc}}\right]$ for neutron stars stars spinning faster than 400 Hz.

Tables

Stage 0	$60$	$73$	$4$	$1.5$	$50$	$0.008$	$56$	full range	full range	all-sky	$1.4\text{×}{10}^{18}$	$\mathrm{13\hspace{0.17em}442\hspace{0.17em}787}$
Stage 1	$60$	$73$	$2$	$0.9$	$450$	$0.002$	$39$	$\mathrm{1\hspace{0.17em}000}$	$1$	$5$	$\mathrm{13\hspace{0.17em}442\hspace{0.17em}787}$	$\mathrm{831\hspace{0.17em}387}$
Bayesian			Prior	$n_{\mathrm{live}}$
Stage 2	$120$	$37$	Jeffreys’	750	$9\text{×}{10}^7$	$4.4\text{×}{10}^5$		$80$	$0.3$	$2$	$\mathrm{831\hspace{0.17em}387}$	$\mathrm{831\hspace{0.17em}387}$
Stage 3	coh.	$1$	GMM	250	$4\text{×}{10}^4$	$8\text{×}{10}^4$		-	-	-	$\mathrm{831\hspace{0.17em}387}$	$2$

1	$\mathrm{848.934\hspace{0.17em}979\hspace{0.17em}05}$	$-3\text{×}{10}^{-10}$	$2:29:34.5244$	$-29:27:08.8573$	$1.3\text{×}{10}^{-08}$	$-3.6\text{×}{10}^{-15}$	$0:00:00.0669$
7	$\mathrm{1\hspace{0.17em}220.425\hspace{0.17em}863\hspace{0.17em}627}$	$-1.12\text{×}{10}^{-09}$	$14:53:42.1490$	$-20:27:02.2693$	$1.8\text{×}{10}^{-09}$	$-1.2\text{×}{10}^{-14}$	$0:00:00.9693$

Equations

$$h(t)=F_{+}(\alpha ,\delta ,\psi ;t)h_{+}(t)+F_{\times }(\alpha ,\delta ,\psi ;t)h_{\times }(t),$$

$$h_{+}(t)=A_{+}\cos \Phi (t)$$

$$h_{\times }(t)=A_{\times }\sin \Phi (t)$$

$$A_{+}$$

$$A_{\times }$$

$$\begin{gathered} \Phi ({\tau }_{\text{SSB}})={\Phi }_0+ \\ 2\pi [f({\tau }_{\text{SSB}}-{\tau }_{0SSB})+\frac{1}{2}\dot{f}({\tau }_{\text{SSB}}-{\tau }_{0,\text{SSB}}{)}^2]. \end{gathered}$$

$$\bar{\mathcal{F}}(\mathbf{\lambda }):=\frac{1}{N_{seg}}\sum _{i=1}^{N_{\text{seg}}}{\mathcal{F}}_i({\mathbf{\lambda }}_i).$$

$$d(m_{\text{sky}})=0.15\sqrt{m_{\text{sky}}}\left[\frac{100\text{Hz}}{f}\right],$$

$$\mathcal{D}=\left[45,65\right][1/\sqrt{\text{Hz}}].$$

$$f_{max,n}(x^{\ell })=n{F}^{n-1}(x^{\ell })f(x^{\ell }),$$

$$\{\begin{array}{ll}2{\bar{\mathcal{F}}}^{\text{cand}}\ge 2{\bar{\mathcal{F}}}^{\mathrm{thr}}& =6.76\\ {\log }_{10}{\hat{\beta }}_{\mathrm{S}/\mathrm{G}\mathrm{L}\mathrm{t}\mathrm{L}}^{\text{cand}}\ge {\log }_{10}{\hat{\beta }}_{\mathrm{S}/\mathrm{G}\mathrm{L}\mathrm{t}\mathrm{L}}^{\mathrm{thr}}& =0.\end{array}$$

$$R^3=\frac{\log Z^3-\log Z^2}{\log Z^2}.$$

$$\begin{array}{rl}{\epsilon }^{\text{UL}}& =1.4\times 10^{-6}\left(\frac{h_0^{\text{UL}}}{1.4\times 10^{-25}}\right)\times \\ & \left(\frac{d}{1\text{kpc}}\right){\left(\frac{\text{170 Hz}}{f}\right)}^2\left(\frac{10^{38}{\text{kg m}}^2}{I}\right).\end{array}$$

$$\alpha^{\textrm{UL}}=0.028\left(\frac{h_{0}^{\textrm{UL}}}{$1\text{\times}{10}^{-24}$}\right)\left(\frac{d}{1\mathrm{kpc}}\right)\left(\frac{$100\text{\,}\mathrm{Hz}$}{f}\right)^{3}.$$

Full text

\savesymbol

tablenum

\restoresymbolSIXtablenum

\NewDocumentCommand\lratio e^ s o ¿\SplitArgument1—m Λ\IfValueT#1​^#1\IfBooleanTF#2\expectarg*\expectvar#4\IfNoValueTF#3\expectarg\expectvar#4\expectarg[#3]\expectvar#4 \NewDocumentCommand\lik e^ s o ¿\SplitArgument1—m L\IfValueT#1​^#1\IfBooleanTF#2\expectarg*\expectvar#4\IfNoValueTF#3\expectarg\expectvar#4\expectarg[#3]\expectvar#4 \NewDocumentCommand\prob e^ s o ¿\SplitArgument1—m p\IfValueT#1​^#1\IfBooleanTF#2\expectarg*\expectvar#4\IfNoValueTF#3\expectarg\expectvar#4\expectarg[#3]\expectvar#4 \NewDocumentCommand\fstatistic e^ s o ¿\SplitArgument1—m F\IfValueT#1​^#1\IfBooleanTF#2\expectarg*\expectvar#4\IfNoValueTF#3\expectarg\expectvar#4\expectarg[#3]\expectvar#4 \NewDocumentCommand\semcohfstatistic e^ s o ¿\SplitArgument1—m N_seg $\bar{\mathcal{F}}$\IfValueT#1​^#1\IfBooleanTF#2\expectarg*\expectvar#4\IfNoValueTF#3\expectarg\expectvar#4\expectarg[#3]\expectvar#4 \NewDocumentCommand\expectvarmm#1\IfValueT#2 \delimsize— #2

1 Introduction

Continuous gravitational waves are expected to be produced by a variety of astrophysical scenarios; from non-axis symmetric rotating neutron stars, destabilization of r-modes (Owen et al., 1998), or from more exotic mechanisms such as emission from axion-like particles surrounding back holes (Zhu et al., 2020). The expected gravitational wave amplitude at Earth is several orders of magnitude smaller than that of previously detected signals from compact binary coalescence events, but because the signal is long-lasting ($\sim$ years) one can integrate the detector data over many months and increase the signal-to-noise ratio significantly.

The most challenging searches for this type of signal are the all-sky surveys, where one looks for a signal from a source that is not known. In these searches the number of waveforms that can be resolved over months of observation time is very large and so the sensitivity of the search is limited by its computational cost.

Our earlier all-sky searches for continuous gravitational waves from isolated neutron stars have investigated the frequency range below 800 Hz (Steltner et al., 2023), with a deep search especially devoted to the high sensitivity range of the instruments, between 30-250 Hz (McGloughlin et al., 2025a). In this paper we expand the search space to significantly higher frequencies: 800.0\text{\,}\mathrm{Hz}$\leq f\leq$1\,686.0\text{\,}\mathrm{Hz}. The resolution in the sky increases with frequency, making this range very hard to cover with high sensitivity. Thanks to the computing power donated by the volunteers of the Einstein@Home project we are able to match the sensitivity of (Steltner et al., 2023) and with this search “smoothly” extend its reach well above the 1400 Hz region, where we expect the signal from the fastest known pulsar (Hessels et al., 2006).

The spin-down range of this search – $\qty{-2.7e-9}{\per}\leq\dot{f}\leq\qty{0.2e-9}{\per}$ – is higher than (Dergachev & Papa, 2025a) but lower than that scoped in (Abbott et al., 2022). This strikes a balance between wanting a high spin-down as high intrinsic spin-down rates can potentially channel more of their lost rotational energy into gravitational wave radiation, and the fact that pulsars typically exhibit spin-downs orders of magnitude smaller than $10^{-9}$ Hz/s.

Finally, the gravitational wave amplitude scales linearly with the degree of deformation of the source, the ellipticity, so that at high frequency less deformed sources are detectable comparable with lower frequencies.

We employ a hybrid follow-up procedure on our signal candidates in which we combine a traditional frequentist grid-based method, used for instance in Steltner et al. (2023), followed by a Bayesian non-deterministic follow-up (Martins et al., 2025). Our results are consistent with a non-detection and we set upper limits on the intrinsic gravitational wave amplitude accordingly.

This paper is organized as follows: In Section 2 we describe the signal model, in Section 3 we outline the generalities of this search, in Section 4 we describe the initial stage of the search run on Einstein@Home, in Section 5 we describe the hybrid follow-up procedure of interesting candidates resulting from the initial search, in Section 6 we present our results, and finally in Section 7 we discuss the astrophysical implications of our search results.

2 The signal

This search targets quasi-monochromatic gravitational wave signals described in (Jaranowski et al., 1998; Cutler & Schutz, 2005). These signals display both frequency and amplitude modulations due to the motion of the Earth. The waveforms take the following form:

[TABLE]

where $F_{+}(\alpha,\delta,\psi;t)$ and $F_{\times}(\alpha,\delta,\psi;t)$ are called the detector beam-pattern functions for both the “$+$” and “$\times$” wave polarizations and $h_{+}$ and $h_{\times}$ are the polarization amplitudes. $F_{+}$ and $F_{\times}$ depend on the inclination and declination of the source $(\alpha,\delta)$, the polarization angle $\psi$ as well as t, the time at the detector. The plus and cross waveforms are

[TABLE]

with amplitudes

[TABLE]

Here, $h_{0}\geq 0$ is the intrinsic gravitational wave amplitude, $0\leq\iota\leq\pi$ is the angle between the angular momentum vector of the neutron star and the line of sight and $\Phi(t)$ is the phase of the gravitational wave at time t.

The phase evolution can be expressed in the solar-system barycenter (SSB) frame as:

[TABLE]

The reference time for the phase evolution is taken as $\tau_{0,\textrm{SSB}}=1246070525.0$.

3 Search Generalities

3.1 The Data

This search uses LIGO public data from the first half of the third observation run (O3a), i.e. between GPS time 1238166018 (Apr 01 15:00:00 GMT 2019) and 1254150018 (Oct 03 15:00:00 GMT 2019).

Following Steltner et al. (2022b) and previous Einstein@Home searches, we remove noise that would degrade the quality of the search results, lines and in the frequency domain and loud glitches in the time-domain.

The input to our searches is in the form of half-hour time-baseline Fourier Transforms (SFTs) as shown in Figure 1.

3.2 Hierarchical Semi-coherent Searches

All-sky searches for continuous gravitational wave signals over a broad frequency range are computationally unfeasible using optimal detection methods. Semi-coherent hierarchical schemes have been suggested in order to achieve the best sensitivities under computational budget constraints (Brady et al., 1998; Brady & Creighton, 2000; Krishnan et al., 2004) and are nowadays commonly used (Abbott et al., 2022; McGloughlin et al., 2025a; Dergachev & Papa, 2025b), to cite recent examples.

A semi-coherent approach entails coherent searches over a baseline $T_{\textrm{coh}}$, shorter compared to the entire observation span, and then the incoherent combination of the results from the various $T_{\textrm{coh}}$ segments, typically performed by summing the values of the coherent detection statistic. We use the semi-coherent Global Correlation Transform (GCT) search proposed by Pletsch & Allen (2009); Pletsch (2010), with the multi-detector $\mathcal{F}$-statistic of (Cutler & Schutz, 2005) as the per-segment coherent detection statistic:

[TABLE]

The parameter space points $\bm{\lambda}_{i}$ that enter in this summation (5) depend on the template parameters $\bm{\lambda}$ of the left-hand-side of Eq. 5. The $\bm{\lambda}_{i}$ are on a grid that is coarser than the grid for $\bm{\lambda}$. This is because the resolution increases with a longer observation baseline. The GCT is an efficient way to pick the set $\{\bm{\lambda}_{i}\}$ to approximate the detection statistic at $\bm{\lambda}$.

3.3 Detection Statistics

For a rotating isolated neutron star the signal waveform is described by the eight parameters defined is Section 2, but the search template waveform only depends on the four phase-evolution parameters $\bm{\lambda}$: the frequency $f$, the spin-down $\dot{f}$ and the two sky coordinates of the source $\alpha$ and $\delta$ (Jaranowski et al., 1998).

The $\mathcal{F}$-statistic is the maximum log-likelihood ratio between the signal hypothesis and a pure Gaussian noise hypothesis. Since it is possible to have large $\mathcal{F}$-statistic values due to non-Gaussian disturbances of the data, we also utilize the more robust detection statistic ${\hat{\beta}}_{{{\mathrm{{S}}}/{\mathrm{{G}}{\mathrm{{L}}}{\mathrm{t{L}}}}}}$ (Keitel, 2016), which tests the signal hypothesis against an expanded noise hypothesis i.e. “G” Gaussian noise or “L” lines or “tL” transient lines.

In this paper we use a grid-based deterministic approach in the first 2 stages (Stages 0-1) of the hierarchy and a stochastic approach for the last two stages (Stages 2-3). Both approaches use the detection statistic of Eq. 5. The first stages, where it is crucial to identify disturbance, are also augmented by the use of the ${\hat{\beta}}_{{{\mathrm{{S}}}/{\mathrm{{G}}{\mathrm{{L}}}{\mathrm{t{L}}}}}}$ statistic.

3.3.1 Search Grids

The Stage 0 and Stage 1 coarse grids in frequency and spin-down are each described by a single parameter, the grid spacing, which is constant over the search range. The sky grid is approximately uniform on the celestial sphere orthogonally projected on the ecliptic plane. The tiling is a hexagonal covering of the unit circle with hexagon edge length

[TABLE]

where $m_{\textrm{sky}}$ is a tunable free parameter that controls the coarseness of the grid.

The fine grid is the same as the coarse grid in frequency and sky, and it is refined in spin-down by a factor $\gamma_{1}$.

Table LABEL:tab:FUtable details the set-ups for the entire hierarchy of searches.

4 Stage-0

4.1 Sensitivity Goal and Search Set-up

Since a signal that is not detected at Stage 0 can never be recovered at later stages, the most consequential choices affecting the sensitivity of the search are those defining its earliest stage, in particular the $T_{\textrm{coh}}$, the grids and the clustering procedure that selects the candidates to follow-up. Ensuring a high detection efficiency at this stage while keeping the computational cost at bay is difficult because the first stage of the search has to investigate the entire parameter space, and it is hence the most computationally challenging of the whole hierarchy.

Our sensitivity goal is to match the sensitivity depth (Behnke et al., 2015) of (Steltner et al., 2023), $\mathcal{D}\approx 56~[{1/\sqrt{\textrm{Hz}}}]$, which is non trivial because we are here probing much higher frequencies. We consider different $T_{\textrm{coh}}$ and grid combinations, and estimate their detection efficiency on a population of target signals.

This target population comprises signals with intrinsic amplitudes $h_{0}$ drawn from a uniform distribution constrained to the targeted sensitivity depth bracket

[TABLE]

We estimate the sensitivity of the search in terms of the detection efficiency as a function of sensitivity depth. The detection criterion is that a signal is detected if it survives the first stage cut. This is of course a necessary condition, and since the follow-up is designed to not lose any signal from the target population that has passed the first cut, it is also sufficient.

As described in Section 4.3, the first cut is made based on detection statistic values of uncorrelated candidates. The de-correlation happens through clustering, which is hard to model in detail. We thus implement a simpler and conservative estimate of the detection efficiency.

We draw a signal at a fixed sensitivity depth from our target population and we imagine performing the Stage 0 all-sky search and clustering in a small volume of parameter space around the signal parameters, on data containing only noise, and on data containing noise+signal. Let’s say that this search would require $N^{small}$ templates.

Let’s indicate with $N^{tot}$ the number of templates covered by the entire search, with $N_{FU}^{tot}$ the total number of candidates that we can follow-up, and their ratio with $k=N_{FU}^{tot}/N^{tot}$. If the entire parameter space volume is divided in small search volumes, then on average each will contribute $k\cdot N^{small}$ follow-up candidates – the loudest ones. We take the small volume to be such that this expected number is 1, and set as a detection criterion that the loudest from the search with the signal has to be louder than the loudest from that from the noise-only search.

So, to estimate the detection efficiency we draw from the distribution $f_{max,n}(x^{\ell})$ of the loudest detection statistic $x^{\ell}$ over $n=N^{small}$ realisations, in two cases, noise-only and noise+signal. We compare the two, and consider the signal recovered if the signal draw is higher than the noise draw. We do this 2000 times and the detection efficiency is simply the ratio of the number of recovered signals to 2000.

The probability density function for the loudest takes the form

[TABLE]

where $x^{\ell}$ is the maximum over $n$ trials of a random variable with probability density function $f(x)$ and cumulative probability function $F(x)$. In our case the random variable is $N_{seg}\bar{\mathcal{F}}$ of Eq. 5, so $f(x)=\chi^{2}_{4N_{seg}}(x,\rho^{2})$ is a chi-square distribution with $4N_{seg}$ degrees of freedom and a non-centrality parameter $\rho^{2}=\rho_{sign}^{2}(1-\mu)$. The mismatch value $\mu$ is drawn from the mismatch distribution of the search grid and $\rho_{sign}^{2}$ is the non-centrality parameter in the case when there is no mismatch between signal and template grid point.

This estimate of the detection efficiency is conservative for two reasons: 1) the loudest is taken on $N^{small}$ templates that are not independent. The first stage cuts are instead applied to de-correlated, hence independent, candidates. This means that in order to be detected in a real search the signal has to be louder than the loudest in a smaller set of trials than the ones simulated in our procedure, hence smaller signals are actually detectable. 2) the clustering does more than just de-correlating (Steltner et al., 2022a) and so the signal candidates selected by clustering are often below the “loudest” level of our simulation (see figure 5 of (Steltner et al., 2023)). These effects result in an underestimate of the detection efficiency by about 15%.

Figure 2 shows the run-time versus sensitivity depth achieved by different set-ups at 90% detection efficiency. The chosen search set-up has an estimated sensitivity depth of $\approx 48~[{1/\sqrt{\textrm{Hz}}}]$. After validation by injection-and-recovery studies with clustering applied, the actual sensitivity depth of our chosen set-up at 90% detection efficiency is measured to be about 16% higher, in line with our sensitivity goal of $56~[{1/\sqrt{\textrm{Hz}}}]$.

4.2 Computational load Distribution on Einstein@Home and Atlas

Of the total computational cost of the Stage 0 search, roughly a fifth was carried out on the Atlas computer cluster at the Max Planck Institute for Gravitational Physics, and the remainder on the Einstein@Home project111https://einsteinathome.org which uses the idle time on individual volunteers computers to conduct large-scale computationally intensive searches (Anderson, 2004; Anderson et al., 2006). Both of these search components are run on Graphical Processing Units (GPUs); The NVIDIA Geforce RTX 2070 SUPER model exclusively on ATLAS while a broad and diverse range of GPU models is used on Einstein@Home.

The total number of templates covered by Stage 0 is $1.4\text{\times}{10}^{18}$. The search is divided into work-units (WUs), where each WU searches approximately 150 sky points, $2\textrm{ Hz}$ in frequency and the entire spin-down range, totaling $7.3\text{\times}{10}^{10}$ templates each. A total of $19.2$ million WUs are required to cover the entire parameter range.

Each WU returns the top $6\,250$ results per $50\text{\,}\mathrm{mHz}$, so in total $250\,000$. The reason for using separate “top-lists” for each $50\text{\,}\mathrm{mHz}$ is so that saturation effects from loud disturbances which produce many high-ranking results are effectively confined.

The average runtime of a WU on an ATLAS GPU is approximately 30 minutes, implying a total runtime cost of 1100 GPU-years. This search ran on Einstein@Home for approximately 9 months.

4.3 Post-processing

The following steps are performed on the WU-results returned by the Stage 0 search:

Banding

All results from the same $50\text{\,}\mathrm{mHz}$ are collated together. These results cover the full spin-down range and the whole sky in that frequency band.

Visual Inspection

All $50\text{\,}\mathrm{mHz}$ results are visually inspected to identify disturbed bands Abbott et al. (2017). We do not set upper limits in these bands because the detection efficiency is hard to assess but seeds from these bands are followed up, in the spirit of not leaving any stone unturned.

Clustering

Nearby templates are not completely independent. Thus disturbances or signals tend to contribute several nearby results, whereas statistical fluctuations rarely form clusters. With our clustering technique (Steltner et al., 2022a) we identify results that are likely due to the same underlying cause, and from there on treat them as a single candidate.

A cluster seed is identified, i.e. a set of signal parameters ${\bm{\lambda}}_{\textrm{seed}}$ and an uncertainty range $\Delta{\bm{\lambda}}$, such that, if the cluster were due to a signal, with high confidence the signal parameters ${\bm{\lambda}}_{s}$ would fall in the interval ${\bm{\lambda}}_{\textrm{seed}}\pm\Delta{\bm{\lambda}}$.

Clustering reduces the $4.8\text{\times}{10}^{12}$ results to a more manageable set of order $13$ million seeds. The clustering parameters are determined based on our target sensitivity goal.

The $13$ million seeds are binned in 50 mHz bands and Figure 3 shows how many seeds are found in each band. There exist large excursions from the mean number of seeds, both in excess and in defect. The large positive excursions either come from bands flagged as disturbed by the visual inspection of the results or from bands that house coherent lines which the detection statistic down-weights, but the clustering recuperates as possible weak candidates. The defect in the number of seeds may arise when there are very loud signals or signal-like artefacts that saturate the returned results and conglomerate in few clusters with a large number of occupants. This is the case of the hardware injection 1 at about 850 Hz: its 50 mHz band displays a very low number of seeds compared to most of the other bands. Between 1460-1530 Hz the cause for the small number of seeds is not completely clear, but it is correlated with a steep slope in the power spectral density of the input data which always happens in the surviving shoulders of lines that were cleaned out based on the known lines list Goetz et al. (2021). We exclude these 50 mHz from the upper limit statement.

5 Follow-Up

All follow-up stages are run exclusively on the CPUs on the ATLAS computing cluster.

5.1 The First Follow-up Stage

The first follow-up aims at decreasing the uncertainty around each candidate seed. This is achieved with the same coherent time baseline as Stage 0, but a finer grid. From each search the loudest ${\hat{\beta}}_{{{\mathrm{{S}}}/{\mathrm{{G}}{\mathrm{{L}}}{\mathrm{t{L}}}}}}$ result is taken as the Stage 1 candidate.

Figure 4 shows the Stage 1 candidates from the search and from the target signal population. The distributions for the candidates and the reference population have not yet separated, and further follow-up is necessary.

Based on the signal population results, we define a veto that identifies candidates that should pass on to the next stage: a candidate survives if

[TABLE]

With these thresholds, approximately $94\%$ of the candidates are vetoed, leaving about $800$-thousand candidates to be examined in the next stage. The thresholds also veto the two dismissed test signals.

We evaluate the uncertainty around each candidate from the signal target population. Two test signals Stage 1 candidates are offset from their true parameter values significantly more than the other signals. We set the search region for the next stage ignoring these two signals and accepting their loss, as their inclusion would incur an unreasonably computational expense in the next stage.

5.2 The Bayesian Follow-up Method

Starting from Stage 2, we employ the Bayesian follow-up method of Martins et al. (2025).

The uncertainty region $\Delta\bm{\lambda}$ from Stage 1 is searched at each subsequent Stage with a stochastic method, which produces an estimate of the posterior probability density of the phase-evolution parameters $\bm{\lambda}\in\bm{\lambda}_{\mathrm{seed}}\pm\Delta\bm{\lambda}$ and of the evidence. The evidence $Z$ measures the support for the presence of a signal. This posterior is used as the prior for the next stage. At the initial Monte Carlo stage we use uninformative priors following Jeffreys’ invariance principle (Jeffreys, 1946; Martins et al., 2025).

We use a nested sampling algorithm (Skilling, 2004) – specifically dynesty (Speagle, 2020) wrapped by the bilby-package (Ashton et al., 2019). Within the algorithm, we use a Bilby-native implementation of a differential-evolution random-walk (“diff”) to propose new samples.

We use $Z^{a}$, where $a=2,3$ identifies the stage, as a test statistic to identify false alarms and significant candidates. Following the description in \textciteMartins:2025jnq, we also define a test-statistic $R^{3}$ based on the change of $\log Z^{3}$ relative to the first Bayesian Stage $2$:

[TABLE]

This statistic has proven effective in identifying noise disturbances, since $Z^{a}$ remains roughly constant for our target signals whereas it tends to decrease when evaluated on longer time baselines (higher “a”) on search candidates.

The method recovers all signals of the target population that survived Stage 1, with a $T_{\textrm{coh}}=$120$$ hours at Stage 2 and with a fully coherent search at Stage 3. The corresponding effective parameter space volumes searched ($N_{\star}$ as listed in Table LABEL:tab:FUtable) agree with our findings in Martins et al. (2025).

5.3 Bayesian Follow-up Results

At Stage 2 the evidences estimated from the reference signal population are not significant compared to the search candidates, which we expect to be noise-dominated (Fig. 5, top panel). At Stage 3, conversely, the bulk of the candidate distribution has clearly separated from the reference signal population and has significantly smaller evidences and smaller $R^{3}$ values (Fig. 5, bottom panel).

Three features of the candidate distribution strike as deviations from the bulk behaviour: 1) Two candidates are compatible with the test-signal population. Both candidates correspond to fake signals added to the data at the detector level, the so-called hardware injections, as described in the next section. 2) A long “tail” in the candidate $R^{3}-\log Z^{3}$ plot, comprising of order $700$ candidates with evidences consistent with the signal reference population, but with $R^{3}\approx-1$, inconsistent with the signals. In other words the evidence at Stage 3 is high, but it was much higher at the previous, less sensitive stage. All of these outliers come from only a few frequency bands. An inspection of the input data and of the results from these bands finds that they are plagued by non-Gaussian noise features in one of the two detectors, with many of these being remnants of insufficient line cleaning caused by too-narrow line-cleaning windows. In most bands, the noise features are isolated and give rise to only a small number of outliers ($\sim$ tens). However, about 300 candidates arise from the same 1500-1501 Hz band, infested by the violin modes (Davis et al., 2021). The visual inspection of the Stage 0 results did not flag this as a disturbed band, and in fact the Stage 1 results from these bands are also largely unremarkable. The first stochastic stage, Stage 2, with much higher effective resolution than the previous stages and not using the line-robust statistic ${\hat{\beta}}_{{{\mathrm{{S}}}/{\mathrm{{G}}{\mathrm{{L}}}{\mathrm{t{L}}}}}}$, produces significant values of the evidence, which however decrease in the last fully coherent stage. This behaviour is fully consistent with a disturbance. The effect of these disturbances is so pervasive in the 1500-1501 Hz band that we cannot confidently evaluate the detection efficiency and hence do not set upper limits (see Section 6.2) for them. 3) An outlier with unremarkable evidence $Z^{3}$, but relatively high $R^{3}$, albeit smaller than observed for the reference population. This candidate stands out because of the very low $Z^{2}$, in fact the top panel of Fig. 5 shows that it has the overall smallest Stage-2 evidence of $\approx 9$. Like for some of the outliers discussed above, this also was caused by insufficiently broad cleaning of a violin mode at $1510.9$ Hz.

These results underline the importance of the identification, characterization, and cleaning of non-Gaussian noise prior to a search (Covas et al., 2018; Davis et al., 2019, 2021; Capote et al., 2025).

6 Results

6.1 Hardware Injection Recovery

The hardware injections are fake signals added to the data by using photon calibrators to physically move the detector mirrors in order to provide a check of the entire detection chain (Biwer et al., 2017). These hardware injections serve as reference signals that provide standard detection benchmarks for any continuous wave search pipeline.

Two hardware injection signals from isolated neutron stars fall in our search range, specifically those with IDs 1 and 7 (LIGO & Virgo, 2022). We recover both, marking the first broad survey to recover hardware injection 7. For this reason, in Figure 6, we show the Stage 3 posteriors for the phase-evolution parameters.

Table LABEL:tab:HIRecoveryMonteCarlo lists the parameters of the hardware injections (transitioned to our reference time) and the offsets between these parameters and the posterior sample that maximized the $\mathcal{F}$-statistic in Stage 3. The results showcase the high resolution of the signal space attainable with continuous gravitational waves. The offsets are larger than those reported in Steltner et al. (2023) which were obtained using both LIGO O3a and O3b data, corresponding to an approximately twice as long time-baseline. The uncertainties are within the expected range for an observation period of half a year, given the examined parameter-range.

6.2 Upper Limits

We set upper limits on the gravitational wave amplitude $h_{0}$ in half-Hz bands using the same methodology as Steltner et al. (2023).

As described in previous sections, the upper limit set in a half-Hz band may not hold in the whole band. The sub-bands excluded from the upper limit statements are given in the supplementary materials (McGloughlin et al., 2025b). We recap here the different situations where this happens:

the $50\text{\,}\mathrm{mHz}$ bands identified as disturbed by the visual inspection process (see Section 4.3). This amounts to about $\approx 0.5\%$ of the searched frequency range. 2. 2.

the $50\text{\,}\mathrm{mHz}$ bands where the number of seeds is lower than the average and associated with steep shoulders of partially-cleaned lines (end of Section 4.3). These make up $\approx 0.8\%$ of the searched frequency range. 3. 3.

the 1 Hz region 1500-1501 Hz described in Section 5.3. 4. 4.

the half-Hz bands that are so contaminated by spectral disturbances that 90% detection efficiency cannot be achieved. This is not due to the contamination per se, but to our frequency-domain line-cleaning procedure of the input data that together with the line also removes any signal. There are 79 such bands.

The upper limits on the gravitational wave amplitude $h_{0}$ can be translated to upper limits on the ellipticity $\varepsilon$ of a source modelled as a triaxial ellipsoid spinning around a principal moment of inertia axis ${I}$ at a distance $d$:

[TABLE]

Figure 8 shows the upper limits on the ellipticity $\varepsilon$ for different distances.

Another possible emission mechanism are r-modes, unstable toroidal fluid oscillations driven by the Coriolis force, emitting at $\approx 4/3$ of the spin-frequency. We translate the upper limits on the gravitational wave amplitude $h_{0}$ in upper limits on the r-mode amplitude (Owen, 2010):

[TABLE]

Figure 9 shows the upper limits on the r-mode amplitude $\alpha$ for different distances $d$.

7 Conclusions

We present our all-sky search results for continuous gravitational waves with frequencies 800.0\text{\,}\mathrm{Hz}$\leq f\leq$1\,686.0\text{\,}\mathrm{Hz} and spin-down $-2.7\times 10^{-9}\text{Hz/s}\leq\dot{f}\leq 0.2\times 10^{-9}$ Hz/s using LIGO O3 public data. We identify the $\approx 13$ million most interesting continuous gravitational wave signal candidates from a total of $1.4\text{\times}{10}^{18}$ templates. We recover all hardware injections within our parameter space, marking the first successful recovery of hardware injection 7.

Searching for such high frequency continuous gravitational wave signals allows one to probe interesting regimes of source ellipticity. Let us consider an extremely optimistic case, in order to determine the extent of ellipticitiy regimes we are possibly probing. For a neutron star as close as $10$ pc (Pagliaro et al., 2023) and spinning as the fastest observed pulsar, at $716\textrm{Hz}$ (Hessels et al., 2006), assuming the canonical value for the sources moment of inertia of $10^{38}\textrm{kg m}^{2}$ our upper limits exclude ellipticities in the $10^{-10}$ region. At the distance of the closest of the closest non-recycled neutron from (Pagliaro et al., 2023), 1.9 kpc, our results exclude ellipticities greater than $\approx 2\times 10^{-7}$ level. These are interesting ranges (Morales & Horowitz, 2024; Johnson-McDaniel & Owen, 2013; Bhattacharyya, 2020).

This is the most sensitive all-sky search performed to date in this parameter space region. Our investigation of multiple search configurations makes us confident that this search setup achieves close to maximum sensitivity at the given computational cost.

This search uses a novel Bayesian method for the follow-up of signal candidates (Martins et al., 2025). Compared to past Einstein@Home searches, the new method significantly reduces the computational cost of the follow-up and the complexity of the search hierarchy.

Acknowledgments

We thank the E@H volunteers, without the support of whom this search could not have happened. We acknowledge the use of Topcat (Taylor, 2005); Bilby (Ashton et al., 2019) and scikit-learn (Pedregosa et al., 2011). This research has made use of data or software obtained from the Gravitational Wave Open Science Center (gwosc.org), a service of LIGO Laboratory, the LIGO Scientific Collaboration, the Virgo Collaboration, and KAGRA.