Black Hole Spin Constraints on the Mass Spectrum and Number of Axion-like Fields

Astrophysical observations of spinning BHs, which span $ 5M_\odot\lesssim M_{\rm BH}\lesssim 5\times 10^8 M_\odot$, can be used to exclude the existence of certain massive bosons via the superradiance phenomenon. In this work, we explore for the first time how these measurements can be used to constrain properties of statistical distributions for the masses of multiple bosonic fields. Quite generally, our methodology excludes $N_{\rm ax}\gtrsim 30$ scalar fields with a range of mass distribution widths and central values spanning many orders of magnitude. We demonstrate this for the specific example of axions in string theory and M-theory, where the mass distributions in certain cases take universal forms. We place upper bounds on $N_{\rm ax}$ for certain scenarios of interest realised approximately as mass distributions in M-theory, including the QCD axion, grand unified theories, and fuzzy dark matter.


I. INTRODUCTION
The Penrose process [1] allows bosonic waves infalling into a Kerr black hole (BH) to emerge with more energy than incident upon entry at the horizon, in exact analogy to other superradiant processes in physics, such as Cherenkov radiation. If the bosons can be confined around the BH by a mirror, then this amplification process continues without limit leading to Press and Teukolsky's "black hole bomb" scenario [2,3]. Massive bosonic fields on a Kerr spacetime possess hydrogenic bound states. In this case the potential barrier provided by the particle mass can play the role of the mirror, leading to a natural realisation of the BH superradiance process for massive bosons in orbits around astrophysical BHs (see Ref. [4] for a review). The historic Laser Interferometry Gravitational-Wave Observatory (LIGO) observations of gravitational waves from the binary coalescence of astrophysical BHs has ushered in a new era of interest in BH physics [5]. Gravitational wave data can be used to infer the mass and spin of the two BHs in the binary. LIGO has the prospects to detect the existence of many hundreds of such events, accurately determining the mass and spin distribution of BHs. The future of BH superradiance constraints derived from LIGO, the growing global network of GW observatories, and future space-based missions, is extremely promising as a probe Certain ranges of σ correspond closely to RMT and M-theory mass spectra, and can also be used to approximate the log-flat spectrum. For 1 σ 20, Nax ≥ 30 is excluded for an extremely wide range of central masses. Constraints neglect axion self-interactions and apply approximately in the limit of large decay constants, fa 10 14 GeV.
The ability to constrain ultralight bosonic fields from BH-scalar condensate systems come in the form of two phenomena. It may be possible to identify the presence of scalar clouds in the vicinity of BHs as emission sources of monochromatic gravitational waves (GWs). The signal frequency, f ∼ µax /π, with boson mass, µ ax could potentially be detected by either ground or space-based GW observatories and proposes to be an exciting methodology to enhance constraints on the mass bounds for bosonic fields. This subject has been extensively discussed in Refs. [11][12][13][14]. The second phenomenon of interest, and the subject of this work, is the spin down of astrophysical BHs. If the superradiance rate is faster than any other astrophysical process affecting the BH mass, M BH , and dimensionless spin, a * , then the BH superradiance process can efficiently reduce these quantities. This occurs when the boson Compton wavelength is of the order of the gravitational radius of the BH. Thus, if a massive boson exists, then astrophysical BHs of particular values in the (M BH , a * ) "Regge plane" (which, according to the no-hair theorems, gives a complete description of spinning BHs) should be absent in observations. The masses and spins of a large number of astrophysical BHs have been measured, often incorporating either Xray reflection spectroscopy or continuum-fitting methods (see Table I for BH parameter measurements and corresponding references). These measurements can be used to probe the possible existence of massive bosons [15,16]. BH superradiance constraints apply to a range of particle physics models, including a possible mass for the graviton or the photon [17] (and indeed to the photon plasma mass near the BH), as well as to exotic particles, such as massive vector (Proca) fields [7], massive spintwo fields [18], and axion-like particles and other massive scalars [11,19,20].
These are powerful and generic exclusions, but they leave many axion models of interest unconstrained. Stellar BHs are too heavy to place constraints on the QCD axion [21][22][23] possessing a decay constant far below the Planck scale [20]. "Fuzzy dark matter (DM)" with µ ax ≈ 10 −22 eV [24][25][26][27][28], which has novel effects on the formation of galaxies, is too light to make predictions about the spin distribution of SMBHs with M BH < 10 9 M that inhabit the centres of galaxies. Finally, the axion mass scale associated to grand unification (GUTs) in M-theory, µ ax ≈ 10 −15 eV [29] is in the "desert" of intermediate mass BHs (IMBHs) which so far have not been observed. There is hope, however, since each of these models is only a small logarithmic distance from the BH superradiance constrained regions, while axion models typically have a spectrum spanning many orders of magnitude [19,30]. All previous studies of BH superradiance constraints on bosons have focused on the range of excluded masses assuming the existence of a single new bosonic field. In the present work we assess, for the first time, what constraints can be drawn on the properties of axion mass distributions from BH superradiance.
Any string or M-theory model that realises one of the models of interest (QCD axion, fuzzy DM, or GUTs) will likely contain a distribution of masses around this value. Even a small spread on a logarithmic scale could lead to strong constraints on the model. The central observation of the present work is that, simply from the statistical overlap between a mass distribution and the BH superradiance bounds, it is possible to place constraints on the allowed mass distributions of axions. Furthermore, these constraints get increasingly more stringent as the number of axions increases, placing upper bounds on the number of axion-like fields. Consider the following toy model. In Ref. [19] it was suggested that axion masses have a log-flat distribution from the Planck scale to the Hubble scale, covering approximately sixty orders of magnitude. The BH superradiance constraints cover approximately four orders of magnitude. Assuming independent and identically distributed draws from the log-flat distribution, this naive model of the axiverse is excluded with probability P = 1 − (56/60) N , which is greater than 95% C.L. if N ax ≥ 44. Clearly, the model with a log-flat prior on the axion mass is excluded by BH superradiance for large numbers of fields. The exclusion is a function of the upper and lower bounds on the mass spectrum. The  The solid black line represents the unity limit for nonrelativistic and relativistic regimes. The dashed line corresponds to α = 0.5, the approximate limit in which the analytical approximation for the instability rate is valid. Dotted lines correspond to frequency ranges for monochromatic gravitational wave emission from the scalar cloud accessible to current and future GW observatories [51][52][53][54][55][56][57].
constraint gets considerably stronger if the upper bound is below the Planck scale, and vanishes if the distribution does not extend below about 10 −11 eV. Such a truncated spectrum, on the other hand, cannot realise many of the models of interest discussed above. Fortunately for phenomenologists, the mass distributions arising from RMT models are not log flat from the Hubble scale to the Planck scale. The log-normal distribution, centred on a particular mean mass, µ ax , and with a variance σ 2 , provides a useful benchmark, covering different types of models. For small σ, it resembles a degenerate spectrum, large σ is approximately log-flat, and intermediate values of σ are statistically similar to eigenvalue distributions found in RMT and M-theory. Fig. 1 summarises our conclusions, showing the allowed number of axionic fields drawn from log-normal distributions as a function of the width and central value.
The structure of this work is as follows: Section II contains a brief review of BH superradiance along with the Regge plane and BH spin measurements, while Section III overviews our models for the axion mass matrix and collects our BH data. In Section IV we present constraints on axion mass spectra from BH mass and spin measurements under a frequentist framework. We first reproduce the known single-field results and then move on to considering mass distributions. We conclude our work in Section V. Further details of our BH superradiance calculations are given in Appendix A. Appendix B describes our statistical methods, which we believe are somewhat novel in this context. Appendix C collects results from Ref. [30] on the axion mass matrix and RMT.

