Simulating the Galactic population of axion clouds around stellar-origin black holes: Gravitational wave signals in the 10-100 kHz band
Jacob R. Sprague Shane L. Larson Zhiyuan Wang Shelby Klomp Andrew Laeuger George Winstone Nancy Aggarwal Andrew A. Geraci Vicky Kalogera


APS/123-QED





^1Center for Interdisciplinary Exploration and Research in Astrophysics (CIERA), Department of Physics and Astronomy, Northwestern University, Evanston, Illinois 60208
^2Center for Fundamental Physics, Department of Physics and Astronomy, Northwestern University, Evanston, Illinois 60208, USA
^3The Division of Physics, Mathematics and Astronomy, California Institute of Technology, Pasadena, California 91125, USA
^4Department of Physics and Astronomy, University of California, Davis, Davis, California 95616, USA



             

§ ABSTRACT 
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, 𝒪(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.





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
 (The LSD Collaboration)
    September 9, 2024
===================================================================================================================================================================================================






§ 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) <cit.><cit.><cit.>. The most notable events included the first NS-NS merger (GW170817) <cit.>, the first highly-asymmetric binary (GW190412) <cit.>, the first merger with an intermediate-mass BH remnant (GW190521) <cit.>, and the first object in the mass gap separating the most massive neutron stars from the lowest-mass BH's (GW190814) <cit.>. The first half of the fourth observing run has already seen a new lower-mass-gap event (GW230529) <cit.>.  
 

Adding 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) <cit.>. The most well-motivated scenario for the origin of this background is the extragalactic population of inspiralling supermassive BH binaries. 
 

