astroparticle physics – black hole physics – elementary particles

1 Introduction↩︎

This paper will introduce a simple yet rigorous statistical framework for constraining ultralight bosons (ULBs) using black hole (BH) superradiance (SR) from observational BH spin estimates. The framework introduces no new physics, using only previously developed models of BHSR. Our approach simplifies the constraint setting procedure to use only a timescales analysis, rather than simulating BH populations and evolution. We argue that this simplified approach is appropriate given the inherent astrophysical uncertainties/systematics, and is user-friendly since the analysis is easy to repeat with modified timescale assumptions. What we develop is a rigorous Bayesian framework to set limits on ULBs that allows to simply and directly use the posterior probabilities on BH model parameters published by observers, and advocate for future publication of BH results in this format. The robust approach to data products will allow for transparent and easy compilation of future large BH datasets. Future improved understandings to the physics of BHSR should be incorporated via modifications to the characteristic timescales involved, and we demonstrate this by comparing the bosenova model to the equilibrium model for ULB self-interactions.

[1]–[3] showed that BHs with large spins can probe the existence of new bosonic degrees of freedom via BHSR. Due to a link between the BH mass and the boson mass, the resulting constraints are relevant for ULBs. Indeed, [3]–[8] placed constraints on ULBs with masses \(\mu\lesssim \SI{e-10}{\eV}\) using BHSR from BH spin measurements.

Ultralight bosons are of particular interest in particle astrophysics at present. The quantum chromodynamics (QCD) axion [9]–[12], an excellent dark matter (DM) candidate [13]–[17], is among this class of particles. There is a far-ranging experimental programme underway [18]–[23] to search for QCD axions or axion-like particles (ALPs), which can also be DM candidates [24]. As we will review later, BHSR in the stellar mass range places constraints on the lower end of the allowed \(\mu\) range of the QCD axion. The ‘fuzzy’ DM model is also in the ULB class, and it is of interest as an alternative to cold DM with distinctive phenomenology on sub-galactic scales [25]–[29], and as the DM candidate with the absolute lowest allowed mass [30], [31]. A combination of cosmological probes and galactic astrophysics probe the fuzzy DM mass from below, while BHSR in the supermassive mass range constrains it from above [18], [32]. String/M-theory also predicts the existence of ULBs covering a wide mass range [1], [33]–[36], and it has recently become possible to use BHSR to rule out certain extra dimensional geometries, thus placing observational constraints on the string landscape [8].

While BH spin estimates are a well-established probe of ULBs, there are a number of issues and subtleties related to the derivation of the constraints just discussed. Extracting the BH mass \(M\) and (dimensionless) spin \(a_*\) from BH data requires an involved analysis [37]–[40], making this a highly specialised task where data sets need to be dealt with on a case-by-case basis. The resulting \((M,a_*)\) estimates can be correlated and, since the strongest constraints come from BHs with close-to-maximal \(|a_*|\), a standard Gaussian approximation may be inadequate. Moreover, the theoretical computation of BHSR-related effects is challenging and depends on the BH evolution. These complications require some compromises between a statistically rigorous analysis and the practical challenges: we wish to include uncertainties while – at the same time – leverage as much information as possible from the data, BH evolution, and BH population parameters.

We argue that a Bayesian framework admits all of these desiderata. We demonstrate the feasibility of our approach based on the X-ray binary BH M33 X-7 and the supermassive BH (SMBH) IRAS 09149-6206. These two sources were selected due to their well-understood nature and the availability of reliable \((M,a_*)\) estimates. Ultimately, however, it is important to recognise that any BH mass and spin measurements – regardless of the BH system or measurement technique utilised – may be used for our method.

We use the \((M,a_*)\) distributions from the two BHs above to obtain constraints on ULB parameters. In particular, these are the first constraints from IRAS 09149-6206. We compare the ULB limits arising from the full Bayesian treatment to other approximate statistical methods to explain the differences which arose in past works, and summarise strengths and weaknesses of each approach. While not the focus of this work, we also discuss the role of different computations of the BHSR rate – both with and without the inclusion of self-interactions.

At this point, let us emphasise that BH spin measurements are not the only way to search for ULBs using BHSR. Moreover, BH spins do not, individually, offer a route to discovery, except for in the case of population statistics [2], [41]. Routes to discovery via BHSR can be pursued via dynamical means from either gravitational waves [42]–[45], binary orbits [46]–[49], or direct signatures of the ULB cloud [50]. These methods require significantly more involved modelling, and thus more inbuilt assumptions compared to constraints from BH spin timescales, but the potential for discovery warrants the additional detailed treatment. Any alternative searches for ULBs through BHSR are complementary to the constraints we set here. In particular, our constraints would add additional information to the fitted ULB parameters in case of a discovery, given that the ULB’s existence has to be consistent with all available data.

The goal of the present work is thus to provide rigorous and transparent BHSR constraints from observational BH spin estimates, which can be used with confidence by the community. Explaining the differences between past works contributes to understanding their underlying assumptions and statistical methodology, and hopefully encourages more groups analysing BH data to make their posterior samples or full analysis pipelines available. To facilitate this process, our code and the associated, redistributed BH data are publicly available [51].

2 Black Hole Evolution including BHSR↩︎

Superradiance (SR) arises due to an instability of the wave equation for a massive scalar field \(\phi\) with mass \(\mu\), \[\Box \phi - \mu^2\phi = 0\, , \label{eq:klein-gordon}\tag{1}\] where we initially neglect self-interactions (introduced in 2.3). The D’Alembertian is computed from the Kerr metric [52], for which we use the ‘mostly positive’ signature. The SR instability arises due to the presence of growing modes for \(\phi\), which can be identified as imaginary eigenvalues of the associated Schrödinger equation. The energy required to create the ‘boson cloud’ is extracted from the mass and spin of the BH. The evolution of the mass and spin are treated quasi-statically, i.e.there is no backreaction on the Kerr spacetime and one simply evolves the parameters in the Klein-Gordon equation 1 .³

The SR rate \(\Gamma\) depends on \(\mu\) and is approximately maximised when the dimensionless ‘gravitational coupling’ of the atom-like, bosonic states around a spinning BH, \[\alpha = G M \mu= \num{0.75} \left(\frac{M}{10\,\mathrm{M}_\odot}\right) \left(\frac{\mu}{\SI{e-11}{\eV}}\right) \, ,\label{eq:alpha}\tag{2}\] is \(\alpha \sim 1\), i.e.it is approximately maximised when the boson Compton wavelength is of order the gravitational radius. Here, \(G\) is Newton’s constant. While we sometimes include the Planck mass \(\mathrm{M}_\text{P}= 1/\sqrt{G} = \SI{1.22e19}{\GeV}\) as a reference scale, we will generally set \(c = G= 1\) in what follows.

Superradiance causes the spin \(a_*\)(and mass \(M\)) of a BH to evolve in time, \[\dot{a}_* = - \GammaM_\text{cloud}/ \omega^R\approx - \GammaN_\text{cloud}\, ,\] where the dimensionless BH spin is \(a_*\equiv a/M = J/M^2\), \(\omega^R\) is the ULB energy (the real part of the ULB frequency), and \(M_\text{cloud}\) and \(N_\text{cloud}\) respectively are the total mass and occupation number of the boson cloud forming around the BH. The latter evolves as \[\dot{N}_\text{cloud} = \GammaN_\text{cloud}+ P_\text{GW}/\omega^R\, , \label{eq:could95growth}\tag{3}\] where \(P_\text{GW}\) is the small power radiated in gravitational waves (GW) from the cloud. We neglect \(P_\text{GW}\) in this work and refer the reader to [54] and [41], [55] for further details on the phenomenology of GW emission from boson clouds.

In the non-relativistic regime, the boson cloud can be thought of as levels in hydrogen atom-like states. The complete system is described by a set of coupled rate equations for the occupation numbers of each level, which grow due to BHSR, along with interactions between the levels due to scattering, and losses due to emissions to infinity and flux across the horizon [7]. An astrophysical BH will also be governed by other processes, having some typical timescale \(\tau_\text{BH}\), which we discuss further in 3.1.

2.1 Computing the free BHSR rates↩︎

Starting from [56], the computation of the BHSR rates has seen a number of advances and improvements. In particular the analytical approximations of [57] provide a straightforward way to compute BHSR rates.⁴ Semi-analytical results, e.g.by [60], allowed a further refinement of the (semi-)analytical computations [46], [61], [62].

It is illuminating to discuss the approximate analytical result. We use the metric of an uncharged, rotating BH with angular momentum \(a\) and mass \(M\) [52], and the coordinate system of [63]. The outer horizon of a Kerr BH is \(r_+/M \equiv 1 + (1 - a_*^2)^{1/2}\).

The bosonic state around the BH can be labelled by ‘quantum numbers’ \(\left|{nlm}\right\rangle\), where \(n \geq 2\) is the principal quantum number,⁵ and has the complex-valued frequency \(\omega_{nlm}= \omega^R_{nlm}+ \mathrm{i}\, \omega^I_{nlm}\). The real part \(\omega^R_{nlm}\) corresponds to the ULB energy and the imaginary part \(\omega^I_{nlm}\), if positive, sets the growth rate of the field amplitude. The BHSR rate is defined as the growth of the particle number and is thus \(\Gamma_{nlm}=2\omega^I_{nlm}\). The latter is the case, and the \(\left|{nlm}\right\rangle\) level is superradiant, if the SR condition is fulfilled, \[\alpha / l \leq 1/2 \, . \label{eq:sr95cond}\tag{4}\] [64], [65] derived general bounds on \(\mu\) for superradiant states, which are also valid for levels with \(n \gg 1\).

