A frequency-resolved atlas of the sky in continuous gravitational waves

Astronomy begun with the observation of light from persistent sources. Gravitational-wave astronomy, on the contrary, started with the observation of signals from the merger of compact stars, which are loud transient events lasting less than a minute [2].

This paper is about a different type of gravitational wave signals: those consistently “ON” over months and even years. Such signals, commonly referred to as continuous gravitational waves, are expected from a variety of astrophysical scenarios and could be detected by the LIGO, Virgo or KAGRA interferometers [3, 4, 5]. Among the scenarios that could give rise to continuous waves detectable by ground-based detectors there are rotating neutron stars emitting due to a sustained quadrupolar deformation, fluid oscillations [6, 7, 8] as well as more exotic scenarios as superradiance from ultralight bosons around spinning black holes [9, 11, 10].

The detection of a continuous gravitational wave will lead to high-precision tests of General Relativity and, depending on the emission mechanism, may probe the interior and physics of neutron stars in an entirely new way or unveil new physics [12]. For this reason there exist a variety of research efforts aiming at detecting such signals [13].

Computational searches combining data collected over weeks or even longer periods are carried out using a broad variety of approaches. This variety arises because the most sensitive methods cannot be used due to computational limitations so different speed-optimized methods have been designed, that sacrifice either breadth of search or sensitivity, or a bit of both. Our loosely coherent algorithms [14, 15, 16] and, in particular, the Falcon search pipeline [17, 18, 19] that is used for the analysis presented here, have performed very well in this trade-off.

In this paper we provide a window on the universe of continuous gravitational wave sources. For the first time ever, every sky location has associated a spectrum of upper limits on the gravitational wave amplitude, as well as a signal-to-noise ratio (SNR) spectrum. These are, nearly-raw, minimally processed results and access to them will enable new and independent searches, as we explain in Section V. Together with the new Einstein@Home results [20], these are the most sensitive results to date.

As it is common in broad surveys for continuous gravitational waves, we follow-up and investigate less than 2% of the results of our search – and such results are marked in the atlas. None of the followups results in a conclusive association with astrophysical sources. The rest of the atlas remains unexplored. Other scientists can use the atlas in combination with electromagnetic and particle data [21, 22, 23, 24, 25], and catalogs [27, 26], in search of significant associations. Many other studies are possible. In the appendix, we provide an example of an interesting search that can be easily performed using only the atlas data.

The atlas contains nearly two billion records. To make it easily accessible, the atlas has been designed so that it can be analyzed on a small personal computer using the MVL library [29, 28]. We provide open-source software (compatible with Linux, Windows and Mac OS) that produces sky maps for wide frequency bands, allows queries of the atlas by sky position and/or frequency, and performs computations using the full scan of the data. Effectively this enables searches within the covered parameter space with the full sensitivity of the Falcon pipeline, but without the need for a large computing cluster. Our hope is to make gravitational wave studies and discoveries accessible to scientists with a modest computing budget.

In the following sections we describe the gravitational wave signals, the analysis method, the search setup (including the follow-up) and the results of the analysis.

The ground-based gravitational wave detectors presently in operation [3, 4, 5] are sensitive in the audio frequency range up to a few kHz, with the lower bound at a few tens of Hz given by the steep “seismic wall”.

The seismic wall at low frequencies precludes the observation of continuous gravitational waves from the decay of the orbit of white dwarf binaries with large separation. For this reason the classical sources of audio-frequency continuous gravitational wave signals are rotating neutron stars with non-axisymmetric deformations. For such sources the gravitationl wave emission is powered by the stars’ rotation.

Theory predicts that the neutron star crust can support quadrupolar deformations of $10^{-6}$ and perhaps even larger for neutron stars with exotic composition, however signals like this have not yet been detected [30, 31, 32].

