of 15
Directional limits on persistent gravitational waves
using data from Advanced LIGO’s first two observing runs
The LIGO Scientific Collaboration and The Virgo Collaboration
(Dated: April 19, 2019)
We perform an unmodeled search for persistent, directional gravitational wave (GW) sources using
data from the first and second observing runs of Advanced LIGO. We do not find evidence for any
GW signals. We place limits on the broadband GW flux emitted at 25 Hz from point sources with
a power law spectrum at
F
α,
Θ
<
(0
.
05
25)
×
10
8
erg cm
2
s
1
Hz
1
and the (normalized) energy
density spectrum in GWs at 25 Hz from extended sources at Ω
α
(Θ)
<
(0
.
19
2
.
89)
×
10
8
sr
1
where
α
is the spectral index of the energy density spectrum. These represent improvements of 2
.
5
3
×
over previous limits. We also consider point sources emitting GWs at a single frequency, targeting
the directions of Sco X-1, SN 1987A, and the Galactic Center. The best upper limits on the strain
amplitude of a potential source in these three directions range from
h
0
<
(3
.
6
4
.
7)
×
10
25
, 1.5
×
better than previous limits set with the same analysis method. We also report on a marginally
significant outlier at 36.06 Hz. This outlier is not consistent with a persistent gravitational-wave
source as its significance diminishes when combining all of the available data.
Introduction
— The stochastic gravitational wave
(GW) background (SGWB) is the superposition of many
sources of GWs in the Universe [1]. Anisotropies in the
SGWB can be generated by spatially extended sources
such as a population of neutron stars in the galactic
plane or a nearby galaxy [2, 3], or from perturbations
in statistically-isotropic backgrounds formed at cosmo-
logical distances such as the compact binary background
[4–8] or the background from cosmic strings [9]. Cross-
correlation based methods have been used to search for
the anisotropic background in previous observing runs
[10, 11] of the initial and Advanced Laser Interferomter
Gravitational-wave Observatory (LIGO) [12], and fu-
ture searches will incorporate data from the Advanced
Virgo [13] detector. Using very similar techniques, one
can also search for point sources with an unknown phase
evolution, which could include rotating neutron stars in
the Galaxy [14–17]. Since a SGWB search is by nature
un-modelled, performing the anisotropic SGWB search
allows us to take an eyes-wide-open approach to explor-
ing the GW sky.
In this paper, we present the results of three com-
plementary searches, which probe different types of
anisotropy.
All of the searches are based on cross-
correlation methods; for a review see [18]. A spherical
harmonic decomposition (SHD) of the GW power on the
sky [11, 19] is optimized to search for extended sources on
the sky with a smooth frequency spectrum. The broad-
band radiometer analysis [14, 15] (BBR) is optimized
for detecting resolvable, persistent point-sources emitting
GWs across a wide frequency band. Finally, the directed
narrowband radiometer (NBR) looks at the frequency
spectrum for three astrophysically interesting directions:
Scorpius X-1 (Sco X-1) [20, 21], Supernova 1987A (SN
1987A) [22, 23], and the Galactic Center [24, 25]. We do
not find a significant detection for any of the searches,
and so we place upper limits on the amplitude of the
anisotropic SGWB, and on point sources with broad and
narrow frequency ranges. Our upper limits improve on
the best results from previous runs [10] by approximately
a factor of 2.5-3 for the broadband searches and a factor
of 1.5 for the narrowband searches. For the narrowband
radiometer search, we find a marginally significant out-
lier in the direction of SN 1987A, when analyzing just
the data from LIGO’s second observing run (O2). Its
significance diminishes, however, when including all of
the available data.
Data
— We analyze strain data from the first (O1) and
second (O2) observing runs of Advanced LIGO’s 4 km
detectors in Hanford, Washington (H1) and Livingston,
Louisiana (L1). The O1 data set used here was collected
from 15:00 UTC on 18 September, 2015 to 16:00 UTC
on 12 January, 2016, while the O2 data set was collected
from 16:00:00 UTC on 30 November, 2016 to 22:00:00
UTC on 25 August, 2017. In O2, linearly coupled noise
was removed from the strain time series at H1 and L1
using Wiener filtering [26–30]. The Virgo (V1) detec-
tor started to collect data from August 2017 but does
not contribute significantly to the sensitivity of SGWB
searches in O2, both because its noise level is much higher
than the LIGO detectors and because it ran for a much
shorter period of time. Therefore, we do not include
Virgo in this analysis. We plan, however, to include Virgo
in the analysis of data from future observation runs.
Our data processing methods follow the procedure
used in O1 [10, 31]. First, we down-sample the strain
time series from 16,384 Hz to 4,096 Hz. We then divide
the data into 192 s, 50% overlapping, Hann-windowed
segments, and apply a cascading 16th order Butterworth
digital high-pass filter with a knee frequency of 11 Hz.
We compute the cross correlation of coincident 192 s seg-
ments at both detectors in the frequency domain, and
then coarse-grain to a frequency resolution of 1
/
32 Hz.
Finally, we optimally combine results from those overlap-
ping time segments to produce the final cross-correlation
estimate [32].
arXiv:1903.08844v2 [gr-qc] 18 Apr 2019
2
In order to account for non-Gaussian features in the
data, we remove segments associated with instrumental
artifacts and hardware injections used for signal valida-
tion [33, 34]. Segments containing known GW signals [35]
are also excluded. Finally, we apply a non-stationarity
cut (see, e.g., [36]) to eliminate segments where the power
spectral density of the noise changes on time scales that
are of the same order as the chosen segment length. In
total these cuts removed 16% of the data, leading to a
total search live-time of 99 days from the O2 run. For
our results where we combine data between the O1 and
O2 observing runs we have a total search livetime of 129
days. In addition, frequency bins associated with known
instrumental artifacts are removed [37]. These frequency
domain cuts discarded 4% of the most sensitive frequency
band for the BBR and SHD searches and 15% of the ob-
serving band for the NBR search. The subtraction of lin-
early coupled noise did not introduce any new frequency
domain cuts.
The broadband searches integrate over frequencies be-
tween 20 and 500 Hz. This range accounts for more than
99% of the sensitivity for the power law spectral models
that we use (see Table 1 of [38]). The narrowband analy-
sis searches over the frequency band from 20 to 1726 Hz
using frequency bins of various sizes depending upon fre-
quency and sky direction. The lower edge of this range
is chosen because of increased noise and non-stationarity
at lower frequencies, while the upper edge of the range
is a product of the filter used to resample the data from
16,384 Hz to 4,096 Hz.
Methods
— The anisotropic SGWB background can
be defined in terms of the dimensionless energy density
gw
(
f,
Θ) per unit frequency
f
and solid angle Θ,
gw
(
f,
Θ) =
f
ρ
c
d
3
ρ
GW
d
f
d
2
Θ
,
(1)
where
ρ
c
= 3
H
2
0
c
2
/
(8
πG
) is the critical energy density
needed to have a spatially flat Universe. We take the
Hubble constant to be
H
0
= 67
.
9 km s
1
Mpc
1
[39]. Fol-
lowing past analyses, we assume that we can factorize
gw
into frequency and sky-direction dependent terms,
gw
(
f,
Θ) =
2
π
2
3
H
2
0
f
3
H
(
f
)
P
(Θ)
.
(2)
This quantity has units of the dimensionless energy den-
sity parameter per steradian. For the radiometer searches
it is useful to define a different representation in terms of
energy flux,
F
(
f,
Θ) =
c
3
π
4
G
f
2
H
(
f
)
P
(Θ)
,
(3)
which has units of erg cm
2
s
1
Hz
1
sr
1
, where
c
is the
speed of light and
G
is Newton’s gravitational constant.
We divide the searches into the
broadband
searches
(SHD and BBR), which produce sky maps where the
flux has been integrated over a broad range of frequen-
cies, and the
narrowband
search (NBR), which looks
at the strain amplitude spectrum in a fixed sky di-
rection. For the broadband searches, we typically as-
sume that the energy spectrum has a power law form,
H
(
f
) = (
f/f
ref
)
α
3
, where
α
=
{
0
,
2
/
3
,
3
}
describes a
range of astrophysical and cosmological models [10], and
f
ref
is a reference frequency which we take to be 25 Hz,
as in [10]. The SHD search looks for sources with a large
angular extent. We express the results in terms of the
spherical harmonic decomposition of Ω
gw
(
f,
Θ) assuming
a power-law in frequency of spectral index
α
. We then
report the energy density in each direction at a reference
frequency of 25 Hz, denoted by Ω
α
(Θ).
For the BBR search, we assume that the angular dis-
tribution of the power is localized in a 1 deg
2
pixel,
P
(Θ) =
P
Θ
0
δ
2
,
Θ
0
). The results of the BBR are then
given in terms of the quantity
F
α,
Θ
0
, which is the flux
evaluated at the reference frequency of 25 Hz, assum-
ing a power law, after integrating over solid angle. The
explicit definitions of
F
α,
Θ
0
and Ω
α
(Θ) are given in the
Technical Supplement.
Finally, the NBR search does not integrate over fre-
quency, and attempts to measure the strain amplitude,
h
0
, of a putative monochromatic source in each fre-
quency bin independently. This includes combining ad-
jacent 0.031 Hz frequency bins together to account for
the Doppler modulation due to the motion of the Earth
around the solar system barycenter and any binary mo-
tion of the source itself [10].
The full description of the methods used to search for
an anisotropic SGWB is presented in the Technical Sup-
plement and in the paper describing the analysis of the
Advanced LIGO O1 data. We follow the notation pre-
sented in that Letter [10].
The searches all generally start by estimating the dirty
map
X
ν
, and its corresponding covariance matrix Γ
μν
,
referred to here as the Fisher matrix [10, 19, 40]. The
dirty map represents an estimate of the GW power as
seen through the detector’s beam matrix.
Given the Fisher matrix Γ
I
μν
and dirty map
X
I
ν
, where
I
labels the observing run, we can form a combined Fisher
matrix and dirty map by summing the results from the
two runs, O1 and O2 [18]
Γ
μν
= Γ
(O1)
μν
+ Γ
(O2)
μν
,
X
μ
=
X
(O1)
μ
+
X
(O2)
μ
.
(4)
From the combined Fisher matrix and dirty map, we can
construct estimators of the power on the sky via:
ˆ
P
μ
=
ν
(
Γ
1
R
)
μν
X
ν
.
(5)
In the above equations,
μ
,
ν
label either pixels (i.e., direc-
tions on the sky) or spherical harmonic components—i.e.,
μ
(
lm
), depending on which basis is used to represent
3
All-sky (broadband) Results
Max SNR (%
p
-value)
Upper limit ranges
O1 Upper limit ranges
α
gw
H
(
f
)
BBR
SHD
BBR (
×
10
8
) SHD (
×
10
8
)
BBR (
×
10
8
) SHD (
×
10
8
)
0
constant
f
3
3.09 (9)
2.98 (9)
4.4 – 25
0.78 – 2.90
15 – 65
3.2 – 8.7
2
/
3
f
2
/
3
f
7
/
3
3.09 (20)
2.61 (31)
2.3 – 14
0.64 –2.47
7.9 – 39
2.5 – 6.7
3
f
3
constant
3.27 (66)
3.57 (27)
0.05 – 0.33
0.19 – 1.1
0.14 – 1.1
0.5 – 3.1
TABLE I. Search information for BBR and SHD. On the left side of the table we show the value of the power-law spectral
index,
α
, and the scaling of Ω
gw
and
H
(
f
) with frequency. To the right we show results for the broadband radiometer (BBR)
and spherical harmonic decomposition (SHD) searches for the combined O1 and O2 analysis, as well as the results from O1 for
comparison. We show the maximum SNR across all sky positions for each spectral index, as well as an estimated
p
-value. We
also show the range of 95% upper limits on energy flux set by the BBR search across the whole sky [erg cm
2
s
1
Hz
1
] and
the SHD range of upper limits on normalized energy density across the whole sky [sr
1
]. These limits use data from both O1
and O2. The median improvement across the sky compared to limits set in O1 is 2.6-2.7 for the BBR search and 2.8-3 for the
SHD search, depending on power-law spectral index.
the sky maps. The subscript ‘R’ on the Fisher matrix
means that regularization has been applied (e.g., singu-
lar value decomposition) in order to perform the matrix
inversion [10].
We can also construct an estimate of the angular power
spectrum,
C
l
, for the SGWB from the estimate of the
spherical harmonics coefficients,
ˆ
P
lm
. The
C
l
’s describe
the angular scale of the structure found in the clean
maps [19]
ˆ
C
l
=
(
2
π
2
f
3
ref
3
H
2
0
)
2
1
1 + 2
l
l
m
=
l
[
|
ˆ
P
lm
|
2
1
R
)
lm,lm
]
.
(6)
We have also used theoretical models for the SGWB
from compact binaries [4] and from Nambu-Goto cosmic
strings [9] to check our assumption that the SGWB en-
ergy density Ω
gw
(
f,
Θ) can be factorised into a spectral
shape term and an angular power term. We find that
both models predict
C
l
’s that follow the appropriate fre-
quency power laws across the frequency range in which
the LIGO stochastic searches are most sensitive, thereby
supporting this assumption.
Broadband radiometer and spherical harmonic decom-
position results
— The sky maps for the BBR search are
shown in Figure 1, and for the SHD search in Figure 2.
Converting maps from the spherical harmonics basis (i.e.
μ
= (
lm
)) to the pixel basis is discussed in detail in [19].
Each column indicates a different value of the spectral in-
dex,
α
. The top row shows a map of the signal-to-noise
ratio (SNR) for each sky direction. The SNR sky maps
are consistent with Gaussian noise (see the
p
-values given
in Table I). Consequently, we place upper limits on the
amount of GW power in each pixel using the methods
outlined in [41]. The bottom rows of Figures 1 and 2
show maps of these upper limits for the BBR and SHD
analyses, respectively. The minimum and maximum 95%
confidence upper limits across all pixels for both the BBR
and SHD searches are shown in Table I. These limits rep-
resent a median improvement across the sky of 2.6-2.7 for
the BBR search and 2.8-3 for the SHD search, depending
on the power-law spectral index,
α
.
Limits on angular power spectra
— We also use the
maps from the SHD analysis to set upper limits on the
angular power spectrum components,
C
l
. The upper lim-
its are shown for three spectral indices in Figure 3. The
upper limit for
α
= 2
/
3 can be compared with theoretical
predictions in the literature for the SGWB from compact
binaries [4–6]. In particular, the calculation in Refs. [4, 5]
gives
C
1
/
2
l
3
×
10
11
sr
1
for 1
l
4 (the calculation
in Ref. [6] gives values that are
10
×
smaller). Similarly,
the upper limit for
α
= 0 can be compared with predic-
tions for the SGWB from Nambu-Goto cosmic strings in
Ref. [9], using the same models for the string network as
in Ref. [42]. Assuming the isotropic component of the
cosmic string SGWB is consistent with the upper lim-
its set by LIGO’s second observing run [38], the dipole
(
l
= 1) can be as large as
C
1
/
2
1
10
10
sr
1
, though
the values for higher multipoles
l >
1 are many orders of
magnitude smaller. These predictions are therefore con-
sistent with the upper limits obtained here, and present
an important target for future observing runs.
It has also been recently shown [43] that the finite
sampling of the galaxy distribution and the compact bi-
nary coalescence event rate induce a shot noise in the
anisotropies of the astrophysical GW background, lead-
ing to a scale-invariant bias term in the angular power
spectrum. Such a bias will dominate over the true cosmo-
logical power spectrum, which to be recovered will need
either sufficiently long observing times or subtraction of
the foreground.
Narrowband radiometer results
— The narrowband ra-
diometer search estimates the strain amplitude,
h
0
, of
a potential source of GWs in three different directions.
The maximum SNR across the frequency band and an
estimate of the significance of that SNR for each direc-
tion are shown in Table II. The uncertainty on the fre-
quency for the SNR reported in Table II is a reflection
of the original (uncombined) frequency bin width. The
ephemeris for Scorpius X-1 has been updated since the
publication of [10], and so the search presented below as-
4
FIG. 1. Broadband radiometer maps illustrating a search for point-like sources. The top row shows maps of SNR, while the
bottom row shows maps of the upper limits at 95% confidence on energy flux
F
α,
Θ
0
[erg cm
2
s
1
Hz
1
]. Three different
power-law indices,
α
= 0
,
2
/
3 and 3, are represented from left to right. The
p
-values associated with the maximum SNR are
(from left to right)
p
= 9%
, p
= 20%
, p
= 66% (see Table I).
FIG. 2. All-sky maps reconstructed from a spherical harmonic decomposition. This search is optimized for extended sources,
and the plots above show SNR (top) and upper limits at 95% confidence on the energy density of the SGWB Ω
α
[ sr
1
]
(bottom). Results for three different power-law spectral indices,
α
= 0
,
2
/
3 and 3 are shown from left to right. These three
different sets of maps have an
l
max
of 3, 4, and 16 respectively. The
p
-values associated with the maximum SNR are (from left
to right)
p
= 9%
, p
= 31%
, p
= 27% (see Table I).
sumes a projected semi-major axis,
a
0
, in the center of
the range presented by [44].
In the direction of Sco X-1 and the Galactic Center,
the maximum SNR is consistent with what one expects
from Gaussian noise. In the direction of SN 1987A, there
is a frequency bin with a 1-sided, single-direction
p
-value
1.7% at 181.8 Hz. This
p
-value includes a trials factor for
the number of the number of frequency bins in the anal-
ysis. Under the assumption that we search over three
independent directions, an extra trials factor would be
applied and this
p
-value rises to 5%. Therefore, we find
no compelling evidence for GWs from the analysis that
combines frequency bins together. We set 95% upper lim-
its on the strain amplitude of a putative sinusoidal grav-
itational wave signal,
h
0
, in each individual frequency
bin, taking into account any Doppler modulation in the
signal as well as marginalizing over inclination angle and
polarization angle of the source [10]. These limits, along
with the 1
σ
sensitivity of the search, are shown in Fig-
ure 4. To avoid reporting our best limits from downward
fluctuations of noise, we take a running median over each
1 Hz frequency band and report the best limit on
h
0
and
the frequency band of that limit in Table II.
The best limits on Sco-X1 set in this paper are higher
than the best limit set in O1 using a model-based cross-
correlation method [20], and are now lower than those set
using hidden Markov model tracking [21]. The torque-
balance limit, set by assuming that torque due to ac-
5
FIG. 3. Upper limits on
C
l
’s at 95% confidence for the SHD
analyses for
α
= 0 (top, black triangles),
α
= 2
/
3 (middle, red
circles) and
α
= 3 (bottom, blue squares). These represent
an improvement in upper limits over O1 of 2.5 – 3 depending
on spectral index,
α
, and
l
.
cretion is equal to the braking torque due to GW emis-
sion, is still around a factor of 5 lower than the limits
set in this paper. The best limits on
h
0
in the direc-
tion of the Galactic Center and SN 1987A are generally
higher than previous modeled searches for isolated neu-
tron stars [22, 24, 25, 45, 46]. It is important to note that
the search presented in this paper is inherently unmod-
eled, meaning it makes no assumption about the phase
evolution of a potential signal past timescales of 192 s.
Outlier at 36.06
Hz
in the O2 data
— In the process of
performing the narrowband radiometer search, a natural
intermediate step of the analysis is to look directly at the
0.03125 Hz bins for the O2 data, before combining with
O1 and before combining over adjacent bins to account
for Doppler modulation. We call these “sub-bins”. For
this intermediate data product, the maximum SNRs for
the Galactic Center, Sco X-1, and SN 1987A are 4.6, 4.3,
and 5.3, respectively. These first two values correspond
to
p
-values greater than 5%, consistent with Gaussian
noise. But for SN 1987A, the maximum SNR of 5
.
3 at
36.0625 Hz has a corresponding
p
-value of 0
.
27%, or 3
σ
,
which is marginally significant.
Assuming that the maximum SNR is due to a pul-
sar which is spinning down due to GW emission, we can
relate the observed strain
h
0
= 7
.
3
×
10
25
(assuming cir-
cular polarization) at
f
= 36
.
06 Hz to other parameters
describing the pulsar:
h
0
=
4
πG
c
4
I
z
f
2
r
,
̇
f
=
G
5
πc
5