2.1.1 Non-relativistic approximation↩︎

To first order in the non-relativistic approximation (NRA), we have [2] \[\begin{align} \omega^R_{nlm}&= \mu\left[1 - \frac{\alpha^2}{2n^2} + \mathcal{O}(\alpha^4)\right] \tag{5} \, , \\ \quad \omega^I_{nlm}&\approx \mu\left(ma_*- 2 r_+\mu\right) \, C_{nl} \, \Pi_{lm} \, \alpha^{4l+4}\tag{6} \\ \text{with} \quad C_{nl} &= \frac{ 2^{4l+2} \, (n+l)! }{ n^{2l+4} \, (n-l-1)! } \left[ \frac{l!}{(2l)! (2l+1)!} \right]^2 \\ \text{and} \quad \Pi_{lm} &= \prod_{k=1}^l \left[ k^2 \, \mathlarger(1 - a_*^2\mathlarger) + \mathlarger(ma_*- 2 r_+\mu\mathlarger)^2 \right] \, . \end{align}\] Terms at higher orders, which also introduce a dependence on the \(l\) and \(m\) quantum numbers, are e.g.discussed by [61].

2.1.2 Higher-order corrections and continued fraction method↩︎

We can further refine the NRA result from 2.1.1. For instance, the Python package SuperRad [55] uses the qnm package [66] to solve the radial Teukolsky equation using Leaver’s continued fraction method [67] in order to compute the relativistic BHSR rates. The numerical solution in the relativistic regime is valid at large \(\alpha\) and is matched consistently onto the Newtonian approximation at small \(\alpha\), found from solving for hydrogen atom-like bound states and including terms up to \(\mathcal{O}(\alpha^5)\) [61]. SuperRad is fast and also has good control over the numerical accuracy of the methods, with a relative error on the SR timescale estimated to be better than 1%. It also includes gravitational wave emission.

However, the limitation of SuperRad is that it includes only modes with \(l = m = 1, 2\), does not include self-interactions, and will not return rates at the highest \(a_*\) values. We thus use our own implementation of the CFM, following [60].

1 shows the results of our implementation for the \(\left|211\right\rangle\) rate \(\Gamma_{211}\)(orange line), which compares well with the result of SuperRad(dotted, black line). To facilitate the numerical root-finding for the complex-valued root \(\omega_{211}\), we also compute the next-to-leading-order (NLO) corrections [62] as a starting point, which we show as a dashed, red line for comparison. The computed rates differ at most by an order of magnitude, which, however, should be kept in mind when comparing different results in the literature. This applies in particular to the NRA (dashed-dotted, blue line).

Since the present work was first drafted, a study by [68] appeared, which addresses the relativistic corrections in more detail. For the dominant \(\left|211\right\rangle\) rate, the corrections are shown to be of order 20%. It would be possible to update the rates used in our public code to include these corrections, but at the present accuracy we believe it would have a small effect. This conclusion is validated by the agreement between the limits we set with our Bayesian simplified timescale comparison and the frequentist full trajectory analysis of [68]. ⁶

2.2 Evolution of the free field↩︎

Due to 3 , the SR instability causes the occupation number to grow exponentially, \(N_\text{cloud}\propto \mathrm{e}^{\Gamma_{{nlm}} t}\), where \(\Gamma_{{nlm}}\) denotes the “free” SR rate \(\omega^I_{nlm}\) from 6 . Eventually, however, the SR condition in 4 ceases to be valid and the cloud growth stops – even in the absence of self-interactions. By this time, the cloud will have extracted some spin \(\Delta a_*\) from the BH, which can be related to the cloud occupation number \(N_{\Delta}\) at the time when an amount of spin \(\Delta a_*\) has been extracted, assuming negligible initial occupation number of the cloud [3], \[N_{\Delta}\approx \num{e76} \left(\frac{1}{m}\right) \left(\frac{\Deltaa_*}{0.1}\right) \left(\frac{M}{10\,\mathrm{M}_\odot}\right)^2 . \label{eq:nmax}\tag{7}\] 7 allows us to estimate whether or not a given cloud occupation level can be significantly reduce the observable spin of a BH. In particular, if we consider the exponential growth of the cloud for the \(\left|nlm\right\rangle\) level, some BH timescale \(\tau_\text{BH}\)(further discussed in 3.1), and the SR timescale \(\tau_\text{SR}= \Gamma_{{nlm}}^{-1}\), we can only obtain constraints on ULB properties if

• \(\left|nlm\right\rangle\) is superradiant according to 4 ,

• SR is fast enough to spin down the BH, \(\tau_\text{SR}< \tau_\text{BH}/\ln(N_{\Delta})\).

Note that the typical value of \(\Deltaa_*\) adopted in 7 is of order the size of the observational errors on \(a_*\). Since the condition on \(\tau_\text{SR}\) only logarithmically depends on \(N_{\Delta}\), and thus \(\Deltaa_*\), even an order of magnitude difference in \(N_{\Delta}\) would not significantly change the derived limits. Still, since the choice of \(\Deltaa_*\) is somewhat arbitrary, the BH evolution should be treated in more detail. This treatment goes beyond the present analysis, where we only compare timescales.

2.3 Self-interactions↩︎

The effect of boson self-interactions has been considered in detail for the case of the dimensionless coupling \(\lambda\) from a \(\phi^4\) terms in the Lagrangian, \[\lambda \equiv \frac{\mu^2}{4!} f^{-2} \, ,\label{eq:self95coupling}\tag{8}\] where \(f\) is the ULB decay constant, which is typically connected to the energy scale of new physics leading to the ULB’s existence (such as symmetry breaking or extra dimensions), while \(f^{-1}\) typically sets the scale of the ULB’s couplings to Standard Model particles. For a generic pseudo-Goldstone boson, 8 arises from Taylor expanding a generic potential with respect to the field periodicity \(f\).

For a given value of \(\mu\), experiments and astrophysical searches for ULBs typically place upper limits on the couplings or \(f^{-1}\), while BHSR places lower limits on \(f^{-1}\), making the methodologies highly complementary. In the case of fuzzy DM, the relic abundance scales with \(f^2\) [16]–[18], with future cosmological probes reaching \(f\sim \SI{e16}{\GeV}\) [69]–[71], which overlaps with the BHSR sensitivity.

For QCD axions, the term in 8 arises from a Taylor expansion of the scalar potential in chiral perturbation theory, and [72] found that \(\lambda = \num{-0.346(22)}\). This prediction is possible since QCD axion models link \(\mu\) and \(f\), \[\mu= \frac{\chi_\text{top}^2}{f} = \SI{5.691(51)}{\nano\eV} \left(\frac{\SI{e15}{\GeV}}{f}\right) \, , \label{eq:qcd95axion95mass}\tag{9}\] where the QCD topological susceptibility \(\chi_\text{top}\) has been computed at NNLO [72], [73]. Note that these relations may not hold for a generic ALP. As we now discuss, self-interactions affect the SR rates.

2.3.1 The equilibrium regime↩︎

Initially, the ULBs will be near the vacuum and the rate \(\Gamma_{{nlm}}\) is the one for free bosons, as computed in 2.1. If the SR rate is faster than the inverse of the relevant BH timescale \(\tau_\text{BH}^{-1}\), the boson cloud can grow (see 2.2). Its subsequent evolution is determined by the strength of ULB self-interactions, for which there are two possibilities.

The first possibility is that the cloud enters an equilibrium between BHSR losses and level transitions caused by self-interactions and gravity [7]. The BHSR rate in the equilibrium regime is effectively reduced compared to the vacuum rate by a factor of \(\eta_\text{eq} < 1\) for the equilibrium number density of bosons in the most superradiant level. In particular, the spin-down timescale in the equilibrium regime is given by \[\tau_\text{eq}= \frac{1}{\eta_\text{eq} \Gamma_{{nlm}}} \, , \label{eq:tau95eq}\tag{10}\] where the factor \(\eta_\text{eq}\) can be computed following the procedure of [7], who consider (transitions between) the \(\left|211\right\rangle\) and \(\left|322\right\rangle\) levels, \[\begin{align} \eta_\text{eq} &= \frac{2\sqrt{3}}{3} \frac{\sqrt{\Gamma_{211}\Gamma_{322\times322}^{211\times\infty}}}{\Gamma_{211\times211}^{322\times\text{BH}}} \, , \label{eq:mod95neq} \\ \Gamma_{322\times322}^{211\times\infty} &\approx \num{1.1e-8} \alpha^8 \muf^{-4} \, , \\ \Gamma_{211\times211}^{322\times\text{BH}} &\approx \num{4.3e-7} \Big(1 + \sqrt{1 - a_*^2}\Big) \, \alpha^{11} \muf^{-4} \, , \end{align}\tag{11}\] where \(\alpha\) is assumed to be sufficiently small and where \(f\) is the ULB decay constant in 8 .⁷

