All-sky search for gravitational wave emission from scalar boson clouds
around spinning black holes in LIGO O3 data
R. Abbott
etal.
*
(The LIGO Scientific Collaboration, the Virgo Collaboration, and the KAGRA Collaboration)
(Received 30 November 2021; accepted 13 April 2022; published 9 May 2022)
This paper describes the first all-sky search for long-duration, quasimonochromatic gravitational-wave
signals emitted by ultralight scalar boson clouds around spinning black holes using data from the third
observing run of Advanced LIGO. We analyze the frequency range from 20 to 610 Hz, over a small
frequency derivative range around zero, and use multiple frequency resolutions to be robust towards
possible signal frequency wanderings. Outliers from this search are followed up using two different
methods, one more suitable for nearly monochromatic signals, and the other more robust towards frequency
fluctuations. We do not find any evidence for such signals and set upper limits on the signal strain
amplitude, the most stringent being
≈
10
−
25
at around 130 Hz. We interpret these upper limits as both an
“
exclusion region
”
in the boson mass/black hole mass plane and the maximum detectable distance for a
given boson mass, based on an assumption of the age of the black hole/boson cloud system.
DOI:
10.1103/PhysRevD.105.102001
I. INTRODUCTION
Theories of beyond standard model physics allow for the
production of ultralight bosons that could constitute a
portion or all of dark matter
[1
–
6]
. If such ultralight bosons
exist, they could appear around rotating black holes due to
quantum fluctuations
[7]
. They would then scatter off and
extract angular momentum from these black holes, and
form macroscopic clouds, i.e. boson condensates, through a
superradiance process
[7,8]
. The energy structure of the
clouds resembles that of a hydrogen atom, earning boson
clouds the name
“
gravitational atoms
”
in the sky. After the
black hole spin decreases below a threshold, the super-
radiance process stops and the cloud depletes over time
mainly through the annihilation of bosons into gravitons,
which generates quasimonochromatic, long-duration gravi-
tational waves when the self-interaction for bosons is weak.
The signal would also have a small, positive frequency
variation, known as a
“
spin-up,
”
that would arise from the
classical contraction of the cloud over time as it loses mass
[7]
. The values of the spin-up depend on whether we
consider scalar
[9]
, vector
[10,11]
, or tensor
[12]
boson
clouds, but in this search we consider scalar boson
cloud signals with very small spin-ups, of maximum
O
ð
10
−
12
Þ
Hz
=
s. Recently, the effect on gravitational-wave
emission from boson self-interactions has been studied for
the scalar case
[13]
, showing that the signal can be
significantly affected, including the magnitude of the
spin-up, as the self-interaction increases.
The gravitational-wave signal frequency depends pri-
marily on the mass of the boson, and weakly on the mass
and spin of the black hole. Recently, exclusion regions of
the black hole/scalar boson mass parameter space were
calculated using upper limits from an all-sky search for
spinning neutron stars in LIGO/Virgo data from the second
observing run (O2)
[14]
. Furthermore, a search for boson
clouds from Cygnus X-1 was performed with a hidden
Markov model (HMM) method
[15]
on the same dataset,
and disfavored particular boson masses for this object
[16]
.
It has been suggested
[17]
that current methods for all-sky
searches, such as the one used in
[14]
, could suffer a
significant sensitivity penalty in the extreme case of a
population of
10
8
black holes in the Milky Way, due to the
superposition of many signals in a small frequency range
(as a consequence of the signal
’
s frequency dependence
at leading order on the boson mass). A detailed study
[18]
has quantified this effect, showing that the actual average
loss of sensitivity is at most of about 15% for a signal
“
density
”
of 1
–
3 signals per frequency bin, while it
rapidly reduces, and becomes negligible, for both lower
and higher densities.
In addition to quasimonochromatic gravitational waves
emitted by boson clouds around isolated black holes, it is
also possible to probe the existence of boson clouds in
binary black hole coalescences. One avenue requires
measurements of spin: in principle, boson clouds will
extract mass and spin from black holes, resulting in low
spin values for black holes compared to those in a universe
without boson clouds; in practice, individual black hole
spins are hard to measure
[19,20]
, and the current spin
distribution of black holes depends on both the mass of
*
Full author list given at the end of the article.
PHYSICAL REVIEW D
105,
102001 (2022)
2470-0010
=
2022
=
105(10)
=
102001(28)
102001-1
© 2022 American Physical Society
boson clouds and the distribution of spins when the black
holes were born
[21,22]
. It is therefore interesting to
combine spin measurements from various black hole
mergers, via hierarchical Bayesian inference, to obtain
constraints on boson cloud/spin interactions
[23,24]
.
Additionally, the presence of a boson cloud will induce
a change in the waveforms used in compact binary
searches, e.g. new multipole moments and tidal effects
(parametrized as Love numbers)
[25,26]
, and such
differences in binary black hole mergers may actually be
detectable depending on the boson mass and field strength
[27]
. It has also been proposed to combine multiband
observations of black hole mergers and boson cloud signals
when future detectors are online
[28]
, and to search for a
stochastic gravitational-wave background from scalar and
vector boson cloud signals
[29
–
32]
. Complementary meth-
ods
[33,34]
and searches
[35,36]
for ultralight scalar and
vector bosons have also been developed that use the
gravitational-wave interferometers as particle detectors,
which provide additional constraints on the boson mass.
Ultralight bosons can therefore have different effects on
different types of gravitational-wave signals.
This paper presents results from the first all-sky search
tailored to gravitational waves from depleting scalar boson
clouds around black holes using the Advanced LIGO
[37]
data in the third observing run (O3). Although the expected
gravitational-wave signal is almost monochromatic, it
could come from anywhere in the sky; thus, the Doppler
modulations from the relative motion of the Earth and the
source, of
O
ð
10
−
4
f
Þ
, where
f
is the frequency of the signal,
are
a priori
unknown. A fully coherent search, in which we
take a single Fourier transform of the whole data set and
look for peaks in the power spectrum after demodulating
the data, is therefore impossible to perform in each sky
direction because of limited computational power. Indeed,
the number of sky positions to search over scales with the
square of both the Fourier transform length and the
frequency, and would amount to
O
ð
10
20
Þ
at high frequen-
cies
[38]
. This means that we must employ semicoherent
methods, in which we break the data into smaller chunks in
time that are analyzed coherently, and then combined
incoherently, to keep the computational cost under control
while retaining the desired sensitivity
[39
–
41]
. Moreover,
we adopt a multiresolution approach, equivalent to con-
sidering several different Fourier transform lengths
[42]
,
in order to be sensitive towards possible unpredicted
frequency fluctuations of the gravitational-wave signal.
This paper is organized as follows: in Sec.
II
, we describe
our model for a gravitational-wave signal resulting from a
depleting scalar boson cloud around a rotating black hole;
in Sec.
III
, we explain our method to search for such
signals; in Sec.
IV
, we provide information on the datasets
we use in our analysis; in Sec.
V
, we present the results of
the search and our upper limits, including their astrophysi-
cal interpretations in terms of constraints in the boson/black
hole mass plane and on the maximum detectable distance;
finally, in Sec.
VI
, we discuss prospects for future work.
Appendices contain details on the follow-up of selected
outliers.
II. THE SIGNAL
An isolated black hole is born with a defining mass, spin
and charge
[43]
. Ultralight boson particles around a black
hole will scatter off it and extract energy and momentum
from it if the so-called
superradiant
condition is satisfied,
i.e.
ω
<m
Ω
[7,44]
where
ω
is the angular frequency of
the boson (related to the boson mass linearly),
m
is the
azimuthal quantum number in the hydrogenlike gravita-
tional atom, and
Ω
is the angular frequency of the rotating
black hole
’
s outer horizon. Because the bosons have a finite
mass, they become gravitationally bound to the black hole,
which allows for successive scatterings and a build-up of a
macroscopic cloud, extracting up to about 10% of the black
hole mass
[7]
. This process is maximized when the
particle
’
s reduced Compton wavelength is comparable to
the black hole radius:
ℏ
c
m
b
≃
GM
BH
c
2
;
ð
1
Þ
where
m
b
is the boson mass-energy,
M
BH
is the black
hole mass,
ℏ
is the reduced Planck
’
s constant, and
c
is the
speed of light. The typical time scale of the superradiance
phase is
[9]
τ
inst
≈
27
M
BH
10
M
⊙
α
0
.
1
−
9
1
χ
i
days
;
ð
2
Þ
where
χ
i
is the black hole initial dimensionless spin, and
α
¼
GM
BH
c
3
m
b
ℏ
ð
3
Þ
is the fine-structure constant in the gravitational atom.
Once the superradiant condition is no longer satisfied,
the growth stops, and the cloud begins to give off energy
primarily via annihilation of particles in the form of
gravitational waves
[15]
, over a typically much longer
timescale which, for
m
¼
1
and in the limit
α
≪
0
.
1
,is
given by
τ
gw
≈
6
.
5
×
10
4
M
BH
10
M
⊙
α
0
.
1
−
15
1
χ
i
years
:
ð
4
Þ
For scalar bosons, these clouds have a much shorter growth
timescale compared to the timescale for them to deplete
[9]
;
thus, if they exist now, they are more likely to be depleting
rather than growing. The gravitational-wave emission is at a
frequency
f
gw
[14]
:
R. ABBOTT
et al.
PHYS. REV. D
105,
102001 (2022)
102001-2
f
gw
≃
483
Hz
m
b
10
−
12
eV
×
1
−
7
×
10
−
4
M
BH
10
M
⊙
m
b
10
−
12
eV
2
:
ð
5
Þ
In fact, this frequency slowly increases in time due to the
cloud self-gravity and particle self-interaction, quantified
by the boson self-interaction constant
F
b
[13]
.Asin
[13]
,
we consider a scalar boson that, in addition to a nonzero
mass, has a quartic interaction with an extremely small
coupling
λ
, and characterize this coupling in terms of the
parameter
F
b
≔
m
b
=
ffiffiffi
λ
p
. In the context of an axionlike
particle,
F
b
would roughly correspond to the symmetry
breaking energy scale.
Specifically, if
F
b
≳
2
×
10
18
GeV, the spin-up is domi-
nated by the cloud self-gravity and is given by
_
f
gw
≈
7
×
10
−
15
m
b
10
−
12
eV
2
α
0
.
1
17
Hz
=
s
:
ð
6
Þ
With this condition on
F
b
, our search probes energies at the
Planck scale, and is sensitive to the QCD axion with a mass
of
∼
10
−
13
–
10
−
12
eV
[13]
. However, for smaller
F
b
a spin-
up term, due to the cloud self-gravity and change of the
black hole mass, of the form
_
f
≈
10
−
10
m
b
10
−
12
eV
2
α
0
.
1
17
10
17
GeV
F
b
4
Hz
=
s
ð
7
Þ
is dominant in the intermediate self-interaction regime,
for
8
.
5
×
10
16
ð
α
=
0
.
1
Þ
GeV
≲
F
b
≲
2
×
10
18
GeV (a term
with a similar scaling appears due to energy level
transitions). In this regime, the signal duration signifi-
cantly shortens, thus reducing the chance of detection.
In the strong self-interaction regime, when
F
b
≲
8
.
5
×
10
16
ð
α
=
0
.
1
Þ
GeV, the spin-up can be parametrized as
_
f
≈
10
−
10
m
b
10
−
12
eV
2
α
0
.
1
17
10
17
GeV
F
b
6
Hz
=
s
;
ð
8
Þ
and the signal strength rapidly decreases with increasing
interaction, making the detection of annihilation signals
basically impossible for current detectors. In Sec.
VC
,we
will briefly discuss the implication of our search results in
relation to the value of
F
b
. Moreover, at the analysis level
we do not want to exclude the possibility that the emitted
signal frequency evolves in a more complicated manner
than we model, e.g. there are random fluctuations. In
Sec.
III D
, we describe a computationally cheap method
to deal with such frequency variations.
The signal we are considering have lower amplitudes
than those from compact binary coalescences
[17]
, with
initial values (for
α
≪
0
.
1
)
h
0
≈
6
×
10
−
24
M
BH
10
M
⊙
α
0
.
1
7
1
kpc
D
ð
χ
i
−
χ
c
Þ
;
ð
9
Þ
where
χ
c
≈
4
α
4
α
2
þ
1
;
ð
10
Þ
is the black hole spin when the superradiance is saturated,
and
D
is the source distance. The emitted signal amplitude
decreases in time as the clouds are depleted, according to
h
ð
t
Þ¼
h
0
1
þ
t
τ
gw
:
ð
11
Þ
The equations above are analytic approximations to
the true behavior of a gravitational wave originating from
boson clouds. When we consider
α
∼
0
.
1
, we must account
for the difference between the numerical and analytic
solutions, which is
∼
3
in energy at the highest
α
considered
here,
α
∼
0
.
15
We obtain this factor of
∼
3
comparing the
black and red curves in Fig. 2 of
[9]
at
α
∼
0
.
15
.
III. SEARCH METHOD
The semi-coherent method we employ is based on the
band sampled data (BSD) framework
[45]
, which is being
used in an increasing number of continuous-wave
searches, due to its flexibility and computational effi-
ciency
[46
–
48]
. BSD files are data structures that store a
reduced analytic strain-calibrated time series in 10-Hz/
one-month chunks. To construct BSD files, we take a
Fourier transform of a chunk of time-series strain data,
extract a 10-Hz band, keep only the positive frequency
components, and inverse Fourier transform to obtain the
reduced analytic signal.
The main steps of the analysis pipeline are schematically
shown in Fig.
1
and are summarized below.
FIG. 1. Scheme of the search pipeline.
ALL-SKY SEARCH FOR GRAVITATIONAL WAVE EMISSION
...
PHYS. REV. D
105,
102001 (2022)
102001-3
A. Peakmap construction
Our analysis starts with a set of BSD files covering the
whole third observing run and frequencies between 20 and
610 Hz. The maximum frequency is chosen to ensure the
computational cost of the analysis is reasonable, and is
consistent with the fact that the most relevant part of the
accessible parameter space of the black hole-boson system
corresponds to the intermediate frequencies. Each of these
time series is divided into chunks of duration
T
FFT
¼
1
Ω
rot
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
c
2
f
10
R
rot
r
;
ð
12
Þ
where
f
10
is the maximum frequency of each correspond-
ing 10-Hz band,
Ω
rot
¼
2
π
=
86164
.
09
rad
=
s is the Earth
sidereal angular frequency, and
R
rot
is the Earth rotational
radius, which we conservatively take as that corresponding
to the detector at the lowest latitude (
R
rot
¼
5
.
5
×
10
6
km
for the Livingston detector). This
T
FFT
length guarantees
that, if we take the Fourier transform of the chunk, the
power spread of any possible continuous-wave signal,
caused by the motion of the Earth with respect to the
source, is fully contained in one frequency bin with a width
of
δ
f
¼
1
=T
FFT
in each chunk. In Fig.
2
, the chunk
duration is shown as a function of the frequency. For each
of these chunks, with at least a 50% of nonzero values, the
square modulus of the Fourier transform is computed (by
means of the FFT algorithm) and divided by an estimation
of the average spectrum over the same time window. The
resulting quantity has a mean value of one independently of
the frequency
—
for this reason it is called the
equalized
spectrum
—
and takes values significantly different from
one when narrow frequency features are present in the data.
We select prominent peaks in the equalized spectra, defined
as time-frequency pairs that correspond to local maxima
and have a magnitude above a threshold of
θ
¼
2
.
5
[38]
.
The collection of the these time-frequency peaks forms the
so-called
peakmap
, see e.g.
[38]
for more details.
B. Doppler effect correction
We build a suitable grid in the sky such that, after the
Doppler effect correction for a given sky direction, any
residual Doppler effect always produces an error in
frequency within half a frequency bin
[38]
. For each sky
direction in the grid, the peaks in the peakmap are shifted in
order to compensate for the time-dependent Doppler
modulation, according to the relation
f
0
;
t
k
¼
f
obs
;
t
k
1
þ
⃗
v
t
k
·
ˆ
n
c
;
ð
13
Þ
where
f
obs
;
t
k
is an observed frequency peak at the time
t
k
,
⃗
v
t
k
is the detector
’
s velocity, and
ˆ
n
is the unit vector
identifying the sky direction.
As we have described in Sec.
II
, the signal has
an intrinsic spin-up, associated with the cloud deple-
tion. Although we do not apply an explicit correction
for the spin-up, the analysis retains most of the
signal power
—
with a maximum sensitivity loss less than
a few percent
1
—
as long as the corresponding frequ-
ency variation during the full observation time
T
obs
is confined within
1
=
2
bin of the frequency evolution
rate
δ
_
f
, given by
δ
_
f
¼
1
2
T
obs
T
FFT
:
ð
14
Þ
The
δ
_
f
value is plotted as a function of the frequency
in Fig.
3
.
C. Peakmap projection
Once we have applied the Doppler correction, we project
the peakmap onto the frequency axis and select the most
significant frequencies for further analysis. Indeed, if a
monochromatic signal comes from a given direction, we
expect its associated peaks to be aligned in frequency and
appear as a prominent peak in the projected peakmap, when
the signal is strong enough. The statistics of the peakmap
are described in detail in
[38]
, so here we only provide a
brief review. In the absence of signals, the probability
p
0
of
selecting a peak depends on the noise distribution. In the
case of Gaussian noise, it is given by
p
0
¼
e
−
θ
−
e
−
2
θ
þ
1
3
e
−
3
θ
;
ð
15
Þ
FIG. 2. Chunk duration
T
FFT
as a function of frequency. The
T
FFT
is fixed within each 10-Hz band.
1
In principle, the maximum sensitivity loss would be of about
36%, when the corrected signal frequency falls in the middle
among two consecutive frequency bins. In practice, however,
the application of the moving average on the peakmaps, see
Sec.
III D
, allows us to recover part of the lost signal power,
reducing the maximum loss to about 19% (when applying a
moving average of two bins, i.e. a window with
W
¼
2
).
R. ABBOTT
et al.
PHYS. REV. D
105,
102001 (2022)
102001-4
where
θ
¼
2
.
5
is the threshold we apply to construct the
peakmap. On the other hand, if a signal with spectral
amplitude (in units of equalized noise)
λ
¼
4
j
̃
h
ð
f
Þj
2
T
FFT
S
n
ð
f
Þ
;
ð
16
Þ
where
̃
h
ð
f
Þ
is the signal Fourier transform, and
S
n
ð
f
Þ
is the
detector noise power spectral density, is present, the
corresponding probability
p
λ
of selecting a signal peak,
for weak signals with respect to the noise level, is given by
p
λ
¼
p
0
þ
λ
2
θ
ð
e
−
θ
−
2
e
−
2
θ
þ
e
−
3
θ
Þ
:
ð
17
Þ
The statistic of the peakmap projection is, for well-behaved
noise, a binomial with expectation value
μ
¼
Np
0
,
where
N
is the number of FFTs, and standard deviation
σ
¼
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
Np
0
ð
1
−
p
0
Þ
p
. In the presence of a signal, the
expected number of peaks at the signal frequency (after
the Doppler correction) is
Np
λ
. It is then useful to introduce
the critical ratio (CR):
CR
¼
n
−
μ
σ
;
ð
18
Þ
where
n
is the actual number of peaks in a given frequency
bin, which, in the limit of large
N
(in practice, greater than a
∼
30
), closely follows a Gaussian distribution with zero
mean and unity standard deviation.
D. Peakmap moving averages
The signal emitted by a boson cloud is not exactly
monochromatic, due to the presence of a spin-up, see
Sec.
II
. Moreover, further unpredicted frequency fluctua-
tions could be present, due to our limited knowledge of the
physical processes actually taking place. A signal with
varying frequency can be characterized by its
coherence
time
, defined as the minimum amount of time over which
the signal frequency variation exceeds the search frequency
bin. For the sake of robustness, in order to deal with spin-
up values which cause a frequency variation larger than
half a frequency bin during the observing time
T
obs
(see
Sec.
III C
) and
—
especially
—
to take into account possi-
ble uncertainties in the signal morphology, we apply a
sequence of moving averages over a window width
W
varying from two to ten frequency bins, starting from the
Doppler corrected peakmap projection. This procedure can
be shown to be equivalent
—
and computationally much
more efficient
—
to build peakmap projections with an
equivalent FFT duration
≈
T
FFT
=W
. While these moving
averages reduce the sensitivity to signals with a character-
istic coherence time larger than
T
FFT
, e.g. purely mono-
chromatic signals, they provide better sensitivity to signals
with a typical coherence time comparable to the equivalent
FFT duration
[42]
. At the same time, they allow us to
partially recover sensitivity to spin-ups larger than those
shown in Fig.
3
.
E. Candidate selection and coincidences
We use the CR as a detection statistic to identify
potentially interesting outliers. Outliers are uniformly
selected in the parameter space by taking the two points
with the highest CR in each peakmap projection, for each
sky position in every 0.05-Hz frequency band, independent
of the actual CR values, as long as it is above a minimum
value of 3.8.
In the next step, we identify coincident pairs among
outliers found in the Hanford and Livingston data, corre-
sponding to the same
W
, by using a dimensionless
coincidence window of 3. In other words, two outliers,
one from Hanford and the other from Livingston, are
considered coincident if the dimensionless distance
d
¼
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
Δ
f
δ
f
2
þ
Δ
λ
δλ
2
þ
Δ
β
δβ
2
s
≤
3
;
ð
19
Þ
where
Δ
f;
Δ
λ
;
Δ
β
are the differences between parameter
values of the outliers found in Hanford (H1) and Livingston
(L1), and
δ
f;
δλ
;
δβ
are the corresponding bin sizes.
2
The coincidence window of 3 is chosen according to the
studies carried out using simulated signals and is consistent
with the choice in standard all-sky searches for neutron
stars
[49,50]
.
FIG. 3. Bin size of the frequency time derivative,
δ
_
f
, a function
of the frequency, given the full observing time of O3 and
T
FFT
chosen for each 10-Hz band. The search by default covers a
frequency time derivative range of
_
f
∈
½
−
δ
_
f=
2
;
δ
_
f=
2
. Such a
range contributes to set constraints on the portion of the boson
mass/black hole mass explored in the search.
2
The sky position bin sizes in general are different for two
outliers identified at different sky positions, so the average among
the corresponding bin sizes is actually used.
ALL-SKY SEARCH FOR GRAVITATIONAL WAVE EMISSION
...
PHYS. REV. D
105,
102001 (2022)
102001-5
F. Postprocessing
Coincident outliers are subject to a sequence of post-
processing steps in order to discard those incompatible with
an astrophysical signal.
The first veto is to check if some of the outliers are
compatible with known narrow-band instrumental disturb-
ances, also known as
noise lines
. An outlier is considered as
caused by a noise line, and thus discarded, if its Doppler
band, defined by
Δ
f
dopp
¼
10
−
4
f
, where
f
is the outlier
frequency, intersects the noise line.
The next is a consistency veto, in which a pair of
coincident outliers are discarded if their CR is not com-
patible with the detectors
’
noise level at the corresponding
frequency. Specifically, we veto those coincident outliers
for which the CR in one detector is more than a factor
of 5 different than that in the other detector, after being
weighted by the detector
’
s power spectral density
[14,50]
.
This choice, also used in standard all-sky searches for
continuous waves from neutron stars, is motivated by the
desire to safely eliminate an outlier only if the discrepancy
is really significant and the CR is large in at least one of the
two detectors (i.e. a real astrophysical signal is not expected
to behave that way).
The previous steps are all applied to coincident outliers
separately for each value of the moving average window
width
W
. In case of outliers coincident across different
W
values, only the one with the highest CR is kept.
3
G. Follow-up
Outliers with
W
¼
1
, the case in which we do not apply
moving averages, could be due to monochromatic signals
or signals with very small fluctuations in frequency, while
outliers with
W>
1
are more compatible with signals
characterized by a larger degree of frequency fluctuation.
These two sets of candidates are then followed up with
different methods, one based on the Frequency-Hough
algorithm and used e.g. in
[49]
, and the other using the
Viterbi method, see e.g.
[16]
. Both methods are briefly
explained in Appendix
B
.
IV. DATA
We use data collected by the Advanced LIGO gravita-
tional-wave detectors over the O3 run, which lasted from
April 1, 2019 to March 27, 2020, with a one-month break in
data collection in October 2019. The duty cycles of the two
detectors, H1 and L1, are
∼
76%
and
∼
77%
, respectively,
during O3.
In the case of a detection, the calibration uncertainties in
the detector strain data stream could impact the estimates of
the boson properties. Without a detection, these uncertain-
ties could affect the estimated instrument sensitivity and
inferred upper limits on astrophysical source proper-
ties. The analysis uses the
“
C01
”
version of calibrated
data,
4
which has estimated maximum uncertainties (68%
confidence interval) of
∼
7%
–
11%
in magnitude and
∼
4
deg
–
9
deg in phase, over the first/second halves of
O3
[51,52]
. The frequency-dependent uncertainties vary
over the course of a run but do not change by large values.
The time-dependent variations could lead to errors in
opposite directions at a given frequency during different
time periods. In addition, the calibration errors and uncer-
tainties are different at the two LIGO sites. By integrating
the data over the whole run, we expect that the impact
from the time-dependent, frequency-dependent calibration
uncertainties cancels out to some extent, and the overall
impact on the inferred upper limits is within a level of
∼
2%
.
Thus we do not explicitly consider the calibration uncer-
tainties in our analysis.
Due to the presence of a large number of transient
noise artifacts, a
gating
procedure
[53,54]
has also been
applied to the LIGO data. This procedure applies an inverse
Tukey window to the LIGO data at times when the root-
mean-square value of the whitened strain channel in the
band of 25
–
50 or 70
–
110 Hz exceeds a certain threshold.
Only 0.4% and
∼
1%
of data is removed for Hanford
and Livingston, respectively, and the improvement in data
quality from applying such gating is significant, as seen in
the stochastic and continuous gravitational-wave analyses
in O3
[55]
.
V. RESULTS
In this section, we present the results of the analysis
described in Sec.
III
. With no candidates surviving the
follow-up, we present the upper limits on the signal strain
amplitude in Sec.
VB
, and the astrophysical implications
in Sec.
VC
.
Coincident outliers produced by the main search are first
pruned in order to reduce their number to a manageable
level. A further reduction is based on a close case-by-case
inspection of the strongest outliers, as described in the
following.
A. Outlier vetoes
First, we remove outliers due to
hardware injections
,
which simulate the gravitational wave signals expected
from spinning neutron stars, added for testing our analysis
methods by physically moving the mirrors as if a signal had
arrived. See Appendix
A
for more details.
Moreover, we find several bunches of outliers which
are not compatible with astrophysical signals, and thus
are also discarded and not followed up. In these cases, the
3
This restriction reduces the computational cost with only a
small loss in sensitivity.
4
This is equivalent to the publicly available frame data
titled
“
HOFT_CLEAN_SUB60HZ_C01_AR
”
at
https://www
.gw-openscience.org/O3/
.
R. ABBOTT
et al.
PHYS. REV. D
105,
102001 (2022)
102001-6
time-frequency peakmaps and spectra show the presence
of broadband disturbances or lines. Examples of such
disturbances are shown in Fig.
4
, where we plot on the
left the time-frequency peakmap in the frequency range
of 28.2
–
28.35 Hz over the whole run and on the right
the corresponding projection onto the frequency axis. A
wide and strong transient disturbance responsible for many
outliers with frequency in the range of 28.279
–
28.283 Hz is
clearly visible.
In addition, we find an excess of outliers at several
multiples of 10 Hz in both H1 and L1 data, starting from
230 Hz. These are artifacts produced by the procedure
to build BSDs.
5
As it is extremely unlikely that the real
astrophysical signals are coincident with multiples of 10 Hz,
we veto all candidates within a range of
10
−
4
Hz around
each multiple of 10 Hz.
Finally, we apply a final threshold on the CR in order to
further reduce the number of outliers to a small enough
number to permit follow up. In particular, we use different
thresholds CR
thr
¼
4
.
1
and 4.4 for outliers with
W
¼
1
and
W>
1
, respectively. This procedure leaves us with 9449
outliers with
W
¼
1
and 7816 outliers with
W>
1
.
The choice of threshold does not affect the upper limits
presented later, since those are based on the outliers with
the highest CR.
Such outliers are followed-up with two different meth-
ods: a method based on the Frequency-Hough algorithm if
the outliers are characterized by
W
¼
1
(Appendix
B1
),
and a method based on a HMM tracking scheme (the
Viterbi algorithm) if the outliers are characterized by
W>
1
(Appendix
B2
). No potential CW candidate
remains after the follow-up. Details of the follow-up
procedure and results are presented in Appendix
C
.
B. Upper limits
Having concluded that none of the outliers is compatible
with an astrophysical signal, we compute 95% CL upper
limits placed on the signal strain amplitude.
The limits are computed, for every 1-Hz sub-band in the
range of 20
–
610 Hz, with the analytic formula of the
Frequency-Hough sensitivity, given by Eq. (67) in
[38]
,
where the relevant quantities, i.e. an estimate of the noise
power spectral density
S
n
and the CR, are obtained from the
O3 data used in this search. Specifically, for each 1-Hz
subband, we take the maximum CR of the outliers,
separately for Hanford and Livingston detectors, and the
detectors
’
averaged noise spectra in the same subband.
We then compute the two 95% confidence-level upper
limits using the equation, and take the worse (i.e. the
higher) among the two.
This procedure has been validated by injecting simu-
lated signals into LIGO and Virgo O2 and O3 data
[56]
.In
particular, this procedure has been shown to produce
conservative upper limits with respect to those obtained
with a much more computationally demanding injection
campaign. It has been also verified that the upper limits
obtained with a classical procedure based on injections are
always above those based on the same analytical formula
but using the minimum CR in each 1-Hz subband. The
two curves, based, respectively, on the highest and the
smallest CR, define a belt containing both a more stringent
upper limit estimation and the sensitivity estimation of our
search. When we will discuss the astrophysical implica-
tions of the search, we will always refer to the upper limits
obtained using Eq. (67) in
[38]
, which are conservative
with respect to limits derived from injections. This same
procedure is being used for the O3 standard all-sky search
for continuous waves from spinning neutron stars
[57]
.
Equation (67) of
[38]
gives the minimum detectable
amplitude at a given confidence level, and has been
derived by marginalizing over sky pos
ition and polariza-
tion of the source.
FIG. 4. Time-frequency peakmap (left) and corresponding projection (right) in the frequency band of 28.2
–
28.35 Hz over the full run.
Several instrumental artifacts are visible, e.g. strong lines with varying amplitudes (
∼
28
.
22
Hz), weaker lines with sudden frequency
variation (
∼
28
.
33
Hz), and broad transient disturbances (vertical bright lines). In particular, the features around 28.28 Hz in the second
half of O3 are responsible for a big bunch of outliers with frequency in the range of 28.279
–
28.283 Hz.
5
The artifacts are a well-known consequence of the 10 Hz
baseline used to build BSD files. They could be avoided by
constructing the BSDs with interlaced frequency bands.
ALL-SKY SEARCH FOR GRAVITATIONAL WAVE EMISSION
...
PHYS. REV. D
105,
102001 (2022)
102001-7
The resulting upper limits are shown in Fig.
5
,as
circles connected by a dashed line. The minimum value
is
1
.
04
×
10
−
25
at 130.5 Hz. In the same figure, the dots
connected by the dashed line correspond to the lower
bound computed using the minimum CR. A comparison
with previous all-sky searches, see e.g. Fig. 4 of
[58]
or
Fig. 6 of
[59]
, shows that our conservative results are better
than most of the past searches, including the recent early
O3 analysis
[58]
. In particular, our minimum upper limit
value improves upon that in
[58]
by
∼
30%
. The main
factors that contribute to the improvement relative to
[58]
are (1) the use of the full O3 C01 gated data, (2) the use of
longer FFTs at least in a portion of the frequency band, and
(3) the restricted spin-up/down range that could impact the
maximum CR and thus the upper limit in each subband. On
the other hand, the all-sky search described in
[58]
, as well
as most of the other past searches, cover a significantly
larger parameter space in both the spin-down/up and
frequency ranges, and thus are sensitive to signals that
are not considered in this analysis. An O2 all-sky search
reported in
[60]
produces upper limits slightly better than
ours, by about 5
–
10% on average (due to the use of a much
longer coherent integration time of days), but obtained over
a smaller parameter space and more limited in scope, being
focused on low-ellipticity sources. See Sec.
VI
for a more
detailed discussion.
We want to stress two important points here. First, our
conservative procedure to compute upper limits, based on
the maximum CR in each subband, produces strain values
that are typically larger than the minimum detectable
amplitudes at the same frequency. That is, the search
sensitivity is expected to be better than the upper limits.
Second, this is the first all-sky search that is optimized for
frequency wandering signals, both at the level of outlier
selection and follow-up. This allows us to achieve a better
sensitivity to such signals compared to that achievable
using the maximum possible FFT duration (which is the
best choice for nearly monochromatic signals).
C. Astrophysical implications
The upper limits presented above can be translated into
physical constraints on the source properties. We interpret
the results in two different ways. First, similar to what has
been done in
[14]
, we compute the exclusion regions in the
boson mass and black hole mass plane, assuming fixed
values of the other relevant parameters, namely the distance
to the system
D
, the initial dimensionless spin of the black
hole
χ
i
, and a time equal to the age of the cloud
t
¼
t
age
.
Indeed, from Eqs.
(9)
and
(11)
, once
D
,
χ
i
, and
t
age
are
fixed, the signal strain amplitude depends only on the
boson mass
m
b
and black hole mass
M
BH
. Thus for each
pair of
ð
m
b
;M
BH
Þ
, we can exclude the presence of a source
with those assumed parameters, if it would have produced
a signal whose strain is larger than the upper limit at a given
frequency in Fig.
5
. In our search, the
α
values probed
roughly fall in the range of [0.02, 0.13]. Equations
(4)
and
(9)
do not hold in the range of
α
≳
0
.
1
, as shown in
Fig. 2 of
[9]
. Specifically, we look at the Brito
þ
curve
derived analytically compared to the black curve obtained
from numerical simulations, and see that they differ in
energy by a factor of 3 at the largest
α
considered in this
search,
α
0.15. This implies that we underestimate
τ
gw
up to
a factor of 3 and overestimate
h
0
up to a factor of
ffiffiffi
3
p
(at
α
≃
0
.
13
). Thus we correct Eqs.
(4)
and
(9)
by these factors
such that the resulting exclusion regions are slightly
conservative in the whole parameter space. In Fig.
6
, the
exclusion regions are plotted for
D
¼
1
kpc (left) and
D
¼
15
kpc (right), assuming a high initial spin of
χ
i
¼
0
.
9
.For
each distance assumption, three different values of the
black hole age are considered. In Fig.
7
, the exclusion
regions are shown for the same parameters as before, except
that a lower initial spin value,
χ
i
¼
0
.
5
is assumed. In both
cases, as expected, the constrained region is smaller when
the source is assumed to be at a farther distance as, in this
case, the signal amplitude at the detector is smaller. For any
given black hole mass, these results improve upon the
constraints obtained in Advanced LIGO O2 data
[14]
for lower boson masses, while they are slightly less con-
straining for higher boson masses. This is because for a given
black hole mass, higher boson masses correspond to higher
signal spin-up [see Eq.
(6)
],andwesimplydonotcovera
large spin-up range. On the other hand, the constraints
described in
[14]
were obtained from the results of a search
not specifically designed for boson clouds. In particular,
restricting the exploration to the parameter space relevant for
the expected boson cloud signals significantly reduces the
trial factor with a consequent implicit gain in the sensitivity
with respect to standard wider parameter space searches.
A second type of interpretation, as shown in Fig.
8
,is
represented by the maximum distance
D
max
at which we
FIG. 5. 95% C.L. upper limits on the signal strain amplitude
obtained from this search (circles, connected by the dashed line).
Values are given in each 1-Hz subband. No value is given for the
subbands where no outlier survives. The dots connected by a
dashed line define the lower bound obtained using the minimum
CR in each subband.
R. ABBOTT
et al.
PHYS. REV. D
105,
102001 (2022)
102001-8
FIG. 6. Exclusion regions in the boson mass (
m
b
) and black hole mass (
M
BH
) plane for an assumed distance of
D
¼
1
kpc (left) and
D
¼
15
kpc (right), and an initial black hole dimensionless spin
χ
i
¼
0
.
9
. For
D
¼
1
kpc, three possible values of the black hole age,
t
age
¼
10
3
;
10
6
;
10
8
years, are considered; for
D
¼
15
kpc,
t
age
¼
10
3
;
10
4
.
5
;
10
6
years are considered.
FIG. 7. Same as Fig.
6
but for black hole initial spin
χ
i
¼
0
.
5
. The assumed distance is
D
¼
1
kpc (left), and
D
¼
15
kpc (right).
FIG. 8. Maximum distance at which at least 5% of a simulated population of black holes with a boson cloud would produce a
gravitational-wave signal with strain amplitude larger than the upper limit in the detectors. The left plot refers to a maximum black hole
mass of
50
M
⊙
, while the right plot to a maximum mass of
100
M
⊙
. The different colored markers correspond to different system ages,
ranging from
10
3
years to
10
7
years, as indicated in the legend. The alignment of points for different ages at the smallest boson masses
(and distances) is the result of a discretization effect due to the finite size grid used in distance.
ALL-SKY SEARCH FOR GRAVITATIONAL WAVE EMISSION
...
PHYS. REV. D
105,
102001 (2022)
102001-9