A constraint on the likehood of existence of rapidly spinning neutron stars with deformations much larger than $10^{-8}$ comes from the lifetime of sources: due to emission of the gravitational waves the spin rate of the neutron star decreases in time. Since the energy loss is proportional to the sixth power of the rotation rate, highly deformed (highly emitting) high frequency sources are less likely to exist in our galaxy because they quickly become slowly rotating neutron stars [33]. In fact, known radio pulsars rotating faster than $250$ Hz ($500$ Hz gravitational wave frequency) have frequency derivative less than $2\times 10^{-13}$ Hz/s [26] (Figure 1), corresponding to maximum deformations of the order of $10^{-8}$ and below.

Examination of properties of known pulsars suggests that a population might exist whose spin evolution is governed by gravitational wave emission, with typical ellipticities in the range of $\approx 10^{-9}$ to $10^{-8}$ [34] (see Figure 1).

The expected signals are nearly monochromatic, with small frequency drifts in time. The frequency $f(t)$ of such signals is well described by Taylor polynomials with few components [18]. For our search the second order Taylor polynomial is sufficient:

We call these “IT2” type signals, where “I” indicates an isolated source (no binary modulation), “T” indicates Taylor polynomial model and “2” gives degree of the polynomial [18].

The signal at the detector is more complex as it presents frequency modulations due to the Doppler effect induced by the relative motion between the source and the receiver. For isolated neutron stars, which are the target of this paper, the largest Doppler shifts are due to Earth orbital motion [35, 36], and are on the order of $10^{-4}~{}f_{0}$, with $f_{0}$ being the signal frequency. Additionally the signal also presents amplitude modulations due to the non-uniform beam-pattern response of the detector across the sky.

The first step of the analysis consists in heterodyning the data at a frequency close to the target signal frequencies and in extracting a narrow band from the heterodyned data. This step is repeated over the different frequency sub-bands which compound the broad search band and it is efficiently accomplished in the Fourier domain.

After this transformation the data for the band of interest is a complex time series $\left\{z(t)\right\}$, down-sampled with respect to the original data. In this representation, the gravitational wave signal has the form

where $h_{0}$ is the intrinsic amplitude of the signal at the detector and the amplitude response $A(t)$ depends on the relative orientation of the source and the detector and hence changes as the Earth rotates. The phase evolution term $\Phi(t)$ at the detector is due to the Doppler-shifted frequency evolution of Eq. 1, it is independent of orientation but varies with sky location, frequency and frequency derivatives [36, 35].

The input data is then the sum of the unknown signal and noise $n(t)$:

The term $n(t)$ is a combination of fundamental noise sources and “technical” noise sources. The former include shot noise and radiation pressure on the instrument mirrors, and can be considered Gaussian. The technical noise sources are due to accidental couplings between noise sources external and internal to the detector and the gravitational wave output channels. These couplings are highly dependent on frequency and are non-stationary. Since a lot of effort has gone into reducing these couplings, there are many “clean” bands where these couplings are small enough to be practically negligible. Even in these clean bands, however, the $n(t)$ is in general not a stationary process – the reason being that the beam power in the detector arms varies in time affecting both the level of fundamental noise sources and the gain of the detector.

Next we consider data with the amplitude response and phase evolution terms taken out:

$\Phi_{0}(t)$ is an assumed phase evolution, in general distinct from the true evolution $\Phi(t)$ when we search for unknown signals. We assume that the amplitude evolution for a signal with phase $\Phi_{0}(t)$ is practically equal to $A(t)$ of the signal with phase evolution $\Phi(t)$ which is correct for signals from sky locations within $1^{\circ}$ of each other and the same gravitational wave polarization. A generic sky position can also be considered leading to the slightly different expression [15] actually used in the search, but here we make this simplifying assumption in order to illustrate the main idea.

We construct a linear time-domain filter ${\mathcal{L}}$, say a band-pass filter, with wide enough bandwidth to leave invariant signals of the type

and with as narrow bandwidth as desired to reject the maximum possible amount of noise.

The power $P=\left\|\mathcal{L}\tilde{z}\right\|^{2}$ can then be used to establish a limit on the strength $h_{0}$ of the gravitational wave signal in the data. This is best achieved based on excess power measurements, i.e. by considering a set of power measurements for a set of phase evolution models $\Phi_{0}$ and by setting an upper limit under the assumption that only a single signal is present [37].