Given values for \((M,a_*)\), we may exclude all ULB models for which, in addition to the conditions in 2.2,⁸ the reduced SR rate in equilibrium is still faster than the typical BH timescale. This can be encoded in the condition \(\tau_\text{eq}< \tau_\text{BH}\).

2.3.2 Bosenovae↩︎

The other possibility for the evolution is that a bosenova occurs [2], [3], [74]–[76]. This happens when the occupancy of the cloud reaches a critical value \(N_\text{crit}\), given by \[N_\text{crit}\approx \num{e78} \, c_0 \, \left(\frac{M}{10\,\mathrm{M}_\odot}\right)^2 \left(\frac{f}{\mathrm{M}_\text{P}}\right)^2 \, \frac{n^4}{\alpha^3} \, , \label{eq:nbose}\tag{12}\] where \(c_0 \approx 5\) [3]. Such an expression can be derived by comparing the order of magnitude size of different terms in the action [8]. Thus, in a single superradiant timescale \(\tau_\text{SR}\) that would otherwise lead to a maximum occupancy \(N_{\Delta}\), as given by 7 , a number of bosenovae \(N_\text{crit}/N_{\Delta}\) occur. Since in each cycle the maximum occupancy is \(N_\text{crit}\), the overall spin-down rate is reduced. This can be interpreted as an increase in the relevant timescale [3]: \[\tau_\text{BN}= \frac{1}{\Gamma_{{nlm}}} \frac{N_{\Delta}}{N_\text{crit}} \ln(N_\text{crit}) \, . \label{eq:tau95bn}\tag{13}\]

Given values for \((M,a_*)\), we may exclude all ULB models for which, in addition to the conditions in 2.2, the reduced SR rate from successive bosenovae is fast enough to reduce the spin. This can be written as \(\tau_\text{BN}< \tau_\text{BH}\).

2.3.3 Bosenova or Equilibrium?↩︎

In both the bosenova and equilibrium scenarios, bosonic self-interactions slow down BHSR. In the case of a bosenova, the cloud is prevented from reaching a high occupancy and exponential growth over the entire SR timescale and instead goes through many short periods where a smaller amount of spin is extracted. Stronger self-interactions in this case decrease the value of \(N_\text{crit}\) and make each cycle shorter, thus extracting less spin overall. In the equilibrium case, stronger self-interactions lead to more rapid boson scattering in superradiant levels, consequently reducing their equilibrium number density, and thus their ability to extract spin from the BH.

[7] approximated the full superradiant evolution by including the most relevant levels and found that an equilibrium is typically reached before a bosenova occurs. This finding was also confirmed by [77], who included relativistic effects. More concretely, for \(\alpha \lesssim 0.2\), [7] found that a two-level description is sufficient to model the BH evolution, the inclusion of higher levels does not affect their results, and a bosenova does not occur. For \(\alpha \lesssim 0.3\), they found an equilibrium for the three-level system \(\left|211\right\rangle \leftrightarrow \left|322\right\rangle \leftrightarrow \left|411\right\rangle\) and argued that a bosenova is also not likely in this case. For \(\alpha \gtrsim 0.3\), higher SR levels become important and solving the related system of rate equations becomes increasingly challenging. Since the focus of the present work is to improve the statistical framework, we leave such computations for future work.

In summary, the bosenova scenario appears to be excluded for \(\alpha \lesssim 0.2\). In this regime, we may use the modified rates based on the results from [7]. For larger \(\alpha\), the modified equilibrium and bosenova rates become less robust, introducing a systematic uncertainty from the assumptions underlying the related computations. In order to compare the possible impact of the two scenarios on the ULB constraints, our results in 4 will use the SR rates as computed in the previous sections for all values of \(\alpha\). We stress again that the rates computed in the previous sections are only robust for values of \(\alpha \lesssim 0.2\) and that the bosenova scenario is strongly disfavoured in that regime.

In contrast to the above, [75] performed simulations of the full Klein-Gordon equation in 3+1 dimensions on a fixed Kerr spacetime with a cosine, i.e.axion-like potential, for the scalar field. This method differs from the rate equation methods, and treats axion self-interactions beyond the quartic approximation. [75] provided evidence that a bosenova occurs, rather than saturation. [78] performed full numerical relativity simulations including back-reaction for binary black holes with a complex self-interacting scalar field, and also found evidence of a bosenova.

Finally, as already noted, BHSR for scalars has not been studied in fully non-linear \((3+1)\)-dimensional numerical relativity due to the challenging timescales involved. Given this fact, and the presence of historic disagreement on the bosenova issue, we leave open the possibility of bosenova and equilibrium and present results for both in our analysis. Comparing both bosenova and equilibrium models demonstrates how our method is easily adaptable to changes in physics assumptions and allows easy comparison of the effect of different theoretical models on ULB constraints.

2.4 Regge trajectories↩︎

As visualised by [3], the BH evolution will be dominated by the fastest SR rate for which the SR condition in 4 is met. This is typically the \(\left|211\right\rangle\) level but, as the BH evolves and its mass and spin decrease, the SR condition will eventually cease to be fulfilled. Once this happens, ULBs from higher levels can still drain the BH’s mass and spin – as long as they meet the SR condition and the SR rate is fast enough to make the change observable. Eventually, there will be no such superradiant levels left, and the existence of ULBs makes a striking prediction: any given BH should not have a spin higher than a certain critical value \(a_*^\text{crit}\), which depends on \(M\), \(\tau_\text{BH}\), \(\mu\), and \(f\). The line \(a_*^\text{crit}(M,\tau_\text{BH},\mu,f)\) is referred to as ‘Regge trajectory’.

If we knew the true values of \((M,a_*,\tau_\text{BH})\), we could exclude the existence of any ULBs whose parameters predict \(a_*^\text{crit} < a_*\). Black holes on or below Regge trajectories might simply have been ‘born’ with a lower spin, and this lack of knowledge about the BH’s history means that they are compatible with the existence of the corresponding ULBs. In reality, we can only estimate the BH parameters, and thus only produce a statistical estimate for the exclusion probability of a given ULB (see 4).

3 Black hole data↩︎

We select two BH data sets as examples to demonstrate the methodology proposed in this work. One BH is the high-mass eclipsing X-ray binary (XRB) M33 X-7, which has frequently been studied in the SR literature. The other one is the SMBH IRAS 09149-6206, which possesses well-defined, independent mass and spin posteriors. It is also the first time that this SMBH has been used to infer ULB constraints. The spin of the XRB was inferred via the thermal fit continuum method, while the spin of the SMBH was inferred from X-ray reflection spectroscopy. We discuss these two techniques in [sec:xray_binaries] [sec:smbh], respectively, and refer the reader to [39] and [40] for detailed, recent reviews.

3.1 BH timescales↩︎

Black holes can acquire angular momentum (be ‘spun up’) due to accretion. The characteristic timescale for a BH to acquire a given mass through accretion at the Eddington limit is given by the Salpeter time, \(\tau_\text{Salp}\approx \SI{4.5e7}{\yr}\), assuming a canonical radiative efficiency of 10%. The Salpeter time is also the appropriate characteristic timescale for BH spin-up from accretion, and thus often one simply chooses \(\tau_\text{BH}= \tau_\text{Salp}\). Although super-Eddington accretion is possible, this approach is conservative. X-ray binaries accrete at sub-Eddington luminosities even in outburst [79], with only the rarer sub-population of ultraluminous X-ray sources [80], [81] producing super-Eddington luminosities.

Indeed, the high spin measurements in XRBs present a puzzle: it is impossible for the XRB to acquire near-maximal spin via accretion over its lifetime, thus requiring high ‘natal’ spin [82]. For active galactic nuclei (AGNs), the vast majority of radiatively efficient systems accrete with Eddington fractions of \(\lambda_\text{Edd} \equiv L_\text{bol} / L_\text{Edd} \sim \numrange{0.01}{1}\) [83], [84], where \(L_\text{bol}\) and \(L_\text{Edd}\) are the bolometric and Eddington luminosities, respectively. Indeed, super-Eddington accretion is generally found to be rare in both observational [85], [86] and theoretical/simulation [87]–[89] studies, especially at low redshift [90]–[93].

An approximate current Eddington fraction can be calculated from the estimated bolometric luminosities and BH masses for each object, with reports of \(\lambda_\text{Edd} \approx 0.1\) for M33 X-7[82] and \(\lambda_\text{Edd} \approx 0.2\) for IRAS 09149-6206[94]. However, accreting BHs are variable sources, so to accurately estimate \(\tau_\text{BH}\) for individual sources would require detailed knowledge of the accretion history over tens of  . For the SMBH case (IRAS 09149-6206), we assume here that \(\tau_\text{BH}= \tau_\text{Salp}/ \bar{\lambda}_\text{Edd}\), where \(\bar{\lambda}_\text{Edd}\) can be thought of as an appropriately weighted time-averaged Eddington ratio. We set \(\bar{\lambda}_\text{Edd} = 0.1\), which is likely to be fairly conservative, but note the inherent astrophysical uncertainty associated with this value. For the stellar mass BH M33 X-7, we instead use the system age (see section 3.2.3). If \(\tau_\text{SR}= \Gamma^{-1} < \tau_\text{BH}/\ln N_\Delta\), then SR dominates the evolution of the BH and \(a_*\) will decrease. Observing a BH in a stable configuration for which we can determine \(a_*\) and \(\tau_\text{BH}\) then allows us to exclude ULBs models that would predict \(\tau_\text{SR}< \tau_\text{BH}\).

