Characterizing the continuous gravitational-wave signal from boson clouds around Galactic isolated black holes

Ultralight bosons can form large clouds around stellar-mass black holes via the superradiance instability and produce continuous gravitational wave signals with frequencies within the range of LIGO and Virgo. Unlike continuous gravitational waves from neutron stars, boson annihilation signals are clustered in frequency and have very small positive intrinsic frequency derivatives. We characterize the expected spin-0 boson annihilation ensemble signal from synthetic populations of isolated Galactic black holes. We explore how the different population parameters affect the resulting signal and consider its detectability by recent searches for continuous gravitational waves. A population of $10^8$ black holes with masses up to $ 30~\mathrm{M}_\odot$ and a flat dimensionless initial spin distribution between zero and unity produces up to a thousand signals loud enough to be in principle detected by these searches. For a more moderately spinning population the number of signals drops by about an order of magnitude, still yielding up to a hundred detectable signals for some boson masses. A non-detection of annihilation signals at frequencies between 100 and 1200 Hz disfavors the existence of scalar bosons with rest energies between $2\times10^{-13}$ and $2.5\times10^{-12}$ eV. However, due to the high signal density between 200 and 300 Hz, we urge caution when interpreting a null result for bosons between $4$ and $6\times 10^{-13}$ eV. Finally, we emphasize that the ensemble signal morphology would be the primary, perhaps sole, indicator that a continuous wave signal has a boson annihilation origin.


I. INTRODUCTION
In the last few years, transient and rapidly evolving gravitational waves have been observed from the mergers of stellar-mass compact objects [1]. Persistent and slowly evolving sources of gravitational waves have been predicted as well and have yet to be detected. These continuous gravitational waves (CWs) are expected to be much weaker than the transient events and have durations much longer than the typical observation time; searches for CWs integrate over long periods of time to extract signal from background.
Canonical sources of CWs include neutron stars with mass deformations or internal fluid oscillation modes (see [2] for a recent review). Scenarios involving more exotic emitters of CWs are also being explored and can provide evidence for -or disfavor the presence of -new physics beyond the Standard Model of particle physics [3][4][5][6][7]. A particularly well-motivated target is the axion, proposed to solve the strong-CP problem in particle physics [8][9][10]; axions and axion-like particles are also promising dark matter candidates (see e.g. [11] for a review).
Bosons such as axions or axion-like particles can form "clouds" with enormous occupation numbers around rotating black holes, and do so rapidly on astrophysical timescales when the black hole size is similar to the boson's Compton wavelength [3,4]. These boson-black hole systems can be thought of as gravitational "atoms," and within this scenario, boson annihilations and level transitions source monochromatic gravitational-wave emission [3,4]. Annihilations in particular produce signals that may be loud enough to be detected in standard CW searches on data from the current generation of groundbased interferometers [12][13][14][15]. Bosons with masses of ∼ 10 −13 to 4 × 10 −12 eV produce annihilation signals with frequencies between ∼ 50 and 2000 Hz, in the LIGO and Virgo band.
The boson cloud forms by extracting energy and an-arXiv:2003.03359v1 [gr-qc] 6 Mar 2020 gular momentum from the black hole, which loses a significant fraction of its natal spin [3,4]. Measurements of the spins of several old, highly rotating black holes in Xray binaries [16,17] have been used to disfavor the range of masses 6 × 10 −13 eV µ b 2 × 10 −11 eV [12,18]. However, since the spin measurements come with a set of systematic uncertainties, it is very important to undertake complementary searches. Direct CW searches for the annihilation signal in LIGO and Virgo data are an independent test in the boson mass range disfavored by rapidly rotating black holes, and may be able to extend the reach to lower masses, as well as discover a new particle.
In order to disfavor a range of boson masses or characterize a potential CW signal, it is crucial to consider the ensemble signal from the population of sources as a whole. The properties of the annihilation signal depend strongly on the mass, spin, distance, and age of each black hole, and so the properties of the ensemble signal depend on the properties of the black hole population. Moreover, to leading order, the annihilation signal frequency is set by the boson rest mass, and so the emission from all boson clouds falls within a small frequency range. This is in contrast to CWs from neutron stars, which are expected to span a broad range of frequencies depending on the rotation rates of the individual stars. The clustering of signals in a narrow frequency band could reduce the effectiveness of CW search methods -which are optimized for weak, isolated signals -in identifying the annihilations, and may recommend the use of other methods entirely (e.g., [15,19]).
In this paper, we study the expected boson annihilation signal from the population of isolated black holes in the Galaxy, of which there are expected to be up to ∼ 10 8 . In order to take all effects into account, we use simulated populations of 10 8 black holes and calculate the expected signal from all of the systems for bosons with energies between 1 × 10 −13 and 4 × 10 −12 eV, corresponding approximately to gravitational-wave frequencies between 50 and 2000 Hz. We investigate the detectability of the ensemble signal in current LIGO data by broad continuous gravitational wave surveys [20][21][22][23], and its dependence on black hole population assumptions.
We review the theory of boson cloud formation and CW emission in Section II and present the Galactic isolated black hole population in Section III. We explore the resultant ensemble signal from the entire population in Section IV, study its detectability in Section V, compare our results to current literature in Section VI, and summarize in Section VII.

II. SIGNAL MODEL
"Gravitational atoms"-macroscopic, gravitationallybound states of ultralight bosons around astrophysical black holes-form rapidly via the superradiance instability. The bosons subsequently annihilate, sourcing coher-ent, monochromatic, and long-lasting gravitational waves [3,4,12]. The frequency of the signal at leading order is given by twice the boson's rest energy, with µ b = m b c 2 , m b the boson mass, c the speed of light, and h the Planck constant. In this work, we focus on signals from gravitationally interacting scalar (spin-0) bosons around stellar-mass black holes; we use "black hole" to exclusively refer to stellar-mass black holes throughout the text. The growth and annihilation timescales of vector (spin-1) bosons are shorter [24][25][26][27][28] and require a separate analysis. We only consider annihilation signals from the fastest-growing bound state and limit our analysis to isolated black holes without external effects such as accretion or binary companions. We summarize the necessary background below and provide further details in App. A; see, e.g., [4,12] for more details and [29] for a review.

A. Cloud formation
Black hole superradiance is a purely kinematic process [30] whereby a wave scattering off a rapidly rotating black hole increases in amplitude by extracting some of the black hole's angular momentum [31,32]. The gravitational potential of the black hole enables massive particle bound states, and the amplitudes of bound states that satisfy the "superradiance condition" increase [33][34][35][36]. For bosons, the growth is exponential and results in a "cloud" with a macroscopic number of particles all occupying the same state. The initial seed can be a vacuum fluctuation, so the process need not rely on an existing abundance of bosons in the black hole's vicinity.
will grow if the initial spin of the black hole, χ i , satisfies the superradiance condition, α 1 2 The inequality is equivalent to the angular velocity of the black hole exceeding the angular velocity of the cloud. For a given α, this is true when the initial spin of the black hole χ i is greater than the critical spin χ c , See Eq. (A8) for gravitational potential energy corrections and dependence on bound state quantum numbers. For α ∈ [0.03, 0.2], χ c ∈ [0.1, 0.7] (Fig. 2). The (n, , m) = (0, 1, 1) level has the largest annihilation power and the fastest growth time, with e-folding time for the number of particles τ inst given by at leading order in α, see also Eq. (A7). The instability timescale is very sensitive to the boson and black hole masses and systems with lighter bosons and lighter black holes take longer to form. The cloud extracts angular momentum from the black hole until χ χ c , at which point the cloud ceases to grow. The black hole loses mass as well as spin in this process, so the values of α and χ c decrease slightly as the cloud grows; we take this evolution into account as described in App. A. At this point, the masses of the black hole and of the cloud are approximately where M BH is the black hole mass at the onset of the cloud formation.
For the systems of interest, M cloud ≈ 0.1-5% M BH , which corresponds to particles in the cloud. The timescale to fully populate the level is then ln(N )τ inst ∼ 180 τ inst .