A. Scalar Fields on Kerr Background
The action for N real scalar fields Ψ i with masses µ i takes the form where ∇ µ is the covariant derivative on the spacetime with metric g. The metric is assumed to be the Kerr metric for a spinning BH. A review of the Kerr geometry is given in Appendix A 1. This geometry is taken as a background. The superradiant process leads to time dependence of the BH mass and spin, but the structure of the metric does not change due to backreaction. It is known for single field superradiance that the backreaction of the scalar condensate on the Kerr geometry is small. This is because, although the cloud can obtain a large mass, it is distributed over a large volume compared to the BH, leading to low scalar energy density (and thus a low source of curvature) in the cloud [58]. Concerns that backreaction is a more severe problem with large numbers of fields as opposed to dealing with a single field can be alleviated considering the properties of the scalar cloud. The gravitational backreaction is a function of M S /M BH , where M S is the total mass in the scalar cloud. There is a maximum value of M S independent of the number of axion fields, which is determined by the BH mass at the initial spin, M BH (a * ), and the irreducible mass after all the spin has been extracted, M BH (a * = 0). N ax fields cannot extract any more total mass than a single field, and for resonant modes the cloud size is of the same order of magnitude for all the fields, therefore gravitational backreaction is not enhanced to a greater severity than the single field case. Non-linearities coming from axion interactions, on the other hand, can increase with the number of fields. We discuss this briefly later.
Thus, neglecting the self-interactions, each field Ψ i evolves independently on the fixed background. In this separable limit, the total rate of the superradiant process is given simply by the sum of the single field rates: Solutions of the single field Klein-Gordon equation are discussed in detail in Appendix A 2, and the superradiance phenomenon for multiple fields is described in terms of these.

B. Superradiance
Astrophysical BHs with a mass M BH and spin J = aM BH will spin down via superradiant instabilities extracting energy and angular momentum [4,11], forming very large gravitationally bound states comprising of a scalar cloud containing exponentially large axion population numbers. Axions bound in this way with a BH form a gravitational atom, where superradiant instabilities are found to be strongest when the Compton wavelength of the field, λ ax =h/µaxc is comparable to the Schwarzschild radius of the BH, r s = 2GMBH /c 2 .
The condition for mode amplification of the scalar field requires the angular velocity of the BH horizon to exceed the angular phase velocity of the wave mode, defining the superradiance condition (see Fig. 3) where m is the spherical harmonic quantum number. The effective angular velocity of the BH as a function of the dimensionless rotation spin parameter is where a * is defined in region 0 ≤ |a * | < 1 as , in Boyer-Lindquist coordinates. The gravitational radius of the BH is, In parts of the following we shall work in units c = = G = 1 such that r g ≡ M BH . The Kerr-Klein-Gordon system admits quasi-bound states with complex eigenfrequencies where {ω R , ω I } ∈ R. Kerr BHs present a critical frequency for superradiant scattering with m representing the angular momentum about the BH spin axis. This defines the stability thresholds for the scalar modes: For values of ω nlm satisfying 0 < ω nlm < ω c the imag- Timescale ratios for the superradiance rates for an axion with mass µax = 10 −11.5 eV compared with a typical BH astrophysical timescale, here taken to be τ Salpeter (Eq. (29)). Each cusp represents the analytical limit beyond which Eq. (3) is satisfied. The limit to the right of the cusp (sold line) represents the ratio defining the nature of the timescales where superradiance is apparent. The red volume defines the limit in the two dimensional BH mass/spin parameter space where superradiance occurs within the defined astrophysical timescale used to map the Regge plane isocontour limits.
inary component is positive defining the superradiant regime. Scalar modes in the presence of the Kerr BH spacetime with scalar mass, µ ax contain a natural confinement mechanism in the limit where they are bounded from escaping via their potential (Eq. (A15)). Modes satisfying these conditions will grow exponentially over time identifying the presence of an instability in the Kerr spacetime. When ω nlm = ω c the imaginary component of the frequency drops out allowing for the formation of bound states or scalar clouds. Aside from regions within a significant proximity to the BH the gravitational potential is ∝ 1 /r where the spherically symmetric properties of the potential to leading order allow for a separation of variables of the field evolution in the background reproducing a Schrödinger type wave-equation (see Section A 2). The equation for the separated radial wave function (Eq. (A10)) is the equivalent to that of the Scalar Coulomb, thereby presenting hydrogenic wavefunctions. To leading order the energy levels for the bound states are well approximated by the spectrum of the hydrogen atom in the non-relativistic limit. When the superradaiance condition is saturated the eigenfrequencies take the approximate form, The orbitals around the BH are indexed by the overtone (n), orbital multi-pole (l) and azumutal (m) quantum numbers satisfy l ≤ n − 1 and |m| ≤ l forming discrete sets, {n,l,m} used to quantise the superradiant behaviour. Superradiance requires evolving modes to corotate with the BH which satisfy, m > 0. Details of the methodology used to determine the approximated eigenspectrum are given in Appendix A. The dimensionless coupling of the gravitational BH-scalar condensate system is, in our choice of units. Fig. 2 presents the coupling strength for potential regions of the axion mass parameter space open to investigation for BH masses spanning the stellar and supermassive limits.

C. Superradiance Rates
The evolution of the axion field is defined by the characteristic eigenfrequencies corresponding to the instability timescales for the unstable modes of the system. The nature of scalar instabilities is well researched covering both the frequency [59][60][61] and time domains [62]. In the frequency regime in order to extract valid quasibound state instability rates, Γ nlm , which depend on the wavefunction near the horizon, either one of two approaches can be implemented. The superradiance rates are defined as the small imaginary component of the energy of the free field solution on the Kerr background. Analysing the region of the parameter space where α ∼ 1, solutions for the unstable modes can be found using a numerical analysis of the wave equation (see Appendix A 4) [61,63,64]. When α surpasses unity WKB methods are formulated to evaluate the rate, presenting an exponential suppression proportional to α where Γ nlm ∝ e −3.7α [11,60].
It has been shown it is possible to find analytical solutions to approximate the instability rate, incorporating matching techniques between different regimes of validity as a a function of α. For a particular bound state if the superradiance condition is satisfied then providing that the instability rate is quicker than relevant astrophysical timescales, wave modes will extract energy and angular momentum from the BH. It has been shown in the α 1 regime known as the "small mass approximation" the evolution of the superradiant instability can be analytically described via a matched asymptotic expansion. This solution was initially derived by Detweiler to solve the Klein-Gordon equation of the scalar field perturbation [59]. Comparing the large r behaviour of the near-region solution with the small r behaviour of the far-region solution yields the allowed values of the small imaginary component of the frequency ω I . The instability rate in the small mass approximation is defined as where, It can be seen from Eq. (14) the superradiance rates for scalar fields scale approximately as which is maximised close to the superradiance boundary.
In Fig. 4 we present the superradiance rates for a range of modes and spins as a function of the axion/BH coupling, µ ax M BH . The fastest growing mode occurs for Γ 011 with the superradiance rates exponentially suppressed for higher values of l. The maximum superradiance rates are found by fixing the values of l and m such that, l = m where m determines the ability to satisfy the superradiance condition in Eq. (3) (See right panel of Fig. 4). The value of Γ nlm has a limited dependance on the overtone mode, n. When the BH possesses significant spin higher order overtone modes for larger values of l = m can present greater superradiance rates as compared to the fundamental overtone mode. Analytically this is apparent for l = m = 4 (see right panel of Fig. 4) where it has also been shown to occur for l = m = 3 considering numerical solutions [65].