For completeness, note that the longest relevant timescale is the age of the observable Universe, which is of order the Hubble time \(\tau_\text{H} = 1/H_0 = \SI{14.5}{\giga\yr}\) [95]. Using \(\tau_\text{H}\) would result in the strongest possible limits, and it is thus a useful timescale for investigating the addressable ULB parameter space. While the ‘age’ of a number of SMBHs – or, perhaps more accurately, the time since they accumulated most of their mass – could be of order \(\tau_\text{H}\), as also suggested by studies with the James Webb Space Telescope [96], [97], the value of \(\tau_\text{BH}\) is rather set by accretion and other processes in the BH evolution, as discussed above.

3.2 X-ray binary BHs and M33 X-7↩︎

X-ray binaries consist of a stellar mass BH accreting matter from a companion star. XRBs are classified as low-mass (LMXB) or high-mass (HMXB), depending on whether the companion is lower or higher mass than the BH, respectively. The accretion process emits electromagnetic radiation, which can be detected using X-ray telescopes. Analysing the data from such measurements yields estimates for \((M,a_*)\) with, in some cases, relatively high statistical precision but potentially large systematic uncertainties due to the modelling of the system.

3.2.1 Mass estimates in stellar-mass BHs↩︎

As reviewed by [38], the most robust method to estimate BH masses in XRBs obtains a dynamical mass measurement, using Kepler’s laws of motion and the presence of a stellar companion. This method is more direct and robust than the methods used for the majority of BHs because the dynamical impact of the BH on its companion can be measured. From spectroscopy, radial velocities of the emission or absorption lines from the companion star can be obtained, which allows a measurement of the orbital period \(P_\text{orb}\) and the radial velocity semi-amplitude \(K_c\). The so-called ‘mass function’ \(F(M)\) can then be defined, which explicitly relates the two radial velocity parameters and \(M\), \[F(M) = \frac{K_c^3 P_\text{orb}}{2\pi G} = \frac{M \sin^3(\theta_i)}{(1+q)^2} \, ,\] where \(\theta_i\) is the binary inclination and \(q\) is the mass ratio. For known inclination and \(q\), one can estimate \(M\), whereas for unknown inclination \(F(M)\) provides a firm lower limit on the mass (since \(1 + q > 1\) and \(\sin \theta_i \leq 1\)).

Dynamical mass measurements are easier to obtain for LMXBs, because the systems typically go through outburst cycles and have long periods of very low accretion activity (quiescence). During quiescent periods, optical spectroscopy can be used to accurately characterise the radial velocity curve of the companion star, as has been carried out for about 20 LMXBs [98], [99]. In HMXBs, the accretion state phenomenology is more complex; of the over 100 known HMXBs, many of which host neutron stars, only a handful of dynamical mass measurements have been possible for BH systems. These include famous sources such as Cygnus X-1, and the source studied here, M33 X-7. In HMXBs hosting BHs, the X-ray source is persistent and true quiescent periods do not occur. However, as HMXB companion stars are, by definition, rather massive, they outshine the accretion disc in the optical band, making spectroscopy of the donor star possible. There are still various systematic effects that can complicate the dynamical BH mass measurements, the main three factors being line formation in (or contamination by) the strong stellar winds of the companion, uncertainties in the donor star mass, and unknown Roche lobe filling factors [38]. Fortunately, M33 X-7 is an eclipsing system, which mitigates some of the effects mentioned above. In addition, it would be possible in future to incorporate uncertainties on, e.g., the donor mass, Roche lobe filling factor, and inclination as additional nuisance parameters into the analysis.

3.2.2 XRB Spin Estimates from the continuum fitting method↩︎

In addition to the X-ray reflection method (see 3.3.2), the continuum-fitting method has widely been used to estimate BH spins in XRBs. This relies on the influence of BH spin on the inner edge of the accretion disc, as higher (prograde) spins lead to an innermost stable circular orbit (ISCO) lying closer to the BH. The closer the ISCO radius (\(R_\text{ISCO}\)) is to the event horizon [40], the more energy will be released per unit accreted matter. The maximum effective temperature of the accretion disc corresponds to that at the ISCO and is set by \(T_\text{eff,max}^4 \propto \lambda_\text{Edd} M R_\text{ISCO}^{-3}\) in the thin-disc model. For a canonical 10% radiative efficiency, \(T_\text{eff,max}\) peaks in the soft X-ray band (i.e.below \(\SI{2}{\keV}\)) for an accreting stellar BH of \(M = 10\,\mathrm{M}_\odot\). In contrast, \(T_\text{eff,max}\) peaks in the extreme UV band for an accreting SMBH with \(M = 10^6\,\mathrm{M}_\odot\) for the same canonical radiative efficiency. The interpretation of observed X-ray spectra of XRBs accreting at moderate rates relative to the Eddington limit (\(\lambda_\text{Edd} = 0.01-0.3\)), paired with suitable accretion disc models, has provided spin constraints for around 10 XRBs. We refer to [40], [100], [101] for relevant reviews.

Generally, the accretion disc model assumed in continuum-fitting spin studies is the geometrically-thin, optically-thick, steady-state accretion disc model due to [102]. This model is an extension to the also well-known Shakura–Sunyaev alpha thin-disc model [103] to the relativistic regime. In the Novikov–Thorne model, mass loss due to winds in the inner disc is assumed to be negligible, and, as in other thin-disc models, the disc acquires non-molecular viscosity due to internal stresses and heat is dissipated locally. These two conditions have fairly sound theoretical and empirical motivation, given that thin-disc models can describe the putative disc component of the observed X-ray spectrum in XRBs during soft, thermally dominated accretion states. Furthermore, fully relativistic continuum-fitting models generally adopt a zero-torque boundary condition at the ISCO.

Continuum-fitting models predict local thermal spectra by combining the predicted thermal disc spectrum – subject to an accretion disc model – with the temperature-dependent colour correction factor. Although the temperature dependence of the latter can be inferred from sophisticated radiative transfer calculations [104], in practise, some continuum-fitting models implement an approximate temperature-dependent colour correction factor, such as that described by [105].

Spin estimates in XRBs from continuum fitting have been achieved with high statistical precision. However, the systematic errors emerging from the assumptions and approximations made by these models must be accounted for. First, as mentioned above, the thin-disc solution is only expected to be appropriate for thermally-dominated discs for bolometric luminosities \(L_\text{bol} \lesssim 0.3 \, L_\text{Edd}\) [106]. Second, even in this moderate luminosity regime, radiation could be emitted from within the ISCO, which is expected to induce an uncertainty in the recovered spin from continuum fitting of \(\sigma_{a_*} = 0.05\) [107]. Third, if the aforementioned zero-torque boundary condition at the ISCO is not satisfied, e.g.due to magnetic fields threading the inner disc, further systematic errors of \(\sigma_{a_*} \leq 0.1\) would be induced [108], [109].

Moreover, astrophysical parameters with inherent uncertainties – such as the inclination to the angular momentum vector of the binary, the BH mass, and the distance and accretion geometry of the XRB system – are also expected to contribute to the overall systematic uncertainty on the inferred BH spin. These nuances have led to variations in the analysis pipelines employed by different research groups amongst different data sets, potentially leading to further systematics from modelling and different statistical methods. Nevertheless, a consistent analysis of multiple XRBs would in principle be possible, thanks to publicly available tools such as the Xspec spectral-fitting package [110].

3.2.3 M33 X-7 posteriors↩︎

In this work, we use the spin posterior of the XRB system M33 X-7 of [82], inferred using a fully relativistic continuum-fitting model. The analysis in [82] was applied to a set of Chandra and XMM-Newton observations of the binary in the thermally-dominated state and applied to epochs where \(L_\text{bol} < 0.3 \, L_\text{Edd}\). Based on previously determined system inclination and BH mass values, [82] constrained the BH spin value to \(a_*= \num{0.84(5)}\). The authors noted that the assumptions made by their fully relativistic continuum-fitting model, added to the assumption that the angular momentum vector of the binary is approximately aligned with the BH spin vector, are expected to induce additional systematic uncertainties.

In 2, we show the parameter sampling distributions from bootstrapping (black dots) for the BH system M33 X-7, which we obtained by digitising Fig. 3 of [82] [111]. Clearly, the \(2\sigma\) contour of an uncorrelated Gaussian (dashed, orange line) does not provide a good description of the distribution, and including correlations (blue line) is more appropriate. We also show the ‘box’ defined by the uncorrelated \(2\sigma\) errors, which captures all possible underlying correlations at the price of discarding the information contained therein. We will discuss the relative importance of these effects on the ULB limits in 5.1.

Note that, in other cases, the resulting \((M,a_*)\) distribution may be highly non-Gaussian. For instance, the somewhat uniformly distributed posterior distribution of GRS 1124-683 [112] provides an extreme example.