B. Gravitational wave emission
We study gravitational-wave signals emitted from the resultant bound state: two bosons within the cloud can annihilate into a single graviton in the presence of the black hole gravitational field. "Transition" signals can also occur if multiple states are simultaneously populated, but are less promising at current sensitivities [12]. We assume the bosons only interact via gravity 2 . The presence of self-interactions can change the evolution and limit the size of the cloud [37,38], or potentially cause the collapse of the cloud in a "bosenova" [4,39,40]. The signal from the bosons in a single macroscopic bound state is coherent, monochromatic, and evolves slowly, thus providing an ideal target for CW searches. For lighter bosons and black holes, the signal properties are essentially unchanged over a Hubble time.
The frequency f GW of the emitted gravitational wave is given by The leading-order gravitational-wave frequency f 0 GW is proportional to the boson rest energy µ b (Eq. (9)). Top: Systems with heavier black holes (orange curves) produce GW signals with lower frequencies than do systems with lighter black holes (blue curves) due to the larger gravitational binding energy. Systems with larger spins (dashed curves) result in slightly higher frequencies than systems with smaller spins (solid curves) due to a positive spin-orbit energy. The two solid curves -corresponding to χi = 0.5 -stop at intermediate boson masses; for heavier bosons, χc > 0.5, and the systems never form. The contribution from the cloud self-energy is of order 10 −3 or less of the black hole's gravitational potential energy.
where the corrections to the leading order frequency, Eq. (1), are due to the gravitational potential of the black hole [41,42], see Fig. 3, and the cloud itself, The power emitted in gravitational waves for α 1 is given by [43], The corrections at larger α are significant, and the power in this regime has been computed numerically [14,44] (see Eq. (A16)). The gravitational waves are produced approximately in a background defined by the final black hole mass, and we use the final value of α to evaluate the power, strain, and timescale expressions. The characteristic strain of the annihilation signal is largest when the cloud first reaches its maximum occu-  14)) ranges from longer than the Hubble time for lighter bosons to less than 100 years for heavier bosons. The longer signals from lighter bosons are weaker compared to shorter signals from heavier bosons (Eq. (13)). pation number. At leading order, where h 0 is the maximal intrinsic gravitational wave amplitude (Eqs. (A17),(A18)). The strain h 0,peak thus rapidly increases for larger values of µ b , M BH (Fig. 4 [20,21,45,46]. As the bosons annihilate and the cloud is depleted, the signal strength h 0 (t) decreases from its peak as [4] h 0 (t) = h 0,peak where τ GW = M peak cloud c 2 /P peak GW is the signal "half-time", The self-binding energy of the cloud decreases with decreasing cloud mass (Eq. (10)), giving a small positive first frequency derivative (also see Eq. (A20)), In the relevant part of parameter space, the frequency driftḟ gw 10 −12 Hz/s ( Figure 22) is smaller than what typical all-sky CW searches can resolve [45,46] and is overall a minor effect in our regime. The frequency drift is discussed further in Sec. IV A.
C. Next-fastest-growing bound state The second-fastest-growing bound state with (n = 0, = m = 2) has parametrically longer growth and annihilation timescales (App. A). It also corresponds to a lower critical spin and thus will reduce the black hole spin below the first level's critical value χ c . The black hole then begins absorbing the first level bosons [4,47], and the "spin-up" rate of the first level absorption balances the spin-down rate caused by the growth of the second level. Ultimately, enough of the first level is absorbed by the black hole that the annihilations effectively cease, and CW emission shuts off [4]. We take this effect into account by setting h 0 (t) = 0 when the second level fully populates, ∼ ln(N ) τ 022 inst after the black hole formation, where The cutoff time in the signal becomes important at α 0.1, reducing the number of signals. The second level also produces CWs, and given the lower critical spin, clouds can form in systems with smaller χ i or larger boson masses. However, the strain is significantly weaker than for the first level; as described in App. A, h 011 0,peak /h 022 0,peak ∼ 90/α 2 [44]. In addition, the emission from the (n = 0, = m = 2) level is no longer dominantly quadrupolar [44], which would require an adhoc CW search. For these reasons, we do not consider the gravitational-wave emission from the second-fastestgrowing state in our study, although it would be interesting to do so in the future.

III. THE BLACK HOLE POPULATION
For a given boson mass, the annihilation signal properties depend on properties of the black hole: the distance d, velocity v, initial M BH , and initial spin χ i (Eqs. (9), (13)), and age τ BH . The total number of black holes in the Galaxy is estimated to be around 10 8 (e.g., [48,49]); based on black hole formation rates of ∼ 0.1 − 0.9 per Milky-Way type galaxy per century, the total number typically varies between 10 7 and 10 8 [50][51][52][53]. The number of isolated black holes is unknown, but is expected to be the same order of magnitude as the total number of black holes (e.g., [54,55]). We therefore take 10 8 as the benchmark number of isolated black holes in the Galaxy; as the total number is uncertain, our results can be rescaled to draw conclusions on different total populations. To manage computation time, we simulate the positions and velocities of 10 6 black holes and perform a "resampling" to scale up to 10 8 isolated black holes in the Galaxy (Sec. III A).

A. Simulating positions and velocities
The spatial and velocity distribution of the current population of isolated Galactic black holes is not known. However, since the black holes formed from massive stars at an earlier time, we approximate the black holes' spatial distribution at birth with the current stellar distribution and evolve their trajectories, including initial bulk velocities and natal kicks, though the Galactic gravitational potential using the procedure described in Tsuna et al. [56]. The Galactic gravitational potential is divided into the bulge, disc, and halo components; the black holes are born in the bulge and disc according to their respective birth rates in the regions, and their positions and velocities are evolved through the Galactic gravitational potential. (See [56] and references therein for more details.) Direct observation of isolated Galactic black holes, such as by radio observations [56,57], would significantly reduce the uncertainties in their properties.
The disk is defined by a cylindrical coordinate system with radial coordinate r, height z, and azimuth φ, while the halo is defined by a spherical coordinate system with radial coordinate R. The bulge is described by both cylindrical and spherical coordinate systems. All coordinates are defined with the Galactic Center at the origin. The black hole birth rate per unit volume in the bulge (subscript "b") and disk ("d") are given the following prescriptions [58,59]: where R b = 120 pc and r d = 2.15 kpc;ρ d is uniform along z for |z| < 75 pc. The proportionality factors in Eqs. (18) and (19) are determined by requiring that the total number of black holes born be 10 6 , with 15% born in the bulge and 85% in the disc [58].ρ d is taken to be constant over time, whileρ b is nonzero and constant between 10 and 8 Gyrs ago [60]. When a black hole is first born, it is given an initial velocity composed of a bulk velocity as well as an individual kick. The bulk velocity in the disk is a piecewise function of distance from the Galactic center, in the φ direction [61]: r < 0.2 225 + 15.625(r − 1.8) 2 0.2 < r < 1.8 225 + 3.75(r − 1.8) 1.8 < r < 5.8 240 r > 5.8 (20) where r is in kpc. The bulk velocity in the bulge is taken from a Maxwell-Boltzmann distribution with a mean of 130 km s −1 [62]. Any asymmetries in the supernova explosion can impart a "kick" to the newly formed compact object. The natal kicks of Galactic neutron stars are known to be as large as hundreds of km/s [63]. Since black holes are many times heavier than neutron stars, relatively smaller natal kicks are expected, depending on the kick mechanism [64,65]. Recent studies of known Galactic stellarmass black holes find that their properties are consistent with much smaller natal kicks than the ones attributed to neutron stars [66,67] (although [68] finds similar natal kick distributions for black holes and neutron stars).
The black holes' initial kicks in our simulation follow a Maxwell-Boltzmann distribution with average 3D velocity of 50 km/s, the smallest average kick velocity used in [56]. This means that black holes tend to stay in the Galactic component in which they were born, and few black holes migrate to high Galactic latitudes. We also consider the impact of larger kick velocities in App. B, finding that natal kicks of 100 km/s give reduced signal numbers due to black holes with larger kicks having larger distance d from the solar system, on average; 1/d ∼ (7.8 kpc) −1 for the 50 km/s population and 1/d ∼ (8.3 kpc) −1 for the 100 km/s population, corresponding to the signals from the latter population being about 10% weaker on average. See Appendix B for further discussion.
After a black hole is born in the simulation, its trajectory is determined by the Galactic gravitational potential. The gravitational potential of the bulge, disc, and halo (subscript "h") are defined as follows [61]: where R ≡ √ r 2 + z 2 , and the constants are,   Fig. 6, as a function of their distance from the Sun. The distribution is symmetric around zero km/s and is widest at 8 kpc, corresponding to the distance of the Galactic Center, as the bulk velocity in the bulge is uniformly distributed in direction. In contrast, since the bulk velocity in the disk is in the φ direction within the Galactic plane, the radial velocity of black holes in the disk (i.e., not near 8 kpc) has smaller magnitude. The scatter in velocities at large distances from the Sun is due to the black holes that have wandered into the halo.
as determined by fits to observed Galactic properties (e.g., rotation curves) in [61]. With this procedure, we obtain the final positions and velocities of 10 6 black holes with respect to the Galactic Center, but require a full population of 10 8 black holes with positions and velocities as measured from the Solar System Barycenter (SSB). Since our simulations assume axial symmetry (i.e., we do not include any structure like the spiral arms), we are free to choose the position of the SSB to be anywhere 8.12 kpc away from the center in the disk, with a velocity prescribed by Eq. (20). We therefore randomly select 100 random points from the circle of radius 8.12 kpc, as illustrated in Figure 6. Thus, we are able to build up a sample of 10 8 black holes at minimal computational cost increase, and for these black holes, we compute the distances d and radial and tangential velocities v rad and v tan as measured from the SSB.
We require that the individual black holes be independent; to check that the distribution is not oversampled, we estimate that on average, the distance between one SSB position assignment and the next is 500 pc. In comparison, the average separation between black holes is O(10 pc), assuming the 10 8 black holes are uniformly distributed in the Galactic plane. The average separation between black holes is therefore over an order of magnitude smaller than the average separation between SSB position assignments, and so the local black hole populations that are "seen" at given SSB position assignments are largely independent of one another.
While this procedure effectively reproduces the spatial and velocity distribution of 10 8 black holes using only 10 6 simulated objects, it does not accurately reproduce the age distribution of young black holes. Given the assumption of a constant formation rate in the last 8 Gyrs, the 10 6 population underestimates the number of black holes with age 10 4 years or less compared to the number that would be produced by a full simulation of 10 8 objects. This effect is relevant for heavier bosons, which have short radiation timescales compared to the typical age of a black hole, and will be discussed further in Sec. IV C.

B. Choice of mass distribution
For the mass distribution, we use the Salpeter function, dN/dM ∝ M −2.35 [69], an empirically determined function that has been shown to apply to Galactic stars, especially those more massive than the Sun [70]. We choose this distribution for simplicity, as the true mass distribution for isolated black holes is unknown. In applying the Salpeter function to black holes we implicitly ignore effects such as mass loss due to stellar winds, which has a larger effect on heavier stars. To better approximate the true mass distribution, we introduce minimum and maximum black hole masses for the population based on both observational and theoretical arguments. We use the same mass distribution for both the bulge and the disk [71], and assume it is unchanged over the lifetime of the Galaxy [70].
Two populations of stellar mass black holes have been detected. All known Galactic black holes reside in binaries, and most are observable in X-rays due to accreting material from their stellar companions; these black holes tend to have masses 20M (with the possible controversial exception of LB-1 [72,73]). In contrast, the mass distribution of black holes in merging binaries detected by LIGO and Virgo is consistent with a mass cutoff of approximately M BH ∼ 40M and a power-law index of around 2 [74][75][76]; however, it is difficult to extrapolate from the properties of the gravitationally detected black holes, as these are found in other Galaxies with unknown star formation histories and metallicities. Indeed, studies suggest that low metallicity environments are required to produce these more massive black holes [77,78]. In contrast, most of the black holes in our simulation are at low Galactic latitudes and close to the Galactic center ( Figure 8), where stars are metal-rich [79,80], although the metallicity could have been different in the past [81].
With these considerations in mind, we take [5M , 20M ] as the standard distribution of black hole masses (Table I). This range includes the known Galactic black holes. A maximum mass of 20M is also consistent with simulations, which suggest that the heaviest black hole that can form at Galactic metallicities is around 20M [78].
Our results are not sensitive to the choice of minimum black hole mass for all but the heaviest boson masses; the loudest signals come from the heavier black holes, and bosons heavy enough to produce large strains around light black holes source few observable signals. See App. C for more discussion on the effect of decreasing the minimum black hole mass.
The choice of maximum black hole mass more directly impacts our simulated signals, as both the signal strength and duration are highly sensitive to the black hole mass. Since the loudest signals come from systems with the heaviest black holes, our choice of 20M is conservative. To evaluate the dependence of our conclusions on the maximum mass, we also consider heavy distributions with a maximum mass of 30M (Table I).

C. Spin magnitude and orientation
The strain observed at the detector, h, will in general be smaller than the characteristic strain, h 0 , by a factor that depends on the geometry of the detector and the source. Relevant parameters are the inclination angle ι (or, more directly, cos ι) and the polarization angle θ, defined by the orientation of the spin axis of the black hole. We assume that the black hole spin axis is primarily determined by the angular momentum axis of the progenitor star. Since there is no known favored spin direction of stars in the Galaxy, we choose cos ι and θ to be uniformly distributed in [−1, 1] and [−π/4, π/4], respectively. We take φ 0 (the phase of the gravitational wave at a chosen reference time) to be uniformly distributed in [0, 2π].
The distribution of black hole spin magnitudes at birth is not well understood; in particular, there is no observational data for the spins of isolated Galactic black holes. Future analysis on the detectability of the boson ensemble signal in binaries would be interesting, as many high-spin channels are thought to occur in binary systems. Some 1D stellar evolution models indicate that a majority of black holes tend to be born with minimal spins (e.g. [82]); this is also compatible with the measure-ments of spins of pulsars extrapolated to their natal spins [16,83]. Models also predict a potentially small fraction (order 10 −3 -10 −2 ) of black holes to be born with very high spins (χ > 0.9) through the chemically homogenous channel, thought to be associated with gamma-ray bursts [84,85].
On the other hand, there are measurements of black hole spins in binary systems. Black holes in binaries observed to merge in LIGO have tended to low spins, although the error on measurements of initial spins is quite large; the inferred spin magnitudes of 90% of black holes are found to be below 0.6 +0. 24 −0.28 or 0.8 +0.15 −0.24 in the aligned or isotropic spin scenarios, respectively [75], and typical spin magnitudes have a preferred low central value of ∼ 0.2 ± 0.2 [76,86]. Black holes in X-ray binary systems have a range of spin values, from low to χ > 0.9, with order half with spin above 0.5 [16,17,87]. While accretion can contribute to the spin-up of a black hole, the stellar companions in most of these systems are not long lived or not massive enough to substantially increase the black hole spin from its natal spin [16,17,88].
Given these uncertainties and range of observations, we take a distribution uniform in the range [0, 1] as a standard and χ i,max = 0.5, 0.3 as moderate and pessimistic spin distributions, respectively (Table I). We also discuss our conclusions in the case in which a subpopulation of black holes has moderate or high spin.

D. Frequency derivatives
Any astrophysical signal has a small apparentḟ due to the proper motion of the source as observed at the detector. For distances that are large compared to the distance the black hole travels over the observational period (true for all of the systems we are interested in, which are at distances of at least a few pc and velocities of hundreds of km/s), this givesḟ app ≈ v 2 tan f /cd where v tan is the black hole velocity tangential to the line of sight. At the highest frequencies,ḟ app < 10 −19 Hz/s for the black hole population.
As discussed in Section II B, the decreasing mass of the boson cloud also causes a small spin-up. For a given M BH , the magnitude ofḟ is largest for large values of µ b . For the detectable systems in our search, the frequency derivatives fall in the range 10 −14 Hz/s ḟ < 6 × 10 −13 Hz/s for both the standard and heavy black hole populations ( Figure 22).
Theseḟ values are smaller than the minimum |ḟ | sensitivity of a typical semi-coherent all-sky search. In general, an all-sky search for the boson annihilation signal does not need to include frequency derivatives. On the other hand, follow-up searches of interesting signal candidates may need to takeḟ into account and could detect a small positive frequency drift, pointing toward the signal origin being a boson cloud.

E. Systems outside the Milky Way
We have focused on the signal from Galactic black holes, as they produce the loudest signals. We can also ask whether the black holes in the closest neighboring Galaxy, Andromeda, will produce a detectable signal. The distance to Andromeda is around 770 kpc [89]. The signal strains would drop by almost three orders of magnitude compared to what is plotted in Figure 4, so that systems in the top right corner of the plot (roughly speaking, µ b > 1.4 × 10 −12 eV and M BH > 15M would produce detectable signals. These systems correspond to α 0.17 ( Figure 1) which have timescales of τ GW < 10 4 years ( Figure 5). It is therefore unlikely for annihilation signals from Andromeda to be detectable by the current generation of detectors.
Closer to home, the Large (LMC) and Small Magellanic Clouds (SMC) are potential sources of signals. The contribution from either of these two dwarf galaxies would be an additional over-density of signals in a particular part of the sky and within a small range of Doppler shifts, since the black holes are traveling at the respective galaxy's bulk velocity (68 km/s for the LMC and 5 km/s for the SMC [90]). The LMC and SMC are an order of magnitude further than the Galactic Center [91,92], so the peak strains are an order of magnitude below the peak of the Milky Way distribution. Since the expected signal is small and the star formation histories of the LMC and SMC are very uncertain, we do not include them in our simulations.

IV. THE ENSEMBLE SIGNAL
In this Section we consider the ensemble of signals produced by annihilations of bosons with masses in the range 1 × 10 −13 to 4 × 10 −12 eV, in clouds around Galactic isolated black holes. To facilitate reproducibility and further studies, we provide the full set of simulated signal populations in an online repository 3 .
We focus here on the standard black hole population (Table I).
A. Signal simulation procedure Each of the 10 8 black holes is tagged with a unique position, velocity, and age. For a given run (i.e., choice of µ b ), we randomly assign a mass and spin to each black hole as described in Sec. III B and Sec. III C and summarized in Table I. For each boson-black hole system, we compare the initial spin χ i to the critical spin χ c (Eq. (4)) to determine whether the cloud forms. If χ i > χ c , we calculate the instability timescale τ inst (Eq. (5)) which defines the e-folding time for the cloud to grow, and take ln(N )τ inst ≈ 180τ inst as the timescale for the cloud to reach its maximum size (discussed in Sec. II A), keeping systems for which this time is less than the black hole age.
To check whether the second level has fully formed and therefore caused the CW emission from the first level to cease, we also calculate τ 022 inst for the (n = 0, = m = 2) level; for systems with black holes older than ln(N 022 )τ 022 inst , we set h 0 = 0 since the first level is no longer radiating.
For systems that pass these checks, we calculate h 0,peak and τ GW (Eqs. (A17), (A19)), and use the time evolution of h 0 (Eq. (14)) to determine its current value. The relevant time t is the time since cloud formation minus the time d/c required for the GW emission to travel the distance d to Earth; if the emission since cloud formation has not yet reached Earth, we set h 0 = 0. We compute the frequency f GW at cloud formation, Eq. (9), and apply the frequency drift (Eq. (16)) over the lifetime of the cloud, which is the black hole age minus ln(N )τ inst . Finally, we determine the apparent signal frequency by applying a Doppler shift where f GW,obs is the observed gravitational-wave frequency due to the source's radial velocity as measured at the SSB. The maximum value of |v rad /c| for our black hole population is 0.0025, and 90% of black holes have values of |v rad /c| < 6 × 10 −4 .

B. Ensemble signal properties
Ensemble signal shape: The source-frame frequency of the GW signal, f GW,source , has a small dependence on M BH through α, resulting in different 'potential energy' corrections for different black hole masses. For a fixed boson mass, increasing the black hole mass produces a signal with a lower frequency (Eq. (10)) and a higher peak strain (Eq. (13)). Thus, in a given ensemble, the signals with lower frequencies have higher strains at cloud formation. The signal half-time τ GW decreases rapidly with increasing black hole mass (Eq. (15)). If τ GW is shorter than the black hole age τ BH , the signal strain today is suppressed by τ GW /τ BH ; increasing the black hole mass above the value that gives τ GW ∼ τ BH increases the peak strain but decreases the strain at the current time. The current strain as a function of frequency for a fixed black hole age and distance as seen in the source frame is shown in the orange curves of Fig. 9.
The maximum frequency in an ensemble in the source frame is f 0 GW , while the minimum value depends on multiple factors: there are no signals below a certain source frequency if a) there are no heavier black holes due to the imposed maximum black hole mass cutoff; b) heavier black holes result in prohibitively large critical spins χ c  (Table I), and boson mass (left to right) µ b = 2 × 10 −13 , 4 × 10 −13 , 8 × 10 −13 , and 20 × 10 −13 eV. Top row: Each point represents the annihilation signal from an individual system, illustrating the range of strains and frequencies of the signals. The dashed horizontal red line corresponds to the approximate upper limit at that frequency reported by recent all-sky CW searches [20,21,46,93]. The thick dashed vertical red line marks the values of f 0 GW , and the dotted red line shows fGW calculated using the mode of the MBH distribution ( Figure 11). Second row: Distribution of strains as a function of distance and black hole mass. Smaller, whiter circles correspond to lighter black holes (minimum 5M ) and larger bluer circles correspond to heavier black holes (maximum 20M ); the differences in circle sizes are exaggerated for visual effect and do not scale linearly with MBH. The horizontal dashed red line again represents current search upper limits.
(Eq. (4)) hence systems with heavier black holes never form; or c) heavier black holes have short superradiance times for the = m = 2 level, τ 022 inst (Eq. (A7)), cutting off the signal before today. The ensemble-signal shape as seen in the source frame is shown in the gray dots in The line-of-sight BH velocity produces an additional Doppler shift for each source, proportional to v rad /c, which "smears" the ensemble signal's distribution in frequency (blue dots in Fig. 9). Nearby black holes produce FIG. 11. Each of the four panels shows the mass MBH and initial spin χi of the black holes that produce the ensemble signals for the four boson masses from Figure 10, assuming the standard black hole population (Table I). In each panel, the scatterplot shows MBH and χi for the systems that are detectable with the current search sensitivities (blue circles with orange outlines) and a factor of 2 improvement in sensitivity (blue circles). The histograms show the distributions of MBH (top) and χi (right) for these sets of systems. The black holes are assigned spins uniformly distributed in [0, 1] and only systems satisfying Eqn. (4) will form clouds. The upper-left set of plots (µ b = 2 × 10 −13 eV) contains no signals at these values of h0, suggesting that a larger increase in sensitivity is necessary to detect the annihilation signals produced by the lightest bosons.
stronger signals than farther ones, while farther black holes tend to have larger line-of-sight velocities (Fig. 7). This means that weaker (more distant) signals tend to have larger relative Doppler shifts. The result is that, at a fixed boson mass, the ensemble of signal frequencies observed at the detector (f GW,obs ) is shaped roughly like a triangle, with a wider base and a central 'peak'. In Fig. 10 (top row) we show the ensemble signals from the standard population of 10 8 black holes at a few illustrative boson masses.
The spread of f GW,obs is greater for heavier bosons for several reasons. First, heavier bosons produce GW emission with larger values of f 0 GW , and as both the Doppler shift (Eq. (24)) and potential energy shift (Eq. (10)) terms are proportional to f 0 GW , higher frequencies result in a larger absolute frequency spread, ∆f GW . In addition, increasing the boson mass from 2 × 10 −13 eV to 8 × 10 −13 eV increases the relative frequency spread, ∆f GW /f 0 GW for two reasons: 1) heavier bosons form systems with a larger range of black hole masses (Fig. 11), leading to a larger range of potential energy corrections (Eq. (10)), and 2) heavier bosons produce detectable signals at larger distances (Fig. 10, bottom row), resulting in larger Doppler corrections. Increasing the boson mass further starts to reduce these effects; the relative frequency spread decreases while the absolute frequency spread continues to increase. In all cases, the signal frequencies lie in the range f ∈ f 0 GW ×{1−0.01, 1+0.0025}, Each of the four panels shows the age τBH, and distance d of the black holes that produce the ensemble signals for the four boson masses from Figure 10, assuming the standard black hole population (Table I). In each panel, the scatterplot shows τBH and d for the systems that are detectable with the current search sensitivities (blue circles with orange outlines) and a factor of 2 improvement in sensitivity (blue circles). The histograms show the distributions of τBH (top) and d (right) for these sets of systems.
where the upper end of the range is given by the maximum Doppler shift in the ensemble and the lower end of the range is given by the largest negative Doppler shift plus the largest observed potential energy correction with α ∼ 0.25.