D. Superradiant Evolution
Sequential to the formational phase of a BH, superradiant evolution can begin via quantum fluctuations in the vacuum where each of the quantised superradiant levels begin to grow exponentially with their corresponding superradiance rates. The fastest-growing level which satisfies the superradiance condition always dominates the initial superradiant evolution until it has extracted enough spin so that the superradiance condition is no longer satisfied. Once the scalar cloud has extracted the maximal spin for the dominant mode the system can be be considered as a (quasi)-stationary hairy BH for astrophysical purposes. The BH energy loss through mass reduction is minimal compared to the shift in angular momentum due to the extend of the scalar cloud. Once the growth of the dominant level has stopped the BH will spend a significant portion of its lifetime on a Regge trajectory (dashed lines in Fig. 5) separating higher mode instability bounds. This can be seen from the basic intuition that as the higher modes of the BH begin to spin down the BH perturbing it from the Regge trajectory the negative component of the eigenfrequency for the previous mode dominates the evolution, spinning up the BH. This process is apparent until a significant portion of the scalar density in the cloud is reduced from the previously dominant level. At this point the BH traverses the Regge plane towards the successive superradiant boundary, repeating the process until the timescales considered are to large for superradiance to occur.
If non-linearities are taken into account level mixing can increase the time spent on the superradiance condition boundary via perturbations of the gravitational potential around the BH. Dissipation of the scalar cloud can occur through processes such as the annihilation of axions into gravitons or unbound axions [11,20]. In general the scalar cloud becomes maximally occupied before annihilation processes begin in the non-relativistic limit. Further complications to the trajectory evolution of the BH could come from the bosenova phenomena, introducing intermediate stages comprising of bursts of GWs and phases spinning down the BH before the superradiance condition is finally saturated. Given the hierarchy of timescales between the superradiant instability and the GW emission from non-linearities when compared to the dynamical time scale of the BH it is possible to study the systems evolution in the quasi-adiabatic approximation for N ax fields [4,13,58]. The total scalar energy flux from the superradiance process through the horizon is, With a disregard for accretion the evolution of the system is described by the following equations where E S is the energy of the scalar cloud. The scalar cloud extracts mass and spin until reaching the saturation point. The final BH spin is, The final mass of the BH after the phase of superradiant evolution is defined by Eq. (20) were the variations in the defining BH parameters are related by, This defines the final mass of the BH: The true evolution of course is a complicated picture where non-linearities must be accounted for along with the properties of each system. In particular for SMBHs their mass are generally accumulated via accretion which requires very significant perturbations in order to match the evolutionary traits a stellar BH may follow for example in terms of traversing the mass-spin Regge plane.

E. The Regge Plane
A fundamental prediction stemming from superradiant instabilities of bosonic fields is the existence of exclusion regions in the BH Regge plane. Estimates of the instability time scale, τ SR partnered with reliable spin measurements for BHs, can be used to impose stringent constraints on the allowed masses of ultralight bosons. These bounds on the parameters of ultralight bosons follow from the requirement that in principle an astrophysical spinning BH should be stable over its lifetime. A superradiant instability time scale which acts faster than core processes such as accretion form observational thresholds on the expected regions of the two-dimensional mass-spin parameter space BHs should fall in. Following the process of superradiant evolution a large number of BH ob-servations should trace out the superradiance condition boundaries, mapping the Regge trajectories given the existence of as yet unidentified fields. For axions the shape of the gaps in the Regge plane are extremely sensitive to variations in the superradiant growth rate with the scalar mass. A BH therefore should be excluded from observational measurements given the existence of an ultralight boson if it's spin is measured above the relevant level curves for different orbital states of the quantised modes for the field. The bounds for bosonic fields with spin are wider than those for axion-like particles and so the potentially large systematic errors in BH spin measurements could act as a current restriction to this approach for spin-0 fields. The axion mass window which can be probed is fixed by the heaviest supermassive BHs with accurate recorded spin measurements along with a lower bound defined by the lightest measured stellar mass BHs.
The current lower and upper bounds on BH masses from X-ray spectroscopy and emission data covers the approximate region which defines the relevant axion mass window as, The isocontours defining the exclusion bounds are a function of the instability timescale and the boson mass. As the axion mass decreases the instability exclusion contours reduce in size. This corresponds to tighter instability regions which require larger spins for more massive BHs. Taking into account accretion and GW emissions can also slightly reduce the bounds in the Regge plane [58]. The timescales associated to the astrophysical pro-  29)) and τ Hubble (Yellow, Eq. (28)). The red/black data points denote mass and spin estimates of the stellar BHs from X-ray/BBH sources presented in Tabel I.
cesses of relevance alter when considering different compact object systems. For rapidly spinning BH candidates in X-ray binary systems or binary BH (BBH) mergers identified as detectable GW sources by LIGO more accurate constraints can be imposed when considering the typical timescales associated to a binary systems lifetime as other astrophysical processes such as accretion are sub-leading in this regard. A typical lower bound approximation for the lifetime of the binary system is given as for the most accurate constraints. The most conservative limits come from exclusion regions constructed using the Hubble time, τ H ∼ 10 10 yrs .
As opposed to stellar binary objects the relevant timescales for AGN in order for superradiance to maximally grow the scalar cloud for each quantised level come from accretion models. A statistical analysis of the exclusion limits over the whole BH mass region defined in Eq. (25) requires us to use a characteristic timescale derived from accretion considerations. The time scale for mass growth increases exponentially with an e-folding time given by a fraction 1/f Edd of the Salpeter time scale, where f Edd is the Eddington ratio for mass accretion. The accretion time scale is estimated using the Salpeter time for a BH radiating at it Eddington limit where σ T is the Thompson cross section and m P is the proton mass [66]. In order to model the accretion time the following parameters can be introduced [13] τ Salpeter = 4.5 × 10 8 yrs where η, the thin-disk radiative efficiency is a function of the spin related to a specific energy at the innermost stable circular orbit (ISCO). We select a typical value for the efficiency, η = 0.1 and the most conservative value of f Edd = 1 to model the effects of accretion. This fixes the superradiant instability timescale as τ SR = 45 Myrs. Increasing the bounds on f Edd allows for more optimistic models incorporating potential periods of super-Eddington accretion. A redefinition of f Edd holds the same equivalence as considering a subpopulation of degenerate mass fields (see Section. IV B) or considering different astrophysical processes to define the superradiance timescale. Such considerations are a limitation in the "logistics" of encapsulating the behaviour of FIG. 6. Isocontour exclusion bounds with calculated total exclusion probabilities in the BH mass-spin Regge plane from superradiant instabilities with a single axion field with mass, µax spanning the limits in Eq. (26). The shaded regions represent instability thresholds shorter than the time scale τ Salpeter in Eq. (29) for each value of the dominant orbital/azimuthal quantum numbers, l = m = 1 to 5. The blue data points are mass/spin estimates of stellar X-ray and BBH systems. The orange points correspond to mass/spin estimates of SMBHs from X-ray reflection spectroscopy. The exclusion probability function (black line) is calculated using the statistical model in Appendix B using the BHs compiled in Table I and is given as a function of the axion mass spanning both the stellar and supermassive regimes.
the total BH spectrum and as such we follow the most conservative limit defined above.
An individual treatment of the instability timescales derived from the properties of the accretion disc stability for each BH candidate can be used to tighten constraints of the field mass exclusions [67]. As the timescale limits for the superradiant instability are increased the limits for each mode, m will begin to saturate to the limits set by the boundaries of the superradiance condition. This effect is most prominent for higher order modes in the spin axis of the Regge plane allowing for enhancements in the potential to constrain ultralight bosons using observations of BHs with spins a moderate fraction of the extremal limit.
In Fig. 5 this is shown in the example exclusion window for a fixed axion mass of µ ax = 10 −11.5 eV in the stellar BH parameter space for each of the instability timescales in Eq. (27), Eq. (28) and Eq. (29). As the considered timescale increases the saturation of the mode bounds in the limit of the superradiance condition sees the greatest enhancement for l = m = 5. The red data points are the X-ray binary system BHs from Table I. The black data points are the primary and secondary sources involved in the BBH coalescence events (GW150914,GW151226 and GW170104) for several LIGO detections. Extremal BHs such as NGC 4051 impose constraints on each of the l = m = 1, 2 and 3 modes demonstrating the ability of well defined rapidly spinning BHs to constrain significant portions of the axion mass parameter space. An axion mass of µ ax ≈ 10 −11.5 eV is therefore tightly constrained by known X-ray binary sources as shown in both  Table I. The primary axis presents the Regge exclusion bounds for an instability time scale τ SR = 45 Myrs as a function of the axion mass, µ ax . The blue/orange data points are the stellar/SM BHs in Table I. The secondary axis displays the probability exclusion function formulated from the statistical model in Appendix B across the total BH mass range. The function "well" corresponds to the absence of any well defined IMBH candidates. Well defined mass and spin measurements for BHs covering the approximate region 10 2 M − 10 6 M could fill the currently inaccessible portion of the parameter space and probe interesting masses for axions associated to GUT and supersymmetric models in string/M-theory. The most promising realisation of detecting BHs in this space comes from the proposed space based gravitational wave observatories such as the Laser Interferometer Space Antenna (LISA) (see Fig. 2).