Finally, the launch of the Laser Interferometer Space 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 extreme-mass-ratio-inspirals (EMRI's), are all highly anticipated LISA sources <cit.>. 
 

These 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 <cit.>.   
 

One 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érot 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 driven <cit.> <cit.> <cit.>. 
 

As 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 <cit.>. The result is a macroscopic cloud of bosons all living in the same state – commonly known as a  `gravitational atom' or `axion cloud' <cit.>. 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 μ = 𝒪(10^-11 - 10^-10) eV.  
 

This scenario can be realized with physics beyond-the-Standard-Model (BSM). For example, a large number of ultralight fields may occur as a result of the compactification of extra dimensions <cit.>. One of these may be the QCD axion – the pseudoscalar boson proposed to solve the strong-CP problem <cit.><cit.><cit.>. 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 <cit.>,

    μ≈ 6× 10^-10 eV(10^16 GeV/f_a)

where 10^16 GeV≡Λ_GUT is the grand unification (GUT) scale. An axion of mass 𝒪(10^-10) eV corresponds to f_a being at the GUT scale. However, as we will see in Sec. <ref>, signals in the LSD band are only expected up to ≈ 32 kHz, corresponding to a 6.6×10^-11 eV boson. 
 

At boson masses 𝒪(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 <cit.>, and BH's in the 1 - 4 M_⊙ range might be formed dynamically in binary neutron star mergers <cit.>, accretion-induced collapsing neutron stars <cit.>, or supernovae with unusually high fallback <cit.> <cit.>. 
 

The 1 - 4 M_⊙ range of BH masses is gradually being populated by microlensing candidates, X-ray binary candidates, and GW events such as GW190814 <cit.> and GW230529 <cit.>. 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.  
 
 
As 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 by 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.  
 

In 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.


§ 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 <cit.>. The scalar field then obeys the Klein-Gordon equation <cit.>,

    [g^μν∇_μ∇_ν - m_*^2]Φ(x, t) = 0

where 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. 
 

In Boyer-Lindquist coordinates, and with the ansatz

    Φ(x, t) = e^-i ω te^i m ϕ S(θ) R(r)

the 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,

    ω = ω_R + iω_I

with the consequence that bound-states must either grow or decay:

    e^-i ω t = e^-i (ω_R + iω_I) t = e^-i ω_R te^ω_I t


    Φ(x, t) = e^ω_I t[e^-i ω_R te^i m ϕ S(θ) R(r)]

For ω_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 <cit.>,

    ω_R < mΩ_H

This 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. 
 

It 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:

    α= r_g/λ_c = GM/c^2μ/ħ c = GMμ/ħ c^3

The `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 <cit.>:

    ħω_R = μ[1 - α^2/2n^2 + (2l - 3n + 1/l + 1/2 - 1/8)α^4/n^4 + 
    2mχα^5/n^3l(l + 1/2)(l + 1) + ...]

The 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 dimensionless product α. This motivates the introduction of a dimensionless eigenfrequency ξ = ξ_R + iξ_I:

    ξ_R≡ħω_R/μ      ξ_I≡ħω_I/μ

Once 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. 
 

From a practical point of view, this also means the superradiance condition

    ξ_R < mχ/2α[1 + √(1 - χ^2)]≡ξ_crit

becomes a tool for rapidly determining, for a given parameter set {μ, M, χ}, which states are superradiant. 
 

Additionally, far from the BH where relativistic effects are negligible, the radial equation reduces to (r measured in units of λ_c) <cit.>

    [-1/2r^2d/dr(r^2d/dr) - α/r + l(l + 1)/2r^2 + 1 - ξ^2/2]R(r) = 0

which is the radial Schrödinger equation for a non-relativistic particle in a Coulomb potential V(r) = α/r – hence the moniker `gravitational atom'. The α^2/(2n^2) term in Eq. <ref> is precisely the `hydrogen atom' solution to Eq. <ref>. The complete small-α solutions to the full Klein-Gordon equation have been computed order-by-order using the method of matched asymptotic expansions <cit.>,



    ξ_R = 1 - α^2/2n^2 + [2l - 3n + 1/l + 1/2 - 1/8]α^4/n^4 +                  
    2mχα^5/n^3l(l + 1/2)(l + 1) + ...


    ξ_I = 2(1 + √(1 - χ^2))[mχ/2α(1 + √(1 - χ^2)) - ξ_R]α^4l + 5  
    ×2^4l + 1(n + l)!/n^2l + 4(n - l - 1)!(l!/(2l)!(2l + 1))^2                      
    ×∏_j = 1^l[j^2(1 - χ^2) + [mχ - 2αξ_R(1 + √(1 - χ^2))]]

In general, ξ is a function of {n, l, m, α, χ}. For fixed α and χ, the superradiance rate is largest when m = l = n - 1. We consider only such bound states in our simulations of the Galactic axion cloud population, reducing the parameter list to {n, χ, α}. 
 

Since our fiducial model of the Galactic BH population will assume M ≥ 5 M_⊙, as well as boson masses 𝒪(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 infinite-series ansatz obey a three-term recurrence relation whose solution is equivalent to the solution of a corresponding non-linear continued-fraction equation <cit.> <cit.>. 
 

Denoting the peak mass of the cloud as M_c, the cloud's growth timescale is given by <cit.>

    τ_c≡τ_nlmln N = τ_nlmln(M_cc^2/μ)

with N the number of bosons in the cloud, and τ_nlm the reciprocal of the superradiance rate,

    τ_nlm≡1/Γ_nlm,     Γ_nlm≡ 2ω_I

τ_nlm is the e-folding timescale, and we follow the authors of <cit.> in taking τ_c as the time to fully grow the bound-state. The factor of two in Γ_nlm occurs because the cloud's density is proportional to the 00-component of the stress-energy, ρ∝ T^0_0∝exp(2ω_It). 
 

As 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 <cit.>. 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 

    M_c≡ M_i - M_f

and the hole's final mass & spin (M_f, χ_f) are given by <cit.><cit.>

    M_f = M_i[m^3 -√(m^6 - 16m^2ξ_R^2α_i^2(m - ξ_Rα_iχ_i)^2)/8ξ_R^2α_i^2(m - ξ_Rα_iχ_i)]


    χ_f = (M_i/M_f)^2(χ_i - m M_c/ξ_Rα_i M_i)

Since 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. 
 

The 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. <ref>) 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. <ref>, <ref>, and <ref>.


§ 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 downward-transitions between bound-states <cit.>. 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. 
 

Since 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 <cit.> <cit.>. Following <cit.>, the GW signal seen by a detector with perpendicular arms takes the general form

    h(t) = F_+(t)a_+cos[ϕ(t)] + F_(t)a_sin[ϕ(t)]

where F_+(t) and F_(t) are the detector's angular pattern functions. The amplitudes a_+ / are expanded in terms of spheroidal harmonics with spin-weight s = -2,

    a_+ /  = -∑_l≥ 2l h_0^(l)[_-2S_l,m,ω±_-2S_l,-m,-ω]
 
where ω = 2ω_R is the GW angular frequency, (l, m) refer to the scalar bound-state, and (l, m) refer to the GW modes, with l≥ 2l and m = 2m. For each mode, there is a polarization-independent characteristic amplitude h_0^(l) <cit.>:

    h_0^(l) = c^4/GM_c/M_f1/2 π^2 M_f f^2 d𝒜_lm(α_i, χ_i)

where f is the GW frequency, d is the source distance, and the 𝒜_lm(α, χ) are dimensionless numerical factors which measure how much energy is carried by each mode. The corresponding luminosity in each mode is given by

    E_GW(l, m, ω) = c^5/4π G(c^3/GM_fω)^2M_c^2/M_f^2𝒜_lm^2(α_i, χ_i)

In principle, the coefficients 𝒜_lm must be computed numerically by solving the Teukolsky equation governing linear perturbations of the Kerr metric. The authors of <cit.> express E_GW in the form

    E_GW = c^5/GM_c^2/M_f^2d E/dt

and invoke an analytic solution for dE/dt which is formally valid for α≪ l, and which remains a good approximation up to α∼ l <cit.>:

    dE/dt = 16^l + 1 l(2l - 1) Γ^2(2l - 1) Γ^2(n + l + 1)α_f^4l + 10/n^4l + 8(l + 1)Γ^4(l + 1)Γ(4l + 3) Γ^2(n - l)

where Γ 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),

    α_f = α_iM_f/M_i

Comparing Eqs. <ref> and <ref>, we see that 𝒜_lm∝√(dE/dt), allowing us to express h_0^(l) directly in terms of dE/dt. Restricting ourselves to the dominant mode m = l = 2l, 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 (2l) henceforth),

    h_0^(2l)(d) = GM_c/c^2d2√(π)M_i/ξ_Rα_iM_f√(d E/dt)