2 also shows the Regge trajectory (solid, grey line and shading; cf.2.4) for the non-interacting \(\left|211\right\rangle\) level and an ULB mass of \(\mu= \SI{0.15e-12}{\eV}\). Since the estimated age \(\tau_\text{age}= \SI{3e6}{\yr}\) of M33 X-7[113] is below \(0.1\,\tau_\text{Salp}\), we choose \(\tau_\text{BH}= \tau_\text{age}\). To illustrate the effect of self-interactions, in this case assuming the equilibrium regime and \(f^{-1}= \SI{4e-16}{\GeV^{-1}}\), we also include the modified Regge trajectory as a dotted, grey line. We can exclude ULBs if most of the \((M,a_*)\) samples can be found above the Regge trajectory.

3.2.4 Effects of the companion star↩︎

The companion of M33 X-7 has a mass of \(M_2 \approx 70\,\mathrm{M}_\odot\) [114] and an orbital period of \(T = \SI{3.45}{\day}\) [115]. Using the equations from Appendix I of [7], we verified that the companion does not have a significant effect on the exclusion regions for the ULB parameters. ⁹ For other systems, however, the companion might play a more important role. In such cases, our proposed statistical framework can include and propagate the related uncertainties on the companion parameters, which would be particularly relevant for the large uncertainty of about 10% on \(M_2\) [114] (see 6 for other such nuisance parameters).

3.3 Supermassive BHs and IRAS 09149-6206↩︎

Supermassive black holes are believed to reside at the centres of all galaxies, and a compelling way to study them observationally is by probing the activity of their associated active galactic nuclei (AGN). AGNs are nuclear regions located at the centre of about \(\numrange{1}{10}\%\) of galaxies. Their integrated luminosity exceeds the stellar luminosity of their host galaxies as a consequence of the accretion onto the SMBH from the surrounding accretion disc. Given the high-quality observations needed, obtaining spin and mass estimates for SMBHs is rather challenging. Here, we summarise how relativistic X-ray reflection and infrared interferometry can respectively be used to constrain the spin and mass of SMBHs in AGNs.

3.3.1 SMBH mass estimates↩︎

The mass of Sgr A*, the closest SMBH to Earth, can be measured dynamically from the motion of stars within the BH gravitational potential [116], [117]. However, in more distant AGNs, the BH’s sphere of influence cannot be resolved and its mass is often estimated from broad emission lines: under the assumption that the broad-line region (BLR) gas is virialised, the width of the line, \(\Delta v\) is related to the BH mass through the relation \(M = \mathfrak{f}\, R_\text{BLR} (\Delta v)^2 / G\), where \(\mathfrak{f}\) is the so-called virial factor, \(\mathfrak{f}\sim \mathcal{O}(1)\), that accounts for BLR geometry and projection effects. Using this virial estimator requires an estimate for the line formation radius \(R_\text{BLR}\), which ideally would be estimated from reverberation lags [118], [119]. If \(R_\text{BLR}\) is not known, then it is usually estimated from the luminosity-size relation, \(R_\text{BLR} \propto L^{1/2}\) – where \(L\) is an appropriately chosen proxy for the bolometric or ionising luminosity – calibrated with the low-redshift empirical relation from reverberation mapped AGNs [120]. There are numerous other ways of estimating BH mass, which we cannot describe exhaustively, each with their own characteristic uncertainties and systematics; notable methods include the use of scaling relations between \(M\) and the stellar velocity dispersion [121], [122] or bulge mass [123], and the observation of megamasers [124].

Recently, the GRAVITY interferometric beam-combiner [125] on the Very Large Telescope has made it possible to estimate \(R_\text{BLR}\) and model the BLR geometry from the differential phases in infra-red interferometry. These measurements are extremely useful for BH masses and have allowed more accurate estimates for a number of luminous AGNs [126], [127]. Additionally, while subject to various assumptions, the modelling provides a well-defined posterior probability distribution for \(M\) for our analysis. One such object observed is IRAS 09149-6206, which is our chosen SMBH case study (see 3.3.3). Our use of the mass posterior from GRAVITY observations mitigates some, but not all, of the systematic uncertainties associated with BH mass estimates.

3.3.2 SMBH Spin estimates from X-ray reflection spectroscopy↩︎

In this work, we adopt the BH spin posterior distribution for IRAS 09149-6206 inferred by [128], who used the X-ray reflection method. In the absence of Compton-thick absorption and of spectral degeneracies induced by warm absorbers and ultra-fast outflows, the detailed modelling of reflected features in X-ray spectrum of AGNs provides a powerful pathway to probing the strong gravity regime in their central SMBHs. We note that several alternative observational techniques have been used to generally probe spins of BHs accreting at lower rates or hosting strong radio-jets [129]–[133]. We refer to [134] for a discussion on BHSR from the spin estimate of M87* from the interpretation of Event Horizon Telescope data; and to [135] for joint constraints on the ultralight boson mass and spin distribution at the time of binary BH formation from the Gravitational-Wave Transient Catalog (GWTC-2).

The disc reflection spectrum arises from the reprocessing of the primary X-ray continuum from the inner and outer regions of the ionised accretion disc, as well as from distant neutral and non-relativistic material surrounding the accretion disc. The primary continuum is thought to originate from the Compton up-scattering of seed photons in the innermost regions of the disc by the electrons in the hot X-ray corona. A fraction of the direct coronal emission towards the observer will back-scatter on the disc and be re-radiated away, with the non-thermal emission of coronal electrons governing the high-energy cutoff tail of the broadband spectral energy distribution (usually in the 100 keV–300 keV range) in type-1 AGNs.

Amongst other features in the X-ray reflection spectrum, the most prominent one is the K\(\alpha\) line, as iron is the most abundant element in the disc and has a high fluorescent yield [136]. Generally, the K\(\alpha\) line receives contributions from reprocessing from the inner regions of the disc, as well as cold, non-relativistic matter surrounding the immediate vicinity of the AGN. A narrow K\(\alpha\) reflection feature (centred at 6.4 keV in the rest frame) often arises due to the reflection from cold matter surrounding the outer accretion disc. However, the most important feature for BH spin inference is the broad K\(\alpha\) feature, centred around 6.4 keV–6.97 keV rest-frame energy, which likely indicates the presence of reprocessed emission from the innermost regions of the accretion disc.

Given that this reprocessing is thought to take place in close proximity to the BH event horizon, the reflected emission is subject to standard Doppler and general relativity effects in the disc, with the latter including light bending and gravitational redshift. The spectral imprints emanating from these effects, especially a distinct red wing to the line, will be more dramatic and noticeable the closer the emission takes place to the event horizon [137]. In addition, relativistic reflection from the innermost regions of the disc will also induce a forest of reflection features down to the soft X-ray band (below 2 keV). Finally, another characteristic feature of the disc reflection spectrum in AGNs is the Compton (reflection) hump, which arises in the range of 20 keV–40 keV due to photo-electric absorption of low-energy photons and the multiple reprocessing of high-energy (coronal) electrons [138].

Current relativistic X-ray reflection models can provide a powerful probe of the BH spin if the inner accretion disc extends down to \(R_\text{ISCO}\). Under this condition and assuming an accretion disc model, these models can then be used to make predictions of reflected spectra parameterised by \(a_*\), given assumptions about:

• the geometry of the corona,

• the disc density at the mid-plane and/or disc ionisation state,

• and the metric of spacetime around the central BH.

State-of-the-art X-ray reflection models [139]–[141] have provided spin constraints for about 20 stellar-mass BHs and 50 SMBHs [40]. The simplifications made by these reflection models, and their possible translation into systematic biases on the inferred BH spin were thoroughly reviewed by [40] and [39]. We also note that, although the majority of current spin estimates drawn from X-ray reflection spectroscopy rely on the detection and interpretation of the broadened K\(\alpha\) line, several SMBH spins in AGNs have been inferred via the reflection interpretation of the soft excess – notably, in Narrow-line Seyfert 1 galaxies [142], [143].

3.3.3 IRAS 09149-6206posterior samples↩︎

In this work, we make use of spin posteriors from a specific analysis of IRAS 09149-6206 by [128], who analysed the broadband X-ray data from the Swift, XMM-Newton, and NuSTAR telescopes. As is commonly the case, this AGN is found to have a complex X-ray spectrum with signatures of X-ray reflection from the inner accretion disc as well as both neutral and ionised absorption by more distant circumnuclear matter along the line of sight. [128] construct a 23-parameter spectral model to describe the primary continuum reflection from a BH accretion disc, which then passes through multiple layers of absorption [128]. Using that model, and the multi-satellite dataset, the authors computed the the posterior distribution of the parameters with the Goodman–Weare algorithm.