F. Black Hole Spin Measurements from Binary Systems and Active Galactic Nuclei
The identification of compact systems has seen a steady increase over the past decades with a number of X-ray binary sources and active galactic nuclei (AGN) now providing well defined measurements for the masses and spins of these systems. Currently the main sources of error for catalogued BHs comes from the systematic errors when modelling the emission of the accreting disc of the system. Both stellar BH and SMBH measurements come from analysing the X-ray spectrum of the accretion disk for identified compact sources. Assuming that General Relativity holds true as a valid description of the spacetime region outside the BH horizon and the ISCO of the accretion disk possesses a monotonic function potential then estimates on the spin of BHs can be made. In principle most BH candidates with well defined parameter estimates come from either thermal continuum fitting of the inner accretion disk or inner disk reflection modelling in order to determine the size of the ISCO. Further to this BH spin data has recently been collected via the observations made in several BBH mergers by LIGO [68][69][70]. Currently such observations contain large errors on both the mass and spin of the BHs when compared to existing X-ray binary system records. The resultant BHs formed from such astrophysical events cannot be included in considerations of constraining the masses of bosons given their timescale for observation is less than typical instability timescales by definition in the process of identification. Generally though future generation ground based detectors are still expected to produce large error measurements on BHs identified in this way and so impose a strong limitation on the accuracy of measurements used for constraints. Improvements in observatory sensitivity with space operated missions such as LISA [71] will open up the potential for a large catalogue of accurate BH measurements capable of probing a large potion of the cosmologically significant sector for axion-like fields. A large exclusion in the fully accessible space could also lead to tight constraints on how the axion population or sub-populations may be distributed when seeking realisations of desirable models in the context of cosmology.
We restrict ourselves to considering only BHs with detailed mass and spin errors. Each BH chosen for our analysis therefore has upper and lower bounds on both their mass and spin with well defined quoted uncertainties. In Table I we present all the stellar BHs and SMBHs used to constrain our axion distributions in Section IV along with their associated references. For a review of compiled stellar BH data see Refs. [72,73] and for SMBHs see Refs. [74,75].

III. THE AXION MASS SPECTRUM
The generic multi axion Lagrangian is: where θ i are the dimensionless axion fields, K ij is the kinetic matrix with mass dimension two, U is a general periodic instanton potential with charge matrix Q and phases, δ. Expanding the potential to the mass term only and diagonalising Eq. (31) can be reduced to the simple form: The spectrum of the model is given by the mass eigenvalues, {µ i }, which can be determined after expanding the instanton potential to quadratic order and obtaining a mass matrix, M ij . Diagonalising these matrices following the methodology in Appendix C 1 gives the mass eigenstates of a spectrum of fields. The canonically normalised dimensionful mass eigenstate fields, φ i , are defined in Eq. (C7) from the eigenvalues of the kinetic matrix. Adopting random matrix models for K ij and M ij it is possible to determine the distribution of {µ i } for various models. This process is reviewed in Appendix C and covered extensively in Ref. [30]. In the following we will consider just two simple models for the mass eigenvalues.
The first follows the celebrated Marčhenko-Pastur law for the eigenvalues of white Wishart matrices [45]. The spectrum is thought to describe Type-IIB string theory models of inflation with large numbers of axions [46]. The limiting distribution as the matrix size goes to infinity is given by on the compact interval where γ + and γ − are defined as, The expectation value of µ 2 ax isμ 2 ax and the shape parameter 0 < β M ≤ 1 determines the spread of the distribution, with large β M giving larger spreads as shown in the right panel of Fig. 7. The distribution for random realisations with finite N ax is shown in the left panel of Fig. 7, and is well fit by the limiting law.
Our second model for the mass eigenvalues follows from the "M-theory axiverse" [29]. In this case the mass eigenvalues follow an approximately log-normal distribution [30], as shown in Fig. 8. 1 In this example, the spread is controlled by the shaping parameter β M = N ax /N inst which takes values 0 < β M ≤ 1. Increasing β M leads to a larger mean and smaller variance. There are five parameters in the model of Ref. [30] in total, but here we use a two-parameter approximate fit: The mean of the log-normal distribution can be related to the expectation value of the 3-cycle volumes in the G 2 1 The naive expectation of log-flat eigenvalues turns out not to be realised for large numbers of fields after applying rotations to the canonical basis.
manifold, V X (see Eq. (C21)), and as in the above example, the variance, σ 2 , can be controlled by the number of instantons in the potential sum. The variance of the log-normal distribution is dimensionless, and so should take on some O(1) value. In Ref. [30] we typically found σ 1. Axion self-interactions can also play an important role in BH superradiance. In principle, by expanding the instanton potential to higher orders our RMT approach could lead to a distribution for the quartic interaction tensor: We are unaware of any study of the distribution of λ ijkl in RMT, and thus the treatment of interactions is beyond the scope of the present work. For sparse charge matrices the flavour changing, non-diagonal, entries in λ ijkl will be rare. If the attractive self-interactions are too strong then the superradiant cloud collapses via a bosenova before it can extract large amounts of spin from the BH. Superradiance can also be shut off by non-linear level mixing, or affected by axion emission due to annihilations [11,16,20]. The interaction tensor can used to calculate these rates, e.g. for axion emission via the φφφ → φ process. The level-mixing will be enhanced if the λ ijkl are non-diagonal and allow scattering of axions of different flavours. Decays from one flavour into another will have a similar effect of additional cooling of the cloud as the axion photon coupling considered in Ref. [11]. The Bosenova critical size, N Bosenova , could also become smaller in such a case due to the increased phase space for the scattering. How these and other non-linear effects compete with the basic increase of the BH superradiance rate and increased probability of mass outliers at large N ax is unclear.
Using the single instanton, dilute gas potential, V (φ) = µ 2 f 2 a [1−cos(φ/f a )] for a single field, it can be shown that the ratio of emission via the quartic interaction compared to graviton emission due to annihilations is given by [11]: The overall strength of the interactions, and their importance relative to gravity, is controlled by the axion decay constants, f a . The f a distributions for multiple fields derived from RMT can be computed (see e.g. Ref. [30]). Distributions with a high probability of small decay constants will have non-linearities dominated by self-interactions, while a for high probability of large decay constants the pure-gravity results can be used. Since we consider BH superradiance dominated by gravity, our results should be understood to apply strictly to distributions dominated by large f a . Taking the single- Limits determined using the SMBHs given in Table I. Right panel: Limits Determined using stellar mass BHs given in Table I. field results of Ref. [20] as a guide, this should be for f a 10 14−16 GeV. In the context of string models, our results should apply well to small volume compactifications [107], whereas self-interactions will play an important role in the Large Volume Scenario [33].