Maximum signal strength:
The power emitted by the boson cloud in gravitational waves is set by the total energy of the cloud, and the timescale τ GW over which the annihilations take place, P ∼ M cloud /τ GW , providing an estimate of the maximum strain in the detector today: Here we have chosen near-maximal parameters for the cloud and black hole mass, as well as the lightest axion mass that produces observable signals; note that M cloud and τ GW are not independent parameters and cannot reach the optimal values for most (µ b , M BH ) combinations. In addition, given the finite birth rate of black holes and their spatial distribution, the nearby black holes tend to be older, while the small number of young black holes are on average farther away (Fig. 12). We also see that black holes younger than 10 5 years contribute to the potentially detectable signal for heavier bosons, but the strain is reduced due to the higher signal frequency and farther typical black hole distance. Thus it is very unlikely to find strains larger than h 0 ∼ 10 −23 . The typical distribution from black holes at 1 kpc as well as the analytical strain amplitudes for different black hole age and spin parameter choices are shown in Fig. 9.
C. Properties of black holes in potentially detectable systems The detectability of a signal depends on the exact search involved. As a guideline, we take the approximate search sensitivities of recent all-sky searches in Advanced LIGO data [20][21][22] at the relevant frequencies.
Black hole distances: By construction, all the signals in the ensemble arise from black holes within the Milky Way (Fig. 12). The number of black holes increases rapidly as a function of distance, and in general, many distant black holes produce signals below the typical search sensitivity (Fig. 10).
If the maximum black hole mass M BH,max = 20M , then for the lighter bosons, the signals are very weak and are only observable from short distances (e.g., less than 2 kpc for µ b = 4 × 10 −13 eV). Increasing the boson mass increases the signal power, and at µ b = 8 × 10 −13 eV, signals at 10 kpc are observable, although most signals are still produced by systems within 2 kpc. If we allow for M BH,max = 30M , systems as far as 15 kpc are detectable for µ b ≥ 8 × 10 −13 eV (Fig. 33).
Black hole masses: We consider black holes with masses between 5M and 20M (30M ) in our standard (heavy) population. For lighter bosons, the signals are weak and the distribution of black holes is therefore peaked toward the heavier black holes (Fig. 11). The signal strain increases with increasing black hole mass, h 0 ∝ M 9 BH (Eq. (13)), and the total number of signals is thus sensitive to the upper cutoff of the mass distribution. This is especially evident for µ b = 2 × 10 −13 eV; while the standard population produces no signals with h 0 > 1.5 × 10 −25 (Fig. 11), the heavy black hole population produces tens of such signals (Fig. 32).
In contrast, for heavier bosons, there are fewer signals from heavier black holes because of prohibitively high critical spins χ c or short instability (τ 022 inst ) and gravitational wave (τ GW ) timescales, as discussed in Section IV B. Because of these effects, for µ b 8×10 −13 eV, the fraction of signals with black holes masses of 20M or above goes to zero (Fig. 11).
The effect of the minimum black hole mass is relatively weak: we find that 5−6M black holes do not contribute any signals above current search sensitivities for boson masses up to 8×10 −13 eV. Black holes as light as 5M can form systems with α > 0.075 for µ b 2 × 10 −12 eV, and  4)). Since χc increases with α, this also results in fewer signals above a given strain for heavier bosons.
comprise ∼20% of the population of detectable signals at µ b ∼ 2 × 10 −12 eV. Thus, decreasing M BH,min would give a slightly larger number of signals for boson masses 2×10 −12 eV. However, the instrument sensitivity is lower at the high frequencies corresponding to these heavier bosons, and few clouds are radiating today (Fig. 14).
Black hole spins: We find that only systems with initial spins of χ i > 0.25 produce strains above current search sensitivities (Figs. 11,13). Systems with χ < 0.25 require α 0.06 for the cloud to form, producing signals that are too weak to observe. For µ b ≥ 8 × 10 −13 eV, there is a trend toward higher black hole spins as the black hole mass increases (Fig. 11) due to the critical spin increasing for larger α.
We show the number of signals with h 0 > 10 −25 as a function of χ i,max for six example boson masses in Fig. 13. The lightest boson mass produces few to no signals with h 0 > 10 −25 ; for the heavier five boson masses, the number of signals with h 0 > 10 −25 drops by at least an order of magnitude in reducing the maximum spin from χ i,max = 1.0 to χ i,max = 0.5 and two orders of magnitude when χ i,max decreases to 0.3. For µ b = 2 × 10 −12 eV and 3.5 × 10 −12 eV, there are no signals at or below χ i,max = 0.3 and 0.5, respectively, since these critical spins would require black holes with mass below 5M to form clouds.
Black hole ages: Given the approximately constantin-time black hole formation rate, young black holes are proportionally less common than old ones. For the lighter bosons, e.g. µ b = 4 × 10 −13 eV, detectable signals arise from heavier black holes (Fig. 11) which are relatively rare; since black holes that are both heavy and young are doubly uncommon, the black holes that produce detectable signals for this boson mass are predominantly 10 8 -10 10 years old, Fig. 12. Bosons with µ b = 8 × 10 −13 eV produce detectable signals around black holes with masses across the entire range, and have intermediate signal times, so these black holes span the entire range of ages. Finally, for the heaviest bosons, e.g. µ b = 2 × 10 −12 eV, the shortness of the gravitationalwave timescale (Eq. (15)) and the instability timescale of the second level (Eq. (17)) means that very old systems have stopped radiating, so that only young black holes produce detectable signals.
As discussed in Sec. III A, our method of producing 10 8 black holes results in fewer young black holes than the population would statistically contain. Specifically, on average we expect O(10) black holes of τ BH ∼ 10 3 years, and O(1) black holes of τ BH ∼ 10 2 years. On the other hand our example population overestimates the number of 10 4 year black holes, as the initial population included an upward fluctuation of the number of these black holes. Thus our population includes 200 black holes that are younger than 10 4 years old, compared to an average expectation of 62. This disproportionately affects heavier bosons; at µ b = 2 × 10 −12 eV, of order half of the signals that are detectable with current search sensitivities are produced by ∼ 10 4 year old black holes. Therefore, our sampling procedure introduces an order-1 uncertainty in the number of signals produced by bosons with µ b ≥ 2 × 10 −12 eV. Figure 15 shows the total number of signals with characteristic strains above given values of h 0 ∈ [10 −26 , 10 −23 ] for six reference boson masses. For most values of h 0 , the number of signals above a given h 0 is highest for µ b = 8 × 10 −13 eV (f GW ≈ 400 Hz).

