of 28
All-sky search for gravitational wave emission from scalar boson clouds around
spinning black holes in LIGO O3 data
The LIGO Scientific Collaboration, The Virgo Collaboration, and The KAGRA Collaboration
(Dated: December 1, 2021)
This paper describes the first all-sky search for long-duration, quasi-monochromatic 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 Hz 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.
I. INTRODUCTION
Theories of beyond standard model physics allow for
the production of ultralight bosons that could consti-
tute a portion or all of dark matter [1]. If such ul-
tralight bosons exist, they could appear around rotat-
ing black holes due to quantum fluctuations [2]. 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 [2].
The energy structure of the clouds resembles that of a
hydrogen atom, earning boson clouds the name “grav-
itational atoms” in the sky. After the black hole spin
decreases below a threshold, the superradiance process
stops and the cloud depletes over time mainly through
annihilation of bosons into gravitons, which generates
quasi-monochromatic, long-duration gravitational 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 [2]. The values of the spin-up depend on whether
we consider scalar [3], vector [4, 5], or tensor [6] 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 emis-
sion from boson self-interactions has been studied for the
scalar case [7], showing that the signal can be signifi-
cantly 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) [8]. Furthermore, a search
Full author list given at the end of the article.
for boson clouds from Cygnus X-1 was performed with
a hidden Markov model (HMM) method [9] on the same
data set, and disfavored particular boson masses for this
object [10]. It has been suggested [11] that current meth-
ods for all-sky searches, such as the one used in [8], 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 depen-
dence at leading order on the boson mass). A detailed
study [12] has quantified this effect, showing that the ac-
tual 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 quasi-monochromatic 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 re-
quires 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 [13, 14], and the
current spin distribution of black holes depends on both
the mass of boson clouds and the distribution of spins
when the black holes were born [15, 16]. It is therefore
interesting to hierarchically combine spin measurements
from various black hole mergers to obtain constraints on
boson cloud/spin interactions [17, 18]. 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 (parameterized as
Love numbers) [19, 20], and such differences in binary
black hole mergers may actually be detectable depend-
ing on the boson mass and field strength [21]. It has
also been proposed to combine multi-band observations
of black hole mergers and boson cloud signals when fu-
ture detectors are online [22], and to search for a stochas-
tic gravitational-wave background from scalar and vec-
tor boson cloud signals [23–25]. Complementary meth-
arXiv:2111.15507v1 [astro-ph.HE] 30 Nov 2021
2
ods [26, 27] and searches [28, 29] 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 [30] data in the third observing run (O3). Al-
though 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 frequencies [31]. This means
that we must employ semi-coherent 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 retain-
ing the desired sensitivity [32, 33]. Moreover, we adopt a
multi-resolution approach, equivalent to considering sev-
eral different Fourier Transform lengths [34], in order to
be sensitive towards possible unpredicted frequency fluc-
tuations of the gravitational-wave signal.
This paper is organized as follows: in section II, we
describe our model for a gravitational-wave signal result-
ing from a depleting scalar boson cloud around a rotating
black hole; in section III, we explain our method to search
for such signals; in section IV, we provide information on
the data sets we use in our analysis; in section V, we
present the results of the search and our upper limits,
including their astrophysical interpretations in terms of
constraints in the boson/black hole mass plane and on
the maximum detectable distance; finally, in section 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 [35]. Ultralight boson particles around a
black hole will scatter off it and extract energy and mo-
mentum from it if the so-called
superradiant
condition is
satisfied, i.e.
ω < m
Ω [2, 36] where
ω
is the angular fre-
quency of the boson (related to the boson mass linearly),
m
is the azimuthal quantum number in the hydrogen-
like gravitational 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 scat-
terings and a build-up of a macroscopic cloud, extracting
up to about 10% of the black hole mass [2]. This pro-
cess 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 [3]
τ
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 en-
ergy primarily via annihilation of particles in the form of
gravitational waves [9], over a typically much longer time
scale 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 time scale compared to the time scale for them to
deplete [3]; 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
[8]:
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 by an
amount that depends on the boson self-interaction con-
stant
F
b
[7]. As in [7], we consider a scalar boson that,
in addition to a non-zero mass, has a quartic interaction
with an extremely small coupling
λ
, and characterize this
coupling in terms of the parameter
F
b
:=
m
b
/
λ
. In the
context of an axion-like particle,
F
b
would roughly cor-
respond to the symmetry breaking energy scale.
Specifically, if
F
b
&
2
×
10
18
GeV, the spin-up is dom-
inated by the depletion of the cloud through particle an-
nihilation 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 ax-
ion with a mass of
10
13
–10
12
eV [7]. However, if
3
the bosons have a non-negligible self-interaction (i.e. for
smaller
F
b
), two other spin-up terms, one due to energy
level transitions and another due to the variation of the
self-interaction energy, appear. The former is given by
̇
f
tr
10
10
(
m
b
10
12
eV
)
2
(
α
0
.
1
)
17
(
10
17
GeV
F
b
)
4
Hz
/
s
,
(7)
which is dominant in the intermediate self-interaction
regime, for 8
.
5
×
10
16
(
α/
0
.
1) GeV
.
F
b
.
2
×
10
18
GeV.
In this regime, the signal duration significantly 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
se
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. V C, we
will briefly discuss the implication of our search results
in relation to the value of
F
b
.
These signals have lower amplitudes than those from
compact binary coalescences [11], 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 satu-
rated, 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 ac-
count for the difference between the numerical and ana-
lytic 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 [3] at
α
0
.
15.
III. SEARCH METHOD
The semi-coherent method we employ is based on the
Band Sampled Data (BSD) framework [37], which is be-
ing used in an increasing number of continuous-wave
searches, due to its flexibility and computational effi-
ciency [38–40]. 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 schemati-
cally shown in Fig. 1 and are summarized below.
FIG. 1. Scheme of the search pipeline.
A. Peakmap construction
Our analysis starts with a set of BSD files covering
the whole third observing run and frequencies between
20 Hz and 610 Hz. The maximum frequency is chosen to
ensure the computational cost of the analysis is reason-
able, and is consistent with the fact that the most rele-
vant part of the accessible parameter space of the black
hole-boson system corresponds to the intermediate fre-
quencies. Each of these time series is divided into chunks
of duration
T
FFT
=
1
rot
c
2
f
10
R
rot
,
(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 rota-
tional radius, which we conservatively take as that cor-
responding 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 Trans-
form 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, the square mod-
ulus of the Fourier Transform is computed (by means of
4
FIG. 2. Chunk duration
T
FFT
as a function of frequency. The
T
FFT
is fixed within each 10-Hz band.
the FFT algorithm) and divided by an estimation of the
average spectrum over the same time window. The re-
sulting 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 spec-
tra, defined as time-frequency pairs that correspond to
local maxima and have a magnitude above a threshold of
θ
= 2
.
5 [31]. The collection of the these time-frequency
peaks forms the so-called
peakmap
, see e.g. [31] 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 [31]. 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 in-
trinsic spin-up, associated with the cloud depletion. Al-
though 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
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 section
- as long as the corresponding frequency 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.
FIG. 3. Bin size of the frequency time derivative,
δ
̇
f
, a func-
tion of the frequency, given the full observing time of O3 and
T
FFT
chosen for each 10-Hz band. The search by default cov-
ers 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.
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. In-
deed, if a monochromatic signal comes from a given di-
rection, we expect its associated peaks to be aligned in
frequency and appear as a prominent peak in the pro-
jected peakmap, when the signal is strong enough. The
statistics of the peakmap are described in detail in [31],
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)
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
|
̃
h
(
f
)
|
2
T
FFT
S
n
(
f
)
,
(16)
III D, allows us to recover part of the lost signal power, reduc-
ing the maximum loss to about 19% (when applying a moving
average of two bins, i.e. a window with W=2).
5
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 de-
viation
σ
=
Np
0
(1
p
0
). In the presence of a signal,
the expected number of peaks at the signal frequency (af-
ter 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 fre-
quency bin, which, in the limit of large N (in practice,
greater than a
30), closely follows a Gaussian distribu-
tion with zero mean and unity standard deviation.
D. Peakmap moving averages
The signal emitted by a boson cloud is not exactly
monochromatic. In order to deal with the uncertainty in
the signal morphology, we apply a sequence of moving av-
erages (MA) 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 effective FFT dura-
tion
T
FFT
/W
. While these moving averages reduce the
sensitivity to signals with a characteristic coherence time
larger than
T
FFT
, e.g. purely monochromatic signals,
they provide better sensitivity to signals with a typical
coherence time comparable to the effective FFT duration
[34].
E. Candidate selection and coincidences
We use the CR as a detection statistic to identify po-
tentially interesting outliers. Outliers are uniformly se-
lected 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, indepen-
dent 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, cor-
responding to the same
W
, by using a dimensionless co-
incidence window of 3. In other words, two outliers, one
from Hanford and the other from Livingston, are consid-
ered coincident if the dimensionless distance
d
=
(
f
δf
)
2
+
(
λ
δλ
)
2
+
(
β
δβ
)
2
3
(19)
where ∆
f,
λ,
β
are the differences between param-
eter 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 accord-
ing to the studies carried out using simulated signals and
is consistent with the choice in standard all-sky searches
for neutron stars [41, 42].
F. Post-processing
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 dis-
turbances, also known as
noise lines
. An outlier is con-
sidered 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 coin-
cident outliers are discarded if their CR is not compatible
with the detectors’ noise level at the corresponding fre-
quency. 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 be-
ing normalized by the detector’s power spectral density
[8, 42]. 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 dis-
crepancy 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 MA 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.
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.
3
This restriction reduces the computational cost with only a small
loss in sensitivity.
6
These two sets of candidates are then followed up with
different methods, one based on the Frequency-Hough
algorithm, and used e.g. in [41], and the other using the
Viterbi method, see e.g. [10]. Both methods are briefly
explained in Appendix B.
IV. DATA
We use data collected by the Advanced LIGO
gravitational-wave detectors over the O3 run, which
lasted from April 1st, 2019 to March 27th, 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 uncertain-
ties in the detector strain data stream could impact the
estimates of the boson properties. Without a detec-
tion, these uncertainties could affect the estimated in-
strument sensitivity and inferred upper limits on astro-
physical source properties. The analysis uses the “C01”
version of calibrated data, which has estimated maxi-
mum uncertainties (68% confidence interval) of
7%/
11% in magnitude and
4 deg/
9 deg in phase, over
the first/second halves of O3 [43, 44]. 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 ad-
dition, the calibration errors and uncertainties are dif-
ferent 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 uncer-
tainties cancels out to some extent, and the overall im-
pact 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 [45, 46] 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 chan-
nel in the band of 25–50 Hz or 70–110 Hz exceeds a
certain threshold. Only 0
.
4% and
1% of data is re-
moved 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 [47].
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. V B, and the astrophysical implications
in Sec. V C.
Coincident outliers produced by the main search are
first pruned in order to reduce their number to a man-
ageable level. A further reduction is based on a close
case-by-case inspection of the strongest outliers, as de-
scribed in the following.
A. Outlier vetos
First, we remove outliers due to
hardware injections
,
which simulate the gravitational wave signals expected
from spinning neutron stars, added for testing our anal-
ysis 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
time-frequency peakmaps and spectra show the presence
of broadband disturbances or lines, Examples of such dis-
turbances 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
4
. 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 or-
der 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 algo-
rithm if the outliers are characterized by
W
= 1 (Ap-
pendix B 1), and a method based on a HMM tracking
scheme (the Viterbi algorithm) if the outliers are charac-
terized by
W >
1 (Appendix B 2). No potential CW can-
didate remains after the follow-up. Details of the follow-
up procedure and results are presented in Appendix C.
4
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.
7
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.
B. Upper limits
Having concluded that none of the outliers is compat-
ible with an astrophysical signal, we compute 95% confi-
dence level (C.L.) 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 Eqn. (67) in [31],
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 [48].
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 for-
mula 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 es-
timation of our search. When we will discuss the astro-
physical implications of the search, we will always refer to
the upper limits obtained using Eqn. (67) in [31], which
are conservative with respect to limits derived from in-
jections. This same procedure is being used for the O3
standard all-sky search for continuous waves from spin-
ning neutron stars [49].
The resulting upper limits, marginalized over the sky
and the source polarization parameters, are shown in
Fig. 5, as circles connected by a dashed line. The mini-
mum 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 com-
parison with previous all-sky searches, see e.g. Fig. 4 of
[50], shows that our conservative results are better than
most of the past searches, including the recent early O3
analysis [50]. In particular, our minimum upper limit
value improves upon that in [50] by
30%. The main
factors that contribute to the improvement 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 [50],
as well as most of the other past searches, cover a signif-
icantly 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 [51] produces upper limits slightly bet-
ter 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 de-
tectable amplitudes at the same frequency. That is, the
search sensitivity is expected to be better than the up-
per 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 dura-
tion (which is the best choice for nearly monochromatic
signals).
8
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.
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 [8], 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 dimen-
sionless 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 figure 5. In
our search, the
α
values probed roughly fall in the range
of [0.02, 0.13]. Eqs. (4) and (9) do not hold in the range
of
α
&
0
.
1, as shown in Fig. 2 of [3]. 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 over-
estimate
h
0
up to a factor of
3 (at
α
'
0
.
13). Thus we
correct Eqs. (4) and (9) by these factors such that the re-
sulting 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 [8]
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)), and we simply do not
cover a large spin-up range. On the other hand, the con-
straints described in [8] were obtained from the results
of a search not specifically designed for boson clouds.
In particular, restricting the exploration to the param-
eter 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
can exclude, as a function of the boson mass, the pres-
ence of an emitting system of a given age
t
age
, assuming a