Marginalising over all other parameters of this global model, the inferred black hole spin is \(a_*= 0.94^{+0.02}_{-0.07}\). The principal ‘nuisance’ parameters that have bearing on the spin are the inclination of the accretion disc, the iron abundance of the disc, the ionisation state of the surface layers of the disc, and the geometry of the X-ray-emitting corona (characterised as a ‘lamppost’ height, also known as the coronal geometry). For typical AGN parameters, the inclination of the disc is readily determined by the high-energy ‘edge’ of the iron line profile, resulting in little degeneracy with the inferred spin (that results from the shape of the low-energy ‘tail’ of the iron line). The iron abundance principally influences the relative strength of the iron line and the Compton reflection hump with secondary impacts on the inferred spin, but is well constrained in datasets that cover both features (such as is the case here). The ionisation state of the accretion disc photosphere strongly influences the shape of the low-energy (\(E < \SI{4}{\keV}\)) reflection spectrum and shifts the energy of the iron line but, again, is well constrained by datasets that cover the 0.5 keV–50 keV band. It is the lamppost height that has the greatest degeneracy with the inferred spin; this is included in the quoted errors on the spin measurement [128].

Our mass posteriors come from the BLR modelling in the analysis of interferometric data by the [94], [125]. This mass estimate is consistent with the estimate from timing by [128].

We show the IRAS 09149-6206 posterior samples for \(M\) [94], [125] and \(a_*\)[128] in 3, assuming a BH timescale of \(\tau_\text{BH}= \tau_\text{Salp}/0.1 = \SI{4.5e8}{\yr}\). Note that, due to the two independent measurements for \(M\) and \(a_*\), there are no correlations between the two parameters in this case. As for M33 X-7, we again include the Regge trajectories for non-interacting (solid, grey line) and self-interacting (dotted, grey line) ULBs.

Since IRAS 09149-6206 has a spin close to the theoretical maximum value, the ‘\(2\sigma\) region’ of the \(a_*\) distribution surpasses that value. Instead of a regular Gaussian, one could use a truncated Gaussian to implement this physical upper limit on \(a_*\).

We also checked that, as is presumably the case for many SMBHs, the \(M\) distribution is better described by a Gaussian in \(\ln(M/\mathrm{M}_\odot)\) than in \(M\). The mean value of \(\ln(M/\mathrm{M}_\odot)\) corresponds to the median for the log-normal distribution of \(M\), and Gaussian error propagation, \(\sigma_M/M \approx |\sigma_{\ln(M/\mathrm{M}_\odot)}|\), is not necessarily a good approximation. Converting quoted errors on \(\ln(M/\mathrm{M}_\odot)\) into an error on \(M\) is thus potentially problematic – even more so when considering the ‘\(2\sigma\) region’. We will see the effect of this in 5.1.

3.4 Gravitational-wave binary BH mergers↩︎

Gravitational-wave (GW) data can be used to estimate the mass and spin of binary black hole (BBH) mergers. In principle, these data can thus be used to constrain ULB properties.

However, current detectors only measure a mass-weighted combination of the spins along the orbital angular momentum [144]. This complicates the determination of the individual spins at the time of the merger, while the unobservable pre-merger history also introduces additional challenges [7], [145]. These caveats apply to ULB constraints from BHB mergers based on the inferred pre-merger spins alone, as e.g.presented by [5], [6], [146].

To ameliorate these issues, [135], [145], [147] propose using a Bayesian hierarchical analysis to account for the pre-merger history using hyperprior distributions. [147] in particular investigate the influence of prior choices and consider an average of possible BBH lifetimes in the range of  × 10⁶– × 10¹⁰. Since the constraints are mainly driven by a few of the observed BBH systems, the authors conclude that more data is needed for more robust conclusions once a broad range of pre-merger scenarios is included in the analysis. We thus leave the including of GW data from BBH mergers for future work.

4 Statistical framework↩︎

With the BHSR rates from [sec:equilibrium] [sec:bosenovae], and the posterior distributions for the BH parameters from [sec:xray_binaries] [sec:smbh], we can now introduce our procedure for computing exclusion regions for ULB models.

As discussed in 2.4, we can exclude ULB models¹⁰ \(\boldsymbol{\alpha} = (\mu,f^{-1})\) for which the true values of \(\boldsymbol{\beta} = (M,a_*,\tau_\text{BH},\dots)\) lie above the Regge trajectory, i.e.when \(a_*> a_*^\text{crit}(M,\tau_\text{BH},\dots,\boldsymbol{\alpha})\). The ellipsis in \(\boldsymbol{\beta}\) is to indicate that there may be more parameters of the BH environment that have an influence on the Regge slope. For this work, however, we will not consider any such parameters and assume \(\tau_\text{BH}\) to be a fixed value. Also note that any other BH nuisance parameters, such as the inclination angle \(i\) or virial factor \(\mathfrak{f}\), can already be marginalised out in the original BH analysis.

In terms of probabilities, we can write the condition on \(a_*\) as a Heaviside function, \(p(a_*| M,\boldsymbol{\alpha}) = \Theta(a_*^\text{crit} - a_*)\), i.e.\(p(a_*| M,\boldsymbol{\alpha}) = 0\) if the BH lies above the Regge slope, and otherwise \(p(a_*| M,\boldsymbol{\alpha}) = 1\). In particular, the BH’s parameters only depend on the ULB parameters when SR is active and the the condition on \(a_*\) applies, such that the conditional probability \(p(\boldsymbol{\beta}|\boldsymbol{\alpha})\) factorises as \(p(\boldsymbol{\beta}|\boldsymbol{\alpha}) = p(\boldsymbol{\beta}) \, p(a_*| M,\boldsymbol{\alpha})\).

While we do not know the true values of \(\boldsymbol{\beta}\), we can obtain a sampling distribution from the BH data \(\text{D}\) to compute the posterior \(p(\boldsymbol{\alpha}|\text{D})\). Using the law of conditional probabilities and Bayes’ theorem, we find that \[\begin{align} p(\boldsymbol{\alpha} | \text{D}) &= \int \! p(\boldsymbol{\alpha}, \boldsymbol{\beta} | \text{D}) \, \mathrm{d}\boldsymbol{\beta} = \int \! p(\boldsymbol{\beta} | \text{D}) \, p(\boldsymbol{\alpha} | \boldsymbol{\beta}, \cancel{\text{D}}) \, \mathrm{d}\boldsymbol{\beta} \\ &= \int \! p(\boldsymbol{\beta} | \text{D}) \, \frac{p(\boldsymbol{\beta} | \boldsymbol{\alpha}) \, p(\boldsymbol{\alpha})}{p(\boldsymbol{\beta})} \, \mathrm{d}\boldsymbol{\beta} \\ &= p(\boldsymbol{\alpha}) \int \! p(\boldsymbol{\beta} | \text{D}) \, p(a_*| M, \boldsymbol{\alpha}) \, \mathrm{d}\boldsymbol{\beta} \tag{14} \\ &\approx \frac{p(\boldsymbol{\alpha})}{N} \sum_{i=1}^{N} p(a_*^i | M^i, \boldsymbol{\alpha}) \, ,\tag{15} \end{align}\] where we use a slash to indicate the removal of the redundant dependence on the data \(\text{D}\) and where, in the last step, we assume that the BH analysis provides \(N\) equally-weighted samples \(\boldsymbol{\beta}^i = (M^i,a_*^i)\) from the posterior \(p(\boldsymbol{\beta} | \text{D})\), which we can use to compute the integral in 14 via Monte Carlo integration.

Since \(p(\boldsymbol{\alpha} | \text{D}) \propto p(\boldsymbol{\alpha}) \, p(\text{D}| \boldsymbol{\alpha})\), we may re-interpret the sum in 15 as a marginal likelihood. Furthermore, the computation of \(p(a_*| M,\boldsymbol{\alpha})\) requires a computation of the BHSR rates and assumptions about the BH evolution. As we will discuss in 5.2, both our computations and assumptions could be refined further. Even then, the framework presented here is general enough to allow for such modifications without having to change to the framework itself.

5 Results and Discussion↩︎

Our main results can be found in 4, where we show the posterior distribution for \((\mu,f^{-1})\) for IRAS 09149-6206(left panel) and M33 X-7(right panel).

We assume log-uniform prior distributions for \(\mu\) and \(f^{-1}\), limited by the ranges shown in 4. The ranges approximately correspond to the relevant (observable) parameter space, except that we restrict \(f< \mathrm{M}_\text{P}\). The latter is expected in string theory, and more generally if the effective field theory is free of corrections from quantum gravity [148].

We use emcee [149] with 48 walkers to generate, after a burn-in run, 60000 samples per chain i.e.a total of around  × 10⁷ samples (total runtime around 600 CPUh). Finally, we compute the 95% credible regions (CRs) at highest posterior density (HPD) for the ULB parameters.

We can see from 4 that the equilibrium regime and the bosenova condition lead to similar disfavoured regions of parameter space in \(\mu\). In particular, the lower \(\mu\) cutoff is set by the condition to achieve a significant cloud occupation while, at larger \(\mu\), the limits are cut off by the SR condition. It has been noted before, and we will see this explicitly in 5.1, that the range of excluded \(\mu\) may be extended by including higher occupation levels with \(n > 2\). However, here we only show results considering the \(\left|211\right\rangle\) since our computed SR rates are most reliable for this case. To further emphasise this, we also show the \(\mu\) value for which we mostly, i.e.for at least 95% of the \(M\) samples for a given BH, obtain \(\alpha < 0.2\). The ULB constraints are most robust in this regime where, moreover, the occurrence of bosenovae is strongly disfavoured (see 2.3.3).