D. Number and density of signals as a function of GW amplitude
For the "standard" black hole population (M BH,max = 20M ), annihilation signals from the lightest bosons are unlikely to be detectable by current CW searches: for µ b = 2 × 10 −13 eV (corresponding to f GW ∼ 100 Hz) few or no systems produce GW emission with intrinsic amplitude above 10 −25 . The loudest signal in our example has a strain a factor of a few below the current sensitivity thresholds (Fig. 10).
At fixed black hole mass, increasing the boson mass increases the strain as µ 9 b . Thus, increasing the boson mass by a factor of 2 to µ b = 4 × 10 −13 eV (f GW ≈ 200 Hz), greatly increases the number of signals that would be detectable by current searches (h 0 > 2 × 10 −25 near this frequency).
Although the peak strain increases with boson mass, fewer signals are expected at heavier boson masses: Fig. 14 illustrates that both the percentage of systems with fully formed clouds and the percentage of systems that are still radiating CWs today decrease with boson mass. The further effect of the signal half-time is not included in Fig. 14 as it depends on the search strain sensitivity and the black hole's age. Fig. 15 shows that   The number of signals with intrinsic amplitude above a given h0 value. The signal number is largest at µ b = 8×10 −13 eV (fGW ≈ 400 Hz) for the standard black hole population (Table I). At higher frequencies, the number of signals decreases because both the signal half-time (Eq. 15) and the signal cutoff time due to the formation of the second level (Eq. (17)) decrease with increasing frequency. The contours show Gaussian statistical uncertainties, an underestimate at small signal number.
for bosons with µ b > 8 × 10 −13 eV, the number of sig-  (Table I). For heavier bosons, the overall density decreases; the number of signals above a given h0 drops due to the combination of effects described in Fig. 15 while the frequency range spanned by the ensemble signal continues to increase. The contours show the uncertainty √ N /∆f , where N is the number of signals and ∆f is the frequency range spanned by the N signals. Note that this does not take into account the additional uncertainty from ∆f , which is itself correlated with N .  17. The maximum density of signals above a given h0 is an order of magnitude larger than the mean density (Figure 16). Frequency ranges with sufficiently high signal densities can cause complications for standard searches for CW signals (Sec. V). The contours show the uncertainty √ N /∆f , where N is the number of signals and ∆f = 0.01 Hz is the size of the frequency bin over which the maximum density is calculated.
nals above a given h 0 begins to decrease with increasing boson mass. Figures 16 and 17 show the mean and maximum den-sity of signals in frequency space above a characteristic strain h 0 , for the same six reference bosons as in Figure 15. Both the mean and maximum signal density are largest for µ b = 5 × 10 −13 eV. This is in contrast with the fact that the number of signals is typically larger for 8 × 10 −13 eV (Fig. 15). Signals from heavier bosons are spread over a larger frequency range as described in Section IV B, thereby decreasing the signal density.