The corresponding GW frequency is given by

    f = ω/2π = 1/2π2μ/ħξ_R≡ f_0ξ_R


    f = f_0ξ_R

where we've introduced the zeroth-order frequency f_0 = ω_0/2π, ω_0≡ 2μ/ħ. 
 

It 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. <ref>, Eq. <ref>). 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. 
 

Eq. <ref> 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 <cit.> (r, in kpc, is the cylindrical radial distance from the Galactic center):

    v_ϕ(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

Denoting the source-frame frequency as f_s, the non-relativistic Doppler-shifted frequency we observe is

    f_obs = (1 - v_r/c)f_s

where 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.  
 

When a cloud finishes growing, it emits GW's whose initial amplitude h_0 is given by Eq. <ref>. As the cloud dissipates, the amplitude decreases as <cit.>

    h(t) = h_0/1 + t/τ_GW

where τ_GW is the time for h to drop to half its initial value.


§  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.  
 

Our knowledge of the stellar-origin BH mass distribution is informed by mass measurements in X-ray binary systems <cit.> <cit.> <cit.>, microlensing events <cit.>, and astrometry <cit.>, as well as through modelling of the complex physics of core-collapse supernovae <cit.>. 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 <cit.>. Not coincidentally, the massive stars which produce BH remnants are also characterized by a power-law distribution, ψ(M)dM ∝ M^-2.35dM – the `Salpeter' function. We will assume M_BH to be Salpeter-distributed on the interval 5 - 20 M_⊙.  
 

BH spins have been measured in several X-ray binaries <cit.>, 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]. 
 

The 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(810^9)] and x ∼ U[3, 10], respectively <cit.>. For the bulge, we assign each BH an age 10^x yr, with x ∼ U[9, log_10(1310^9)] <cit.>.  
 
 
We assume black holes are distributed in space according to the mass profiles of the disks & bulge described in Ref. <cit.>. Both disks have the same axisymmetric form, with the corresponding scale lengths, scale heights, and surface densities quoted in Table <ref>:

    ρ_disk(r, z, ϕ) = Σ_d, 0/2 z_d e^-z/z_d e^-r/R_d

The bulge is also axisymmetric, with the corresponding parameters also given in Table <ref>:

    ρ_b = ρ_b,0/(1 + r'/r_0)^αe^-(r'/r_cut)^2


    r' ≡√(r^2 + (z/q)^2)

































We 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/∑_iM_i, i ∈ {thin, thick, bulge}. The disk masses are obtained by integrating ρ_disk, with the radial integral cut-off at to 25 kpc, and the vertical integral cut-off at 3 scale heights. This gives 3.9710^10 M_⊙ and 1.510^10 M_⊙ for the thin and thick disks, respectively. We take the bulge mass to be 8.910^9 M_⊙, the value quoted in <cit.>. 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.


§  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. 
 

At 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. <ref>, determines which bound-state grows first.

    ξ_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)


Since ξ_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.  
 
 
The 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.999410^7 yrs – still plenty of time left for further superradiance.  
 
 
The 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.910^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. <ref>). The signal observed today was emitted d/c = 3300 yrs ago, so the corresponding amplitude h(t) = 6.910^-27 (Eq. <ref> with t = t_r - d/c). 
 

Our 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. 
 

Our 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. 
 

Those BH's with an extant cloud are assigned a location in the Milky Way (Eqs. <ref> and <ref>). 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 clouds 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: 


    For a given μ, M_BH, χ, and BH age, find the lowest superradiant value of n. → If τ_c > τ_uni, the BH is removed from the simulation. 

    → 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).  

    → 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.  

    → 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. 



§  GW'S FROM THE AXION CLOUD POPULATION




 §.§  Cloud populations


The total number of stellar-origin black holes has been estimated to be 𝒪(10^8) from the Milky Way's supernova rate of 𝒪(1) century^-1 <cit.>, and from population-synthesis estimates <cit.>. 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 <cit.>. We have simulated the axion cloud population for μ = (3, 3.5, 4, 4.5, 5, 5.5, 6, 6.5)×10^-11 eV.  
 

The 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.  
 

For 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 sample mean & sample standard deviation of N_c over the 20 trials, and then multiplying them by 20 and √(20), respectively. 
 

An ensemble of GW signals from axion clouds is a scatter plot in the h vs. f plane, as in Figs. <ref>, <ref>, and <ref>. The distribution of amplitudes and frequencies is not random, but consists of well-defined bands corresponding to the various occupied bound-states. 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.  
 

Also reflecting this difficulty is the rapid decline in the number of clouds N_c with increasing boson mass μ (Fig. <ref>). For μ=310^-11 eV, N_c = (9.323±0.007)×10^5, while at μ=6.510^-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. 
 

In the introduction (Sec. <ref>), we noted a potential connection between the QCD axion and the GUT scale Λ_GUT (Eq. <ref>): An axion of mass 𝒪(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 𝒪(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 𝒪(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 𝒪(10^-11) eV bosons.





































 §.§  Resolvable signals


The standard result for coherent detection of a continuous monochromatic signal, h(t) = h_0cos (ω t), is that the signal-to-noise ratio (SNR) ρ grows as the square-root of the coherent integration time T_coh <cit.>

    ρ = h_0√(T_coh)/√(S_n(f))

where √(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. Although the LSD is an Earth-bound detector for which the observed signal is modulated by the Earth's daily (diurnal) rotation and orbital motion, we, instead, 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). This scenario isolates the intrinsic sensitivity of the LSD to continuous-wave signals from incidental factors, such as the Earth's rotational and orbital periods.  
 

Taking 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 <cit.>, 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. <ref>). 
 

The `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. <ref> and <ref>) of the source. Thus, for a given set of extant clouds, the number of individually-resolved 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 `re-shuffling'. The mean & standard deviation are then computed as

    N_res = 1/N_reshuffle∑_i = 1^N_reshuffleN_res, i


    σ_res = 1/N_reshuffle - 1∑_i = 1^N_reshuffle(N_res, i - N_res)^2

With 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 𝒪(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, 𝒪(1) resolvable signals appear at μ = 3× 10^-11 eV, while a 1-m instrument does not have the required sensitivity to detect individual sources. 
 

In the event a continuous monochromatic signal is 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 <cit.>. 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. <ref>. 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.












    

 §.§  Unresolved signals




For all boson masses, the majority of GW signals have amplitudes less than 10^-23, with the weakest having h = 𝒪(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. <ref>), we bin the cloud amplitudes (with bin width δ f = 10^-2f_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 an rms amplitude, defined as follows, with each bin. 
 

We 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. <ref>. All axion clouds emit monochromatic signals,

    h_i(t) = h_0,icos(2π f t + ϕ_i)

where 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 time-averaged power associated with that bin. The square-root of the power represents an effective amplitude h_eff of the confusion foreground in the bin,

    h_eff= √(1/T_c∫_0^T_cdt [∑_ih_i(t) ]^2)

We then create a new bin with center frequency f_c|_new and width δ f|_new,

    f_c|_new = f_c|_old + 10^-2δ f|_old


    δ f|_new = 10^-2f_c|_new

and 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 adjacent 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. <ref>).  
 

A 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. <ref>, is then compared to a threshold ρ_t. We continue to require ρ_t = 10. The numerator and denominator of Eq. <ref> (h_eff√(T_coh) and √(S_n(f_c)), respectively) are shown separately in Fig. <ref>, and their ratio (the SNR) is shown in Figs. <ref>, <ref>, and <ref> for the 1-, 10-, and 100-m instruments, respectively. 
 

We 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 = 𝒪(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. <ref>). A 100-m instrument could detect the foreground with large ρ_eff if μ∈ (3-6)×10^-11 eV (Fig. <ref>). In the range (3-3.5)×10^-11 eV, the peak value of ρ_eff is 𝒪(10^3), and remains 𝒪(10^2) up to 5× 10^-11 eV. 
 









§ 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. 
 

The 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.   
 
 
The 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 𝒪(1) resolvable signals occur in our simulation at μ = 3×10^-11 eV.    
 

Meanwhile, 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. 
 

Finally, 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 searches. A full template search will require us to pixelate the sky, and it is an open question what is the smallest angular size Δθ for which the number of templates in a blind all-sky search is not prohibitively large: We do not want the time required for data analysis to become larger than the four-month observation period. Transverse proper motions will also need to be accounted for if, over the observation time, a source moves out of the pixel on the sky it was in initially. 
 

Diurnal and annual modulations of the GW frequency & amplitude must also be included in our detection scheme. The angular dependence of the detector sensitivity, encoded in the detector's pattern functions, produces a daily modulation of the signal amplitude. The rotational and orbital motions of the Earth produce time-varying Doppler shifts. The Doppler modulation can be corrected, but doing so requires knowledge of the source position to a precision determined by the observation time and GW frequency. 
 

In this work, we have used the SNR as a baseline detection statistic. We have not yet determined what are the most appropriate detection statistics & criteria for the monochromatic signals from individual clouds. Moreover, 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. Additionally, 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.






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 gravitational-wave 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 PHY-2111544, 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 University 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. 


*



§ 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, `superradiance', 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. 
 