From the relative posterior densities, it appears that the constraints on IRAS 09149-6206 are weaker than for M33 X-7, which we will explore in more detail in 5.1. This highlights the relevance of a precise determination of SMBH masses to obtain meaningful constraints. While only BHs with a high spin can give strong limits on ULB parameters, the connection between \(M\) and \(\mu\) via \(\alpha = M \mu\) implies that we can only locate these constraints if \(M\) is measured precisely enough. Large uncertainties on \(M\) prevent us from doing so and, although it is clear that such a BH would give strong constraints on ULBs, an improvement in the BH mass measurement is required to tell us about the corresponding \(\mu\) range. The primary and secondary criteria for suitable BH candidates for obtaining ULB constraints should thus be a high \(a_*\) and precisely determined \(M\), respectively.

Note that the adopted priors implicitly assume the existence of ULBs within their possible parameter range, and that we can always define credible regions of the posterior distribution to exclude part of the ULB parameter space. Such constraints are only meaningful if the prior is sensible – i.e.in particular if such ULBs exist. In order to assess the probability of ULBs’ existence, one has to perform a model comparison, e.g.using Bayes’ factors. Here, we are interested in setting limits, and thus do not perform this additional step.

More generally, the choice of prior on \((\mu,f^{-1})\) will necessarily affect the exclusion regions. [150] analysed this in more detail for axion models and found that the prior dependence can be very strong. Still, at least with regards to changing the range of the log-uniform priors, we do not expected extreme deviations as for e.g.the choice between log-uniform and uniform priors. In particular, when combining data from multiple BHs, the overall prior range will increase when considering the ‘sensitive \(\mu\) region’ for all BHs, making a joint fit increasingly conservative.

We also want to stress that the marginal distributions on \(\mu\) and \(f^{-1}\) will not be very constraining:¹¹ there is practically no limit on \(f^{-1}\) alone due to the available \(\mu\) parameter space with high posterior density outside of the region related to the inferred \(M\) range. Similarly, the region of \(\mu\) that can be constrained when marginalising \(f^{-1}\) is much narrower than suggested by the two-dimensional exclusion region. When considering self-interactions, it is important to show the two-dimensional contour instead of assuming that the excluded one-dimensional \(\mu\) range will correspond to the range suggested by the two-dimensional contour at low \(f^{-1}\). This would even be more so the case in a frequentist framework.

5.1 Comparison with other works↩︎

Surveying the literature, we identify two main strategies to obtain ULB constraints, which we refer to as the ‘box method’ and ’ Gaussian distance method’, respectively.

The ‘box method’.

[3] applied the following checks for BH parameter estimates \((\hat{M},\hat{a}_*)\) with associated uncertainties \((\sigma_{\hat{M}},\sigma_{\hat{a}})\) and \(\tau_\text{BH} = \min(\tau_\text{age}, \SI{40}{\mega\yr})\). A given ULB model can be excluded if at least one \(\left|nlm\right\rangle\) with \(l \leq 5\) (i.e.\(n \leq 6\)) and all \(M \in [\hat{M} - 2\sigma_{\hat{M}}, \hat{M} + 2\sigma_{\hat{M}}]\) and \(a = \hat{a} - 2\sigma_{\hat{a}}\) fulfil the SR and bosenova conditions (see 2.3.2). Comparing to 2 3, these conditions demand that the rectangle defined by the mass and spin measurement \(\pm 2 \sigma\) in each direction is contained within the envelope of the Regge trajectory, and that bosenovae do not spoil the spindown of the BHs. [7] also used the ‘box method’, but based on the equilibrium scenario and conditions (see 2.3.1).

The ‘Gaussian distance method’.

[5], [6] proposed to interpret \((\hat{M}, \hat{a}_*)\) estimates and their errors as uncorrelated Gaussian measurements. They then defined ‘effective errors’ with respect to the Regge trajectory, which allowed them to compute a distance measure between the uncorrelated Gaussian distribution of the measurement and the area above the Regge trajectory. Their goal was to compute the weighted overlap between Gaussian distribution and the Regge trajectory “in a numerically efficient manner”, as visualised in [5]. In this sense, their underlying logic is similar to ours with the difference being their assumption about the uncorrelated Gaussian distribution and their computational scheme which will, as they pointed out, only approximately give the desired result when considering levels with \(n > 2\) [5].

Comparison with our method.

To facilitate a direct comparison with the literature in 5, we first ignore self-interactions – i.e.take \(f\to \infty\) or, equivalently, set \(\lambda = 0\) – and compute the posterior probability at fixed values of \(\mu\).¹² Instead of implementing the ‘Gaussian distance method’, we will use our method, but with an uncorrelated Gaussian instead of the full posterior distributions. This is equivalent to what [5] wanted to compute. For illustrative purposes, for the uncorrelated Gaussian, we choose estimators in \(M\) for both BHs, even though the underlying distribution of the SMBH is better described by a Gaussian in \(\log(M/\mathrm{M}_\odot)\). Finally, note that the ‘box method’ gives a binary outcome; the probability is either \(P = 0\) or \(P = 1\).

Consider the left panel of 5 for IRAS 09149-6206. The ‘box method’ finds \(P = 1\) for all ULB masses and, indeed, the lines for the other methods and assumptions also remain above the exclusion threshold. However, thanks to our probabilistic interpretation, our method could yield constraints by combining the data from multiple BHs – even when they individually would not constraint ULB parameters. Further note that the deviations between the full and the approximated results mainly come from the approximate Gaussian error propagation, as the data themselves are uncorrelated. 5 also demonstrates explicitly that the inclusion of more levels can extend the bounds to higher ULB masses. This is because the higher \(n\), and thus higher \(l\), allow us to fulfil the SR condition in 4 for larger \(\mu\). For these constraints, we only consider the SR rates \(\Gamma_{{nlm}}\) of the individual levels with no regards to the BH evolution and time required to fill these levels.

The case of M33 X-7 is shown in the right panel of 5. Given the clearly visible correlations in the \((M,a_*)\) distribution in 2, it is somewhat surprising that neglecting them appears to be an acceptable approximation. This might only be a somewhat lucky accident for M33 X-7, but at least in this case the approximation seems to be inconsequential. The inclusion of levels with \(n \leq 6\) again extends the ULB mass range of the constraints.

6 presents a direct comparison with ULB constraints from M33 X-7 in the \((\mu,f^{-1})\) plane. [6] used the bosenova prescription for self-interactions, a Gaussian approximation for the BH errors, and included levels up to \(n = 6\) (dashed, purple line). The result broadly agrees with our implementation (dashed, black line) of the bosenova prescription in terms of \(f^{-1}\), while we observe that the inclusion of higher SR levels extends the limits from [6] to higher values of \(\mu\). Our results disagree with [6] at the low mass end by a factor of \(\mathcal{O}(2)\) on \(\mu\). This is due to [6] [5], [8] neglecting the large \(\ln(N_{\Delta})\) factor required to extract significant spin, and the effect of the companion. [7] used the ‘box’ statistical method, and the equilibrium model for self-interactions (solid, red line). As mentioned before, our result (solid, black line) gives rise to a slightly wider exclusion region on \(\mu\), due to the overly conservative nature of the box method. Note that our implementation of the box method (dotted, black line) closely reproduces the results by [7] for both \(\mu\) and \(f^{-1}\).

Advantages of our approach include that we do not assume Gaussianity, which is not always justified – in particular for the treatment of upper limits on \(a_*\). It is important to stress again that, for SMBHs, Gaussianity only appears to be justified for \(\ln(M/\mathrm{M}_\odot)\) rather than \(M\). While past work also ignored potential correlations between \(M\) and \(a_*\), we do not find this to be a major issue in our examples. Our approach also turns out to be computationally feasible with Monte Carlo integration, alleviating concerns of [5].

As a final remark, note that, in some sense, the ‘box method’ is conservative as it accounts for all possible correlations and at least some non-Gaussian features of the \((M,a_*)\) distribution. This naturally leads to weaker exclusion limits, i.e.smaller exclusion regions for \(\mu\), as can be seen in 5. However, the \((M,a_*)\) distribution contains additional information, and discarding it seems unnecessary. Moreover, the ‘\(2\sigma\) interval’ relies on a (possibly ill-defined) extrapolation from the ‘\(1\sigma\) error’ of \(\hat{M}\) and \(\hat{a}_*\).¹³ Due to this, it is not straightforward to apply the ‘box method’ to multiple BHs or to compare to limits obtained at a different exclusion level.

5.2 Limitations and possible extensions↩︎

As explored before, more powerful limits on ULB masses can be derived using higher SR levels. While we can compute the free SR rates for \(n > 2\), the most reliable limits on the self-coupling are currently only possible for \(n = 2\). One important extension of existing work is to compute the equilibrium rates when higher levels are relevant. Furthermore, one may include more interactions, such as stimulated decay in the case of the axion-photon coupling \(g_\gamma \propto f^{-1}\).