V. SIGNAL DETECTABILITY
We now consider the signal detectability in the context of the results of current CW searches with LIGO data, using several black hole mass and spin distributions (Table I ).

A. Number of detectable signals
For each boson mass in the range µ b ∈ [1 × 10 −13 , 4 × 10 −12 ] eV we calculate the number of detectable signals for the five black hole populations (Table I) at the corresponding frequencies f GW,obs ∈ [50,2000] Hz. For a signal to be detectable, we require that its intrinsic strain h 0 is greater than the 95% upper limit value in the {f GW,obs } frequency bin. Finally, we multiply the total number of detectable signals for each boson mass by 95%.
We use the 95% upper limits for the low-ḟ CW search in the 20-600 Hz range [20,21] and the upper limits for the CW search in the 10-2048 Hz range [22]. For ease of reading, we refer to the searches [20,21] and [22] by their nicknames "Falcon" and "Freq Hough," respectively. For bosons with µ b < 1.2 × 10 −12 eV we use the near-zero spin-down upper limits from the Falcon searches and for µ b ≥ 1.2 × 10 −12 eV we use the upper limits from the very broad spin-down Freq Hough search.
For each boson mass, we calculate the number of detectable signals for the five black hole populations (Table I) at the corresponding frequencies {f GW,obs }. For a signal to be detectable, we require that its intrinsic strain h 0 is greater than the 95% upper limit value in the {f GW,obs } frequency bin. Finally, we multiply the total number of detectable signals for each boson mass by 95%.
The number of detectable signals for the five black hole populations is shown in Fig. 18 [20][21][22], for bosons in the mass range 1 × 10 −13 eV to 4 × 10 −12 eV. The Falcon search [20,21] provides 95% confidence strict upper limits on CW emission with near-zero frequency derivative between 20 and 600 Hz. These upper limits are valid for any sky position and polarization of the source, even the ones with the most unfavorable coupling to the detectors. In contrast the upper limits from the Frequency Hough search [22] hold for 95% of the entire population. For this reason, even at the same sensitivity, the Frequency Hough upper limits are lower than the Falcon upper limits. We use the Falcon upper limits for bosons with µ b < 1.2 × 10 −12 eV and the Frequency Hough upper limits for bosons with µ b ≥ 1.2 × 10 −12 eV. We consider the standard (black circles) and heavy (blue circles) black hole populations, as well as the moderate (gray squares), heavy moderate (blue squares) and pessimistic (gray triangles) spin populations. The fluctuations are due to stationary lines and other artifacts in the detectors, which decrease the CW search sensitivity in nearby frequency bins.
black hole population (M BH,max = 20M ). For the heavy black hole population, the largest number of detectable signals occurs for µ b = 4 × 10 −13 eV; for the lighter black hole population, the largest number corresponds to µ b = 6 × 10 −13 eV. In the case that M BH,max < 20M , fewer signals would be produced at the lightest boson masses. At a fixed value of α, the strain is linearly proportional to black hole mass, h 0,peak ∝ M BH , so even decreasing the maximum black hole mass by a factor of 1.5 (already in tension with known Galactic black holes [94]) would shift the peak of the detectable signals to bosons a factor of 1.5 heavier.
For µ b > 8 × 10 −13 eV, the two black hole populations produce the same number of detectable signals. For these heavier bosons, systems with heavier black holes form less frequently and have significantly shorter emission lifetimes and do not contribute to the event rates.
For a given black hole mass distribution, the number of detectable signals is most sensitive to the black hole spin distribution, through the superradiance condition of Eq. (4). We consider three different uniform spin distributions: a reference distribution with χ i ∈ [0, 1], a moderate distribution with χ i ∈ [0, 0.5], and a pessimistic distribution with χ i ∈ [0, 0.3].
The reference and moderate spin distributions produce detectable signals for boson masses 2 × 10 −13 eV µ b 2.5 × 10 −12 eV, although the number of detectable signals decreases by an order of magnitude from χ i,max = 1 to 0.5, as fewer black holes have initial spins above the critical spin. For the pessimistic spin distribution, only a handful of signals are detectable, and only for boson masses 4 × 10 −13 eV µ b 1.3 × 10 −12 eV. This, in conjunction with Fig. 13, suggests that if all black holes are born with spins smaller than 0.3, no bosons will produce signals detectable by current CW searches.