In 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.  
 
 
The story for waves runs analogously: An incident wave with amplitude ℐ splits into a part transmitted into the BH (with amplitude 𝒯) and a part which escapes (the reflected wave with amplitude ℛ). If the transmitted wave is counter-rotating, the rotational energy of the BH decreases, leading to an outgoing wave with ℛ > ℐ. 
 

The 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 <cit.><cit.>. As a result, the growth of an axion cloud begins immediately after the birth of a BH. 
 

An 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:

    ds^2 = -[1 - 2GMr/c^2ρ^2]c^2dt^2 - 4GMar sin^2θ/c^2ρ^2cdtdϕ
     + ρ^2/Δdr^2 + ρ^2dθ^2 + [r^2 + a^2 + 2GMa^2rsin^2θ/c^2ρ^2]sin^2θ dϕ^2

where M is the BH mass, J is the BH angular momentum, ρ^2≡ r^2 + a^2cos^2θ, a ≡ J/(Mc) is the Kerr parameter, and Δ≡ r^2 - 2r_gr + 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

    r_± = r_g[1 ±√(1 - χ^2)]

It follows that χ is restricted to the interval 

    0 < χ < 1

The event horizon is located at r = r_+, and the angular velocity of the horizon is

    Ω_H = c χ/2 r_+

