Ultralight scalar fields can experience runaway superradiant amplification near spinning black holes resulting in a macroscopic axion cloud which slowly dissipates via continuous monochromatic gravitational waves. For a particular range of boson masses mathcalO1011  1010 eV an axion cloud will radiate in the 10  100 kHz band of the Levitated Sensor Detector LSD. Using fiducial models of the mass spin and age distributions of stellarorigin black holes we simulate the presentday Milky Way population of these hypothetical objects. As a first step towards assessing the LSDs sensitivity to the resultant ensemble of GW signals we compute the corresponding signaltonoise ratios which build up over a nominal integration time of 107 s assuming the projected sensitivity of the 1m LSD prototype currently under construction as well as for future 10m and 100m concepts. For a 100m cryogenic instrument hundreds of resolvable signals could be expected if the boson mass mu is around 3times1011 eV and this number diminishes with increasing mu up to approx 5.5times1011 eV. The much larger population of unresolved sources will produce a confusion foreground which could be detectable by a 10m instrument if mu in 34.5times1011 eV or by a 100m instrument if mu in 36times1011 eV.
['Astrophysics - High Energy Astrophysical Phenomena', 'General Relativity and Quantum Cosmology']
Simulating the Galactic population of axion clouds around stellarorigin black holes Gravitational wave signals in the 10100 kHz band
{'Simulating the Galactic population of axion clouds around stellar-origin black holes: Gravitational wave signals in the 10 - 100 kHz band': "Jacob R. Sprague 1 , Shane L. Larson 1 , Zhiyuan Wang 2 , Shelby Klomp 2 , Andrew Laeuger 3 , George Winstone 2 , Nancy Aggarwal 4 , Andrew A. Geraci 2 , and Vicky Kalogera 1 ∗  \n(The LSD Collaboration) 1 Center for Interdisciplinary Exploration and Research in Astrophysics (CIERA), Department of Physics and Astronomy, Northwestern University, Evanston, Illinois 60208 2 Center for Fundamental Physics, Department of Physics and Astronomy, Northwestern University, Evanston, Illinois 60208, USA 3 The Division of Physics, Mathematics and Astronomy, California Institute of Technology, Pasadena, California 91125, USA and 4 Department of Physics and Astronomy, University of California, Davis, Davis, California 95616, USA (Dated: October 29, 2024)  \nUltralight scalar fields can experience runaway 'superradiant' amplification near spinning black holes, resulting in a macroscopic 'axion cloud' which slowly dissipates via continuous monochromatic gravitational waves. For a particular range of boson masses, O (10 -11 - 10 -10 ) eV, an axion cloud will radiate in the 10 - 100 kHz band of the Levitated Sensor Detector (LSD). Using fiducial models of the mass, spin, and age distributions of stellar-origin black holes, we simulate the present-day Milky Way population of these hypothetical objects. As a first step towards assessing the LSD's sensitivity to the resultant ensemble of GW signals, we compute the corresponding signal-to-noise ratios which build up over a nominal integration time of 10 7 s, assuming the projected sensitivity of the 1-m LSD prototype currently under construction, as well as for future 10-m and 100-m concepts. For a 100-m cryogenic instrument, hundreds of resolvable signals could be expected if the boson mass µ is around 3 × 10 -11 eV, and this number diminishes with increasing µ up to ≈ 5 . 5 × 10 -11 eV. The much larger population of unresolved sources will produce a confusion foreground which could be detectable by a 10-m instrument if µ ∈ (3 -4 . 5) × 10 -11 eV, or by a 100-m instrument if µ ∈ (3 -6) × 10 -11 eV.", 'I. INTRODUCTION': "The era of gravitational-wave (GW) astronomy is in full-swing. During their first three observing runs, the GW interferometers Advanced LIGO and Advanced Virgo detected 90 compact binary coalescences (CBC) involving neutron stars (NS) and stellar-mass black holes (BH) [1][2][3] . The most notable events included the first NS-NS merger (GW170817) [4], the first highlyasymmetric binary (GW190412) [5], the first merger with an intermediate-mass BH remnant (GW190521) [6], and the first object in the mass gap separating the most massive neutron stars from the lowest-mass BH's (GW190814) [7] . The first half of the fourth observing run has already seen a new lower-mass-gap event (GW230529) [8].  \nAdding to the excitement, evidence for a stochastic background has been reported in the 15-year dataset from the North American Nanohertz Observatory for Gravitational Waves (NANOGrav) [9] . The most wellmotivated scenario for the origin of this background is the extragalactic population of inspiralling supermassive BH binaries.  \nFinally, the launch of the Laser Interferometer  \nSpace Antenna (LISA) in the mid-2030's will open up the millihertz band for exploration. The Galactic population of compact binaries, and the extragalactic population of supermassive BH binaries and extrememass-ratio-inspirals (EMRI's), are all highly anticipated LISA sources [10].  \nThese observatories cover multiple windows in the GW spectrum from the nanohertz up to several hundred Hz. The push to higher frequencies is now underway, with cosmic strings, axion clouds, primordial black hole (PBH) binaries, and early-universe stochastic backgrounds as the main science drivers [11].  \nOne such concept, currently in development at Northwestern University, is the Levitated Sensor Detector (LSD). With sensitivity to GW's at tens to hundreds of kHz, the LSD employs optically-trapped micron-scale disks as GW sensors. The instrument is a Michelson interferometer with two perpendicular 1-meter Fabry-P'erot arm cavities. In each arm, a disk is levitated at an antinode of a standing-wave formed by two counter-propagating beams. The trapped object behaves like a driven damped harmonic oscillator, with the corresponding trap frequency being widely-tunable with laser intensity. The periodic changes in arm-length induced by a GW manifest as a periodic shift in the position of the antinode. If the trap frequency matches the GW frequency, the levitated sensor is resonantly  \nAs a resonant detector, the LSD is well-suited to search for continuous monochromatic signals. A popular scenario involves the interaction between spinning black holes and 'ultralight' bosonic fields - i.e . those with masses several orders-of-magnitude smaller than an electron-volt (eV). Such fields can extract rotational energy from spinning BH's via 'superradiant amplification' of certain bound-states [15]. The result is a macroscopic cloud of bosons all living in the same state - commonly known as a 'gravitational atom' or 'axion cloud' [16]. These oscillating non-axisymmetric clouds generate continuous monochromatic GW's at a frequency primarily determined by the boson's mass. Tens to hundreds of kHz corresponds to µ = O (10 -11 -10 -10 ) eV.  \nThis scenario can be realized with physics beyond-theStandard-Model (BSM). For example, a large number of ultralight fields may occur as a result of the compactification of extra dimensions [17]. One of these may be the QCD axion - the pseudoscalar boson proposed to solve the strong-CP problem [18][19][20]. The axion is a Goldstone boson of a spontaneously-broken global symmetry which acquires a small mass through non-perturbative effects. Its mass, µ , is determined by the energy scale f a associated with the broken symmetry [17],  \nµ ≈ 6 × 10 -10 eV ( 10 16 GeV f a ) (1)  \nwhere 10 16 GeV ≡ Λ GUT is the grand unification (GUT) scale. An axion of mass O (10 -10 ) eV corresponds to f a being at the GUT scale. However, as we will see in Sec . VIA, signals in the LSD band are only expected up to ≈ 32 kHz, corresponding to a 6 . 6 × 10 -11 eV boson.  \nAt boson masses O (10 -11 ) eV, superradiance occurs optimally for BH's with masses between 0 . 1 and a few solar masses. Sub-solar BH's may exist as PBH's [21], and BH's in the 1 -4 M ⊙ range might be formed dynamically in binary neutron star mergers [22], accretion-induced collapsing neutron stars [23], or supernovae with unusually high fallback [24] [25].  \nThe 1 -4 M ⊙ range of BH masses is gradually being populated by microlensing candidates [26] and GW events such as GW190814 [7] and GW230529 [8]. Since the mass distribution for these objects is still unknown, we will limit our attention to stellar-origin BH's with masses between 5 and 20 M ⊙ , typical of BH's in X-ray binary systems.  \nAs a first step towards building the LSD search pipeline, we simulate the Galactic population of axion clouds with 5 -20 M ⊙ BH hosts. The essential data returned by these simulations are the gravitational-wave frequency & dimensionless strain amplitude emitted  \nby each cloud. Together with the LSD's projected sensitivity curve, we estimate the number of resolvable signals, i.e . those whose signal-to-noise ratio (SNR) rises above a given threshold after a coherent observation time T coh = 10 7 s (a little less than four months). We adopt the idealization of a 'freely-floating' detector orbiting the Milky Way at the same radius as the Solar System, but not situated on a rotating planet orbiting a star. In doing so, we neglect the amplitude and frequency modulations induced by the Earth's sidereal rotation and by its orbital motion in the Solar System. Our results establish a baseline from which a more in-depth analysis, including the aforementioned modulations, can be undertaken in future work.  \nIn Sections II & III, we introduce the essential physics of axion clouds and their GW emission. To simulate the population of axion clouds, we require a model of the stellar-origin BH population. The parameters of a black hole - mass, spin, age, and location in the Milky Way - are taken to be independent random variables, and we discuss their distributions in Section IV. The procedure for determining whether a BH of given mass, spin, and age presently hosts an axion cloud is described in Section V. The simulated cloud populations, and the corresponding ensembles of GW signals, are discussed in Section VI. Section VII provides a summary of the results, as well as tasks for future work. Throughout the paper, we adopt the metric signature ( -, + , + , +), and we retain all factors of G , c , and ℏ . We hope our decision not to set physical constants to unity will make this work more accessible to those unaccustomed to the conventions of fundamental physics theory.", 'II. SUPERRADIANT BOUND-STATES': "The creation of macroscopic clouds around spinning black holes can occur for any massive bosonic field. The simplest scenario, and the one we adopt, is that of an electrically-neutral massive scalar field freely propagating in the Kerr spacetime; We denote the BH mass and dimensionless spin by M and χ ≡ Jc/ ( GM 2 ), respectively ( J is the BH angular momentum). We also assume no self-interactions to avoid complications such as the bosenova instability [17]. The scalar field then obeys the Klein-Gordon equation [15],  \n[ g µν ∇ µ ∇ ν -m 2 ∗ ] Φ( ⃗x, t ) = 0 (2)  \nwhere the constant m ∗ has dimensions of inverse length; In the quantum theory of a scalar field, the physical meaning of m ∗ is 1 /λ c , where λ c ≡ ℏ / ( mc ) = ℏ c/µ is the reduced Compton wavelength of the boson, m is the mass of the particle, and µ = mc 2 .  \nIn Boyer-Lindquist coordinates, and with the ansatz  \nΦ( ⃗x, t ) = e -iωt e imϕ S ( θ ) R ( r ) (3)  \nthe Klein-Gordon equation separates into two ordinary differential equations (ODE's) for R ( r ) and S ( θ ). We seek a bound-state solution which is 'in-going' at the event horizon - i.e. a solution which goes to zero at infinity and looks like an in-going-wave at the horizon. The in-going boundary condition causes the eigenfrequency ω to be complex,  \nω = ω R + iω I (4)  \nwith the consequence that bound-states must either grow or decay:  \ne -iωt = e -i ( ω R + iω I ) t = e -iω R t e ω I t  \n= ⇒ Φ( ⃗x, t ) = e ω I t [ e -iω R t e imϕ S ( θ ) R ( r ) ]  \nFor ω I > 0, the field amplitude grows exponentially. A necessary and sufficient condition for the growth of a bound-state with azimuthal number m is that the event horizon's angular speed Ω H (times m ) be faster than the oscillation of the field [15],  \nω R < m Ω H (5)  \nThis requirement is called the 'superradiance condition'. As the field amplitude grows, the BH loses rotational energy, and Ω H decreases until the inequality becomes an equality. At that point, the superradiant growth ceases, and the resultant bound-state slowly dissipates by emitting GW's.  \nIt is conventional to define a dimensionless 'coupling parameter' α as the ratio of the BH's gravitational radius r g to the reduced Compton wavelength λ c of the scalar field:  \nα := r g λ c = GM c 2 µ ℏ c = GMµ ℏ c 3 (6)  \nThe 'weak-coupling' limit, defined by α ≪ 1, corresponds to the Compton wavelength of the boson being much larger than the characteristic size r g of the BH. In this limit, the bound-state energy, given by the real part of ω , can be written in closed-form [27]:  \nℏ ω R = µ [ 1 -α 2 2 n 2 + ( 2 l -3 n +1 l +1 / 2 -1 8 ) α 4 n 4 + 2 mχα 5 n 3 l ( l +1 / 2)( l +1) + ... ] (7)  \nThe small 'fine-structure' corrections beyond the leading 1 depend on the angular momentum of the cloud and the spin of the BH. The quantity in large square brackets depends on the BH & boson masses only through their  \ndimensionless product α . This motivates the introduction of a dimensionless eigenfrequency ξ = ξ R + iξ I :  \nξ R ≡ ℏ ω R µ ξ I ≡ ℏ ω I µ (8)  \nOnce we have computed ξ over a sufficiently large region of the ( α, χ ) parameter space for all bound-states { n, l, m } of interest, we can freely plug-in any BH masses and axion masses of our choosing. For example, taking M = 10 7 M ⊙ and µ = 10 -17 eV, we get α = 0 . 748. The same value is obtained taking M = 10 M ⊙ and µ = 10 -11 eV. The essential consequence is that, for a given BH spin χ , the same set of superradiant bound-states exists for both scenarios.  \nFrom a practical point of view, this also means the superradiance condition  \nξ R < mχ 2 α [ 1 + √ 1 -χ 2 ] ≡ ξ crit (9)  \nbecomes a tool for rapidly determining, for a given parameter set { µ, M, χ } , which states are superradiant.  \nAdditionally, far from the BH where relativistic effects are negligible, the radial equation reduces to ( r measured in units of λ c ) [16]  \n[ -1 2 r 2 d dr ( r 2 d dr ) -α r + l ( l +1) 2 r 2 + 1 -ξ 2 2 ] R ( r ) = 0 (10)  \nwhich is the radial Schrodinger equation for a nonrelativistic particle in a Coulomb potential V ( r ) = α/r - hence the moniker 'gravitational atom'. The α 2 / (2 n 2 ) term in Eq . 7 is precisely the 'hydrogen atom' solution to Eq . 10. The complete smallα solutions to the full KleinGordon equation have been computed order-by-order using the method of matched asymptotic expansions [16],  \nξ R = 1 -α 2 2 n 2 + [ 2 l -3 n +1 l +1 / 2 -1 8 ] α 4 n 4 + 2 mχα 5 n 3 l ( l +1 / 2)( l +1) + ... (11)  \nξ I = 2 ( 1 + √ 1 -χ 2 ) [ mχ 2 α ( 1 + √ 1 -χ 2 ) -ξ R ] α 4 l +5 × 2 4 l +1 ( n + l )! n 2 l +4 ( n -l -1)! ( l ! (2 l )!(2 l +1) ) 2 × l ∏ j =1 [ j 2 (1 -χ 2 ) + [ mχ -2 αξ R ( 1 + √ 1 -χ 2 )] ] (12)  \nIn general, ξ is a function of { n, l, m, α, χ } . For fixed α and χ , the superradiance rate is largest when  \nFIG. 1. The real part ξ R of the n = 5 bound-state, for BH spin χ = 0 . 995, plotted up to the associated maximum superradiant value of the coupling parameter, α max = 2 . 06.  \n<!-- image -->  \nFIG. 2. The ξ R curve in Fig . 1 corresponds to a BH spin of χ = 0 . 995. As mentioned in the main text, there are small spin-dependent corrections to ξ R . Had we plotted several curves in Fig . 1, each corresponding to a different spin, they would lie so close together as to be almost indistinguishable. To spread them out, we plot the fractional deviations between the curve in Fig . 1 and the corresponding curves for a variety of spins. As indicated on the y-axis, the fractional deviations are negative: For a given α , lower BH spin implies lower ξ R ; In turn, this implies a lower GW frequency via Eq . 27.  \n<!-- image -->  \nm = l = n -1. We consider only such bound states in our simulations of the Galactic axion cloud population, reducing the parameter list to { n, χ, α } .  \nSince our fiducial model of the Galactic BH population will assume M ≥ 5 M ⊙ , as well as boson masses O (10 -11 -10 -10 ) eV, the corresponding values of α are always greater than 1, but still of order unity. In this 'intermediate' regime, there are no closed-form solutions for ξ . As detailed in Appendix A, we must resort to the series-solution method for solving the radial Klein-Gordon equation. The coefficients of the infiniteseries ansatz obey a three-term recurrence relation whose solution is equivalent to the solution of a corresponding non-linear continued-fraction equation [15] [28].  \nDenoting the peak mass of the cloud as M c , the cloud's growth timescale is given by [29]  \nτ c ≡ τ nlm ln N = τ nlm ln ( M c c 2 µ ) (13)  \nwith N the number of bosons in the cloud, and τ nlm the reciprocal of the superradiance rate,  \nτ nlm ≡ 1 Γ nlm , Γ nlm ≡ 2 ω I (14)  \nτ nlm is the e-folding timescale, and we follow the authors of [29] in taking τ c as the time to fully grow the bound-state. The factor of two in Γ nlm occurs because  \nthe cloud's density is proportional to the 00-component of the stress-energy, ρ ∝ T 0 0 ∝ exp(2 ω I t ).  \nAs the cloud grows, the BH gradually loses mass and angular momentum. The growth timescales are long enough to permit an adiabatic treatment of the BH's evolution [30]. The metric can be thought of as Kerr with slowly changing M and χ . Denoting the initial BH parameters as ( M i , χ i ), the cloud's mass is  \nM c ≡ M i -M f (15)  \nand the hole's final mass & spin ( M f , χ f ) are given by [30][31]  \nM f = M i [ m 3 -√ m 6 -16 m 2 ξ 2 R α 2 i ( m -ξ R α i χ i ) 2 8 ξ 2 R α 2 i ( m -ξ R α i χ i ) ] (16)  \nχ f = ( M i M f ) 2 ( χ i -mM c ξ R α i M i ) (17)  \nSince our simulation of the Galactic axion cloud population requires us to follow the evolution of each BH-cloud system, - of which there could be millions - we save computation time by relying on these expressions for the final BH parameters.  \nThe final mass & spin become the new parameters ( M f → M i , χ f → χ i ) for determining which bound-state will grow after the present cloud has dissipated. For our simulations, the superradiance condition (Eq . 9) is used to determine, from the set { 1 , 2 , 3 , ... } , the smallest value of m for which superradiance occurs. The final state of the BH-boson system at the cessation of cloud growth is determined by Eqs . 15, 16, and 17.", 'III. GRAVITATIONAL WAVES FROM AXION CLOUDS': "At a particle-physics level, GW production by axion clouds can be understood in terms of two processes: annihilation of two bosons to a single graviton (with the BH absorbing the recoil momentum), and downwardtransitions between bound-states [32]. However, just as superradiance is a purely classical kinematic effect, the GW emission can also be understood classically in terms of the cloud's time-dependent quadrupole moment. That being said, the GW signals considered in this work correspond to the annihilation channel.  \nSince our simulation of the Galactic axion cloud population requires us to compute the GW amplitude for each cloud, - of which there could be millions we save computation time by relying on semi-analytic formulas for the amplitudes [31] [32]. Following [30], the  \nGW signal seen by a detector with perpendicular arms takes the general form  \nh ( t ) = F + ( t ) a + cos[Φ( t )] + F × ( t ) a × sin[Φ( t )] (18)  \nwhere F + ( t ) and F × ( t ) are the detector's angular pattern functions, and the phase Φ( t ) encodes the frequency evolution in the detector frame,  \nΦ( t ) = Φ 0 +2 π ∫ T T 0 f ( T ' ) dT ' = Φ 0 +2 π [ f 0 ( T -T 0 ) + 1 2 ˙ f 0 ( T -T 0 ) 2 + . . . ] (19)  \nwhere ˙ f is the frequency derivative, quantities with a zero-subscript are their values at a reference time T 0 , and T ( t ) is the time at the Solar System barycenter, related to the time at the detector t by the addition of the Rømer, Shapiro, and Einstein delays.  \nThe amplitudes a + / × are expanded in terms of spheroidal harmonics with spin-weight s = -2,  \na + / × = -∑ ˜ l ≥ 2 l h ( ˜ l ) 0 [ -2 S ˜ l, ˜ m, ˜ ω ± -2 S ˜ l, -˜ m, -˜ ω ] (20)  \nwhere ˜ ω = 2 ω R is the GW angular frequency, the parameters ( l, m ) refer to the scalar bound-state, and ( ˜ l, ˜ m ) refer to the GW modes, with ˜ l ≥ 2 l and ˜ m = 2 m . For each mode, there is a polarization-independent characteristic amplitude h ( ˜ l ) 0 [30]:  \nh ( ˜ l ) 0 = c 4 G M c M f 1 2 π 2 M f f 2 d A ˜ l ˜ m ( α i , χ i ) (21)  \nwhere f is the GW frequency, d is the source distance, and the A ˜ l ˜ m ( α, χ ) are dimensionless numerical factors which measure how much energy is carried by each mode. The corresponding luminosity in each mode is given by  \n˙ E GW ( ˜ l, ˜ m, ˜ ω ) = c 5 4 πG ( c 3 GM f ˜ ω ) 2 M 2 c M 2 f A 2 ˜ l ˜ m ( α i , χ i ) (22)  \nIn principle, the coefficients A ˜ l ˜ m must be computed numerically by solving the Teukolsky equation governing linear perturbations of the Kerr metric. The authors of [31] express ˙ E GW in the form  \n˙ E GW = c 5 G M 2 c M 2 f d ˜ E dt (23)  \nand invoke an analytic solution for d ˜ E/dt which is formally valid for α ≪ l , and which remains a good approx-  \nion up to α ∼ l [32]:  \nd ˜ E dt = 16 l +1 l (2 l -1)Γ 2 (2 l -1)Γ 2 ( n + l +1) α 4 l +10 f n 4 l +8 ( l +1)Γ 4 ( l +1)Γ(4 l +3)Γ 2 ( n -l ) (24)  \nwhere Γ is the gamma function, and α f denotes the value of α corresponding to the final mass of the BH (i.e . after the cloud has finished growing),  \nα f = α i M f M i (25)  \nComparing Eqs . 22 and 23, we see that A ˜ l ˜ m ∝ √ d ˜ E/dt , allowing us to express h ( ˜ l ) 0 directly in terms of d ˜ E/dt . Restricting ourselves to the dominant mode ˜ m = ˜ l = 2 l , we obtain a closed-form solution for the characteristic amplitude, which we use without abandon to compute the GW amplitudes of the axion clouds resulting from our simulations (we will drop the superscript (2 l ) henceforth),  \nh (2 l ) 0 ( d ) = GM c c 2 d 2 √ πM i ξ R α i M f √ d ˜ E dt (26)  \nThe corresponding GW frequency is given by  \nf = ˜ ω 2 π = 1 2 π 2 µ ℏ ξ R ≡ f 0 ξ R f = f 0 ξ R (27)  \nwhere we've introduced the zeroth-order frequency f 0 = ω 0 / 2 π , ω 0 ≡ 2 µ/ ℏ .  \nIt is often remarked that the GW frequency is proportional to twice the axion mass, f ∝ 2 µ . We see that this is, indeed, true in the smallα limit by noting that ξ R → 1 as α → 0 (Fig . 1, Eq . 11). The frequency monotonically decreases with increasing α , and for axion clouds in the kHz band, with stellar-mass BH hosts (where α is generically greater than 1), GW frequencies can be upwards of 10% smaller than the nominal value f 0 .  \nEq . 27 gives the frequency as measured in the rest-frame of the axion cloud. For an observer located elsewhere in the Milky Way, the measured signal is Doppler-shifted due to the differential rotation of the Galaxy. We assume all bodies in the Galaxy move in the azimuthal direction, ⃗v = v ϕ ˆ ϕ , and we assume the following Galactic rotation curve [33] ( r , in kpc, is the  \ncylindrical radial distance from the Galactic center):  \nv ϕ ( r ) (km/s) =          265 -1875( r -0 . 2) 2 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 (28)  \nDenoting the source-frame frequency as f s , the nonrelativistic Doppler-shifted frequency we observe is  \nf obs = ( 1 -v r c ) f s (29)  \nwhere v r is the line-of-sight component of the relative velocity between source and observer. v r is defined to be positive when the source and observer are moving away from each other.  \nWhen a cloud finishes growing, it emits GW's whose initial amplitude h 0 is given by Eq . 21. As the cloud dissipates, the amplitude decreases as [29]  \nh ( t ) = h 0 1 + t/τ GW (30)  \nwhere τ GW is the time for h to drop to half its initial value.", 'IV. THE GALACTIC POPULATION OF ISOLATED STELLAR-ORIGIN BLACK HOLES': "With the results of the previous sections in hand, we can follow the 'superradiance history' of any given BH - i.e. we can determine the sequence of scalar field bound-states, their growth & dissipation timescales, the BH mass and spin decrements, and, above all, the GW frequency & amplitude of each successive cloud. To simulate the entire Galactic population of axion clouds, we must assign each BH a mass, spin, age, and location - taken to be independent random variables in accordance with known or assumed distributions.  \nOur knowledge of the stellar-origin BH mass distribution is informed by mass measurements in X-ray binary systems [34] [35] [36], microlensing events [37], and astrometry [38], as well as through modelling of the complex physics of core-collapse supernovae [39]. Known BH's typically have masses between 5 M ⊙ and 20 M ⊙ , and power-law models are favored when fitting the mass function of low-mass X-ray binaries [34]. Not coincidentally, the massive stars which produce BH remnants are also characterized by a power-law distribution, ψ ( M ) dM ∝ M -2 . 35 dM - the 'Salpeter' function. We will assume M BH to be Salpeter-distributed on the interval 5 -20 M ⊙ .  \nBH spins have been measured in several X-ray  \nbinaries [40], but none have been measured for isolated BH's. In the case of binaries, the distribution of spin magnitudes is more-or-less uniform, so we take the BH spin to be uniformly distributed, χ ∼ U [0 , 1].  \nThe stellar content of the Milky Way can be divided into three primary regions - the thin disk, the thick disk, and the central bulge. The age distribution of stellar-origin BH's is tied to the star formation history in each region. As the Milky Way's star formation history is a topic of ongoing research, we take an agnostic approach by assigning each BH an age of 10 x yr, with x uniformly distributed on an interval which varies between the three Galactic regions. For the thin disk and thick disk, we take x ∼ U [3 , log 10 (8 × 10 9 )] and x ∼ U [3 , 10], respectively [41]. For the bulge, we assign each BH an age 10 x yr, with x ∼ U [9 , log 10 (13 × 10 9 )] [42].  \nWe assume black holes are distributed in space according to the mass profiles of the disks & bulge described in Ref . [43]. Both disks have the same axisymmetric form, with the corresponding scale lengths, scale heights, and surface densities quoted in Table I:  \nρ disk ( r, z, ϕ ) = Σ d, 0 2 z d e -| z | /z d e -r/R d (31)  \nThe bulge is also axisymmetric, with the corresponding parameters also given in Table I:  \nρ b = ρ b,0 ( 1 + r ' r 0 ) α e -( r ' /r cut ) 2 (32)  \nr ' ≡ √ r 2 +( z/q ) 2  \nTABLE I. Physical parameters for the empirical stellar-mass distribution of the Milky Way inferred in Ref . [43]. The values for those with an asterisk (*) are the means marginalized over all other parameters in the model. The disk scale heights are the best-fitting values from Ref . [44], and the stellar bulge model is an axisymmetric modification of the result from Ref . [45] in which the assumption of a constant mass-to-light ratio in the bulge was used to convert photometric data into a mass model.  \nWe apportion the BH's among the three Galactic regions according to the fractions f thin , f thick , and f bulge , defined by f i = M i / ∑ i M i , i ∈ { thin, thick, bulge } . The disk masses are obtained by integrating ρ disk , with the radial integral cut-off at 25 kpc, and the vertical integral cut-off at 3 scale heights. This gives 3 . 97 × 10 10 M ⊙ and 1 . 5 × 10 10 M ⊙ for the thin and thick disks, respectively. We take the bulge mass to be 8 . 9 × 10 9 M ⊙ , the value quoted in [43]. The corresponding f i are 62%, 24%, 14%, respectively. We will assume the Galactic population of N BH BH's to be apportioned likewise: 62% in the thin disk, 24% in the thick disk, and 14% in the bulge.", 'V. SIMULATION PROCEDURE': "The simulation is a procedure by which, for a given axion mass, and from an initial population of N BH BH's sprinkled throughout the Milky Way, we determine the number N c of extant axion clouds. Each simulation outputs the physical properties, distances, and the GW frequencies & amplitudes of the N c clouds.  \nAt the outset, each BH is assigned a mass, spin, and age. We will illustrate the procedure with an example, and then summarize the procedure with a flowchart: Taking µ = 4 × 10 -11 eV, consider the evolution of a 5 M ⊙ , χ = 0 . 95 BH with an age of 10 8 yrs. The superradiance condition, Eq . 9, determines  \nwhich bound-state grows first.  \n         ξ R = 1 . 03 ξ crit = 0 . 24 ( m = 1) ξ R = 0 . 76 ξ crit = 0 . 48 ( m = 2) ξ R = 0 . 89 ξ crit = 0 . 73 ( m = 3) ξ R = 0 . 94 ξ crit = 0 . 97 ( m = 4) (33)  \nSince ξ R > ξ crit for m = 1 , 2 , 3, the first superradiant bound-state is n = 5 , l = m = 4, and it grows on a timescale of τ c = 3 . 4 yrs. The BH's mass and spin are decreased to 4 . 94 M ⊙ and 0 . 938, respectively. Once the cloud has finished growing, it dissipates on a timescale τ GW = 0 . 8 yrs. The time from the BH's birth to the cloud's dissipation is only τ c + τ GW = 4 . 2 yrs, leaving plenty of time for new clouds to develop. We denote by t r the time remaining to the present. In this case, t r = 10 8 -4 . 2 ≈ 10 8 yrs.  \nThe next bound-state is n = 6 with τ c = 3536 yrs. The BH's mass and spin are decreased to 4 . 67 M ⊙ and 0 . 84, respectively. Once the cloud has finished growing, it dissipates on a timescale τ GW = 2445 yrs. At this point, t r = 9 . 9994 × 10 7 yrs - still plenty of time left for further superradiance.  \nThe next (and final) bound-state is n = 7 with τ c = 6 × 10 7 yrs. The BH's mass and spin are decreased to 4 . 5 M ⊙ and 0 . 74, respectively. Once the cloud has finished growing, t r = 3 . 9 × 10 7 yrs remain. The dissipation timescale τ GW = 7 × 10 7 yrs. Since τ GW > t r , the n = 7 cloud is still present today. It has an initial mass M c = 0 . 16 M ⊙ , and it radiates at f = 18 . 9 kHz. Placing the source at d = 1 kpc (for example), the initial strain amplitude h 0 = 10 -26 (Eq . 26). The signal observed today was emitted d/c = 3300 yrs ago, so the corresponding amplitude h ( t ) = 6 . 9 × 10 -27 (Eq . 30 with t = t r -d/c ).  \nOur simulation of the Galactic cloud population consists of applying the foregoing procedure to each of the BH's in the galaxy. If a given BH only permits a bound-state whose growth timescale is greater than the age of the universe ( τ c > τ uni = 1 . 38 × 10 10 yr), the host BH is removed from the simulation.  \nOur criterion for whether a given cloud is still present today is τ GW > t r . For each black hole, there are only two final options: Either a cloud has finished growing and is still present today, or a cloud is growing on a timescale greater than the age of the universe.  \nThose BH's with an extant cloud are assigned a location in the Milky Way (Eqs . 31 and 32). Earth is assigned to an arbitrary, but fixed, point on the circle of radius 8 . 3 kpc in the Galactic midplane. For a cloud located at distance d , we check the inequality ct r > d to determine if there has been enough time for GW's to propagate to Earth since the cloud formed. Those  \nclouds for which d > ct r are presently unobservable, and we retain only those clouds for which ct r > d . We summarize this section with the following flowchart:  \nFor a given µ , M BH , χ , and BH age, find the lowest superradiant value of n .  \n→ If τ c > τ uni , the BH is removed from the simulation.  \n→ Otherwise, the dissipation timescale τ GW determines whether a new cloud will start growing in accordance with τ GW > t r (cloud still present) or τ GW < t r (cloud has dissipated, and a new cloud begins growing).  \n→ Repeat the previous steps until one of two possibilities is obtained: a.) A cloud is growing with τ c > the age of the universe, or b.) A cloud is still present & radiating GW's today.  \n→ If the cloud hasn't dissipated yet, assign it a random position, and compute the GW strain at Earth's location only if the travel-time inequality ct r > d is true.", 'A. Cloud populations': "The total number of stellar-origin black holes has been estimated to be O (10 8 ) from the Milky Way's supernova rate of O (1) century -1 [46], and from populationsynthesis estimates [47]. We take N BH = 10 8 , bearing in mind that the true number could be larger by a factor of a few, or even another order-of-magnitude [48]. We have simulated the axion cloud population for µ = (3 , 3 . 5 , 4 , 4 . 5 , 5 , 5 . 5 , 6 , 6 . 5) × 10 -11 eV.  \nThe output of a simulation is a collection of all extant BH-cloud systems in the Milky Way. Those BH's which have experienced the growth of a single cloud are described by a list comprising the BH age, the initial and final values of the BH mass & spin, the bound-state { n, l, m, ξ R , ξ I } , the cloud's properties mass M c , growth timescale τ c , and dissipation timescale τ GW , - the source distance d , and the GW frequency and amplitude ( f, h ) . BH's which have experienced the growth of multiple bound-states are each characterized by a set of such lists, one per bound-state. The GW frequency and amplitude are only computed for the extant cloud, all previous bound-states having already dissipated.  \nFor a given axion mass, the number of extant clouds is a random variable whose mean and standard deviation are estimated by performing twenty simulations with 5 × 10 6 BH's per simulation, computing the  \nsample mean & sample standard deviation of N c over the 20 trials, and then multiplying them by 20 and √ 20, respectively.  \nAn ensemble of GW signals from axion clouds is a scatter plot in the h vs. f plane, as in Figs . 3, 4, and 5 . The distribution of amplitudes and frequencies is not random, but consists of well-defined bands corresponding to the various occupied boundstates. The lowest bound-state resulting from our simulations is n = 6, reflecting the general difficulty for stellar-mass BH's to produce clouds in the LSD band.  \nAlso reflecting this difficulty is the rapid decline in the number of clouds N c with increasing boson mass µ (Fig . 6) . For µ = 3 × 10 -11 eV, N c = (9 . 323 ± 0 . 007) × 10 5 , while at µ = 6 . 5 × 10 -11 eV, the number has dropped to 130 ± 10. N c goes to zero around 6 . 6 × 10 -11 eV, corresponding to a nominal upper limit of ≈ 32 kHz for signals expected in the LSD band. Higher-frequency signals could occur from BH's with M BH < 5 M ⊙ , especially in light of the recent discoveries of lower-mass-gap objects.  \nIn all cases, the distribution of GW frequencies occurs below the nominal value f 0 ∝ 2 µ due to the positive scaling of gravitational redshift with BH mass. This interpretation is confirmed by a plot of the sourceframe GW frequencies f GW,s and initial BH masses M i for all extant clouds in a given simulation (Fig . 7, with µ = 3 × 10 -11 eV). For each scalar bound-state, there is a tight relationship, with more massive BH's producing lower-frequency clouds.  \nIn the introduction (Sec . I), we noted a potential connection between the QCD axion and the GUT scale Λ GUT (Eq . 1): An axion of mass O (10 -10 ) eV corresponds to f a ≈ Λ GUT . If the solution to the strong-CP problem is tied to GUT phenomenology, then discovery of an O (10 -10 ) eV axion would be an exciting, albeit indirect, form of evidence for grand unification. The number of clouds in the Milky Way dropping to zero around 6 . 6 × 10 -11 eV would seem to preclude the possibility of detecting an O (10 -10 ) eV axion - and, by extent, of probing GUT-scale physics with the LSD. Lower-mass-gap BH's could produce clouds at higher µ , thereby reviving hopes of finding a GUT-scale axion. Another possibility is that Λ GUT is model-dependent, giving rise to a range of possible values including 10 17 GeV, which corresponds to O (10 -11 ) eV bosons.", 'B. Resolvable signals': "The standard result for coherent detection of a continuous monochromatic signal, h ( t ) = h 0 cos( ωt ), is that the signal-to-noise ratio (SNR) ρ grows as the square-root of  \nthe coherent integration time T coh [11]  \nρ = h 0 √ T coh √ S n ( f ) (34)  \nwhere √ S n ( f ) is the one-sided amplitude spectral density (ASD) of the detector noise (the 'sensitivity curve') evaluated at the GW frequency, and the trapping frequency of the levitated sensor is constant during the entire observation time. The LSD is an Earth-bound detector for which the observed signal, Eq . 18, experiences both amplitude modulation via the angular-dependence of the sensitivity, and phase modulation via the Earth's daily (diurnal) rotation and orbital motion. These corrections require knowledge of both the Earth's position and the source position to high accuracy. In an all-sky blind search for axion clouds, the source position is not known ahead of time, so a realistic search will require a large number of templates corresponding to many sky locations. A coherent search for O (10) kHz GW's over the full observation time ('fully-coherent search') of four months is not computationally-feasible because of the need to take fine steps in parameter space. A 'semi-coherent' search, in which the observation time is divided into several segments, sacrifices sensitivity for a great reduction in the number of templates needed to perform a coherent search on each segment [49].  \nAs a means of setting upper limits on the expected number of resolvable signals, we compute the SNR for the idealized case of a detector freely orbiting the Milky Way at the same radius as the Solar System (i.e . not attached to a planet or star system). There is a positive frequency derivative due to the gradual dissipation of the cloud which, however, is estimated to be too small to be detected [29]. The signal remains perfectly coherent over the full observation time, and a hypothetical search performed at the Solar System barycenter would simply involve Fourier-transforming the data and looking for lines in the power spectrum.  \nTaking T coh = 10 7 s, and with the projected sensitivity curves for the current 1-m LSD prototype, as well as for future 10-m and 100-m versions [13], we compute the corresponding SNR's for all sources in the galaxy. We count those with ρ > ρ t as resolvable, and we adopt the threshold ρ t = 10 (Fig . 8).  \nThe 'loudness' of a signal is determined primarily by the source distance. The distance, in turn, is a random variable determined by the randomly-assigned position vector (Eqs . 31 and 32) of the source. Thus, for a given set of extant clouds, the number of individuallyresolved sources N res will vary each time we re-assign their position vectors. We estimate the mean & standard deviation of N res for a given population of extant clouds by laying them down in the Galaxy N reshuffle = 100 times and counting how many are resolvable in each  \nFIG. 3. The population of axion clouds, with µ = 3 × 10 -11 eV. From a population of 10 8 BH's, the simulation returns N c = (9 . 32 ± 0 . 03) × 10 5 extant clouds. Of these, 8 . 3 × 10 5 (89%) satisfy the GW travel-time condition (Sec . V). Since GW amplitudes & frequencies are, necessarily, only computed for those clouds satisfying the travel-time condition, it should be understood that only those clouds are represented in this figure, as well as in Figs . 4 and 5.  \n<!-- image -->  \nFIG. 4. The population of axion clouds, with µ = 4 × 10 -11 eV. From a population of 10 8 BH's, the simulation returns N c = (1 . 17 ± 0 . 01) × 10 5 extant clouds. Of these, 9 . 5 × 10 4 (81%) satisfy the GW travel-time condition.  \n<!-- image -->  \nFIG. 5. The population of axion clouds, with µ = 6 × 10 -11 eV. From a population of 10 8 BH's, the simulation returns N c = 900 ± 200 extant clouds. Of these, 620 (69%) satisfy the GW travel-time condition. The apparent gap in the n = 9 band at 28 kHz is an artifact due to the small number of clouds.  \n<!-- image -->  \nFIG. 6. The number of extant axion clouds N c drops rapidly with increasing boson mass µ , and it extrapolates to zero around 6 . 6 × 10 -11 eV, corresponding to f ≈ 32 kHz. In light of the lower-mass-gap objects found by LIGO-Virgo, we note that superradiant instabilities with a O (10 -11 ) eV boson are stronger for M BH < 5 M ⊙ than for M BH > 5 M ⊙ . Since the BH's involved in our simulations are of the latter type, the true number of clouds, for any boson mass, could be greater than our estimate by a factor depending on the mass distribution & total number of lower-mass-gap BH's in the Milky Way.  \n<!-- image -->  \n're-shuffling'. The mean & standard deviation are then computed as  \n¯ N res = 1 N reshuffle N reshuffle ∑ i =1 N res ,i (35)  \n¯ σ res = 1 N reshuffle -1 N reshuffle ∑ i =1 ( N res ,i -¯ N res ) 2 (36)  \nWith a 100-m detector, assuming µ = 3 × 10 -11 eV, ¯ N res = 600 with ¯ σ res = 20. In the most pessimistic case ( µ = 5 . 5 × 10 -11 eV), there are only O (1) resolvable signals, and we have not estimated the associated uncertainty. The 10 -26 kHz range is where we expect resolvable signals to be present for a 100-m LSD. For a 10-m instrument, O (1) resolvable signals appear at µ = 3 × 10 -11 eV, while a 1-m instrument does not have the required sensitivity to detect individual sources.  \nIn the event a continuous monochromatic signal  \nFIG. 7. The source-frame GW frequencies f GW,s and initial BH masses M i for the extant clouds in the 3 × 10 -11 eV simulation. The distribution of frequencies below the nominal value f 0 ∝ 2 µ is a result of gravitational redshift.  \n<!-- image -->  \nis detected by the LSD, we will have to answer the question: Is this signal from an 'axion cloud' - a superradiant bound-state of a scalar (spin-0) field - or from a cloud involving a spin-1 ('Proca') field? In general, Proca fields give rise to stronger GW signals than scalar fields [50]. As a result, we would expect resolvable signals from Proca clouds to be found at greater distances than those from scalar clouds. For the 100-m detector, with µ = 3 × 10 -11 eV, the resolvable signals are depicted in terms of their SNR's and source distances in Fig . 9. The vast majority are less than 3 kpc away. Turning this on its head, the detection of a continuous monochromatic signal with an inferred distance significantly greater than 3 kpc could be a potential indicator of a spin-1 field.", 'C. Unresolved signals': "For all boson masses, the majority of GW signals have amplitudes less than 10 -23 , with the weakest having h = O (10 -29 ). The unresolvable signals incoherently combine to form a Galactic confusion foreground which manifests as an excess noise in the detector. As before, we neglect the diurnal and annual modulations of the background, and instead provide a preliminary estimate of the foreground's strength compared to the nominal 1-m, 10-m, and 100-m LSD sensitivity curves. In a strain-frequency plot (e.g . Fig . 3), we bin the cloud amplitudes (with bin width δf = 10 -2 f c , where f c is the center frequency of a given bin, and the factor 10 -2 is the full-width-at-half-maximum (FWHM) of the trapped object's response function around f c ), and we associate  \nan rms amplitude, defined as follows, with each bin.  \nWe start by creating a bin centered on the frequency of the cloud with the smallest GW frequency in a strain-frequency plot, e.g. Fig . 3 . All axion clouds emit monochromatic signals,  \nh i ( t ) = h 0 ,i cos ( 2 πft + ϕ i ) (37)  \nwhere the phases ϕ i are uniformly-distributed between 0 and 2 π , and i runs over all clouds in the bin. The squared sum of all signals in the bin is time-averaged over a period T c = 1 /f c , where f c is the frequency at the center of the bin; The result is a dimensionless timeaveraged power associated with that bin. The squareroot of the power represents an effective amplitude h eff of the confusion foreground in the bin,  \nh eff := √ √ √ √ 1 T c ∫ T c 0 dt [ ∑ i h i ( t ) ] 2 (38)  \nWe then create a new bin with center frequency f c ∣ ∣ new and width δf ∣ ∣ new ,  \nf c ∣ ∣ new = f c ∣ ∣ old +10 -2 δf ∣ ∣ old (39)  \nδf ∣ ∣ new = 10 -2 f c ∣ ∣ new (40)  \nand we compute h eff for this bin. The center frequency is shifted rightwards by a fraction (arbitrarily chosen to be 10 -2 ) of the previous bin width so that  \nFIG. 8. The SNR's ρ of resolvable GW signals from simulated axion cloud populations, for the 100-m LSD. As a function of µ , the number of resolvable signals drops to zero beyond 5 . 5 × 10 -11 eV, thereby constraining the expected frequency range of resolvable sources to f < 26 kHz.  \n<!-- image -->  \nFIG. 9. Scatter plot of the SNR's ρ and source distances d of resolvable signals for the 100-m LSD, taking µ = 3 × 10 -11 eV. There appears to be an upper limit d ≲ 3 kpc, and certainly a hard upper limit d < 4 kpc, suggesting a way to distinguish between spin-0 and spin-1 bosonic clouds: Proca clouds are generally louder GW sources than axion clouds, so a GW signal with d ≫ 3 kpc could be taken as evidence of a spin-1 field.  \n<!-- image -->  \nadjacent bins overlap, ensuring some degree of continuity in h eff vs. f . We continue until we reach the rightmost end of the cloud population. Each bin is then characterized by an ordered pair ( f c , h eff ) (Fig . 10).  \nA preliminary method for estimating the LSD's sensitivity to the confusion foreground is to treat each pair ( f c , h eff ) as if they were the frequency and amplitude of a hypothetical monochromatic signal whose corresponding effective SNR ρ eff , computed via Eq . 34, is then compared to a threshold ρ t . We continue to require ρ t = 10. The numerator and denominator of Eq . 34 ( h eff √ T coh and √ S n ( f c ) , respectively ) are shown separately in Fig . 11, and their ratio (the SNR) is shown  \nin Figs . 12, 13, and 14 for the 1 -, 10 -, and 100 -m instruments, respectively.  \nWe find that a single 1-m LSD does not appear to have the required sensitivity to detect the foreground for any value of µ . A 10-m detector could detect the foreground with ρ eff = O (10) if the axion mass µ ∈ (3 -4) × 10 -11 eV, while in the (4 -4 . 5) × 10 -11 eV range, only the peak of the foreground rises to the threshold, and just barely so (Fig . 13). A 100-m instrument could detect the foreground with large ρ eff if µ ∈ (3 -6) × 10 -11 eV (Fig . 14). In the range (3 -3 . 5) × 10 -11 eV, the peak value of ρ eff is O (10 3 ), and remains O (10 2 ) up to 5 × 10 -11 eV.", 'VII. CONCLUSION': "We have produced Galactic-scale populations of the hypothetical GW sources known as 'axion clouds' with the axion mass chosen to correspond to frequencies in the 10 -100 kHz band. By computing superradiant bound-states up to n = 9, we have accounted for nearly all clouds with growth timescales less than the age of the universe.  \nThe largest number of clouds occurs for the lightest boson mass capable of producing GW's at the frequencies of interest. This was to be expected, as superradiance occurs more readily for small α ∝ µM BH . For a BH of mass M ≥ 5 M ⊙ , the smallest value of α is obtained with the smallest allowed boson mass, 3 × 10 -11 eV. In this most optimistic case, the total number of extant clouds is close to 1 million.  \nThe population of axion clouds has been assumed to be spatially-distributed within the Milky Way in the same way as the stellar disks and central bulge. Statistically, some may be near enough that the continuous monochromatic signal can be detected by observing over a long enough period of time, e.g . 10 7 s, such that the SNR rises above a given threshold ρ t ; We have imposed a stringent threshold ρ t = 10, but we leave it for future work to determine the most appropriate threshold for our search pipeline. For a 100-m instrument, several hundred resolvable signals are predicted to occur if µ ≈ 3 × 10 -11 eV, but this number could be upwards of an order-of-magnitude larger if the total number of stellar-origin BH's is also larger than we have assumed (see the comment made at the beginning of Sec . VI). For a 10-m detector, only O (1) resolvable signals occur in our simulation at µ = 3 × 10 -11 eV.  \nMeanwhile, the ensemble of unresolved signals produces a confusion foreground which is estimated to be detectable with potentially large SNR by a 100-m LSD, assuming µ ∈ ( 3 -6 ) × 10 -11 eV, or by a 10-m instrument at moderate SNR, assuming µ ∈ ( 3 -4 . 5 ) × 10 -11 eV.  \nFinally, we note the following limitations of this work, as well as directions for future work: First, since isolated BH's have no EM counterpart, we do not know, ahead of time, the direction to these GW sources. Targeted & directed searches for axion clouds will, therefore, not be possible for isolated BH's, and we must resort to blind all-sky searches. Doppler modulations of the GW frequency can be accounted for by including the source's right ascension & declination as additional parameters in the signal model for Bayesian parameter estimation. To avoid the number of demodulation templates becoming prohibitively large, - we do not  \nwant the time required for data analysis to be greater than the four-month observation period - we will resort to semi-coherent techniques for constructing a detection statistic. A final complication might be that a source has a transverse proper motion large enough to change its sky location within the observation time.  \nIn this work, we have used the SNR as a baseline detection statistic. Since we have not yet developed the full semi-coherent search pipeline, we do not yet know what is the most appropriate detection statistic for continuous monochromatic signals. Another outstanding question pertains to our scan strategy: Given that we will take data at each trapping frequency for 4 months, how finely should the set of frequencies be discretized? At what frequency should we begin our searches? Figs . 8 and 11 both indicate the 14 -20 kHz range to be the most optimistic, in terms of the expected number of resolvable signals and the strength of the confusion foreground, respectively.  \nAdditionally, our treatment of the confusion foreground has not accounted for the intrinsic anisotropy of the signal: The axion clouds will be distributed throughout the disks and bulge of the Milky Way, so the strength of the foreground will vary over the sky in a complicated way. Searches for stochastic signals typically involve an 'excess-power' method, as well as cross-correlation between multiple detectors. Plans to build a second 1-m instrument at UC Davis (in addition to the Northwestern detector) are in development, so while a single 1-m detector might not have the requisite sensitivity, the prospects for a two- or multi-detector scheme are an exciting avenue of future study.", 'ACKNOWLEDGMENTS': "We would like to thank Vedant Dhruv for making public his Mathematica notebook for scalar bound-states in Kerr. We also thank Timothy Kovachy for clarifying issues of numerical precision when using Mathematica 's root-finding routines; and Richard Brito for several clarifying discussions on axion clouds & their gravitationalwave emission. JS, AG, and SL are supported by the W.M. Keck Foundation. AG, GW, and NA are supported in part by NSF grants PHY-2110524 and PHY2111544, the Heising-Simons Foundation, the John Templeton Foundation, and ONR Grant N00014-18-1-2370. NA is partially supported by the CIERA Postdoctoral Fellowship from the Center for Interdisciplinary Exploration and Research in Astrophysics at Northwestern University and the University of California-Davis. SL is also supported by EPSRC International Quantum Technologies Network Grant EP/W02683X/1 and is grateful for EPSRC support through Standard Research Studentship (DTP) EP/R51312X/1. VK is supported by a CIFAR Senior Fellowship and through Northwestern Uni-  \nFIG. 10. The effective amplitudes h eff , as defined in Eq . 38, for the binned GW amplitudes of a simulated population of axion clouds with µ = 3 × 10 -11 eV. The width of each bin is a factor 10 -2 of the central frequency f c , reflecting the FWHM of the detector response when the trap frequency is f c .  \n<!-- image -->  \nFIG. 11. The confusion-limited foreground of unresolved axion clouds, as a function of the boson mass µ , taking T coh = 10 7 s. The 1 -, 10 -, and 100 -m LSD sensitivity curves are depicted by the red, blue, and black dashed curves, respectively.  \n<!-- image -->  \nversity through the D.I. Linzer Distinguished University Professorship. A.L. is supported by the Fannie and John Hertz Foundation. This work used the Quest computing facility at Northwestern.", 'Appendix: Superradiant bound-states': "The creation of an axion cloud corresponds to an instability of the Kerr space-time due to the presence of a massive scalar field. The amplifying mechanism, 'su-  \nFIG. 12. The corresponding SNR ρ eff of the confusion-limited foreground with respect to the 1-m LSD sensitivity curve. The dashed line depicts the detection threshold ρ t = 10 adopted in this paper, and the solid line represents a lower threshold of 5. The foreground does not rise above either threshold.  \n<!-- image -->  \nFIG. 13. The corresponding SNR ρ eff of the confusion-limited foreground with respect to the 10-m LSD sensitivity curve. The dashed line depicts the detection threshold ρ t = 10 adopted in this paper, and the solid line represents a lower threshold of 5. Over an observation time of 10 7 s, the 10-m instrument has the required sensitivity to detect the foreground in the range (3 -4 . 5) × 10 -11 eV, although in the (4 -4 . 5) × 10 -11 eV range, only the peak rises above the threshold, and only barely so. If a lower threshold ρ t = 5 were adopted, the range of boson masses could be extended up to 5 × 10 -11 eV.  \n<!-- image -->  \nFIG. 14. The corresponding SNR ρ eff of the confusion-limited foreground with respect to the 100-m LSD sensitivity curve. The dashed line depicts the detection threshold ρ t = 10 adopted in this paper, and the solid line represents a lower threshold of 5. In the (3 -6) × 10 -11 eV range, the foreground is predicted to be detectable by a 100-m instrument with effective SNR's upwards of a thousand in the (3 -3 . 5) × 10 -11 eV range.  \n<!-- image -->  \nperradiance', is the Penrose process in which rotational energy is extracted by a bosonic wave rather than by a particle. In the process, the Kerr BH loses mass and angular momentum, subject to the condition that its 'irreducible mass' does not decrease.  \nIn the Penrose scenario, a particle travelling through a BH's ergoregion can split in two, one of which falls into the hole, while the other escapes to infinity. If the orbital angular momentum of the infalling particle is of opposite sign to that of the hole, the BH's loses rotational energy to the escaping particle: Energy has been extracted from the ergoregion.  \nThe story for waves runs analogously: An incident wave with amplitude I splits into a part transmitted into the BH (with amplitude T ) and a part which escapes (the reflected wave with amplitude R ). If the transmitted wave is counter-rotating, the rotational energy of the BH decreases, leading to an outgoing wave with R > I .  \nThe novelty of a massive scalar field is that its mass acts like a mirror: Unlike a massless field, a massive field can become trapped in a bound-orbit, leading to continuous extraction of rotational energy. The end result of the runaway amplification is a macroscopic scalar field bound-state - the 'axion cloud'. In an astrophysical context, rather than a wave incoming from infinity, the initial seed for superradiance can be any arbitrary quantum fluctuation in the scalar field, even if the field is in its classical ground state [30][29] . As a result, the growth of an axion cloud begins immediately after the birth of a BH.  \nAn axion cloud's binding energy (which determines the GW frequency) and growth timescale depend on the dynamics of the scalar field. For the scenario we have adopted, the field obeys the Klein-Gordon equation on the Kerr space-time. The Kerr metric describes an axisymmetric, neutral, and rotating black hole:  \nds 2 = -[ 1 -2 GMr c 2 ρ 2 ] c 2 dt 2 -4 GMar sin 2 θ c 2 ρ 2 cdtdϕ + ρ 2 ∆ dr 2 + ρ 2 dθ 2 + [ r 2 + a 2 + 2 GMa 2 r sin 2 θ c 2 ρ 2 ] sin 2 θdϕ 2 (A.1)  \nwhere M is the BH mass, J is the BH angular momentum, ρ 2 ≡ r 2 + a 2 cos 2 θ , a ≡ J/ ( Mc ) is the Kerr parameter, and ∆ ≡ r 2 -2 r g r + a 2 , where we have defined the gravitational radius r g ≡ GM/c 2 . In terms of the dimensionless Kerr parameter, χ ≡ a/r g = Jc/ ( GM 2 ), the inner and outer horizons - the two roots of ∆ = ( r -r + )( r -r -) - are  \nr ± = r g [ 1 ± √ 1 -χ 2 ] (A.2)  \nIt follows that χ is restricted to the interval  \n0 < χ < 1 (A.3)  \nThe event horizon is located at r = r + , and the angular velocity of the horizon is  \nΩ H = cχ 2 r + (A.4)  \nThe scalar field obeys the Klein-Gordon equation,  \n[ ∇ µ ∇ µ -m 2 ∗ ] Φ( ⃗x, t ) = 0 (A.5)  \nwhere ∇ µ is the covariant derivative with respect to the Kerr metric, and, as mentioned in the text, m ∗ has the quantum-mechanical interpretation as the reciprocal of the boson's Compton wavelength. In Boyer-Lindquist coordinates, the Klein-Gordon equation is separable via the ansatz  \nΦ( ⃗x, t ) = Re [ e -iωt e imϕ S ( θ ) R ( r ) ] (A.6)  \nInvoking the identity  \n∇ µ ∇ µ Φ = 1 √ -g ∂ µ [ √ -gg µν ∂ ν Φ ] (A.7)  \n√ -g = ρ 2 sin θ  \nthe Klein-Gordon equation separates into two 2 nd -order linear homogeneous ODE's for R ( r ) and S ( θ ):  \nD θ [ S ]+ [ χ 2 α 2 ( ξ 2 -1 ) cos 2 θ -m 2 sin 2 θ +Λ ] S ( θ ) = 0 (A.8)  \nD r [ R ] + [ α 2 ξ 2 ( r 2 + χ 2 ) 2 -4 χmαξr + m 2 χ 2 -∆ ( α 2 r 2 + χ 2 α 2 ξ 2 +Λ )] R ( r ) = 0 (A.9)  \nD θ ≡ 1 sin θ d dθ [ sin θ d dθ ] , D r ≡ ∆ d dr [ ∆ d dr ] (A.10)  \nWe have expressed the decoupled equations in terms of the dimensionless variables ( χ, α and ξ ) used in the main text . The radial coordinate in A.9 is measured in units of r g .  \nBound-state solutions must go to zero at infinity and be in-going at the event horizon . The in-going condition means that R ( r ) ∝ e -ikr ∗ as r ∗ → -∞ , with r ∗ the Kerr tortoise coordinate which maps the event  \nhorizon to -∞ ,  \ndr ∗ dr = r 2 + a 2 ∆ (A.11)  \nThis means that plane waves at the event horizon ( r ∗ → -∞ ) can only move 'to the left', i.e . into the black hole.  \nThe spectra of both bound-states and BH quasinormal modes can be found via Leaver's continuedfraction method [51] [15] . The radial function R ( r ) is represented by an infinite series,  \nR ( r ) = ( r -r + ) -iσ ( r -r -) iσ + β -1 e qr ∞ ∑ n =0 a n ( r -r + r -r -) n (A.12)  \nσ = α ( 1 + √ 1 -χ 2 )( ξ -ξ crit ) √ 1 -χ 2 (A.13)  \nq = α √ 1 -ξ 2 (A.14)  \nβ = α 2 ( 1 -2 ξ 2 ) q (A.15)  \n(The quantity we denote by β is the same as the quantity denoted by χ in Ref . [15].) With this ansatz, A.9 implies a three-term recurrence relation for the unknown coefficients a n ,  \n{ α n a n +1 + β n a n + γ n a n -1 = 0 , n = 1 , 2 , . . . a 1 = -β 0 α 0 a 0  \n(A.16)  \nwhere the coefficients α n , β n and γ n are defined by  \n     α n = n 2 +( c 0 +1) n + c 0 β n = -2 n 2 +( c 1 +2) + c 3 γ n = n 2 +( c 2 -3) n + c 4 (A.17)  \nand c 0 , c 1 , c 2 , c 3 and c 4 are given by  \nc 0 = 1 -2 iαξ -2 i ( αξ -mχ 2 ) √ 1 -χ 2 (A.18)  \nc 1 = -4 + 4 i [ αξ -iα √ 1 -ξ 2 ( 1 + √ 1 -χ 2 )] + 4 i ( αξ -mχ 2 ) √ 1 -χ 2 -2 [ α 2 ξ 2 + α 2 ( 1 -ξ 2 ) ] α √ 1 -ξ 2 (A.19)  \nc 2 = 3 -2 iαξ -2 [ α 2 ( 1 -ξ 2 ) -α 2 ξ 2 ] α √ 1 -ξ 2 -2 i ( αξ -mχ 2 ) √ 1 -χ 2 (A.20)  \nc 3 = 2 i ( αξ -iα √ 1 -ξ 2 ) 3 α √ 1 -ξ 2 + χ 2 α 2 ( 1 -ξ 2 ) -Λ lm -1 + 2 √ 1 -χ 2 ( αξ -iα √ 1 -ξ 2 ) 2 +2 imχα √ 1 -ξ 2 -( αξ -iα √ 1 -ξ 2 ) 2 α √ 1 -ξ 2 +2 α √ 1 -ξ 2 √ 1 -χ 2 + 2 i √ 1 -χ 2 [ 1 + ( αξ -iα √ 1 -ξ 2 ) 2 α √ 1 -ξ 2 ] [ αξ -mχ 2 ] (A.21)  \nc 4 = ( αξ -iα √ 1 -ξ 2 ) 4 α 2 (1 -ξ 2 ) + 2 iξ ( αξ -iα √ 1 -ξ 2 ) 2 √ 1 -ξ 2 -2 i ( αξ -iα √ 1 -ξ 2 ) 2 ( αξ -mχ 2 ) α √ 1 -ξ 2 √ 1 -χ 2 (A.22)  \nThe series coefficients are related by an infinite continued-fraction [28]  \na n +1 a n = -γ n +1 β n +1 -α n +1 γ n +2 β n +2 -... (A.23)  \nContinued-fractions are commonly written in the slightly-less cumbersome notation  \na n +1 a n = -γ n +1 β n +1 -α n +1 γ n +2 β n +2 -α n +2 γ n +3 β n +3 -. . . (A.24)  \nSince a 1 /a 0 = -β 0 /α 0 , we obtain a condition whose roots are the desired bound-state frequencies:  \nβ 0 -α 0 γ 1 β 1 -α 1 γ 2 β 2 -α 2 γ 3 β 3 -. . . = 0 (A.25)  \nStrictly speaking, the radial and angular eigenvalues, ξ and Λ, must be found simultaneously. Leaver's method can also be applied to A.8 [52], resulting in a continued-fraction condition analogous to A.25. We then have two equations for the two unknowns.  \nConveniently, we can reduce the root-finding problem to merely solving A.25 by using the Mathematica function 'SpheroidalEigenvalue'. With the change of variable  \nz = cos θ , and in terms of the following quantities,  \nγ 2 ≡ χ 2 α 2 (1 -ξ 2 ) (A.26)  \nλ ≡ Λ -γ 2 (A.27)  \nthe angular equation A.8 takes the standard form implemented in Mathematica :  \n(1 -z 2 ) d 2 S dz 2 -2 z dS dz + [ γ 2 (1 -z 2 ) + λ -m 2 1 -z 2 ] S ( z ) = 0 (A.28)  \n'SpheroidalEigenvalue' yields λ , and 'SpheroidalPS' yields S ( z ). The continued-fraction equation A.25, with Λ replaced by γ 2 + 'SpheroidalEigenvalue', can then be solved for ξ with the Mathematica function FindRoot .  \nFor our axion cloud simulations, we have needed to compute ξ for bound-states up to, and including, n = 9. As an example, we have plotted the real & imaginary parts of the n = 8 bound-state in Figs . 15, 16 and 17; Figs . 15 and 16 are analogous to Figs . 1 and 2.  \nFIG. 15. The real part ξ R of the n = 8 bound-state, for BH spin χ = 0 . 995, plotted up to the associated maximum superradiant value of the coupling parameter, α max = 3 . 75.  \n<!-- image -->  \nFIG. 16. The fractional deviations between Fig . 15 and the ξ R curves for various other spins.  \n<!-- image -->  \nFIG. 17. The imaginary part of the n = 8 bound-state eigenfrequency.  \n<!-- image -->  \nFIG. 18. The radial profiles ∣ ∣ Φ( r ; t = 0 , θ = π/ 2 , ϕ = 0) ∣ ∣ of the n = 5, l = m = 4 superradiant bound-state for a BH of spin χ = 0 . 99 at four consecutive values of α . The peak of the axion cloud shifts towards the BH with increasing α . This makes sense by comparison with the hydrogen atom: For the l = n -1 states of the hydrogen atom, the most probable radius r mp is inversely proportional to the electromagnetic fine-structure constant α EM : r mp = n 2 a B ∝ n 2 /α EM , where a B is the Bohr radius. While α EM is actually a constant, the analog for scalar field bound-states in Kerr, α ∝ µM , is different for each BH. As a result, for fixed n , the most probable radius for the l = n -1 bound-states decreases with α .  \n<!-- image -->  \nFIG. 19. The same four radial profiles of Fig . 18, with the addition of the lowest superradiant bound-state for α = 1 . 5 (with the BH spin still χ = 0 . 99). This bound-state is m = 5, and its peak is to the right of the m = 4 states. This is consistent with the hydrogen-atom expectation that r mp increases with n .  \n<!-- image -->  \n- [1] e. a. Abbott, B. P. (LIGO Scientific Collaboration and Virgo Collaboration), Gwtc-1: A gravitational-wave transient catalog of compact binary mergers observed by ligo and virgo during the first and second observing runs, Physical Review X 9 , 031040 (2019).\n- [2] e. a. Abbott, R. (LIGO Scientific Collaboration and Virgo Collaboration), Gwtc-2: Compact binary coalescences observed by ligo and virgo during the first half of the third observing run, Physical Review X 11 , 021053 (2021).\n- [3] The LIGO Scientific Collaboration, the Virgo Collaboration, the KAGRA Collaboration, and e. a. Abbott, R., GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo During the Second Part of the Third Observing Run, arXiv e-prints , arXiv:2111.03606 (2021), arXiv:2111.03606 [gr-qc].\n- [4] R. Abbott et al., LIGO Scientific Collaboration, and Virgo Collaboration, GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral, Physical Review Letters 119 , 161101 (2017), arXiv:1710.05832 [gr-qc].\n- [5] R. Abbott et al., LIGO Scientific Collaboration, and Virgo Collaboration, GW190412: Observation of a binary-black-hole coalescence with asymmetric masses, Physical Review D 102 , 043015 (2020), arXiv:2004.08342 [astro-ph.HE].\n- [6] R. Abbott et al., LIGO Scientific Collaboration, and Virgo Collaboration, GW190521: A Binary Black Hole Merger with a Total Mass of 150 M ⊙ , Physical Review Letters 125 , 101102 (2020), arXiv:2009.01075 [gr-qc].\n- [7] R. Abbott et al., LIGO Scientific Collaboration, and Virgo Collaboration, GW190814: Gravitational Waves from the Coalescence of a 23 Solar Mass Black Hole with a 2.6 Solar Mass Compact Object, The Astrophysical Journal Letters 896 , L44 (2020), arXiv:2006.12611 [astro-ph.HE].\n- [8] The LIGO Scientific Collaboration, the Virgo Collaboration, and the KAGRA Collaboration, Observation of Gravitational Waves from the Coalescence of a 2 . 5 -4 . 5 M ⊙ Compact Object and a Neutron Star, arXiv e-prints , arXiv:2404.04248 (2024), arXiv:2404.04248 [astro-ph.HE].\n- [9] T. N. Collaboration, The nanograv 15 yr data set: Evidence for a gravitational-wave background, The Astrophysical Journal Letters 951 , L8 (2023).\n- [10] P. Amaro-Seoane, J. Andrews, M. Arca Sedda, A. Askar, Q. Baghi, R. Balasov, I. Bartos, S. S. Bavera, J. Bellovary, C. P. L. Berry, E. Berti, S. Bianchi, L. Blecha, S. Blondin, T. Bogdanovi'c, S. Boissier, M. Bonetti, S. Bonoli, E. Bortolas, K. Breivik, P. R. Capelo, L. Caramete, F. Cattorini, M. Charisi, S. Chaty, X. Chen, M. Chru'sli'nska, A. J. K. Chua, R. Church, M. Colpi, D. D'Orazio, C. Danielski, M. B. Davies, P. Dayal, A. De Rosa, A. Derdzinski, K. Destounis, M. Dotti, I. Dutan, I. Dvorkin, G. Fabj, T. Foglizzo, S. Ford, J.-B. Fouvry, A. Franchini, T. Fragos, C. Fryer, M. Gaspari, D. Gerosa, L. Graziani, P. Groot, M. Habouzit, D. Haggard, Z. Haiman, W.-B. Han, A. Istrate, P. H. Johansson, F. M. Khan, T. Kimpson, K. Kokkotas, A. Kong, V. Korol, K. Kremer, T. Kupfer, A. Lamberts, S. Larson, M. Lau, D. Liu, N. Lloyd-  \nRonning, G. Lodato, A. Lupi, C.-P. Ma, T. Maccarone, I. Mandel, A. Mangiagli, M. Mapelli, S. Mathis, L. Mayer, S. McGee, B. McKernan, M. C. Miller, D. F. Mota, M. Mumpower, S. S. Nasim, G. Nelemans, S. Noble, F. Pacucci, F. Panessa, V. Paschalidis, H. Pfister, D. Porquet, J. Quenby, A. Ricarte, F. K. Ropke, J. Regan, S. Rosswog, A. Ruiter, M. Ruiz, J. Runnoe, R. Schneider, J. Schnittman, A. Secunda, A. Sesana, N. Seto, L. Shao, S. Shapiro, C. Sopuerta, N. C. Stone, A. Suvorov, N. Tamanini, T. Tamfal, T. Tauris, K. Temmink, J. Tomsick, S. Toonen, A. Torres-Orjuela, M. Toscani, A. Tsokaros, C. Unal, V. V'azquez-Aceves, R. Valiante, M. van Putten, J. van Roestel, C. Vignali, M. Volonteri, K. Wu, Z. Younsi, S. Yu, S. Zane, L. Zwick, F. Antonini, V. Baibhav, E. Barausse, A. Bonilla Rivera, M. Branchesi, G. Branduardi-Raymont, K. Burdge, S. Chakraborty, J. Cuadra, K. Dage, B. Davis, S. E. de Mink, R. Decarli, D. Doneva, S. Escoffier, P. Gandhi, F. Haardt, C. O. Lousto, S. Nissanke, J. Nordhaus, R. O'Shaughnessy, S. Portegies Zwart, A. Pound, F. Schussler, O. Sergijenko, A. Spallicci, D. Vernieri, and A. Vigna-G'omez, Astrophysics with the Laser Interferometer Space Antenna, Living Reviews in Relativity 26 , 2 (2023), arXiv:2203.06016 [gr-qc].  \n- [11] N. Aggarwal, O. D. Aguiar, A. Bauswein, G. Cella, S. Clesse, A. M. Cruise, V. Domcke, D. G. Figueroa, A. Geraci, M. Goryachev, H. Grote, M. Hindmarsh, F. Muia, N. Mukund, D. Ottaway, M. Peloso, F. Quevedo, A. Ricciardone, J. Steinlechner, S. Steinlechner, S. Sun, M. E. Tobar, F. Torrenti, C. Unal, and G. White, Challenges and opportunities of gravitationalwave searches at MHz to GHz frequencies, Living Reviews in Relativity 24 , 4 (2021), arXiv:2011.12414 [grqc].\n- [12] A. Arvanitaki and A. A. Geraci, Detecting HighFrequency Gravitational Waves with Optically Levitated Sensors, Physical Review Letters 110 , 071105 (2013), arXiv:1207.5320 [gr-qc].\n- [13] N. Aggarwal, G. P. Winstone, M. Teo, M. Baryakhtar, S. L. Larson, V. Kalogera, and A. A. Geraci, Searching for New Physics with a Levitated-Sensor-Based Gravitational-Wave Detector, Physical Review Letters 128 , 111101 (2022).\n- [14] G. Winstone, Z. Wang, S. Klomp, G. R. Felsted, A. Laeuger, C. Gupta, D. Grass, N. Aggarwal, J. Sprague, P. J. Pauzauskie, S. L. Larson, V. Kalogera, A. A. Geraci, and LSD Collaboration, Optical Trapping of High-Aspect-Ratio NaYF Hexagonal Prisms for kHz-MHz Gravitational Wave Detectors, Physical Review Letters 129 , 053604 (2022), arXiv:2204.10843 [physics.optics].\n- [15] S. R. Dolan, Instability of the massive Klein-Gordon field on the Kerr spacetime, Physical Review D 76 , 084001 (2007), arXiv:0705.2880 [gr-qc].\n- [16] D. Baumann, H. S. Chia, J. Stout, and L. ter Haar, The spectra of gravitational atoms, Journal of Cosmology and Astroparticle Physics 2019 (12), 006, arXiv:1908.10370 [gr-qc].\n- [17] A. Arvanitaki and S. Dubovsky, Exploring the string axiverse with precision black hole physics, Physical Review D 83 , 044026 (2011), arXiv:1004.3558 [hep-th].\n- [18] R. D. Peccei and H. R. Quinn, CP conservation in the presence of pseudoparticles, Physical Review Letters 38 , 1440 (1977).\n- [19] S. Weinberg, A new light boson?, Physical Review Letters 40 , 223 (1978).\n- [20] F. Wilczek, Problem of strong p and t invariance in the presence of instantons, Physical Review Letters 40 , 279 (1978).\n- [21] B. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama, Constraints on primordial black holes, Reports on Progress in Physics 84 , 116902 (2021), arXiv:2002.12778 [astroph.CO].\n- [22] E. D. Barr, A. Dutta, P. C. C. Freire, M. Cadelano, T. Gautam, M. Kramer, C. Pallanca, S. M. Ransom, A. Ridolfi, B. W. Stappers, T. M. Tauris, V. Venkatraman Krishnan, N. Wex, M. Bailes, J. Behrend, S. Buchner, M. Burgay, W. Chen, D. J. Champion, C. H. R. Chen, A. Corongiu, M. Geyer, Y. P. Men, P. V. Padmanabh, and A. Possenti, A pulsar in a binary with a compact object in the mass gap between neutron stars and black holes, Science 383 , 275 (2024), arXiv:2401.09872 [astro-ph.HE].\n- [23] J. C. Siegel, I. Kiato, V. Kalogera, C. P. L. Berry, T. J. Maccarone, K. Breivik, J. J. Andrews, S. S. Bavera, A. Dotter, T. Fragos, K. Kovlakas, D. Misra, K. A. Rocha, P. M. Srivastava, M. Sun, Z. Xing, and E. Zapartas, Investigating the Lower Mass Gap with Low-mass X-Ray Binary Population Synthesis, The Astrophysical Journal 954 , 212 (2023), arXiv:2209.06844 [astroph.HE].\n- [24] T. Sukhbold, T. Ertl, S. E. Woosley, J. M. Brown, and H. T. Janka, Core-collapse Supernovae from 9 to 120 Solar Masses Based on Neutrino-powered Explosions, The Astrophysical Journal 821 , 38 (2016), arXiv:1510.04643 [astro-ph.HE].\n- [25] T. Ertl, S. E. Woosley, T. Sukhbold, and H. T. Janka, The Explosion of Helium Stars Evolved with Mass Loss, The Astrophysical Journal 890 , 51 (2020), arXiv:1910.01641 [astro-ph.HE].\n- [26] K. Kruszy'nska, glyph[suppress] L. Wyrzykowski, K. A. Rybicki, M. Maskoli¯unas, E. Bachelet, N. Rattenbury, P. Mr'oz, P. Zieli'nski, K. Howil, Z. Kaczmarek, S. T. Hodgkin, N. Ihanec, I. Gezer, M. Gromadzki, P. Mikoglyph[suppress]lajczyk, A. Stankeviˇci¯ut˙e, V. ˇ Cepas, E. Pakˇstien˙e, K. ˇ Siˇskauskait˙e, J. Zdanaviˇcius, V. Bozza, M. Dominik, R. Figuera Jaimes, A. Fukui, M. Hundertmark, N. Narita, R. Street, Y. Tsapras, M. Bronikowski, M. Jabglyph[suppress]lo'nska, A. Jabglyph[suppress]lonowska, and O. Zi'oglyph[suppress]lkowska, Lens parameters for Gaia18cbf - a long gravitational microlensing event in the Galactic plane, Astronomy and Astrophysics 662 , A59 (2022), arXiv:2111.08337 [astro-ph.SR].\n- [27] D. Baumann, H. S. Chia, and R. A. Porto, Probing ultralight bosons with binary black holes, Physical Review D 99 , 044001 (2019), arXiv:1804.03208 [gr-qc].\n- [28] W. Gautschi, Computational Aspects of Three-Term Recurrence Relations, SIAM Review 9 , 24 (1967).\n- [29] S. J. Zhu, M. Baryakhtar, M. A. Papa, D. Tsuna, N. Kawanaka, and H.-B. Eggenstein, Characterizing the continuous gravitational-wave signal from boson clouds around Galactic isolated black holes, Physical Review D 102 , 063020 (2020), arXiv:2003.03359 [gr-qc].\n- [30] M. Isi, L. Sun, R. Brito, and A. Melatos, Directed searches for gravitational waves from ultralight bosons,  \n- Physical Review D 99 , 084042 (2019), arXiv:1810.03812 [gr-qc].\n- [31] C. Yuan, R. Brito, and V. Cardoso, Probing ultralight dark matter with future ground-based gravitationalwave detectors, Physical Review D 104 , 044011 (2021), arXiv:2106.00021 [gr-qc].\n- [32] H. Yoshino and H. Kodama, Gravitational radiation from an axion cloud around a black hole: Superradiant phase, Progress of Theoretical and Experimental Physics 2014 , 043E02 (2014), arXiv:1312.2326 [gr-qc].\n- [33] A. Irrgang, B. Wilcox, E. Tucker, and L. Schiefelbein, Milky Way mass models for orbit calculations, Astronomy and Astrophysics 549 , A137 (2013), arXiv:1211.4353 [astro-ph.GA].\n- [34] W. M. Farr, N. Sravan, A. Cantrell, L. Kreidberg, C. D. Bailyn, I. Mandel, and V. Kalogera, The Mass Distribution of Stellar-mass Black Holes, The Astrophysical Journal 741 , 103 (2011), arXiv:1011.1459 [astro-ph.GA].\n- [35] M. C. Miller and J. M. Miller, The masses and spins of neutron stars and stellar-mass black holes, Physics Reports 548 , 1 (2015), the masses and spins of neutron stars and stellar-mass black holes.\n- [36] Y. Shao and X.-D. Li, Population Synthesis of Black Hole X-Ray Binaries, The Astrophysical Journal 898 , 143 (2020), arXiv:2006.15961 [astro-ph.HE].\n- [37] K. C. Sahu, J. Anderson, S. Casertano, H. E. Bond, A. Udalski, M. Dominik, A. Calamida, A. Bellini, T. M. Brown, M. Rejkuba, V. Bajaj, N. Kains, H. C. Ferguson, C. L. Fryer, P. Yock, P. Mr'oz, S. Kozglyph[suppress]lowski, P. Pietrukowicz, R. Poleski, J. Skowron, I. Soszy'nski, M. K. Szyma'nski, K. Ulaczyk, glyph[suppress]L. Wyrzykowski, R. K. Barry, D. P. Bennett, I. A. Bond, Y. Hirao, S. I. Silva, I. Kondo, N. Koshimoto, C. Ranc, N. J. Rattenbury, T. Sumi, D. Suzuki, P. J. Tristram, A. Vandorou, J.-P. Beaulieu, J.-B. Marquette, A. Cole, P. Fouqu'e, K. Hill, S. Dieters, C. Coutures, D. Dominis-Prester, C. Bennett, E. Bachelet, J. Menzies, M. Albrow, K. Pollard, A. Gould, J. C. Yee, W. Allen, L. A. Almeida, G. Christie, J. Drummond, A. Gal-Yam, E. Gorbikov, F. Jablonski, C.-U. Lee, D. Maoz, I. Manulis, J. McCormick, T. Natusch, R. W. Pogge, Y. Shvartzvald, U. G. Jørgensen, K. A. Alsubai, M. I. Andersen, V. Bozza, S. C. Novati, M. Burgdorf, T. C. Hinse, M. Hundertmark, T.-O. Husser, E. Kerins, P. Longa-Pe˜na, L. Mancini, M. Penny, S. Rahvar, D. Ricci, S. Sajadian, J. Skottfelt, C. Snodgrass, J. Southworth, J. Tregloan-Reed, J. Wambsganss, O. Wertz, Y. Tsapras, R. A. Street, D. M. Bramich, K. Horne, I. A. Steele, and RoboNet Collaboration, An Isolated Stellar-mass Black Hole Detected through Astrometric Microlensing, The Astrophysical Journal 933 , 83 (2022), arXiv:2201.13296 [astroph.SR].\n- [38] K. El-Badry, H.-W. Rix, E. Quataert, A. W. Howard, H. Isaacson, J. Fuller, K. Hawkins, K. Breivik, K. W. K. Wong, A. C. Rodriguez, C. Conroy, S. Shahaf, T. Mazeh, F. Arenou, K. B. Burdge, D. Bashi, S. Faigler, D. R. Weisz, R. Seeburger, S. Almada Monter, and J. Wojno, A Sun-like star orbiting a black hole, Monthly Notices of the Royal Astronomical Society 518 , 1057 (2023), arXiv:2209.06833 [astro-ph.SR].\n- [39] C. L. Fryer, K. Belczynski, G. Wiktorowicz, M. Dominik, V. Kalogera, and D. E. Holz, Compact Remnant Mass Function: Dependence on the Explosion Mechanism and Metallicity, The Astrophysical Journal 749 , 91 (2012),  \narXiv:1110.1726 [astro-ph.SR].  \n- [40] C. S. Reynolds, Observational Constraints on Black Hole Spin, Annual Review of Astronomy and Astrophysics 59 , 117 (2021), arXiv:2011.08948 [astro-ph.HE].\n- [41] M. Kilic, J. A. Munn, H. C. Harris, T. von Hippel, J. W. Liebert, K. A. Williams, E. Jeffery, and S. DeGennaro, The Ages of the Thin Disk, Thick Disk, and the Halo from Nearby White Dwarfs, The Astrophysical Journal 837 , 162 (2017), arXiv:1702.06984 [astro-ph.SR].\n- [42] A. Helmi, Streams, substructures, and the early history of the milky way, Annual Review of Astronomy and Astrophysics 58 , 205 (2020), https://doi.org/10.1146/annurev-astro-032620-021917.\n- [43] P. J. McMillan, Mass models of the Milky Way, Monthly Notices of the Royal Astronomical Society 414 , 2446 (2011), arXiv:1102.4340 [astro-ph.GA].\n- [44] M. Juri'c, ˇ Z. Ivezi'c, A. Brooks, R. H. Lupton, D. Schlegel, D. Finkbeiner, N. Padmanabhan, N. Bond, B. Sesar, C. M. Rockosi, G. R. Knapp, J. E. Gunn, T. Sumi, D. P. Schneider, J. C. Barentine, H. J. Brewington, J. Brinkmann, M. Fukugita, M. Harvanek, S. J. Kleinman, J. Krzesinski, D. Long, J. Neilsen, Eric H., A. Nitta, S. A. Snedden, and D. G. York, The Milky Way Tomography with SDSS. I. Stellar Number Density Distribution, The Astrophysical Journal 673 , 864 (2008), arXiv:astroph/0510520 [astro-ph].\n- [45] N. Bissantz and O. Gerhard, Spiral arms, bar shape and bulge microlensing in the Milky Way, Monthly Notices of the Royal Astronomical Society 330 , 591 (2002), arXiv:astro-ph/0110368 [astro-ph].  \n- [46] R. Diehl, H. Halloin, K. Kretschmer, G. G. Lichti, V. Schonfelder, A. W. Strong, A. von Kienlin, W. Wang, P. Jean, J. Knodlseder, J.-P. Roques, G. Weidenspointner, S. Schanne, D. H. Hartmann, C. Winkler, and C. Wunderer, Radioactive 26 Al from massive stars in the Galaxy, Nature 439 , 45 (2006), arXiv:astro-ph/0601015 [astro-ph].\n- [47] A. Olejak, K. Belczynski, T. Bulik, and M. Sobolewska, Synthetic catalog of black holes in the Milky Way, Astronomy and Astrophysics 638 , A94 (2020), arXiv:1908.08775 [astro-ph.SR].\n- [48] E. Agol and M. Kamionkowski, X-rays from isolated black holes in the Milky Way, Monthly Notices of the Royal Astronomical Society 334 , 553 (2002), arXiv:astroph/0109539 [astro-ph].\n- [49] K. Riles, Searches for continuous-wave gravitational radiation, Living Reviews in Relativity 26 , 3 (2023), arXiv:2206.06447 [astro-ph.HE].\n- [50] N. Siemonsen and W. E. East, Gravitational wave signatures of ultralight vector bosons from black hole superradiance, Physical Review D 101 , 024019 (2020), arXiv:1910.09476 [gr-qc].\n- [51] E. W. Leaver, An Analytic Representation for the QuasiNormal Modes of Kerr Black Holes, Proceedings of the Royal Society of London Series A 402 , 285 (1985).\n- [52] E. Berti, V. Cardoso, and M. Casals, Eigenvalues and eigenfunctions of spin-weighted spheroidal harmonics in four and higher dimensions, Physical Review D 73 , 024013 (2006), arXiv:gr-qc/0511111 [gr-qc]."}