In our approach there is an inherent trade-off between the maximum phase mismatch of phase evolutions $\left\{\Phi(t)\right\}$ from $\Phi_{0}(t)$ and the narrowness of the filter $\mathcal{L}$. The narrower is the filter, the less noise it lets through and the more sensitive the search is.

The structure of the manifold of phase evolutions $\left\{\Phi(t)\right\}$ used to design the filter is very important but its description goes well beyond the scope of this paper [36].

An important characteristic of the manifold of phase evolutions is the maximum phase variation $V$ over the timescale $\delta$:

When $V(\delta)$ is close to $0$, the filter $\mathcal{L}$ can be a simple averaging of the data. As the variation grows, the filter $\mathcal{L}$ needs to include an explicit - and more computationally expensive - demodulation. As $V(\delta)$ exceeds $\pi$, the exponent $\exp i\left(\Phi(t)-\Phi_{0}(t)\right)$ changes sign. The value of $\delta$ when $V(\delta)=\pi$ is known as coherence length.

By picking $\mathcal{L}$ as projection onto a constant, we obtain a very narrow filter that admits only $\Phi(t)-\Phi_{0}(t)=\textrm{const}$ - this is known as a “fully-coherent” search. At the other extreme, if $\mathcal{L}$ is the identity operator which leaves its input unchanged, the power $P=\left\|\mathcal{L}\tilde{z}\right\|^{2}=\left\|\tilde{z}\right\|^{2}$ discards the phase information and we arrive at a method with the shortest coherence length possible given the bandwidth of the heterodyne procedure used to obtain $z(t)$. The loosely coherent method [14, 15, 16] employed by Falcon uses the filter $\mathcal{L}$ with a tunable coherence length [16].

In practice in order to search for unknown sources anywhere on the sky one needs to repeat the above procedure for many different assumed phase evolutions $\Phi_{0}$, corresponding to different source sky positions and signal frequency and frequency derivatives. We thus arrive at the fundamental difficulty of all-sky searches: the number of different phase evolutions that need to be considered and hence the computational demands, scales with coherence length, which refines the resolution in sky, frequency and frequency-derivatives.

The Falcon pipeline used in this paper employs a loosely coherent algorithm [14, 15, 16] that finds the maximum power $P$ for a set of assumed evolutions $\left\{\Phi_{0}^{i}(t)\right\}$. Operating on a set allows to frame the problem as the optimization of a highly oscillatory function and take advantage of efficiencies of scale.

The typical analysis takes several months to complete on a facility of the scale of the ATLAS computing cluster [38], so with every new analysis we strive to better take advantage of existing data and computational infrastructure. While limits of computational hardware dictate the choice of the algorithm used, compilation of the code, scheduling, threading of the compute jobs on the physical cores and even aspects such as data layout are carefully considered for every run because they can have a large impact on computational speed.

Based on the power $P$, the algorithm establishes upper limits [37] and SNRs for every analyzed portion of the sky and every frequency band. This data set is quite large and in the past was considered impractical to record, so only summary information was presented, obtained by grouping the results by frequency.

For this search, instead, thanks to the use of Mapped Vector Library (MVL) files, we could keep and distribute all the results – upper limits and SNRs, adding up to over 70 GB of data. The MVL technology was originally used to extract large amounts of engineering data describing the pipeline’s operation in real time. Only later it was realized that it is also perfectly suited for the scientific output of the pipeline. The key property of MVL files is that they are optimized for memory mapping. For example, with a notebook with 16 GB RAM and solid state drive, one can open the atlas in its entirety and access any data by simple array subscription, as if the notebook had an order of magnitude more memory than it actually does.