A. Single Field
In this short section, we begin the presentation of our results by computing single field limits to check our methodology is consistent with other results in the literature. Our statistical methods are described in Appendix B, and we calculate the exclusion probability, P ex (µ ax ).
Treating the stellar BHs and SMBHs as a single data set, our results for a single axion field with mass µ ax are shown in Fig. 6, superimposed on the Regge plane with the data. In this combined data set the exclusion probability remains finite over a range of intermediate axion masses due to the large mass errors on the lightest SMBHs. The absence of IMBHs means that the regions with P ex (µ ax ) > 0.68 ("1σ exclusion") do not overlap between the two datasets and they can be considered separately.
The exclusion probability for the stellar BH data set is shown in the right panel of Fig. 9. The high quality of these measurements, and the large number of them, leads to a smooth exclusion probability. At the 95% C.L. the stellar BHs exclude: The exclusion probability for the SMBH data set is shown in the left panel of Fig. 9. The data is of general poorer quality than the stellar data, with certain systems containing significantly large mass errors. It is also much sparser, with fewer SMBHs in the set. The sparseness of the data leads to oscillatory features in the exclusion probability, driven by the shape of the BH superradiance contours for each of the modes, with the exclusions being driven by individual BHs. This causes the probability of exclusion to oscillate between the 95% C.L when transitioning between certain BHs (faded red lines in the left panel of Fig. 9). The largest candidate, Fairall 9 drives the non-monotonic nature of the function at low axion masses. The large mass errors lead to non-zero exclusion probability extending to large axion masses. Taking the outer edge of the 95% C.L. region, the SMBHs exclude: Our exclusions for the stellar BH and SMBH datasets are consistent with the results of Refs. [20,67], after accounting for the differences in the data sets and methodology used. In particular comparing to Ref. [20] our choice to include BBH coalescence events with large masses when partnered with their large uncertainties push the constraints to incorporate lower masses, increasing the lower bound on the axion mass exclusion.

B. Degenerate Masses
We now begin to consider cases with multiple axion masses. The degenerate case is trivial to treat for any number of N ax axions with identical masses, µ ax . Since   29)) for an axion mass µax = 10 −12.75 eV. Large values of Nax effectively correspond to greater superradiance instability timescales considering a single field. Green data points are mass/spin estimates of X-ray binary stellar BH candidates. Blue data points are primary and secondary sources from BBH coalescence detections at LIGO.
the rate is additive in N ax we have: Therefore, setting τ BH Γ tot = 1 is equivalent to the single field case with the timescale rescaled as τ N = N ax τ BH . Thus, for the degenerate case the exclusion probabilities are trivial to compute for any N ax , and they will simply grow wider for increasing N ax corresponding to larger rates such as those shown in Fig. 5. In Fig. 11 we show the effect on the Regge plane with a degenerate population of axions with masses µ ax = 10 −12.75 eV. It is clear that an increase in N ax can lead to an exclusion on µ ax where there was not one in the single field case (purple limits). As the instability thresh-holds sweep through the Regge plane as a function of the axion mass, the wider instability limits possess the ability to "catch" lighter BHs in their exclusion bounds. We present the exclusion probabilities P ex (µ ax ) for various values of log 10 N ax for each regime in the left and right panels of Fig. 10. The contours in Fig. 11 always increase in the direction of smaller M BH , and so the constraints in Fig. 10 only broaden relative to the single field case for smaller axion masses. For SMBHs, where the higher harmonics play a role in the exclusion, the oscillations in the exclusion probability at high mass are also mildly affected. This is shown in the inset of the left panel of Fig. 10. Extremely large values of N ax quench the oscillations from the instability bounds of the higher order modes, saturating the upper bounds on the constraints.
The 95% excluded regions for µ ax for the degenerate case change by less than an order of magnitude compared to the single field case for N ax 10 5 . This shows that the increase in the superradiance rate for multiple fields (i.e. the rate sum in Eq. (2)) can be virtually neglected when computing the exclusion probability, even in the most extreme case of a very large number of degenerate superradiant fields.

C. Mass Distributions
There are two effects on BH superradiance constraints for mass distributions. The first is the effect of rate addition, the second is the effect of an overlap between the mass distribution and the exclusion probability. The results of the previous section show that even for the extreme case of degenerate masses, this effect is virtually negligible in the the exclusion probability for µ ax . Rate addition will be even more negligible for mass distributions with finite width, where off-resonant superradiance rates are exponentially suppressed. This leaves probability overlap as the dominant effect for mass distributions of finite width. With the effect of rate addition neglected, the exclusion probability for a mass distribution is trivial to construct from the exclusion probability for a single mass from the overlap integral. We use the probability that a model is allowed, since this trivially accounts for the combinatorics, and the excluded probability is in turn found trivially from this. Let P al (µ ax |N ax = 1) = 1 − P ex (µ ax |N ax = 1) be the probability that a given axion mass is allowed, assuming just one axion field. We then have that in a given model M with one axion, the probability that some parameters θ are allowed is P al (θ, N ax = 1|M) = dµ ax p(µ ax |θ, M)P al (µ ax |N ax = 1) , (43) where dµ ax p(µ ax |θ, M) is the probability distribution for µ ax in the model. The single axion allowed regions were evaluated numerically in Section IV A and the integral in Contours representing the 95% exclusion for Marčhenko-Pastur axion mass distributions as a function of the distribution shape, βM, and number of fields, Nax, for various distribution mean scales,μax. Regions above the contours are excluded. Large numbers of fields are constrained for a significant region of the probable axion mass space, with Nax ≥ 30 constrained for a wide range of βM over the considered scales.
Eq. (43) can be evaluated numerically given p(µ ax |θ, M). The above trivially generalises to the case of N ax fields: The exclusion probability for N ax fields is then given by P ex (θ, N ax |M) = 1 − P al (θ, N ax |M).

The Marčhenko-Pastur Distribution
The Marčhenko-Pastur (MP) distribution depends on two parameters: a mean mass,μ ax , and a shape parameter, β M . In order to probe the potential of a spectrum of fields scanning the Regge plane analogous to our single field constraints we highlight several interesting configurations. Consider the caseμ ax = 10 −13 eV, shown in Fig. 12, left panel. A single axion at this mass is excluded by the stellar BH data. However, for large spreads, i.e. β M → 1, the mode of the distribution moves to smaller values of the mass (shown in the right panel of Fig. 7), which are not constrained. Eventually at still larger β M the mode moves down to masses excluded by the SMBH data. In Fig. 13 we present the axion mass window open to superradiance as distribution mean scalesμ ax . Increasing N ax makes the exclusion probability grow, and for all β M there is a maximum N ax ≈ 20 above which the model is excluded at better than the 95% C.L for all β M (see Fig. 13). The maximum N ax allowed grows with β M . Now consider the caseμ ax = 10 −15 eV, shown in Fig. 12, right panel. In this case, the mean mass is in between the stellar and SMBH exclusions, and is allowed by the data. Thus, increasing β M now increases the exclusion probability. The non-zero exclusion probability at µ ax = 10 −15 eV coming from the SMBH data lowest mass points with large error causes the exclusion probability to grow as N ax increases even for small β M . Once again, there is a maximum N ax ≈ 50 above which the model is excluded at better than the 95% C.L for all β M (see Fig. 13). The maximum N ax allowed decreases with β M .
Motivated by the peak in the Calabi-Yau distribution along the self-mirror manifold line, Fig. 14 shows constraints onμ ax at fixed β M = 0.5. The excluded region has the same approximate shape as the single field exclusions for small N ax . As N ax → 1 the exclusion limits trace out the constraints for the single field case up to statistical fluctuations about the mean scale. The softer edges to the untouched regions when compared with the single field exclusion bounds come from the non-equidistant logarithmic spread of the mass spectrum about the mean scale when β M = 0.5. Reducing β M relaxes the limits to fully match the single field case in the low N ax limit. Increasing the number of fields, the model is excluded at better than the 95% C.L. for the range of mean masses shown for all N ax 100.