2
I
z
(2
πf
)
5
.
(7)
We use a fiducial value for the moment of inertia
I
z
= 10
39
kg
·
m
2
.
If the source is associated with
SN 1987A, then the distance to Earth is approximately
r
= 51 kpc [47, 48], leading to an ellipticity

= 3
×
10
2
and spin down
̇
f
=
7
.
7
×
10
8
Hz
/
s. But this value
of the spin down parameter is inconsistent with the fact
that the signal is seen in only one frequency bin. For
the signal to remain in a single frequency bin, we need
r
.
1 kpc (corresponding to
̇
f
=
2
.
9
×
10
11
Hz
/
s),
but the ellipticity

= 5
×
10
4
is still much larger than
that predicted for typical pulsars. So the signal does not
appear to be consistent with GW emission from a pul-
sar. In addition, an all-sky search for continuous GWs
has set limits more stringent than our estimate of
h
0
for
this outlier [46].
Using the techniques described in [37], we have not
been able to identify a coherent instrumental witness
channel that would explain this large SNR. But the fact
that the sky direction of the maximum SNR is close to
the equatorial pole is consistent with the behavior of in-
strumental noise lines, since the equatorial poles have no
sidereal-time modulation. The signal appears to turn on
during O2, with the SNR exceeding 1 on March 13th,
2017, as shown in Figure 5, but it does not exhibit any
significant short-term non-stationarity biasing the esti-
mate of the cross correlation. This turn-on feature of the
cumulative SNR is not evidence of a real signal, however,
as we have performed simulations of Gaussian noise con-
ditioned on getting a maximum SNR
5, and have found
examples where a turn-on like this can be produced. In
addition, upon combining O2 and O1 data together, the
SNR of this frequency bin is reduced to 4.7, which cor-
responds to a
p
-value of 10%, which is consistent with
noise.
Conclusions
— We have placed upper limits on the
anisotropic SGWB using three complementary methods.
In each case we do not find conclusive evidence for a GW
signal, and so we place upper limits by combining data
from Advanced LIGO’s first and second observing runs.
A marginal outlier at a frequency of 36.06 Hz was seen
by the narrowband radiometer search in O2 in the di-
rection of SN 1987A; however it does not appear in the
combined O1+O2 data and is not consistent with a per-
sistent signal. We will continue to monitor this particu-
lar frequency bin during the next observing run, taking
advantage of the greater confidence that comes with in-
creased observation periods and more sensitive detectors.
In the future, the anisotropic searches will include data
from Advanced Virgo as well, and can be used to study
specific astrophysical models. Additionally, new algo-
rithms can take advantage of folded data to produce a
wider search of every frequency and sky position [49–52].
Acknowledgements
— The authors gratefully acknowl-
edge 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 construc-
tion and operation of the GEO600 detector. Additional
support for Advanced LIGO was provided by the Aus-
tralian Research Council. The authors gratefully ac-
knowledge the Italian Istituto Nazionale di Fisica Nucle-