The scalar field obeys the Klein-Gordon equation,

    [∇_μ∇^μ - m_*^2]Φ(x, t) = 0

where ∇_μ 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

    Φ(x, t) = Re[ e^-i ω te^i m ϕ S(θ) R(r)]

Invoking the identity

    ∇_μ∇^μΦ = 1/√(-g)∂_μ[√(-g) g^μν∂_νΦ]


    √(-g) = ρ^2sinθ

the Klein-Gordon equation separates into two 2^nd-order linear homogeneous ODE's for R(r) and S(θ):







    𝒟_θ[S] + [χ^2α^2(ξ^2 - 1 )cos^2θ -m^2/sin^2θ + Λ]S(θ) = 0


    𝒟_r[R] + [α^2ξ^2(r^2 + χ^2)^2 - 4χ m αξ r +  m^2χ^2
     - Δ(α^2r^2 + χ^2α^2ξ^2 + Λ) ]R(r) = 0


    𝒟_θ≡1/sinθd/dθ[sinθd/dθ],     𝒟_r≡Δd/dr[Δd/dr]
 

We have expressed the decoupled equations in terms of the dimensionless variables (χ, α and ξ) used in the main text. The radial coordinate in <ref> is measured in units of r_g. 
 

Bound-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 horizon to -∞,

    dr_*/dr = r^2 + a^2/Δ
 