The M-theory Axiverse: the QCD axion, GUTs, and Fuzzy DM
The axion mass spectrum of the M-theory axiverse [29] was computed from RMT models in Ref. [30] and is well described by a log-normal distribution. The 95% excluded region in (σ, N ax ) for the log-normal distribution across a range of central values is shown in Fig. 1. We now derive BH superradiance constraints on three scenarios of interest realised approximately from this simple model for the M-theory axiverse. The M-theory axiverse with GUT scale unification predicts the existence of an axion with which arises from fixing a single modulus to give the correct GUT scale coupling, α GUT = 1/25 arising from a 3-cycle with volume V X = 25 in string units (see Appendix C 3). We model this by fixing the log-normal mean to log 10μax = −15. The fuzzy DM model [24][25][26][27][28] posits that DM composed of axions with mass has certain desirable properties that could lead to its being favoured over standard cold DM by observations of galactic structure. We model this by fixing the lognormal mean to log 10μax = −22. The QCD axion [21-   Each panel corresponds to three unique scenarios which determine the mean of the mass spectrum required to maximise the probability of drawing the desired masses detailed in Section IV C 2. In the limit σ 1 the total probability for Nax = 1 → ∞ converges to zero as the spread crosses the bounds probable by BH spin measurements. The behaviour in the limit σ 1 is determined by the accuracy of the available BH mass/spin measurements.

23] mass is given by
In order to realise the QCD axion in M-theory, some light eigenstate in the "pure M-theory" spectrum should receive its mass dominantly from QCD instantons. Furthermore, the VEV of this field should be not far displaced from θ = 0 to solve the strong-CP problem. These two conditions together require that there is at least one eigenstate in the pure M-theory spectrum with [29]: We model this by fixingμ ax and σ such that µ ax,low is within 95% of the probability at the lower end of the distribution after N ax draws. This fixesμ ax (σ, N ax ) in terms of standard error functions: With the above fixed, one linear combination of axions receives its mass from QCD instantons. Therefore, we remove one axion from the M-theory distribution and replace it with the QCD axion. The probability that the QCD axion in M-theory is allowed based on BH superradiance data is thus: P al (σ, N ax ) = P al (µ ax,QCD |N ax = 1) dµ ax p[µ ax |σ,μ ax (σ, N ax )]P al (µ ax |N ax = 1) Constraints on the distribution parameters of each of these benchmark models are shown in Fig. 15. While none of these models are ruled out for a single axion, in all cases the exclusion probability starts to become significant for non-zero distribution widths and large numbers of fields. In all cases, the maximum allowed value of N ax increases for very large σ. For large σ the distribution is approximately log-flat with respect to the data exclusions, and increasing the width simply reduces the probability of overlap.
The GUT model has a small exclusion probability at zero width due to the large mass errors on the lightest SMBHs (NGC 4051 and MCG-6-30-15 ). The GUT model is excluded at better than the 95% C.L. for all widths σ < O(100) for N ax 100. The fuzzy DM model is excluded at better than the 95% C.L. for all widths 1 σ 10 3 if N ax 100.
The QCD axion model is the least constrained by the data. The mass of the QCD axion with f a 10 17 GeV is not excluded itself by BH superradiance, nor is the light mass µ ax,low required from the M-theory part of the spectrum. There is a small range of intermediate widths where the distribution does overlap the excluded region, excluding 0.2 σ 4 if N ax 100 at 95% C.L. while σ 4 is allowed for all N ax < 1000 considered.

Comment on Fuzzy DM and BH Superradiance
Recently it has been claimed that the global 21cm signal [108], which is strong evidence that the Universe was undergoing reionization at redshift z re ≈ 17, places a lower bound on the fuzzy DM mass of µ ax ≥ 5-8 × 10 −21 eV [109,110]. This result is extremely interesting since, if it is to be believed in its accuracy, it significantly shrinks the gap between fuzzy DM bounds from BH superradiance and structure formation. In the context of the present work, if fuzzy DM is realised from a mass distribution, then respecting the reionization bound and BH superradiance demands an extremely narrow distribution with a small number of light fields. If the gap between fuzzy DM constraints from BH superradiance and reionization is closed, either by the measurement of spins of the most massive SMBHs, or improvements on the lower limit to z re , then fuzzy DM with no self-interactions will be completely excluded. Rescuing fuzzy DM from BH superradiance constraints in such a case would require self interaction strengths corresponding to decay constants f a 10 16 GeV. Low decay constants open the door to new fuzzy DM phenomenology [111][112][113], but may become increasingly hard to realise in small-volume string compactifications.

V. DISCUSSION AND CONCLUSIONS
BH superradiance places strong constraints on the possible existence of light bosonic fields with small selfinteractions, in particular on axion-like fields. Many authors have considered these constraints for the case of a single new light field. The excluded ranges of axion mass are: A model with multiple axions is excluded if just one field lies in these ranges. We have studied this possibility, and used BH superradiance to exclude certain distributions of axion masses. The constraints become more severe with larger numbers of axion-like fields due to the increased probability of drawing an outlier. This allows us to place constraints on the number of axion-like fields, N ax .
Models for axions coming from string theory and Mtheory typically involve many axion-like fields. These fields have their masses determined by microscopic quantities related to the geometry of the compact space. Their masses, however, are expected to follow particular statistical distributions independently of the microscopic details. We have considered various different distributions, log-flat, log-normal, and Marčhenko-Pastur, using BH superradiance to bound both the parameters of the distribution, and, more significantly, the number of light axions within that distribution.
Constraints on N ax from a process such as BH superradiance, which relies only on the existence of the vacuum fluctuations of the given field, are extremely powerful, and could be used in this context to bound the dimensionality of phenomenologically consistent moduli spaces in string/M-theory. Indeed we have seen that the benchmark value of N ax ≈ 30 found in the majority of known Calabi-Yau manifolds can be excluded for a wide range of distribution parameters. Only a small number of fields should obtain masses anywhere in the BH superradiance region from 10 −10 eV µ ax 10 −20 eV, which can be accommodated with a single very wide distribution σ 30, or bimodal distributions containing only very light or relatively heavy axions. Our analysis has neglected axion self-interactions, which shut off BH superradiance if they are strong, and other constraints, for example coming from the relic abundance. It would be interesting in this regard to combine our previous analysis in Ref. [30] with the current analysis and compute, in addition to axion masses, the axion decay constants, relic density, and self-interaction potential. The present work is more model-independent, since it does not rely on any cosmological assumptions, and applies to any model for light scalars with sufficiently small self-interactions. The extended and combined analysis will be the subject of future work.

ACKNOWLEDGMENTS
We acknowledge useful conversations with Bobby Acharya, Katy Clough and Chakrit Pongkitivanichkul. The work of MJS is supported by funding from the UK Science and Technology Facilities Council (STFC). DJEM is supported by the Alexander von Humboldt Foundation and the German Federal Ministry of Education and Research. The 3+1 dimensional spacetime region outside the horizon of a rotating Kerr BH is described by the invariant line element, ds 2 = g αβ dx α dx β which, using the standard Boyer-Lindquist coordinates (t, r, θ, φ) and metric signature [−, +, +, +], takes the form which is invariant under time translations, possessing a Killing vector. The metric functions are defined as, The zero solutions of Eq. (A3) define two horizons, an inner Cauchy horizon at r − with the larger root at r + defining the outer physical event horizon. The characteristic limits of each BH horizon as a function of the dimensionless spin are displayed in the panels of Fig. 16.
As the spin of the BH approaches the extremal limit, a * = 1 the inner and outer horizons coincide. A defining property of Kerr BHs is existence of an a surface external to the outer horizon known as the ergosurface. The ergosurface is defined by the static limit roots, g tt = 0 with the coordinates, As the BH spin approaches the static Schwarzschild solution, a * → 0 the ergosurface and outer horizon coincide. The region between the outer horizon and ergosurface defines the ergoregion. Inside the ergoregion the vector, ξ µ in the time coordinate basis becomes spacelike, ξ µ ξ ν g µν = g tt > 0. This property allows for a Killing energy in the presence of a BH to be negative inside the ergoregion, leading to the superradiant amplification of the infalling waves associated to the bosonic field. The event horizon angular velocity for observers at spacial infinity for the BH is, The dynamics of the linearised massive scalar in the Kerr spacetime are governed by the Klein-Gordon wave equation.