B. Density of signals
The number of detectable signals in a parameter space search cell determines the best method for signal detection. If the total signal is an incoherent superposition of many signals that cannot be resolved individually, then the total signal is a stochastic process that can best be detected with cross-correlation methods [95]. On the other hand, if the signals are sparse enough that they can be individually resolved, the matched filtering and semicoherent techniques that are commonly used to search for CW signals should be applied [20,[96][97][98]. In between these two extremes is a regime with some loud and resolvable signals on a background of weaker and unresolved ones, that may act like confusion noise.
We now want to understand what situation applies here. To address this question we consider the FreqHough [22] and Falcon searches [20,21].
We bin the signal frequency with the resolution of the first stage of the respective CW search. We compute the ensemble signal for a standard black hole population and count the number of signals in each frequency bin with amplitude equal to or greater than the amplitude upper limit from the given CW search, as a function of the signal frequency. For every bin we also determine the mean and maximum signal amplitude over the ensemble signals with frequency in that bin and amplitude greater than the search upper limit amplitude. Figures 19,20 show these quantities for µ b = 4 × 10 −13 eV which is close to the boson mass yielding the highest number of detectable signals per bin. The maximum number of signals in a single FreqHough bin is ∼ 10, whereas for Falcon it is 2. This difference is mostly due to the fact that the FreqHough frequency bin is 14 times larger than the Falcon bin.
Although there are many FreqHough signal-frequency bins that contain multiple signals, these signals are in principle still resolvable by the search because they come from different sky locations. However, the high concentration of detectable signals in many neighboring frequency bins calls for caution in applying standard CW search techniques: Due to large-scale parameter correlations of the detection statistics, sufficiently strong signals produce high detection statistic values even at template parameters far from the signal parameters. This is evident for instance in [99]; in spite of this being a search aimed at the Galactic Center, it still detects very well all LIGO hardware-injected fake signals at sky locations far from the Galactic Center. This means that when the ensemble signal comprises many loud signals in nearby frequency bins, even if each signal comes from a different sky location, it will contribute to the detection statistics at the templates of other signals, giving rise to a sort of  [20,21] than in the FreqHough search [22] is that the frequency bin of FreqHough is significantly larger (by a factor of 14) than the frequency bin of Falcon. The maximum initial spin of the black hole population is taken to be 1.

contamination/confusion noise.
A large signal density has consequences because the complex all-sky CW search pipelines [20][21][22][23]100] comprise several follow-up stages with increasing coherence lengths, with vetoes and thresholds designed to reject noise fluctuations and coherent disturbances while preserving signal candidates. The procedures are carefully calibrated with extensive Monte Carlos on fake isolated signals added to real noise. It is unlikely that settings that have been tuned to isolated signals would be equally effective for the ensemble signals explored in this study, especially in regions of high signal density. We expect FIG. 20. Same plots as those of Fig. 19 but having assumed that the maximum initial spin of the black hole population is 0.5.
that the detection efficiency of searches for which the density and number of detectable signals per bin is larger will be impacted more severely; so, for example, the impact will be greater for the Freq. Hough search [22] than the Falcon search [20,21].
We note that in Fig.19 the frequencies where the signals are the strongest are also the frequencies with the highest density of signals. This is not the case in general. As discussed in Section IV B, the signal strain tends to increase toward lower frequencies; however, the number of heavy black holes decreases due to the black hole mass distribution, and the number of clouds formed decreases as well, so the number of signals per frequency bin can be larger for relatively lighter black holes (and thus relatively higher frequencies). Thus the strongest signals are at or below the frequencies of the densest signals, as shown in Appendix E.
A high density of loud signals in a small frequency range could also have a more subtle negative impact on the ability of standard CW searches to detect signals from the ensemble. The reason is that the average noise level data is "normalized out" of the detection statistic, but if this is not done carefully the danger is there of removing the signal: To address the non-stationary character of the noise, most CW search methods (and all of the ones that we refer to in this paper) whiten the data at an early point in the analysis by dividing the Fast Fourier Transform of the data taken over some short time-baseline with the amplitude spectral density estimated from the data itself. Under the assumption that the data is noise-dominated, with at most a single detectable signal, the noise level is estimated through an average of the amplitude spectral density in frequency: since the signal is concentrated in a frequency bin or two, the "average" over tens of bins is insensitive to its presence. For instance, FreqHough [22] searches use an autoregressive average on the power spectral density over 0.02 Hz [101,102] and the Einstein@Home searches [23,100] use running medians over 0.06 Hz frequency bins.
If the ensemble signal is a collection of many signals distributed over multiple frequency bins loud enough to increase the apparent noise level at those frequencies, then the normalization procedure will down-weight the signals and degrade the detection efficiency with respect to when the signal is concentrated in a single frequency bin. In Figure 21 we plot the amplitude spectral density (ASD) expected for the ensemble signal corresponding to µ b = 7 × 10 −13 eV. There is a clear excess power due to the ensemble signal extending over 1 Hz, which is much larger than any of the frequency ranges used to estimate the noise. This means that the normalization factors would include the signal and would significantly impact the detection ability of the search. More plots like that of Fig. 21 are given in Appendix F, for different combinations of boson mass, maximum black hole mass and maximum initial black hole spin. For 3 ≤ µ b ≤ 7 ×10 −13 eV, this effect is certainly not negligible, except for ensembles with initial maximum black hole spin 0.3 or less and the maximum black hole mass of 20M .