Our search uses data of the LIGO H1 and L1 interferometers, from the publicly released O3a set [39]. We limit our search to the LIGO detectors because no other detector comes within a factor of 2 of their sensitivity. We set O3b data [40] aside as an independent data set for verification of any surviving outlier. This is a standard approach for searches over a very broad parameter space: due to the large trials’ factor, noise alone produces outliers that can only be vetoed through verification on a different data set (see for example [41]). We note that this procedure can be used because the same physical signal is assumed present in both data sets.

The O3a data set is contaminated by loud short-duration excesses of noise [42, 43, 44]. To address this, we first prepare short Hann-windowed Fourier transforms (SFTs) with 0.25 s duration. We then monitor the frequency band 20-36 Hz for power excess that indicates noise contamination and exclude the affected SFTs from the analysis. The remaining SFTs are searched with the Falcon pipeline. The 20-36 Hz band was chosen because it is outside of the search range and heuristically we found it to be an effective “witness band” to identify spectral disturbances.

We search for nearly monochromatic signals of the IT2 type given by Eq. 1 with the signal frequency $f_{0}$ and its first derivative $f_{1}$ defined at GPS epoch $t_{0}=1246070000$ (2019 Jul 2 02:33:02 UTC).

The search targets emission in the $500-1000$ Hz frequency range, covering frequency derivatives $|f_{1}|\leq 5\times 10^{-11}$ Hz/s and $|f_{2}|\leq 4\times 10^{-20}$ Hz/$\textrm{s}^{2}$. Compared to our earlier analyses on O2 data [19] we use a more sensitive data set and we search a spindown range larger by a factor $\gtrsim$ 16. This provides a larger margin ($\times 2.9$ at $1000$ Hz) for Doppler shifts due to stellar motion and energy loss due to magnetic fields (Figure 1).

The analyzed frequency band was chosen to make the first gravitational wave atlas simpler to use, while providing the most useful data. The lower cutoff of 500 Hz was picked to avoid numerous detector artifacts, while the upper cutoff of 1000 Hz was imposed to limit total data size, as the size of the atlas scales as cube of upper frequency.

The pipeline employs 4 stages (Table 1) with coherence length increasing from 12 hours to 6 days. Results exceeding the SNR threshold of one stage pass to the next stage. The upper limits and the atlas are constructed based on the results of the first two stages only. This simplifies the upper limit procedure without significant impact on the sensitivity. The last stage of the pipeline has long enough coherence time that signals that have survived up to that stage can now be identified using the data from the individual interferometers separately. This means that signal-consistency checks can be performed by comparing the results from each detector.

The Falcon pipeline computes upper limits on the intrinsic amplitude $h_{0}$ of continuous gravitational waveforms as a function of signal frequency and source sky position. We provide upper limits [45] for circularly polarized gravitational waves, which are generated by optimally oriented sources, as well as worst-case upper limits computed by maximizing over polarization. Such upper limits hold for sources with the most unfavorable orientation.

The upper limits are part of the frequency-resolved atlas of the sky, which also includes the measured SNR as a function of signal frequency and source sky-position. Specifically, the atlas contains sky-resolved data for every $45$ mHz signal-frequency band. The sky resolution increases proportionally to frequency, and at $1000$ Hz at any location in the sky there is a grid point within a $0.26^{\circ}$ radius. For each sky pixel and frequency bin we also provide the highest SNR value and the corresponding frequency. Our simulations show that if the maximum SNR exceeds 9 at the closest sky grid point to the signal location, the signal frequency is within $0.6$ mHz of the frequency of the maximum, with 95% confidence.

The upper limits are summarized in Figures 2 and 3, by taking the maximum over the sky in 45 mHz frequency bands. This is the traditional way to present continuous wave upper limit results.

The upper limits can be taken as a measure of the sensitivity of our search, and recast in terms of the equatorial ellipticity $\varepsilon$ of sources that the search is able to detect:

where $I_{zz}=10^{38}$ kg m² is the moment of inertia of the star with respect to the principal axis aligned with the rotation axis and $d$ its distance.