The Klein-Gordon Wave Equation
The classical massive scalar field obeys the Klein-Gordon wave equation, The massive Klein-Gordon equation on a Kerr spacetime background allows for a separation of variables with an infinite discrete set of complex eigenfrequencies ω lmn , of the form in Eq. (7). The Klein-Gordon wave equation following a separation of variables is expressed by two coupled ordinary differential equations. Using the Teukolsky formulism [114] the separated ODEs for the radial and angular parts, R lm (r) and S lm (θ) respectively are,   ∆∂ r (∂ r R) + (ω 2 (r 2 + a 2 ) 2 − 4ar g rmω + a 2 m 2 − ∆(µ ax r 2 a 2 ω 2 + l(l + 1))R(r) = 0 .
The first ODE in Eq. (A9) determines the angular component, S lm (θ) of the scalar eigenfunction. The angular solutions of Eq. (A9), S lm are the the spheroidal harmonics which are required to be regular at the pole boundaries, θ = 0 and θ = π. These boundary conditions single out a discrete family {K lm } of angular eigenvalues also known as the coupling constant which characterise the massive scalar. The angular eigenvalues can either be found using an expansion in the limit that aω and aµ ax → 0 where the expansion of K lm , gives Λ lm → l(l + 1) + O(a 2 ω 2 ) in the non rotating limit where, when k = 0, an analytical expression can be extracted. Higher orders of k require numerical solutions. The function inside the sum defines the so called spheroidicity. These angular eigenvalues can also be found via Leavers' continued fraction method (Appendix A 4) or Hughes' spectral decomposition method. A rescaling of the radial function introducing, ψ lm = √ r 2 + a 2 R lm , along with a definition of the Regge-Wheeler tortoise coordinate where, allows for the radial Teukolsky equation (Eq. (A10)) to be expressed in the form of a Schrödinger like wave equation, The effective potential is defined as: We require solutions to Eq. (A10) with boundary conditions defining an outgoing solution tending to zero at spacial infinity and purely incoming waves at the event horizon. In the limit where ωm 1 and µ ax m 1, Eq. (A10) is susceptible to analytic methods. These boundary conditions correspond to modifications of the radial solutions in the limits, where k + ≡ ω − mΩ H . In the low energy limit, ωM BH 1, the radial equation is amenable to the method of matched asymptotics.
It has been shown that analytic solutions for small values of α can be found using approximate solutions at large and small radii in terms of hypergeometric functions, where matching techniques are used at an intermediate radius to obtain the superradiance rates to leading order in α [59]. In this limit analytical methods utilise the fact that the radial mode functions, R lm (r) can be approximated in asymptotic regimes by known analytical functions. For each region the equations can be reduced to the form of a confluent hypergeometric function.
Regions far from the BH outer horizon adhering to r r g whilst ensuring we are in the µ ax M BH 1 regime allow the ODE in Eq. (A10) to be approximated as where the axion momentum is given by, The solutions can be extracted by defining: where the value of δν represents a small complex number which describes the deviation away from the pure hydrogenic spectrum. When the axion momentum satisfies k 2 > 0 we are presented with a series of quasi-bound state solutions. This equation is the same form of the Schrödinger equation which governs the electron in the hydrogen atom. The solution to Eq. (A19) can be expressed as where U is the confluent hypergeometric function of the second kind. In the limit r r g , Eq. (A10) is solved analytically where the approximate solution takes the form The form of the equation in Eq. (A23) presents a solution infalling at the horizon

Numerical Solutions for Bound States
Analytical based methods for approximating the superradiance rates suggest a maximal value for the approximate regime µ ax M BH ∼ 1. In order to probe this region of the parameter space it is required to solve the radial mode function ODEs eigenvalue problem using numerical techniques. See [63,64] for an initial study incorporating Leaver's continued fraction method [115] for numerical calculations and Dolan's work [61] for an extensive study of the expanded parameter space, providing numerical solutions using a three-term recurrence relation and the continued fraction method.
The radial function R(r) is assumed to take the follow-ing form of the the infinite series A substitution of Eq. (A27) into Eq. (A10) obtains the three term relation for the expansion coefficients a n for n > 0, n ∈ N α 0 a 1 + β 0 a 0 = 0 , (A31) α n a n+1 +β n a n + γ n a n−1 = 0 , where, The values of the constants c1, c2, c3 and c 4 are expressed as functions dependant on the parameters, ω, σ, m and angular eigenvalues, Λ lm (Eq. (A11)), where (A41) The three factor recurrence relation can be solved in terms of a continued fraction if we take the assumption that the factor an+1 /an → 0 as n → ∞ obtaining, Rearranging Eq. (A31) to give and substituting in n = 0 into Eq. (A42) gives the condition for the eigenvalue equation for bound state eigenfrequencies of the form in Eq. (7), which can be solved using numerical method techniques.

Appendix B: Statistical Model
We model the BH data in Table I with two dimensional multivariate gaussian distributions for both x = M BH and y = a * . There are N d data points d i comprising the dataset {d i }. For each point in the data set the values of M BH and a * and their associated errors become centred data values (x,ȳ) with errors (σ x , σ y ). We are interested in the probability that a given model, M, is excluded given the data, {d i }: P ex (M|{d i }). Since a single data point in the disallowed region would exclude the model, P ex (M|{d i }) is given by the probability that any single data point is above the BH superradiance isocontour boundaries for each value of l. For a large number of data points, this is a relatively tricky combinatorial problem. However, the probability is normalised such that: Now we can use the binomial theorem (or a simple probability tree) to note that P allowed (M|{d i }) is simply the cumulative probability that all data points simultaneously fluctuate below the isoctontour: and P allowed (M|d i ) is simply the volume of the bivariate Gaussian contained outside the isocontour boundary given by the function y = f (x). To evaluate P allowed (M|d i ) in a numerically efficient manner, we make two simplifying assumptions. Firstly, we assume zero covariance between x and y. Secondly, the error on the two-dimensional data can be evaluated using an effective one dimensional error [116,117]. These two simplifications allow us to use the standard error function to evaluate P allowed (M|d i ), rather than the more numerically expensive integral under the curve.
The shape of the BH superradiance contours y = f (x), which only have support over finite x, requires this procedure to be evaluated in two separate regimes. Where the contour is defined, we use the contour as y = f (x) and evaluate the effective one dimensional error in y, Σ y , as: When the contour is not defined for a given x, we instead use the inverse function x = g(y) and evaluate the effective error in x, Σ x , as: The effective errors are represented visually in Fig. 17.
Since our functions are all given numerically, the inverse function and its derivative are trivial to evaluate given the original function. A complication arises since g(y) is multivalued, taking two values g 1 and g 2 for a single y. We choose to evaluate the derivative g (ȳ) at the nearest part of the contour (i.e. the value g i which minimisesx − g(ȳ)), and evaluate the error function between the two values g 1 and g 2 . This approximation only affects P allowed (M|d i ) for values close to unity, while P ex (M|{d i }) is dominated by the smallest values of P allowed (M|d i ) contained well within the contours where f (x) has support and is single valued.
The use of the effective errors, Eqs. (B3, B4), assumes that, for a given data point, the functions f (x) and g(y) are smooth at the mean value over the range of the errors. When a BH data point with large errors sits close to a cusp in the contours the exclusion probability computed from the effective error is smaller than the true answer. Cusps in the total contour are caused by the meeting of individual contours with different l values, each of which are smooth. A more exact procedure would thus be to compute the probability individually for each l contour, and then compute the cumulative probability from a product over l. This would increase the number of likelihood evaluations by l max × N d , and for speed of computation we do not perform this more accurate calculation. The more accurate calculation would give larger exclusion probabilities (reducing the overall effective size of BH errors), and so the approximate computation is more conservative in the sense that it does not give overly strong exclusions.
The most general form for the multi-axion action for fields below any compactification, moduli stabilisation or PQ symmetry scales is of the form given in Eq. (31). We set the axion field alignment used to determine the diameter of the fundamental domain to only include P = N , where we always possess sufficient instanton contributions N for each axion field, P . We restrict our considerations to non-perturbative terms with trivial charges, Q j,i = 1 N . We therefore only need consider enhancements to the N ax field space diameter defined as the longest distance between vertices in the polytopes defining the field ranges via the pythagorean sum from the N-flation model and kinetic alignment in our models considered in Section C 2. Lattice alignment as well as alignment theories possessing P ≥ N are beyond the scope of this work. See Refs. [50,[118][119][120] for extensive details of axion field alignment.
We begin in the lattice basis, with an axion defined by a single cosine potential possessing a shift symmetry obeying, θ i → θ i + 2π. We diagonalise and canonically normalise the axion field space metric K ij moving to the kinetic basis with the unitary rotation U ij where, (C1) We define the axion decay constants, f a , from the eigenvalues of K ij in the lattice basis in Planck units, We can now define the canonically normalised field as, In the kinetic basis the effective Lagrangian takes the form, where the new mass matrix is defined as, Moving to the mass eigenstate basis with a further unitary rotation, V ij gives, In this basis the mass eigenstate fields are defined as, with the effective Lagrangian,

The Random Matrix Theory Mass Spectrum
A systematic construction of the axion decay constant and mass spectrum in explicit realisations of the string axiverse is an extremely complex and numerically comprehensive task to undertake. In general, considerations need to be made for leading instanton corrections to the superpotential (Eq. (C10)), a calculation of the full scalar potential along with a minimisation of polynomial expressions with potentially many variables when considering realistic numbers of apparent axions or moduli. See Ref. [121] for a detailed discussion of the complexities of the string landscape. The effective field theory approach in Section C 1 can benefit from the simplistic nature of RMT inspired models on the grounds of universality. In these models universality dictates that the distribution of the physical dimensional parameters are characterised by some mean scale and variance. For generic field space metric considerations the kinetic matrix in Eq. (31) can be well described by a matrix belonging a class of matrices of the Wishart form, This formalism can also be extended in the small field approximation to govern the properties of the axion mass matrix. The mean scale of the axion population defines the phenomenological properties of the fields. The spread of the spectrum is controlled by a shaping index β K,M ∈ (0, 1] which has been shown to have theoretical foundations relating to the total dimension of the moduli space [30,46]. Below we present a series of RMT inspired models based on charge quantisation and field space alignment considerations. The defining features of each of our models we consider are presented in Table II. Each model presents a modest hierarchy in complexity regarding the initial basis and determined mass eigenstate spectrum. Our first model is based on N-flation type field alignment [46] acting as our simplest strawman model. The fields in canonical coordinates are defined by the effective field space metric, . The canonical field ranges in this basis are defined as φ i =f a θ i wheref a is a scaling factor introduced to represent the degenerate decay constant scales arising from the diagonal kinetic matrix. The field space diameter is given by the pythagorean sum over the N-Dimensional hyper-rectangle (see Table II). The resulting mass spectra is displayed in Fig. 7.
When beginning in the lattice basis a spectrum of decay constants now scale the initial field sampling. It has been shown the kinetic matrix K ij , could belong to the Gaussian orthogonal Wishart ensemble with i.i.d gaus-  Table II.   TABLE II. RMT models considered in this work and extensively covered in Ref. [30]. Detailed are the initial basis considerations for each effective model along with the relevant sampling procedures and field space diameters. The values of σ represent the approximate spread a population or subpopulation of axions would have in each model. √ N ax f a,max . √ N ax f a,max . sian entries [50,122]. The entries for the sub-matrices Y ij in Eq. (C9) composing each of the kinetic and mass matrices are drawn from normal distributions with zero mean and unit variance. The axion decay constant spectrum is given by the limiting Marčhenko-Pastur law (Fig. 7). In the mass eigenstate basis the mass spectrum is now rotated by the non trivial rotations between the lattice and kinetic basis. Universality dictates a convergent mass spectrum well modelled by a log-normal distribution with its limited variance, σ defined by the bounded nature of the initial Wishart structure. The mass spectrum in this model is presented in Fig. 18a. If the entries of the sub-matrices are not selected as i.i.d gaussian entries and instead selected from a log-uniform distribution, the resulting matrix of the Wishart form will now reside in a class of rank one spiked Wishart matrices. In the original models [123][124][125] the matrices are defined by the class, W R (Σ, M ) where a single element of the covariance matrix deviates from unity inducing a phase transition in the distribution for the largest eigenvalues. In the limit N ax → ∞ a bulk region forms supported by the Marčhenko-Pastur limiting law and one singular eigenvalue is repulsed from the bulk as shown in Fig. 18b. The approximate order of the singular eigenvalue is λ ax ≈ O(N ax ), which defines the enhancement of the diameter of field space by a factor of √ N ax due to the large hierarchy between f a,max and the second largest eigenvalue. This behaviour governs the axion decay constant spectrum in the model (up to canonical normalisation factors) as shown in Fig. 18c. In the limit β M = 1 the mass spectra converges to the white Wishart case. For values of β M < 1 the convergent mass spectrum is well modelled by a log-normal distribution plus two positively and negatively logarithmically repulsed regions enhancing the total spectral width.

The M-Theory Mass Spectrum
It has been shown in M-theory compactified on G 2 manifolds with an absence of fluxes it is possible to stabilise both the moduli and axions in order to realise a spectrum of ultra-light axions in the low energy spectrum of its four-dimensional effective supergravity theory. We follow the explicit realisation of the string axiverse in [29,126,127]. In such models the moduli are stabilised in a non-supersymmetric minima, with all axions pairing up with geometric moduli where all moduli superfields possess PQ symmetries. The superpotential in this model takes the form A k e ib k F k , (C10) with order O(1) constants, A k . The first two terms in Eq. (C10) come from strong gauge dynamics in the hidden sector using up one combination of axions where φ 1 is a holomorphic composite field made of hidden sector matter fields. In general the rest of the fields present in the summation come from non-perturbative physics such as membrane instantons and serve as a fundamental feature of such compactification models. We only need to consider the higher order correctional terms assumed to be generated from membrane instantons where b k = 2πI and I ∈ Z with the gauge kinetic functions, It is always possible to find realistic arguments determining the number of non-perturbative effects as larger than the number of axions, N Inst > N ax giving rise to sufficient independent terms in the superpotential. To consider a spectrum of axions we integrate out the moduli and heavy axion combinations such that the relevant effective superpotential becomes whereΛ i are the associated mass scales for each nonperturbative effect. The potential now takes the following form, In Ref. [30] it was shown an expansion of the periodic potentials to quadratic order reveals the axion mass matrix, where N j i = b i N j i is a rectangular matrix of size (N ax , N ), C r = which in terms of sub-matrix structure following the philosophy of Eq. (C9) gives, This defines the following form for the sample submatrix, where i, j = 1, . . . , N ax and r = 1, . . . , N , A ir is a rectangular matrix of size (N ax , N ) with a normalisation factor 1 /N. In order to define our mass scales of interest in the Mtheory axiverse consider the general form for the superpotential in Eq. (C11). F i represents the gauge kinetic functions which are linear combinations of the moduli superfields, The generalised volume of the corresponding 3-cycles is calculated from, The geometric moduli are stabilised in terms of a single parameter V X which represents the stabilised volume of the three-cycle supporting the hidden sector. In order to realise a GUT in the low energy limit of the theory, at least one of the gauge kinetic functions must give rise to the expected value of the GUT coupling constant, The average value of V X therefore fixes the mass scales of the spectrum of axions appearing in the visible sector (Fig. 8). We parameterise the axion mass distribution in terms of the average value of the three-cycle volume distribution, V X via the relationship which contains the parameters we statistically sample to determine the nature of σ used throughout the basis of this work (see Section III D of Ref. [30] for details of the parameters used). The number of parameters and hierarchy of scales involved in the statistical sampling of the three-cycle volume dictate a large spread in the mass eigenstates covering many decades. The nature of uni-versality ensures a convergence to a normal distribution over these scales [30].