VI. EXISTING LITERATURE AND RESULTS
Several works have made projections of the detectability of the boson annihilation signal at Advanced LIGO and Virgo [12][13][14]. These studies predict tens to thousands of signals detectable with Advanced LIGO/Virgo CW searches, with a range of assumptions about black hole mass and spin distributions in combination with often only idealized search concepts. Other works have searched the gravitational wave data for boson annihilation signals leveraging existing search pipelines and/or results from other astrophysical searches.
all black holes are formed with spin 0.6, and a range of 3-100 M black hole masses in both cases. We find that the number of signals below ∼ 3 × 10 −13 eV falls off very rapidly for both M max,BH = 30M to M max,BH = 20M (Figure 18), indicating that the excluded masses of [22] rely on a speculative population of heavy black holes. Our analysis of the ensemble signal shows that any bound derived from a single spin or distance assumption can result in misleading conclusions. Furthermore, the ensemble signal from bosons with mass between 3 and 4×10 −13 eV would likely be detected with a lower efficiency than that of an isolated signal and hence non-detections from a standard CW search cannot automatically be translated in exclusions for boson signals at the same frequency at the same h 0 level. A cross-correlation search technique is considered in [19]: the boson annihilation ensemble signal from outside the Milky Way is treated as a stochastic background [15]. Based on their null results, [19] exclude the range 2.0-3.8 ×10 −13 eV under an "optimistic" black hole spin distribution (defined as unity maximum dimensionless spin and letting the minimum vary) but not under a "pessimistic" one (zero minimum spin to and letting the maximum vary). The range of black hole masses considered is 3-100M . The amplitude spectral density estimates based on our ensemble signal suggest that their null result could exclude Galactic signals from some boson masses even for maximum initial black hole spins of 0.5, as long as the black hole masses go at least up to 30M .
A way to bypass the uncertainties in the system's parameters is to consider a known black hole, for example the remnant of a merger of two compact objects [13]. For the LIGO black holes the formation time is precisely known, and the spin and mass are measured with good precision but the remnant is too far to yield a signal detectable by the current interferometers [13,103]. Turning to galactic sources, [104] have searched for CW emission from Cygnus X-1, a well known system: a black hole of ∼ 15M , 5 million years old at a distance of 1.86 kpc. The absence of a detection disfavors the existence of bosons in the range 5.8-8.6 × 10 −13 eV. This result however assumes a near-extremal initial black hole spin of 0.99 and implies a low current spin, χ ∈ [0.25, 0.36], which is in tension with continuum and reflection measurements of the spin, χ 0.95 [105,106].

VII. CONCLUSIONS
In this work, we simulate the ensemble annihilation signal from a population of isolated Galactic black holes and predict the signal detectability using recent all-sky CW searches. We study the dependence of the ensemble signal on black hole population properties. We show how population assumptions, in particular on black hole mass and spin distributions, can strongly affect the detectability of the ensemble signal. We propose that the results of future boson annihilation signal searches are interpreted using the ensemble-signal paradigm that we present here. We make the ensemble signal population parameters used in this work publicly available to facilitate future studies.
Our analysis is the first study to characterize the unique shape of the ensemble signal -constructed from a range of distributions of the underlying parametersand suggests a clear way to distinguish between these signals and the signals produced by rotating neutron stars.
We have considered the signals from five different black hole populations. We use two mass distributions: the 'standard' black hole mass distribution (maximum black hole mass M BH,max = 20M ) covers the range of masses for known Galactic black holes [94] and is therefore conservative, while the heavy black hole population (M BH,max = 30M ) allows us to test how the properties of the ensemble signal depend on the heaviest black holes. We also consider three different initial spin distributions: the optimistic spin distribution of χ i ∈ [0, 1], the moderate spin distribution of χ i ∈ [0, 0.5], and the pessimistic spin distribution of χ i ∈ [0, 0.3]. For all populations, we assume a total of 10 8 Galactic black holes. The predicted number typically varies between 10 7 and 10 8 and is highly uncertain, as there are no direct observations; the number of signals scales proportionally to black hole number.
Assuming there are 10 8 isolated black holes in the Galaxy with M BH ∈ [5M , 30M ] and initial spin χ i ∈ [0, 1], bosons with masses ∼2.5-−17 × 10 −13 eV produce 100 or more signals with amplitude above the upper limits in the O1/O2 CW searches, and masses ∼ 4.5 × 10 −13 eV produce over 1000 signals above those upper limits [20,21,46,93]. For a lighter population of black holes (M BH ∈ [5M , 20M ]), the boson mass range yielding signal amplitudes above existing upper limits is not significantly different than the range for the heavier population -the lowest boson mass yielding 100 or more signals increases to ∼ 3.5 × 10 −13 eV -and the highest number of signals is a factor of ∼2 lower.
If we take the maximum black hole spin at birth to be 0.5 rather than 1.0, keeping M BH,max = 20M , we find that the number of detectable signals drops by an order of magnitude. This suggests that black holes with χ i > 0.5 produce 90% of the "detectable" signals, so even if only 10% of black holes are born with χ i > 0.5, we predict tens to hundreds of "detectable" signals for 3 × 10 −13 eV µ b 17 × 10 −13 eV. The number of "detectable" signals drops by at least another order of magnitude if the maximum spin at birth is 0.3, and black holes with χ i between 0.3 and 0.5 produce ≈ 10% of the detectable signals. If instead all black holes are born with spins 0.3, we expect ten or fewer signals across the entire frequency range. Given other uncertainties such as the overall number of Galactic black holes, in this regime even a single detectable signal is not assured.
We caution that for boson masses in the range ≈3-10 × 10 −13 eV, the sensitivity of standard CW search methods has to be reassessed: Search methods that have been carefully tuned for the regime of quiet and very sparse signals should be characterized again on a dense ensemble of loud signals. The interplay between the down-weighting of the signal from the normalization of the data and the overlap in parameter space of the detection statistic from different signals in the ensemble and how these elements factor in a multi-stage follow-up process, needs to be investigated. For boson masses between 4-6 × 10 −13 eV, even for moderately rotating black hole populations, the GW ensemble signal is very prominent, which is a dramatically different regime than the one assumed by standard CW searches.
While the detection of a single particularly loud signal from a nearby black hole could serendipitously occur even with a highly "mismatched" search, ultimately the identification of a boson annihilation signatureas opposed to a CW signal from a compact rotating object -lies in the identification of the ensemble signal. Searches for gravitational-wave signals from boson clouds around black holes are only starting to be explored. This work lays the foundations for this type of detectability study.