From the data modelling side, our method can also easily include additional parameters to propagate their uncertainties into the ULB exclusion regions. One such parameter is the BH timescale \(\tau_\text{BH}\), as discussed in 3.1, which could be included with a physically-informed prior distribution or constrained by related data sets. We also did not consider possible interactions between the boson cloud and the SMBH environment, which could affect the cloud growth. For instance, [151] computed the boson cloud’s perturbative stability in a thin accretion disc or in the presence of a stellar halo and found that their stability estimator worsens as \(\alpha^q\) with \(q \gtrsim 5\).

Additional nuisance parameters can also be included when inferring the mass and spin of a given BH to better reflect their uncertainties. In 3, we provided a summary of various aspects of the BH analysis to highlight the underlying assumptions and possible sources of uncertainties. In the case of X-ray reflection modelling, several parameters, e.g.the inclination and iron abundance of the accretion disc, are well-known to be degenerate with the spin in state-of-the-art reflection models. In addition, in the case of continuum-fitting models used in XRBs, spin inference is affected by intrinsic astrophysical uncertainties, such as the inclination angle, distance, or the BH’s mass. More information on these parameters could be provided by involving more constraints, such as data from current and upcoming X-ray polarimetry missions including the currently operating IXPE observatory.

Moreover, in addition to including more parameters, we could also employ a hierarchical Bayesian analysis, considering both individual BH data as well as BH population parameters and relations. This has already been explored in the context of BBH mergers (see 3.4), which we did not consider in this work. The past work assumed hyperpriors for the initial masses and spins of the merging BH, which we could also employ for SMBHs and XRBs. This would allow us to compute the final BH mass and spin from the full BH evolution including SR, which can be compared to the BH data [152]. Finally, note that various nuisance parameters, such as the virial factor \(\mathfrak{f}\)[153] or a possible \(M\)-\(L\) relationship between the BH mass and luminosity [154] might be exploited to refine the \((M,a_*)\) inference for BHs. While the complexity of our analysis would increase when employing hierarchical model, these are natural extensions of a Bayesian framework, which – apart from refining the accuracy of its inference – could help to improve convergence and gain new insights into BH physics.

Compared to both alternative methods, our Bayesian statistical approach more easily generalises to a joint analysis of multiple BHs and the inclusion of other constraints: the Bayesian logic is naturally extendable while allowing us to include any correlations or additional uncertainties through nuisance parameters where necessary. This is useful when considering BHSR constraints in global fits of axion models, such as the analysis by [150] within the GAMBIT global-fitting framework. Apart from previously used BH data [6], such a global fit could also consider the well-studied XRB GRS 1915+105, whose spin has been estimated from both the X-ray reflection and continuum-fitting methods [155], [156]. Another promising prospect is the inclusion of future data on tidal disruption events near SMBHs from the Vera C. Rubin Observatory Legacy Survey of Space and Time (LSST), which depend on the BH spin and can constrain ULBs in the range of \(\SI{e-20}{\eV} < \mu< \SI{e-18}{\eV}\) [151].

6 Summary and conclusions↩︎

We have presented a universal approach for obtaining constraints on ultralight bosons (ULBs) from black hole (BH) data. We exemplified our methodology with two well-studied BHs: M33 X-7, a stellar-mass BH, and IRAS 09149-6206, a supermassive BH. For each of these, we used posterior samples for the BH’s mass and spin \((M,a_*)\), and computed the posterior probability on the ULB model parameters \((\mu,f^{-1})\). In particular, these are the first ULB constraints derived from IRAS 09149-6206.

Our method introduces a rigorous statistical framework, extending previous approaches and preserving their advantages while, at the same time, giving more accurate or powerful constraints. In particular, our method appears to be the only proposed methodology to self-consistently derive ULB constraints from superradiance (SR) in multiple BHs. We only rely on samples from the \((M,a_*)\) posterior distribution alone and need not reproduce the entire BH data analysis tool chain.

Our approach also allowed us to clarify the statistical nature and reasons for disagreement between previous works (see 5.1). In particular, we improve on previous methods by capturing correlations and non-Gaussianities of the \((M,a_*)\) distribution, providing a more complete treatment of statistical uncertainties from the BH analyses. Still, we emphasise again that BH inference is affected by potentially large systematic uncertainties related to the assumptions and limitations of accretion disc models, BH evolution, bosenova vs equilibrium scenarios, or computation of the SR rates.

Let us summarise the most important results from our work and the comparison to the literature. In this initial study we observe that

• the proposed Bayesian approach automatically includes all information from the data while being easily realisable,

• neglecting correlations between \(M\) and \(a_*\) may be acceptable,

• the Gaussian approximation holds for \(\ln(M/\mathrm{M}_\odot)\) instead of \(M\) for SMBHs,

• the Gaussian approximation is poor when \(|a_*| \sim 1\),

• the ‘box method’ leads to overly conservative limits and does not leverage the full information contained in the BH data,

• inclusion of higher levels up to \(n > 2\) can potentially constrain higher masses by up to an order of magnitude,

• compared to the bosenova, the equilibrium method for self-interactions reduces the maximum value of \(f^{-1}\) probed by any given BH by approximately one order of magnitude.

From the theoretical side, we are still limited by the availability of multi-level equilibrium BH evolution predictions, which would allow us to extend constraints for larger \(\mu\) values. From the data side, a wider applicability of our method is limited by the availability of a (sampled) distribution of \((M,a_*)\), and we thus encourage all groups deriving these constraints to make their likelihoods or posteriors publicly available. Each new data set will at worst linearly increase the computational time required for the sampling and more likely be sublinear due to faster convergence and due to abandoning parameter points that fall below some likelihood threshold after considering part of the data. Overall, the computational effort should be manageable.

Our method has improved on the statistical interpretation of BHSR constraints compared to previous works [5], [7], and we have identified causes of disagreement in the literature due to rates models and statistics. Our conclusions have consequences for BHSR constraints on the string theory landscape. [8] computed constraints on the landscape using a version of the implementation of rates and statistics by [6]. We have verified the statistical approximation of [5], [6] as self-consistent and valid. An error in the low mass end rates of [5], [6] is relatively small and should not affect conclusions about constraints on the landscape, due to the logarithmic nature of the \(\mu\) distribution. Likewise, using the equilibrium rather than bosenova model for self-interactions will lead to only a slight shift in the maximum topological numbers probed by BHSR, due to the relatively small shift in \(f^{-1}\) and slow trends in landscape models on the variation of this parameter. On the other hand, inclusion of higher levels in [8] significantly widens the \(\mu\) range that BHSR can probe. The effect this will have on conclusions relating to the landscape is unclear. [8] used a large number of overlapping BHs in the stellar regime, so the effect from any single BH will be small, and the limits should only be affected by the inclusion of higher levels on the lowest mass BH in the sample. Nonetheless, this underscores the importance of computing full BH Regge trajectories in future improvements of BHSR constraints, in order to leverage the constraining power of higher levels accurately, and the importance of extending our methods to large numbers of BHs if the posterior samples can be made available.

Possible extensions of our Bayesian method are discussed in detail in 5.2, including additional parameters related to the BH timescale, Bayesian hierarchical modelling of BH populations, or a companion in binary systems such as M33 X-7. Assuming hyperprior distributions for the initial BH mass and spin would further allow us to compute the final BH mass and spin from the full SR evolution. This treatment is more refined than simply comparing the associated SR and BH timescales. Similarly, this approach can also incorporate the constraints from binary BH mergers, following [135], [145], [147]. We can also include the physics of direct gravitational wave emission from the boson cloud in our model, which has been considered by [41]–[43], [46], [157], [158].

Moreover, constraints on \(\mu\) and \(f\) from BHSR are highly complementary to experimental and astrophysical probes of ULBs, with particular application to the QCD axion, fuzzy DM, and string theory. Our systematic approach will allow this complementarity to be leveraged in global statistical analyses.

To further facilitate these, and possibly other, future extensions, we made our software code available on Github [51].

Acknowledgements↩︎

We thank Masha Baryakhtar, Marios Galanis, Olivier Simon, and Sam Witte for helpful, substantial input on our draft and discussions related to their ongoing and past works; Lijun Gou and Xueshan Zhao for \((M,a_*)\) samples from [112]; Jinyi Shangguan for \(a_*\) samples from the [94], [125]; Dominic J.Walton for \(a_*\) samples from [128]; and Jakob van den Eijnden, Rob Fender, and Matt Jarvis for helpful discussions on BH masses.

SH has received funding from the European Union’s Horizon Europe research and innovation programme under the Marie Skłodowska-Curie grant agreement No 101065579. DJEM is supported by an Ernest Rutherford Fellowship from the STFC, Grant No. ST/T004037/1 and by a Leverhulme Trust Research Project (RPG-2022-145). JSR acknowledges the support from the Science and Technology Facilities Council (STFC) under grant ST/V50659X/1 (project reference 2442592) and the Institute of Astronomy, as well as from a NASA Astrophysics Data Analysis Program (grant 80NSSC24K0617). JHM acknowledges funding from a Royal Society University Research Fellowship. We acknowledge the CINECA award under the ISCRA initiative, for the availability of high-performance computing resources and support.

We made use of the BibCom [159] and WebPlotDigitizer [160] tools, as well as the Python packages numba, numpy [161], and scipy [162].

Data Availability↩︎

The Python software code and data underlying this article are available on Github at https://github.com/sebhoof/bhsr.

References↩︎