Comprehensive All-sky Search for Periodic Gravitational Waves in the Sixth Science Run LIGO Data
We report on a comprehensive all-sky search for periodic gravitational waves in the frequency band 100-1500 Hz and with a frequency time derivative in the range of Hz/s. Such a signal could be produced by a nearby spinning and slightly non-axisymmetric isolated neutron star in our galaxy. This search uses the data from the Initial LIGO sixth science run and covers a larger parameter space with respect to any past search. A Loosely Coherent detection pipeline was applied to follow up weak outliers in both Gaussian (% recovery rate) and non-Gaussian (% recovery rate) bands. No gravitational wave signals were observed, and upper limits were placed on their strength. Our smallest upper limit on worst-case (linearly polarized) strain amplitude is near 169 Hz, while at the high end of our frequency range we achieve a worst-case upper limit of . Both cases refer to all sky locations and entire range of frequency derivative values.
In this paper we report the results of a comprehensive all-sky search for continuous, nearly monochromatic gravitational waves in data from LIGO’s sixth science (S6) run. The search covered frequencies from 100 Hz through 1500 Hz and frequency derivatives from Hz/s through Hz/s.
A number of searches for periodic gravitational waves have been carried out previously in LIGO data (2); (3); (4); (5); (6); (7); (8); (9); (10); (11), including coherent searches for gravitational radiation from known radio and X-ray pulsars. An Einstein@Home search running on the BOINC infrastructure (12) has performed blind all-sky searches on S4 and S5 data (13); (14); (15).
The results in this paper were produced with the PowerFlux search program. It was first described in (2) together with two other semi-coherent search pipelines (Hough, Stackslide). The sensitivities of all three methods were compared, with PowerFlux showing better results in frequency bands lacking severe spectral artifacts. A subsequent article (4) based on the data from the S5 run featured improved upper limits and a systematic outlier follow-up search based on the Loosely Coherent algorithm (16).
The analysis of the data set from the sixth science run described in this paper has several distinguishing features from previously published results:
A number of upgrades to the detector were made in order to field-test the technology for Advanced LIGO interferometers. This resulted in a factor of about two improvement in intrinsic noise level at high frequencies compared to previously published results (4).
The higher sensitivity allowed us to use less data while still improving upper limits in high frequency bands by % over previously published results. This smaller dataset allowed covering larger parameter space, and comprehensive exploration of high frequency data.
This search improved on previous analyses by partitioning the data in month chunks and looking for signals in any contiguous sequence of these chunks. This enables detections of signals that conform to ideal signal model over only part of the data. Such signals could arise because of a glitch, or because of influence of a long-period companion object.
The upgrades to the detector, while improving sensitivity on average, introduced a large number of detector artifacts, with % of frequency range contaminated by non-Gaussian noise. We addressed this issue by developing a new Universal statistic (17) that provides correct upper limits regardless of the noise distribution of the underlying data, while still showing close to optimal performance for Gaussian data.
We have observed no evidence of gravitational radiation and have established the most sensitive upper limits to date in the frequency band 100-1500 Hz. Our smallest 95% confidence level upper limit on worst-case (linearly polarized) strain amplitude is near 169 Hz, while at the high end of our frequency range we achieve a worst-case upper limit of . Both cases refer to all sky locations and entire range of frequency derivative values.
Ii LIGO interferometers and S6 science run
The LIGO gravitational wave network consists of two observatories, one in Hanford, Washington and the other in Livingston, Louisiana, separated by a 3000 km baseline. During the S6 run each site housed one suspended interferometer with 4 km long arms.
While the sixth science run spanned a months period of data acquisition, this analysis used only data from GPS 951534120 (2010 Mar 02 03:01:45 UTC) through GPS 971619922 (2010 Oct 20 14:25:07 UTC), for which strain sensitivity was best. Since interferometers sporadically fall out of operation (“lose lock”) due to environmental or instrumental disturbances or for scheduled maintenance periods, the dataset was not contiguous. The Hanford interferometer H1 had a duty factor of 53%, while the Livingston interferometer L1 had a duty factor of 51%. The strain sensitivity was not uniform, exhibiting a % daily variation from anthropogenic activity as well as gradual improvement toward the end of the run (18); (19).
Non-stationarity of noise was especially severe at frequencies below Hz, and since the average detector sensitivity for such frequencies was not significantly better than that observed in the longer S5 run (4), this search was restricted to frequencies above Hz.
A detailed description of the instruments and data can be found in (20).
Iii The search for continuous gravitational radiation
In this paper we assume a classical model of a spinning neutron star with a rotating quadrupole moment that produces circularly polarized gravitational radiation along the rotation axis and linearly polarized radiation in the directions perpendicular to the rotation axis. The linear polarization is the worst case as such signals contribute the smallest amount of power to the detector.
The strain signal template is assumed to be
where and characterize the detector responses to signals with “” and “” quadrupolar polarizations (2); (3); (4), the sky location is described by right ascension and declination , the inclination of the source rotation axis to the line of sight is denoted , and the phase evolution of the signal is given by the formula
with being the source frequency and denoting the first frequency derivative (which, when negative, is termed the spindown). We use to denote the time in the Solar System barycenter frame. The initial phase is computed relative to reference time . When expressed as a function of local time of ground-based detectors the equation 2 acquires sky-position-dependent Doppler shift terms. We use to denote the polarization angle of the projected source rotation axis in the sky plane.
The search has two main components. First, the main PowerFlux algorithm (2); (3); (4); (22); (23); (24) was run to establish upper limits and produce lists of outliers with signal-to-noise ratio (SNR) greater than 5. Next, the Loosely Coherent detection pipeline (16); (25); (4) was used to reject or confirm collected outliers.
Both algorithms calculate power for a bank of signal model templates and compute upper limits and signal-to-noise ratios for each template based on comparison to templates with nearby frequencies and the same sky location and spindown. The input time series is broken into % overlapping s long segments which are Hann windowed and Fourier transformed. The resulting short Fourier transforms (SFTs) are arranged into an input matrix with time and frequency dimensions. The power calculation can be expressed as a bilinear form of the input matrix :
Here denotes the detector frame frequency drift due to the effects from both Doppler shifts and the first frequency derivative. The sum is taken over all times corresponding to the midpoint of the short Fourier transform time interval. The kernel includes the contribution of time dependent SFT weights, antenna response, signal polarization parameters and relative phase terms (16); (25).
The main semi-coherent PowerFlux algorithm uses a kernel with main diagonal terms only and is very fast. The Loosely Coherent algorithms increase coherence time while still allowing for controlled deviation in phase (16). This is done by more complicated kernels that increase effective coherence length.
The effective coherence length is captured in a parameter , which describes the amount of phase drift that the kernel allows between SFTs, with corresponding to a fully coherent case, and corresponding to incoherent power sums.
Depending on the terms used, the data from different interferometers can be combined incoherently (such as in stages 0 and 1, see Table 2) or coherently (as used in stages 2, 3 and 4). The coherent combination is more computationally expensive but provides much better parameter estimation.
The upper limits (Figure 1) are reported in terms of the worst-case value of (which applies to linear polarizations with ) and for the most sensitive circular polarization ( or ). As described in the previous paper (4), the pipeline does retain some sensitivity, however, to non-general-relativity GW polarization models, including a longitudinal component, and to slow amplitude evolution.
|First harmonic of violin modes||343.25-343.75 Hz, 347-347.25 Hz|
|Second harmonic of violin modes||686.25-687.5 Hz|
|Third harmonic of violin modes||1031.00-1031.25 Hz|
The 95% confidence level upper limits (see Figure 1) produced in the first stage are based on the overall noise level and largest outlier in strain found for every template in each 0.25 Hz band in the first stage of the pipeline. The Hz bands are analyzed by separate instances of PowerFlux (4). A followup search for detection is carried out for high-SNR outliers found in the first stage. Certain frequency ranges (Table 1) were excluded from the analysis because of gross contamination by detector artifacts.
iii.2 Universal statistics
The detector sensitivity upgrades introduced many artifacts, so that in % of the sensitive frequency range the noise follows non-Gaussian distributions which, in addition, are unknown. As the particular non-Gaussian distribution can vary widely depending on particular detector artifacts, the ideal estimate based on full knowledge of the distribution is not usually available. In the previous analysis (2); (3); (4), the frequency bands where the noise distribution was non-Gaussian were not used to put upper limits. However, in the present case this approach would have resulted in excluding most of the frequency bands below 400 Hz and many above 400 Hz; even though the average strain sensitivity in many of these frequency bands is better than in the past.
To make use of the entire spectrum, we used in this work the Universal statistic algorithm (17) for establishing upper limits. The algorithm is derived from the Markov inequality and shares its independence from the underlying noise distribution. It produces upper limits less than % above optimal in case of Gaussian noise. In non-Gaussian bands it can report values larger than what would be obtained if the distribution were known, but the upper limits are always at least 95% valid. Figure 2 shows results of an injection run performed as described in (4). Correctly established upper limits are above the red line.
|Stage||Instrument sum||Phase coherence||Spindown step||Sky refinement||Frequency refinement||SNR increase|
|Initial/upper limit semi-coherent||NA||NA|
iii.3 Detection pipeline
The detection pipeline used in (4) was extended with additional stages (see Table 2) to winnow the larger number of initial outliers, expected because of non-Gaussian artifacts and larger initial search space. This detection pipeline was also used in the search of the Orion spur (5).
The initial stage (marked 0) scans the entire sky with semi-coherent algorithm that computes weighted sums of powers of s Hann-windowed SFTs. These power sums are then analyzed to identify high-SNR outliers. A separate algorithm uses universal statistics (17) to establish upper limits.
The entire dataset was partitioned into 7 segments of equal length and power sums were produced independently for any contiguous combinations of these stretches. As in (5) the outlier identification was performed independently in each stretch.
High-SNR outliers were subject to a coincidence test.Â For each outlier with in the combined H1 and L1 data, we required there to be outliers in the individual detector data that had , matching the parameters of the combined-detector outlier within a distance of on the sky, mHz in frequency, and Hz/s in spindown.Â However, the combined-detector SNR could not be lower than either single-detector SNR.
The identified outliers using combined data are then passed to followup stage using Loosely Coherent algorithm (16) with progressively tighter phase coherence parameters , and improved determination of frequency, spindown and sky location.
As the initial stage 0 only sums powers it does not use relative phase between interferometers, which results in some degeneracy between sky position, frequency and spindown. The first Loosely Coherent followup stage also combines interferometer powers incoherently, but demands greater temporal coherence (smaller ) within each interferometer, which should boost SNR of viable outiers by at least 20%. Subsequent stages use data coherently providing tighter bounds on outlier location.
The testing of the pipeline was done above Hz and included both Gaussian and non-Gaussian bands. We focused on high frequency performance because preliminary S6 data indicated the sensitivity at low frequencies was unlikely to improve over S5 results due to detector artifacts.
The followup code was tested to recover % of injections % above the upper limit level assuming uniform distribution of injection frequency. (Figure 3). Recovery of signals injected into frequency bands which exhibits non-Gaussian noise was % (Figure 4). Our recovery criterion demanded that an outlier close to the true injection location (within mHz in frequency , Hz/s in spindown and radHz in sky location) be found and successfully pass through all stages of the detection pipeline. As each stage of the pipeline only passes outliers with an increase in SNR, this resulted in an outlier that strongly stood out above the background, with good estimates of the parameters of the underlying signal.
It should be noted that the injection recovery curve in Figure 3 passes slightly below the % level for equal to the upper limit. However, the upper limits are based on power levels measured by stage 0, independent of any follow-up criteria.Â That is, we can say with % confidence that a signal above the upper limit level is inconsistent with the observed power, even though such a (hypothetical) signal might not pass all of our follow-up criteria to be “detected”.Â The main reason that these injections fail to be detected is the different sensitivities of the H1 and L1 detectors.Â When one interferometer is less sensitive sensitive we can still set a good upper limit, but the initial coincidence criteria requires that an outlier be marginally seen in both interferometers. In the previous analysis (4) the interferometers had similar sensitivity and the curve passed through the intersection of the green lines (horizontal axis value of , vertical axis value of %).
iii.4 Gaussian false alarm event rate
The computation of the false alarm rate for the outliers passing all stages of the pipeline is complicated by the fact that most outliers are caused by instrumental artifacts for which we do not know the underlying probability distribution. In principle, one could repeat the analysis many times using non-physical frequency shifts (which would exclude picking up a real signal by accident) in order to obtain estimates of false alarm rate, but this approach is very computationally expensive. Even assuming a perfect Gaussian background, it is difficult to analytically model the pipeline in every detail to obtain an accurate estimate of the false alarm rate, given the gaps in interferometer operations and non-stationary noise.
Instead, following (5), we compute a figure of merit that overestimates the actual Gaussian false alarm event rate. We simplify the problem by assuming that the entire analysis was carried out with the resolution of the very last stage of follow-up and we are merely triggering on the SNR value of the last stage. This is extremely conservative as we ignore the consistency requirements that allow the outlier to proceed from one stage of the pipeline to the next; the actual false alarm rate could be lower.
The SNR of each outlier is computed relative to the Loosely Coherent power sum for 501 frequency bins spaced at Hz intervals (including the outlier) but with all the other signal parameters held constant. The spacing assures that correllations between neighboring sub-bins do not affect the statistics of the power sum.
To simplify computation we assume that we are dealing with a simple distribution with the number of degrees of freedom given by the timebase divided by the coherence length and multiplied by a conservative duty factor reflecting interferometer uptime and the worst-case weights from linearly-polarized signals.
Thus to find the number of degrees of freedom we will use the formula
with the duty factor taken to be and giving the phase coherence parameter of the Loosely Coherent search. The duty factor was chosen to allow for only % interferometer uptime and only one quarter of the data receiving high weights from our procedure, which weights the contribution of data inversely as the square of the estimated noise (22); (23).
Thus we define the outlier figure of merit describing Gaussian false alarm (GFA) event rate as
where defines the number of degrees of freedom as given by equation 4, gives the probability for a distribution with degrees of freedom to exceed , and is the estimated number of templates.
We point out that the GFA is overly conservative when applied to frequency bands with Gaussian noise, but is only loosely applicable to bands with detector artifacts, which can affect both the estimate of the number of degrees of freedom of the underlying distribution and the assumption of uncorrelated underlying noise.
The PowerFlux algorithm and Loosely Coherent method compute power estimates for gravitational waves in a given frequency band for a fixed set of templates. The template parameters usually include frequency, first frequency derivative and sky location.
Since the search target is a rare monochromatic signal, it would contribute excess power to one of the frequency bins after demodulation. The upper limit on the maximum excess relative to the nearby power values can then be established. For this analysis we use a universal statistic (17) that places conservative 95% confidence level upper limits for an arbitrary statistical distribution of noise power. The universal statistic has been designed to provide close to optimal values in the common case of Gaussian distribution.
Most natural sources are expected to have negative first frequency derivative, as the energy lost in gravitational or electromagnetic waves would make the source spin more slowly. The frequency derivative can be positive when the source is affected by a strong slowly-variable Doppler shift, such as due to a long-period orbit.
The large gap in data taking due to installation of Advanced LIGO interferometers provided an opportunity to cover an extended parameter space (Figure 5). With respect to previous searches, we have chosen to explore comprehensively both negative and positive frequency derivatives to avoid missing any unexpected sources in our data.
The upper limits obtained in the search are shown in figure 1. The numerical data for this plot can be obtained separately (21). The upper (yellow) curve shows the upper limits for a worst-case (linear) polarizations when the smallest amount of gravitational energy is projected towards Earth. The lower curve shows upper limits for an optimally oriented source. Because of the day-night variability of the interferometer sensitivity due to anthropogenic noise, the linearly polarized sources are more susceptible to detector artifacts, as the detector response to such sources varies with the same period. The neighborhood of 60 Hz harmonics is shown as circles for worst case upper limits and dots for circular polarization upper limits. Thanks to the use of universal statistic they do represent valid values even if contaminated by human activity.
Each point in figure 1 represents a maximum over the sky: only a small excluded portion of the sky near ecliptic poles that is highly susceptible to detector artifacts, due to stationary frequency evolution produced by the combination of frequency derivative and Doppler shifts. The exclusion procedure is described in (4) and applied to % of the sky over the entire run.
A few frequency bands shown in Table 1 were so contaminated that every SFT was vetoed by data conditioning code and the analysis terminated before reaching universal statistic stage. While the universal statistic could have established upper limits with veto turned off, we opted to simply exclude these bands, as the contamination would raise upper limits to be above physically interesting values.
If one assumes that the source spindown is solely due to emission of gravitational waves, then it is possible to recast upper limits on source amplitude as a limit on source ellipticity. Figure 6 shows the reach of our search under different assumptions on source distance. Superimposed are lines corresponding to sources of different ellipticities.
|Hardware injection ip8|
|Hardware injection ip3,|
|Non Gaussian, disturbed H1 spectrum|
|Hardware injection ip2|
|Non Gaussian, Line in H1, disturbed spectrum in L1|
|Induced by loud hardware injection ip4,|
|Non Gaussian, highly disturbed H1+L1 spectra|
|Highly disturbed H1 spectrum, stationary line area|
|Line in H1 at 566.085 Hz|
|Disturbed H1 and L1 spectrum|
|Hardware injection ip7, sloping H1 and L1 spectra|
|Highly disturbed H1 spectrum, stationary line area|
|Induced by loud hardware injection ip8|
|Lines in H1 and L1, Non Gaussian|
|Highly disturbed H1 spectrum|
|Induced by loud hardware injection ip4,|
|Non Gaussian, highly disturbed H1+L1 spectra|
|Highly disturbed H1 spectrum|
|Line in H1 at 451.5 Hz|
The detection pipeline produced 16 outliers (Table 3). Each outlier is identified by a numerical index. We report SNR, decimal logarithm of Gaussian false alarm rate, as well as frequency, spindown and sky location.
The “Segment” column describes the persistence of the outlier through the data, and specified which contiguous subset of the 7 equal partitions of the timespan contributed most significantly to the outlier: see (5) for details.Â A continuous signal will normally have [0,6] in this column (similar contribution from all 7 segments), or on rare occasions [0,5] or [1,6].Â Any other range is indicative of a non-continuous signal or artifact.
During the S6 run several simulated pulsar signals were injected into the data by applying a small force to the interferometer mirrors. Several outliers were due to such hardware injections (Table 4). The full list of injections including those too weak to be found by an all-sky search can be found in (26). The hardware injection ip3 was exceptionally strong with a clear signature even in non-Gaussian band.
The recovery of the injections gives us confidence that no potential signal were missed. Manual followup has shown non-injection outliers to be caused by pronounced detector artifacts.
We have performed the most sensitive all-sky search to date for continuous gravitational waves in the range 100-1500 Hz. We explored both positive and negative spindowns and placed upper limits on expected and unexpected sources. At the highest frequencies we are sensitive to neutron stars with an equatorial ellipticity as small as and as far away as kpc for favorable spin orientations. The use of a universal statistic allowed us to place upper limits on both Gaussian and non-Gaussian frequency bands. The maximum ellipticity a neutron star can theoretically support is at least according to (27); (28). Our results exclude such maximally deformed pulsars above Hz pulsar rotation frequency ( Hz gravitational frequency) within kpc.
A detection pipeline based on a Loosely Coherent algorithm was applied to outliers from our search. This pipeline was demonstrated to be able to detect simulated signals at the upper limit level for both Gaussian and non-Gaussian bands. Several outliers passed all stages of the coincidence pipeline; their parameters are shown in table 3. However, manual examination revealed no true pulsar signals.
The authors gratefully acknowledge the support of the United States National Science Foundation (NSF) for the construction and operation of the LIGO Laboratory and Advanced LIGO as well as the Science and Technology Facilities Council (STFC) of the United Kingdom, the Max-Planck-Society (MPS), and the State of Niedersachsen/Germany for support of the construction of Advanced LIGO and construction and operation of the GEO600 detector. Additional support for Advanced LIGO was provided by the Australian Research Council. The authors gratefully acknowledge the Italian Istituto Nazionale di Fisica Nucleare (INFN), the French Centre National de la Recherche Scientifique (CNRS) and the Foundation for Fundamental Research on Matter supported by the Netherlands Organisation for Scientific Research, for the construction and operation of the Virgo detector and the creation and support of the EGO consortium. The authors also gratefully acknowledge research support from these agencies as well as by the Council of Scientific and Industrial Research of India, Department of Science and Technology, India, Science & Engineering Research Board (SERB), India, Ministry of Human Resource Development, India, the Spanish Ministerio de Economía y Competitividad, the Conselleria d’Economia i Competitivitat and Conselleria d’Educació, Cultura i Universitats of the Govern de les Illes Balears, the National Science Centre of Poland, the European Commission, the Royal Society, the Scottish Funding Council, the Scottish Universities Physics Alliance, the Hungarian Scientific Research Fund (OTKA), the Lyon Institute of Origins (LIO), the National Research Foundation of Korea, Industry Canada and the Province of Ontario through the Ministry of Economic Development and Innovation, the Natural Science and Engineering Research Council Canada, Canadian Institute for Advanced Research, the Brazilian Ministry of Science, Technology, and Innovation, Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP), Russian Foundation for Basic Research, the Leverhulme Trust, the Research Corporation, Ministry of Science and Technology (MOST), Taiwan and the Kavli Foundation. The authors gratefully acknowledge the support of the NSF, STFC, MPS, INFN, CNRS and the State of Niedersachsen/Germany for provision of computational resources.
This document has been assigned LIGO Laboratory document number LIGO-P1500219-v19.
- B. Abbott et al. (LIGO Scientific Collaboration), Phys. Rev. D 77, 022001 (2008).
- B. P. Abbott et al. (LIGO Scientific Collaboration), Phys. Rev. Lett. 102, 111102 (2009).
- B. Abbott et al. (The LIGO and Virgo Scientific Collaboration), Phys. Rev. D 85, 022001 (2012)
- J. Aasi et al. (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. D 93, 042006 (2016).
- B. Abbott et al. (LIGO Scientific Collaboration), M. Kramer, and A. G. Lyne, Phys. Rev. Lett. 94, 181103 (2005).
- B. Abbott et al. (LIGO Scientific Collaboration), M. Kramer, and A. G. Lyne, Phys. Rev. D 76, 042001 (2007).
- B. Abbott et al. (LIGO Scientific Collaboration), Phys. Rev. D 76, 082001 (2007).
- B. Abbott et al. (LIGO Scientific Collaboration), Astrophys. J. Lett. 683, 45 (2008).
- B. P. Abbott et al. (LIGO Scientific Collaboration and Virgo Collaboration), Astrophys. J. 713, 671 (2010).
- J. Abadie et al. (LIGO Scientific Collaboration), Astrophys. J. 722, 1504 (2010).
- The Einstein@Home project is built upon the BOINC (Berkeley Open Infrastructure for Network Computing) architecture described at http://boinc.berkeley.edu/.
- B. Abbott et al. (LIGO Scientific Collaboration), Phys. Rev. D 79, 022001 (2009).
- B. P. Abbott et al. (LIGO Scientific Collaboration), Phys. Rev. D 80, 042003 (2009).
- B. P. Abbott et al. (LIGO Scientific Collaboration), Phys. Rev. D 87, 042001 (2013).
- V. Dergachev, Class. Quantum Grav. 27, 205017 (2010).
- V. Dergachev, Phys. Rev. D 87, 062001 (2013).
- B. Abbott et al. (LIGO Scientific Collaboration), Rep. Prog. Phys. 72, 076901 (2009).
- J. Abadie et al. (LIGO Scientific Collaboration and Virgo Collaboration), arXiv:1203.2674 [gr-qc], http://arxiv.org/abs/1203.2674
- J. Aasi et al. (LIGO Scientific Collaboration and Virgo Collaboration), arXiv:1410.7764 [gr-qc], http://arxiv.org/abs/1410.7764
- See EPAPS Document No. [number will be inserted by publisher] for numerical values of upper limits.
- V. Dergachev, LIGO technical document LIGO-T050186 (2005), available in https://dcc.ligo.org/
- V. Dergachev, LIGO technical document LIGO-T1000272 (2010), available in https://dcc.ligo.org/
- V. Dergachev and K. Riles, LIGO Technical Document LIGO-T050187 (2005), available in https://dcc.ligo.org/
- V. Dergachev, Phys. Rev. D 85, 062003 (2012)
- J. Aasi et al. (LIGO Scientific Collaboration and Virgo Collaboration), The Astrophysical Journal 813 1 (2015)
- C. J. Horowitz and K. Kadau, Phys. Rev. Lett. 102, 191102 (2009).
- N. K. Johnson-McDaniel and B. J. Owen, Phys. Rev. D 88, 044004 (2013)