Superradiance rate
The instability rate of the boson is ( [36], [18]) where χ the dimensionless black hole spin and C n m (α, χ) 1 is a numerical coefficient which is a function of the bound state quantum numbers as well as (α, χ), We define α ω ≡ αω R /µ, i.e., the α corrected for the potential energy contribution to the boson energy. The fastest-growing level for a light scalar, the n = 0, = m = 1 level, has ω 011 sr 1 24 χα 8 µ b . The rate for = 1 is larger than that for = 2 by a factor of ∼ C 011 /C 022 α −4 ∼ 10 3 α −4 .
The instability e-folding timescale is then given by with ln(N ) = ln(M BH c 2 /µ b ) ∼ 180 instability timescales required to fully saturate the cloud growth. The analytical form is accurate to within a factor of 2 for α 0.25 and within 50% for α 0.2 4 .
The maximal peak strain is related to the power emitted by the axion cloud as, where the two gravitational wave polarizations are given by As the power is emitted from the vicinity of a Kerr black hole, the angular dependence is not exactly as defined in Eq. (A18) but is instead specified by spinweighted spheroidal harmonics, which depend on α and the black hole spin χ [14,108]. However, for α, χ f corresponding to α 0.3 (0.2), the standard quadrupolar emission is an excellent approximation to within 5% (3%) or better in h × and to within 10% (5%) or better in h + . Given that the exact angular power calculation is computationally intensive and the CW pipelines are optimized for strain angular dependence according to Eq. (A18) we neglect the extra effect of spin in our analysis. We also focus only of the˜ =m = 2 gravitational wave mode, which dominates the total power for χ ∼ χ c and α 0.35 [44].
We consider only the gravitational wave emission from the n = 0, = m = 1 cloud; this bound state produces the largest strain: using the values of [44], the strain of the first level is larger than that of the second by h 0 011 /h 0 022 ∼ A 011 /A 022 α −2 ∼ 90. The time evolution of the signal as the cloud depletes through GW radiation is also related to the power emitted and is given by The calculations of the power have been performed in the point particle approximation, i.e. the back-reaction of the cloud on the metric is neglected. During the process of superradiance, the black hole spins down and loses mass to the cloud; thus the emission is taking place approximately in a background defined by the final black hole mass and we use the final, smaller, value of α to evaluate the strain and timescale expressions. The presence of the cloud may be viewed as an additional contribution to the mass in the Kerr metric, in which case the initial value of α could be used as an approximation [107]. Our expression is more conservative (as the final α is smaller and thus the power is reduced) and is correct in the limit when the cloud has and/or depletes to masses much smaller than the black hole mass, as is the case for much of our old, long-lasting signals. The difference between using the initial and final α gives approximately a 50% change in the total power, which is a conservative estimate of the overall uncertainty in the rate.

Frequency drift
As the cloud annihilates to gravitational waves, its gravitational energy decreases, leading to a positive frequency drift, There is a larger, negative frequency drift at early times as the cloud is growing, but we neglect this in our analysis as the strain is small at these times. In Figure 22 we check that theḟ caused by the cloud's decreasing mass is still smaller than can be resolved by the first stage of current all-sky CW searches. This is true for both the standard (M BH,max = 20M ) and heavy (M BH,max = 30M ) populations. In both cases, the maximumḟ from the changing cloud mass is 6 × 10 −13 Hz/s.

Appendix B: Higher natal kicks
We examine the ensemble signals that would be produced by black holes assuming they are born with natal kicks of average 100 km/s. The black holes in the 100 km/s population move somewhat faster on average, but the apparentḟ due to the proper motion is still many orders of magnitude smaller than anything resolvable by CW searches and is not a concern.
Black holes born with faster natal kicks are more likely to have sufficiently high speeds to escape the Galactic bulge and disk and travel further on average. Thus, they are more likely to be found at larger distances from both the Galactic Center and the Earth (Fig. 23). Fig. 24 shows the ratio of the sum of 1/d for all the black holes within a distance d; this is smaller than one and decreases with decreasing distances from Earth until very small distances, which are dominated by small number fluctuations. The ratio of the sum of 1/d approximates Fig. 25, the ratio of the number of signals above a given value of h 0 , which is in general less than one and decreases with increasing h 0 until the loudest signals (h 0 > 10 −24 ), a regime that is dominated by the few closest black holes. The smaller number of signals above a given h 0 for the population with larger natal kicks also produces smaller maximum (Fig. 31) and mean densities (Fig. 30).

Appendix C: Light black hole population
We compare the ensemble signals from the standard black hole population (M BH ∈ [5, 20]M ) with the ensemble signals from a lighter black hole population (M BH ∈ [3, 20]M ). We maintain the distribution shape (Salpeter function) as well as the total number of black holes (10 8 ); reducing the minimum black hole mass therefore reduces the number of black holes of all other masses.
For µ b = 3.5 × 10 −12 eV, the addition of the lighter black holes results in a relatively larger number of signals for signals with h 0 3 × 10 −25 ; a similar but smaller effect is seen for µ b = 2 × 10 −12 eV. This is due to the fact that at these heavy boson masses, lighter black holes can more easily satisfy the critical spin condition while still producing relatively large signals. However, the current CW search upper limits are around 10 −24 near the signal frequencies produced by these heavy bosons, so the addition of the lighter black holes does not affect our estimates of the number of signals detectable by current (This ratio is affected by small number statistics at distances much less than 1 kpc.) Overall, this ratio is smaller than unity and decreasing until 1 kpc, consistent with the fact that the population of black holes with faster kicks in general produces fewer signals above a given h0 until h0 ∼ 10 −24 (Fig. 24). CW searches. The overall effect is that the total number of signals decreases by about a factor of two for almost all boson masses except for the heaviest ones (Fig. 28). In reality, the black hole mass distribution should turn over at the lightest masses [94,109]. Therefore, Fig. 28 can be considered a lower bound on the ratio of the number of signals.  The ratio of mean signal density above a given h0 between the faster and slower kicks populations reflects the number of signals above a given h0 (Fig. 25). In general, the black holes with faster natal kicks produce signals with larger Doppler shifts, thereby increasing the frequency range spanned by the ensemble and decreasing the mean signal density further.

Appendix D: Heavy black hole population
We examine the ensemble signal for a heavy population of black holes, with M BH,max = 30M , and compare the resultant ensemble signals with the signals from the standard population.
We show the ratio of the heavy to the standard population in the number of signals (Fig. 29), mean density of signals (Fig. 30), and maximum density of signals (Fig. 31). The ratios are generally consistent with one. The exceptions are the number and density of signals for the lightest boson mass of 2 × 10 −13 eV, which is a factor of several hundred larger for the heavy population, and   The ratio of the mean signal density above a given h0 for the heavy versus standard black hole populations is approximately unity for the five heavier boson masses. For the lightest mass (µ b = 2 × 10 −13 eV), the mean density is over ten times larger for the heavier black hole population due to the . Appendix F: Contribution to the noise amplitude spectral density from the ensemble signal These plots show the expect amplitude spectral density from the ensemble signal from a Fourier transform of varying time-baseline (T SF T ), under different assumptions on the population parameters, in particular the maximum black home mass and maximum initial spin. The T SF T at the different frequencies is chosen to reflect existing search strategies. The values of these parameters is indicated in the plot title.  [22].