This means that plane waves at the event horizon (r_*→ -∞) can only move `to the left', i.e. into the black hole. 
 
 
The spectra of both bound-states and BH quasi-normal modes  can be found via Leaver's continued-fraction method <cit.> <cit.>. The radial function R(r) is represented by an infinite series,

    R(r) = (r - r_+)^-iσ(r - r_-)^iσ + β - 1e^qr∑_n = 0^∞a_n(r - r_+/r - r_-)^n


    σ = α(1 + √(1 - χ^2))(ξ - ξ_crit)/√(1 - χ^2)


    q = α√(1 - ξ^2)


    β = α^2(1 - 2ξ^2)/q

(The quantity we denote by β is the same as the quantity denoted by χ in Ref. <cit.>.) With this ansatz, <ref> implies a three-term recurrence relation for the unknown coefficients a_n,

    α_na_n + 1 + β_na_n + γ_na_n - 1 = 0,     n = 1, 2, …
    
            a_1 = -β_0/α_0a_0

where the coefficients α_n, β_n and γ_n are defined by

    α_n = n^2 + (c_0 + 1)n + c_0
    β_n = -2n^2 + (c_1 + 2) + c_3
    γ_n = n^2 + (c_2 - 3)n + c_4

and c_0, c_1, c_2, c_3 and c_4 are given by

    c_0 = 1 - 2 i αξ -2 i (αξ -m χ/2)/√(1-χ ^2)


    c_1 = -4 + 4i[αξ - iα√(1-ξ ^2)(1 + √(1 - χ^2))]           
     + 4i (αξ -m χ/2)/√(1 - χ^2) - 2 [α ^2 ξ ^2+α ^2 (1-ξ^2)]/α√(1 - ξ^2)


    c_2 = 3 - 2i αξ - 2 [α^2 (1 - ξ^2)-α^2 ξ^2]/α√(1-ξ^2) - 2 i (αξ -m χ/2)/√(1 - χ^2)


    c_3    = 2i (αξ - iα√(1 -ξ^2))^3/α√(1 - ξ^2) + χ^2α^2(1 - ξ^2) 
        - Λ_lm - 1 + 2√(1 - χ^2)(αξ - iα√(1 - ξ^2))^2 
        + 2imχα√(1 - ξ^2) - (αξ - iα√(1 - ξ^2))^2/α√(1 - ξ^2)
        + 2α√(1 - ξ^2)√(1 -χ^2)
        + 2i/√(1 - χ^2)[1 + (αξ - iα√(1 - ξ^2))^2/α√(1 - ξ^2)][αξ - mχ/2]


    c_4 = (αξ - iα√(1 - ξ^2))^4/α^2 (1 - ξ^2) + 2iξ(αξ - iα√(1 - ξ^2))^2/√(1 - ξ^2)
     
        - 2i (αξ - iα√(1 - ξ^2))^2 (αξ
       -mχ/2)/α√(1 - ξ^2)√(1 - χ^2)
 

The series coefficients are related by an infinite continued-fraction <cit.>

    a_n + 1/a_n = -γ_n + 1/β_n + 1 - α_n + 1γ_n + 2/β_n + 2 - …

Continued-fractions are commonly written in the slightly-less cumbersome notation

    a_n + 1/a_n = -γ_n + 1/β_n + 1 -α_n + 1γ_n + 2/β_n + 2-α_n + 2γ_n + 3/β_n + 3-…


Since a_1/a_0 = -β_0/α_0, we obtain a condition whose roots are the desired bound-state frequencies:

    β_0 - α_0γ_1/β_1 -α_1γ_2/β_2 -α_2γ_3/β_3 -… = 0
 

Strictly speaking, the radial and angular eigenvalues, ξ and Λ, must be found simultaneously. Leaver's method can also be applied to <ref> <cit.>, resulting in a continued-fraction condition analogous to <ref>. We then have two equations for the two unknowns. 
 

Conveniently, we can reduce the root-finding problem to merely solving <ref> by using the Mathematica function `SpheroidalEigenvalue'. With the change of variable z = cosθ, and in terms of the following quantities,

    γ^2≡χ^2α^2(1 - ξ^2)


    λ≡Λ - γ^2

the angular equation <ref> takes the standard form implemented in Mathematica:

    (1 - z^2)d^2S/dz^2 - 2zdS/dz + [γ^2(1 - z^2) + λ - m^2/1 - z^2]S(z) = 0
 

`SpheroidalEigenvalue' yields λ, and `SpheroidalPS' yields S(z). The continued-fraction equation <ref>, with Λ replaced by γ^2 + `SpheroidalEigenvalue', can then be solved for ξ with the Mathematica function FindRoot.  
 

For 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. <ref>, <ref> and <ref>; Figs. <ref> and <ref> are analogous to Figs. <ref> and <ref>.
