This search is sensitive to sources with ellipticity of $10^{-8}$ up to 150 pc away, improving on our O2 results by a factor of $\approx$ 1.5 and on the full O3 results of [46] and [10] by a factor of $\approx$ 1.7 and 1.5, respectively. Signals sourced by ellipticities of $10^{-7}$ would have been detected from sources within a distance of 500 pc and, depending on orientation, as large as 1.5 kpc, at the highest frequencies. For lowest frequencies the distance is up to 200 pc, and as large as 600 pc for optimally oriented sources.

The results in [46] cover a much broader spin-down range than our search. They exclude neutron stars with ellipticities larger than $10^{-6}$ within a few kiloparsecs, motivating a search like ours for lower ellipticity sources. We accomplish this with a spindown range that covers ellipticities as high as $10^{-7}$, while improving the sensitivity of the search.

The atlas is made possible due to our rigorous analytical upper limit procedure [37]. The analytical formula is required to efficiently produce billions of individual upper limits, and the statistical rigor is a necessity to ensure that any subset of results is individually valid regardless of any peculiarities in the distribution of underlying data. In contrast, many searches [46, 10, 13, 20, 47] use Monte-Carlo simulations that establish upper limits over the whole sky. This method is impractical to use on many sky locations.

The most sensitive population-average all-sky upper limits in [46] did not use Monte-Carlo simulations, but rather an analytical formula. While the authors provide a scaling prescription to derive upper limits at arbitrary sky positions, its validity is not obvious. In particular, it is not clear why the scaling factor derived from average antenna beam pattern function over the population should be representative of the 95% percentile population amplitudes. This confusion may lead to an underestimate of the upper limits. Even for the all-sky case the formula is heuristic, at best. It was derived under the assumption of independent identically distributed Gaussian noise and neglects the variation of the noise level in real data. The deviations from Gaussianity are also ignored which precludes from establishing rigorous upper limits.

A counterpart of our atlas is produced by the radiometer search for the gravitational wave stochastic background [48]. Both the radiometer and this search are affected by loud continuous signals. However, this search is specifically tuned to continuous signals as it, for example, tracks the signal Doppler shifts due to the Earth motion over months of observation. This makes it significantly more sensitive (we estimate by at least a factor of $10$) to continuous waves – and also more computationally expensive – than the radiometer search.

All-sky surveys for continuous gravitational waves are extremely computationally intensive: the results of this search cost in excess of 36 million CPU-hours. Our atlas [45] will allow others to benefit from this investment: astronomical observations of pulsations from a certain location in the sky may be checked against high SNR occurrences at the same location and consistent frequencies; the atlas can immediately be scanned across the whole frequency range at an interesting sky position, for example that of a supernova remnant (as we demonstrate in the appendix), in search of a significant result; with an estimate of the distance of a gravitational wave source, upper limits on its ellipticity can be readily computed; models of neutron star populations can be directly convolved with the atlas upper limits to make predictions of their detectability, and the same holds for primordial black hole populations emitting continuous gravitational waves through their orbital energy loss [49], or for tens of solar mass black hole populations sourcing continuous waves through super-radiance [9]. The fine granularity of our atlas (Figures 4 and 5) produces high-number statistics and, in case of a coincidence with some other observation, this will greatly aid in the estimation of the chance probability of any finding.

Five outliers survive all stages (Table 2). Three correspond to simulated signals [50] injected via radiation pressure on the detectors’ test masses (Table 3).

One is induced by a large instrumental artifact in H1 at 993.300 Hz. This instrumental artifact is quite loud, but is not included in the known line list [39]. The reanalysis using H1 data alone yields an SNR of 16.3, while more sensitive L1 data gives an SNR of only 9.2.

One outlier with SNR just above threshold does not correspond to any identifiable detector artifact, however a reanalysis of this outlier using public LIGO data from the O3a and O3b [40] runs combined, yields a decrease in the SNR – from 16 on O3a data to 13 using O3a+b – with a 6 day coherence length. This is not consistent with what is expected from the IT2 continuous signal model (Eq.  1) assumed in this search. Thus the outlier is either an astrophysical signal that deviates significatly from the IT2 model, or the outlier is due to terrestrial contamination.

*