Black Hole Superradiance Signatures of Ultralight Vectors

TL;DR

This paper analyzes how ultralight vector particles can extract energy from rotating black holes through superradiance, leading to observable gravitational wave signals and constraining particle masses using astrophysical data.

Contribution

It provides analytical calculations of superradiant growth rates for vectors around black holes and explores observational strategies for detecting these particles via gravitational waves.

Findings

Superradiance times for vectors can be as short as a second around stellar black holes.

Measurements of black hole spins constrain vector particle masses in specific ranges.

Gravitational wave detectors could observe signals from vector particle annihilations.

Abstract

The process of superradiance can extract angular momentum and energy from astrophysical black holes (BHs) to populate gravitationally-bound states with an exponentially large number of light bosons. We analytically calculate superradiant growth rates for vectors around rotating BHs in the regime where the vector Compton wavelength is much larger than the BH size. Spin-1 bound states have superradiance times as short as a second around stellar BHs, growing up to a thou- sand times faster than their spin-0 counterparts. The fast rates allow us to use measurements of rapidly spinning BHs in X-ray binaries to exclude a wide range of masses for weakly-coupled spin-1 particles, $5 \times 10^{-14} - 2 \times 10^{-11}$ eV; lighter masses in the range $6 \times 10^{-20} - 2 \times 10^{-17}$ eV start to be constrained by supermassive BH spin measurements at a lower level of confidence. We also explore routes to detection of new vector particles possible with the advent of gravitational wave (GW) astronomy. The LIGO-Virgo collaboration could discover hints of a new light vector particle in statistical analyses of masses and spins of merging BHs. Vector annihilations source continuous monochromatic gravitational radiation which could be observed by current GW observatories. At design sensitivity, Advanced LIGO may measure up to thousands of annihilation signals from within the Milky Way, while hundreds of BHs born in binary mergers across the observable universe may superradiate vector bound states and become new beacons of monochromatic gravitational waves.

Tables

Table 1. Table 1: Stellar-mass BHs that set limits on vectors as shown in Figure  5 (data compiled in McClintock  et al. ( 2014 ) unless otherwise specified). Mass and spin errors are quoted at 1 ​ σ 1 𝜎 1\sigma and 2 ​ σ 2 𝜎 2\sigma , respectively  Steiner and McClintock ( tion ) . We use the lower continuum-fitting spin values for GRO J1655-40 Reynolds ( 2014 ) , and τ b ​ h = τ Salpeter / 10 subscript 𝜏 𝑏 ℎ subscript 𝜏 Salpeter 10 \tau_{bh}=\tau_{\mathrm{Salpeter}}/10 for GRS1915+105 to set a conservative limit.

#	Object	Mass ($M_{\odot }$)	Spin	Age (yrs)	Period (days)	$M_{\mathrm{comp}.\mathrm{star}}$ ($M_{\odot }$)
1	GRS1915+105	${12.4}_{-1.8}^{+2.0}$ Reid et al. (2014a)	$>0.95$ Reid et al. (2014b); Steiner and McClintock (tion)	$3-5\times {10}^9$ Dhawan et al. (2007)	33.85 Steeghs et al. (2013)	0.47 $\pm$ 0.27 Steeghs et al. (2013)
2	Cyg X-1	$14.8\pm 1.0$	$>0.99$ Gou et al. (2014)	$4.8-7.6\times {10}^6$ Wong et al. (2012)	5.599829 Gou et al. (2011)	17.8 Gou et al. (2011)
3	GRO J1655-40	$6.3\pm 0.5$	${0.72}_{-.24}^{+.16}$ Steiner and McClintock (tion)	$3.4-10\times {10}^8$ Willems et al. (2005)	2.622 Willems et al. (2005)	2.3 - 4 Willems et al. (2005)
4	LMC X-1	$10.91\pm 1.4$	${0.92}_{-.18}^{+.06}$ Gou et al. (2009)	$5-6\times {10}^6$ Orosz et al. (2009)	3.9092 Ruhlen et al. (2011)	31.79$\pm$3.48 Ruhlen et al. (2011)
5	M33 X-7	$15.65\pm 1.45$	${0.84}_{-.10}^{+.10}$Steiner and McClintock (tion)	$2-3\times {10}^6$ Orosz et al. (2007)	3.4530 Pietsch et al. (2006)	$\gtrsim 20$ Pietsch et al. (2006)

Table 2. Table 2: Supermassive BHs that set limits on vectors as shown in Fig.  5 (compiled in Reynolds ( 2013 , 2014 ) ; our analysis excludes BHs without a mass error estimate). The mass and spin errors are quoted at 1 ​ σ 1 𝜎 1\sigma and 90 % percent 90 90\% confidence interval, respectively.

#	Object	Mass (${10}^6{M}_{\odot }$)	Spin
1	Ark 120	$150\pm 19$	${0.64}_{-0.11}^{+0.19}$
2	Mrk 79	$52.4\pm 14.4$	$0.7\pm 0.1$
3	NGC 3783	$29.8\pm 5.4$	$>0.88$
4	Mrk 110	$25.1\pm 6.1$	$>0.89$
5	Mrk 335	$14.2\pm 3.7$	${0.83}_{-0.13}^{+0.09}$
6	MCG-6-30-15	${2.9}_{-0.16}^{+0.18}$	$>0.98$
7	NGC 4051	$1.91\pm 0.78$	$>0.99$

Table 3. Table 3: ‘Flat-space’ results for the leading- α 𝛼 \alpha gravitational wave emission power/ ( N m 2 ​ G r g 4 ) superscript subscript 𝑁 𝑚 2 𝐺 superscript subscript 𝑟 𝑔 4 (N_{m}^{2}\frac{G}{r_{g}^{4}}) for vector bound states. For given j 𝑗 j and ℓ ℓ \ell , and we take m = j 𝑚 𝑗 m=j , and n = 0 𝑛 0 n=0 .

$j$	$\ell =j-1$	$\ell =j$	$\ell =j+1$
1	$\frac{32{\alpha }^{12}}{5}$	$\frac{{\alpha }^{16}}{2560}$	$\left(5\times {10}^{-8}\right){\alpha }^{20}$
2	$\frac{{\alpha }^{16}}{126}$	$\left(1\times {10}^{-6}\right){\alpha }^{20}$	$\left(1\times {10}^{-11}\right){\alpha }^{24}$
3	$\left(6\times {10}^{-6}\right){\alpha }^{20}$	$\left(7\times {10}^{-10}\right){\alpha }^{24}$	$\left(4\times {10}^{-15}\right){\alpha }^{28}$
4	$\left(2\times {10}^{-9}\right){\alpha }^{24}$	$\left(2\times {10}^{-13}\right){\alpha }^{28}$	$\left(2\times {10}^{-13}\right){\alpha }^{32}$
5	$\left(4\times {10}^{-13}\right){\alpha }^{28}$	$\left(2\times {10}^{-17}\right){\alpha }^{32}$	$\left(6\times {10}^{-23}\right){\alpha }^{36}$

Table 4. Table 4: Parameters of astrophysical distributions used to estimate number of events.

Quantity	Optimistic	Realistic	Pessimistic
$M_0$ ($M_{\odot }$)	7.9	4.7	4.7
$P\left(a_{\ast }\right)$	$90\%$ above 0.8	$30\%$ above 0.8	uniform
BHFR (per century)	0.9	0.38	0.08

Equations

$$\frac{\omega }{m}<{\Omega }_H,$$

$${\mathcal{L}}_{\mathcal{A}}=-\frac{1}{4}{F}_{\mu \nu }{F}^{\mu \nu }-\frac{1}{2}{\mu }^2{A}_{\mu }{A}^{\mu }.$$

$$D_{\mu }{F}^{\mu \nu }={\mu }^2{A}^{\nu },$$

$$A^{\mu }(t,x)=\frac{1}{\sqrt{2\mu }}\left({\Psi }^{\mu }(x)e^{-i\omega t}+c.c.\right).$$

$$({\omega }^2-{\mu }^2){\Psi }^{\nu }\simeq -{\nabla }^2{\Psi }^{\nu }+{\omega }^2(1+g^{00}){\Psi }^{\nu }.$$

$$(\omega -\mu ){\Psi }^{\nu }\simeq -\frac{{\nabla }^2}{2\mu }{\Psi }^{\nu }+\frac{\alpha }{r}{\Psi }^{\nu },$$

$${\Psi }_i=R^{n\ell }(r)Y_i^{\ell ,jm}(\theta ,\varphi ),$$

$$\omega \simeq \mu \left(1-\frac{{\alpha }^2}{2(n+\ell +1{)}^2}\right).$$

$$A_i=\frac{1}{\sqrt{\pi }{\mu }^{1/2}{a}^{3/2}}{e}^{-r/a}\left(\begin{array}{c}-\cos \omega t\\ -\sin \omega t\\ 0\end{array}\right),$$

$$A_0=\frac{1}{\sqrt{\pi }{\mu }^{3/2}{a}^{5/2}}{e}^{-r/a}\sin \theta \cos (\omega t-\varphi ),$$

$$\Gamma \sim {\alpha }^{2j+2\ell +5}(m{\Omega }_H-\omega ),$$

$${\Gamma }_{\mathrm{scalar}}\sim {\alpha }^{4\ell +5}(m{\Omega }_H-\omega ).$$

$$\Gamma \simeq 4a_{\ast }{\alpha }^6\mu ,$$

$$N_m\simeq \frac{G{M}_{BH}^2\Delta a_{\ast }}{m}\simeq \frac{10^{77}}{m}{\left(\frac{M_{BH}}{10M_{\odot }}\right)}^2\frac{\Delta a_{\ast }}{0.1};$$

$$\rho \sim \frac{\alpha M}{a^3}\sim \frac{{\alpha }^7}{G^3{M}^2}\sim {\left(20\text{MeV}\right)}^4{\left(\frac{10M_{\odot }}{M}\right)}^2{\left(\frac{\alpha }{0.2}\right)}^7;$$

$${\Gamma }_{\text{sr}}{\tau }_{\text{bh}}\ge \log N_{\text{m}}.$$

$$A_{\mu }=\frac{1}{\sqrt{2\mu }}\sum _i\sqrt{N_i}\left({\Psi }_{\mu }^i{e}^{-i\omega t}+c.c.\right),$$

$$T_{\mu \nu }\supset {\mu }^2{A}_{\mu }{A}_{\nu }+\cdots =\frac{1}{2\mu }\sum _{i,j}\sqrt{N_i{N}_j}\left(\underset{\text{‘annihilations’}}{\underbrace{{\Psi }_{\mu }^i{\Psi }_{\nu }^j{e}^{-i({\omega }_i+{\omega }_j)t}}}+\underset{\text{‘transitions’}}{\underbrace{{\Psi }_{\mu }^i{{\Psi }_{\nu }^j}^{\ast }{e}^{-({\omega }_i-{\omega }_j)t}}}+h.c.\right)+\dots ,$$

$$P\sim 60{\alpha }^{12}\frac{G{N}^2}{r_g^4}.$$

$$\frac{dP}{d\Omega }\propto (1+6{\cos }^2{\theta }_k+{\cos }^4{\theta }_k),$$

$$f_{\mathrm{GW}}$$

$$\simeq 1.3\text{ kHz}\frac{\alpha }{0.2}\frac{10M_{\odot }}{M}.$$

$$N(t)=\frac{N_0}{1+{\Gamma }_{\text{ann}}{N}_{0}t},$$

$$t_{\mathrm{sig}}\simeq (N_m{\Gamma }_{\mathrm{ann}}{)}^{-1}\simeq 80\text{ sec }{\left(\frac{0.2}{\alpha }\right)}^{11}\left(\frac{M}{10M_{\odot }}\right).$$

$$\dot{f}\simeq 10^{-8}Hz/sec\times \left(\frac{\mu }{10^{13}\mathrm{eV}}\right){\left(\frac{\alpha }{0.2}\right)}^3\frac{250\mathrm{hr}}{t_{\mathrm{sig}}}.$$

$$h_{\text{det}}=\mathrm{SNR}{C}_{\text{tf}}\frac{\sqrt{S_h}}{N_{\text{seg}}^{1/4}{T}_{\text{coh}}^{1/2}},$$

$$h={\left(\frac{4GP}{r^2{\omega }_{\mathrm{GW}}^2}\right)}^{\frac{1}{2}},$$

$${\mathbb{P}}_{\mathrm{abs}}\simeq {\left(\frac{(j-s)!(j+s)!}{(2j)!(2j+1)!!}\right)}^2\prod _{n=1}^j\left[1+{\left(\frac{\omega -m\Omega }{n\kappa }\right)}^2\right]2\left(\frac{\omega -m\Omega }{\kappa }\right){\left(\frac{A_H\kappa }{2\pi }\omega \right)}^{2j+1},$$

$${\Psi }_i\simeq {\Psi }_i^H$$

$${\Psi }_0\simeq {\Psi }_0^H$$

$$Y_i^{\ell ,jm}(\theta ,\varphi )=$$

Full text

I Introduction

Light, weakly-coupled new particles are a feature of many beyond Standard Model (SM) physics scenarios Essig et al. (2013); Arvanitaki et al. (2010a, b), and depending on how they couple to the SM, there are a plethora of ways in which they can be searched for experimentally. At an absolute minimum, any new particles must interact gravitationally. However, the very high value of the Planck scale results in very small rates for the purely-gravitational production of new particles.

One way to get around very small production rates is to take advantage of coherent enhancement. As in a laser, it can lead to exponential amplification of a very small initial population (including vacuum fluctuations). For light bosonic particles, such a setup can be realized around a spinning black hole (BH), where gravitational bound states of the bosonic field can extract energy and angular momentum from the BH through ‘rotational superradiance’, the wave analogue of the Penrose process Zeldovich (1971); Misner (1972); Starobinskii (1973). This process results in the exponential growth of some bound states, eventually spinning down the BH and creating a very high occupation number bosonic ‘cloud’ around it.

Using BH superradiance as a tool to search for new light bosons was proposed in Arvanitaki et al. (2010a), and spin-0 superradiance dynamics and phenomenology have been studied in many papers Ternov et al. (1978); Zouros and Eardley (1979); Detweiler (1980); Dolan (2007); Arvanitaki and Dubovsky (2011); Yoshino and Kodama (2014); Arvanitaki et al. (2015); Brito et al. (2015a, b); Arvanitaki et al. (2017). While there are superradiant (growing) bound states for any scalar mass, the growth rates are very suppressed for Compton wavelengths significantly larger or smaller than the BH size. Observationally, the two main signatures are the lack of rapidly-spinning BHs, and the monochromatic gravitational waves (GWs) sourced by the bosonic cloud itself Arvanitaki and Dubovsky (2011); Arvanitaki et al. (2015). Precise measurements of BH spins through X-ray observations constrain sufficiently weakly coupled spin-0 particles over a wide mass range, while near-term gravitational wave observations have the potential to discover such particles at other masses Arvanitaki and Dubovsky (2011); Arvanitaki et al. (2015, 2017).

In this paper, we study the theory of spin-1 superradiance around BHs, and use the associated astrophysical phenomenology to search for and constrain gravitationally-coupled massive vector particles below $10^{-11}$ eV. It is technically natural for vectors to have small masses, and there are a number of string theory constructions which result in very light and weakly coupled spin-1 particles Goodsell et al. (2009); Camara et al. (2011). Computationally, the extra degrees of freedom in vector fields, as compared to scalar fields, make solving for the behavior of bound states more challenging. However, for vector Compton wavelengths larger than the BH size, the bound states are ‘non-relativistic’ and hydrogen-like, and we analytically compute the growth and decay rates of such states. Our calculation applies around Kerr BHs of arbitrary spin, which is significant since astrophysical signatures from superradiance are dominated by fast-spinning BHs. The dynamics of light vector fields around slowly rotating and/or Schwarzschild BHs have been studied in a number of previous works, including a combination of analytic and numerical techniques Rosa and Dolan (2012); Pani et al. (2012a, b) and EFT operator matching Endlich and Penco (2016). Our results agree with numerical computations in the close-to-Schwarzschild limit, and explain the scaling of superradiance rates with vector mass, as found in the numerical computations of Rosa and Dolan (2012). Extrapolated to more relativistic bound states, they are also in approximate agreement with time-domain numerical results East (2017); Witek et al. (2013). These results are summarized in Section II.3 and discussed in detail in Appendix A.

The fact that a vector bound state can carry spin angular momentum means that bound states which carry no ‘orbital’ angular momentum can still grow by superradiant extraction of BH spin. The hydrogenic wavefunctions of such bound states have larger density near the BH than superradiant light scalar bound states, since the latter must carry orbital angular momentum. As a result, vector bound state growth rates can be parametrically larger than those for scalars. Consequently, vector superradiance can be fast enough to be astrophysically relevant over a wider range of vector and BH masses. In Section III, we derive constraints on a range of vector masses from existing measurements of rapidly rotating BHs, and describe future signals that may be observed in statistical analyses of BH spin-mass distributions in LIGO data.

In Section IV, we compute approximate GW emission rates from the non-relativistic vector bound states, with further discussion in Appendix C. We estimate the expected number of sources that GW observatories such as Advanced LIGO could observe in continuous wave searches for coherent, monochromatic GWs. The higher density of the zero orbital angular momentum vector bound state near the BH, where the potential is fast-varying, results in a larger high-momentum component able to emit GWs, and gives a parametrically larger emission rate than superradiant scalar bound states. As a result, GW signals from vector bound states could be visible over cosmological distances, from across a large fraction of the observable universe.

In this paper, we consider light vectors with purely gravitational couplings. If a vector has other couplings to itself or to SM matter, our results apply if they are small enough such that gravitational interactions with the BH remain dominant. We leave the question of how large such couplings would have to be to change this story, and whether such large couplings could have other observational consequences, to future work.

II Vector bound states and superradiance

II.1 Spin-1 superradiance

Superradiant scattering of massless spin-1 particles (e.g. the photon) in a Kerr background was calculated by Teukolsky and Press (1974). As for spin-0 particles, the condition for a mode to be amplified is for the angular velocity of the BH horizon to be larger than the angular phase velocity of the wave mode,

[TABLE]

where $\omega$ is the frequency of the mode, $m$ its angular momentum about the BH spin axis, and $\Omega_{H}$ the angular velocity of the BH horizon. This condition on the energy to angular momentum ratio of the wave can be derived from the black hole area theorem; superradiant modes are those whose emission increases the BH area, so their emission is allowed classically Bekenstein (1973).

Rotational superradiance of spin-1 particles can occur around systems other than black holes. Historically, the first example worked out was Zeldovich’s calculation Zeldovich (1971) for superradiant scattering of EM waves from a conducting cylinder. In general, whenever an object is able to absorb waves, and the superradiance condition in Eq. 1 holds (with $\Omega_{H}$ taken to be the angular velocity of the object), it will be thermodynamically favored for the object to spin down by emitting superradiant wave modes. For the gravitationally-coupled vectors we are considering, this means that spinning objects other than BHs, such as neutron stars, planets, etc., will also superradiate. The reason why we focus on BHs is that the amplification rates for superradiant vector modes around other objects will generally be many, many orders of magnitude lower. This is because other objects absorb gravitationally coupled vectors much less efficiently than black holes, and spin much slower than a near-extremal BH of comparable size.

Since there are no gravitational bound states of a massless vector around a BH, superradiant amplification only occurs during a single ‘pass’, and amplifies the incoming energy by at most $4.4\%$ Teukolsky and Press (1974). The ‘black hole bomb’ thought experiment of Press and Teukolsky (1972) evaded this limitation by placing a mirror around the BH. A superradiant mode confined by the mirror can then be amplified exponentially, eventually extracting an $\mathcal{O}(1)$ fraction of the BH’s spin. The addition of a vector mass effectively allows the BH’s gravitational potential to act as a ‘mirror’, confining superradiant modes around the black hole and enabling their exponential growth. The rest of this section will investigate the dynamics of these gravitationally-bound states.

II.2 Bound states

We consider a free massive spin-1 field $A_{\mu}$ with Lagrangian111We use the $-+++$ signature.

[TABLE]

In this paper, we assume a Stuckelberg-type mass for the vector, i.e. there are no additional light degrees of freedom; we comment briefly on the problems associated with a low-scale Higgs mechanism in Section II.6. $A^{\mu}$ obeys the Proca equation of motion,

[TABLE]

where $F_{\mu\nu}\equiv\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}$ is the usual field strength. This implies the ‘Lorentz condition’ $D_{\mu}A^{\mu}=0$. To find the bound states of a massive vector around a BH, we need to solve the Proca equation in the Kerr metric (as we will see, the BH horizon means that there are generally only metastable ‘quasi-bound’ states Page (1976)). In general, this has to be done numerically. However, for vectors with Compton wavelength large compared to the size of the black hole, we expect the $r_{g}/r$ part of the metric to be most important, where $r_{g}\equiv GM_{\rm BH}$. In these circumstances, the bound states are ‘non-relativistic’ and hydrogen-like.

Non-relativistic bound states will oscillate with frequency $\omega\simeq\mu$, where we write

[TABLE]

Quasi-bound states will have complex $\omega$ (when we compare $\omega$ to real quantities, e.g. inequalities such as $\omega<m\Omega_{H}$, we implicitly mean the real part of $\omega$). For $r\gg r_{g}$, we will assume that $\Psi^{\mu}$ varies slowly on scales $\mu^{-1}$, i.e. that its momentum components are non-relativistic, and also that the metric is close to flat. Then, keeping only the terms that are not suppressed by small momenta, the Proca equation becomes

[TABLE]

In the Kerr metric at $r\gg r_{g}$, $g^{00}\simeq-(1-2r_{g}/r)$ so we have a Schrödinger-type equation describing motion in a $1/r$ potential,

[TABLE]

where $\alpha\equiv r_{g}\mu$.

Using the Lorentz condition, $\partial_{t}A_{0}\simeq\partial_{i}A_{i}$, we can solve for $\Psi_{0}$ in terms of $\Psi_{i}$, and find all of the bound states by solving the $\Psi_{i}$ Schrödinger equation. Since the $1/r$ part of the potential is spherically symmetric, we can (at leading order) separate $\Psi_{i}$ into radial and angular functions,

[TABLE]

where the angular functions are ‘pure-orbital vector spherical harmonics’ Thorne (1980), which are eigenfunctions of the orbital angular momentum operator: $-r^{2}\nabla^{2}Y_{i}^{\ell,jm}=\ell(\ell+1)Y_{i}^{\ell,jm}$. The $j$ label corresponds to the total angular momentum of the state, while $m$ is the total angular momentum along the $z$ axis, and $\ell\in\{j-1,j,j+1\}$ is the ‘orbital angular momentum’ (see App. A).222Our notation differs from some existing literature: Rosa and Dolan (2012) mis-identifies $\ell$ as the total angular momentum of the bound states, and Rosa and Dolan (2012); Pani et al. (2012a, b) have $\ell$ and $j$ swapped compared to our usage. Our notation is the same as Endlich and Penco (2016). The basis used in Rosa and Dolan (2012) is the ‘pure-spin vector spherical harmonics’ of Thorne (1980); the linear combinations of these found in Rosa and Dolan (2012) for non-relativistic bound states at long distances are precisely the pure-orbital harmonics.

The equation for the radial wavefunction $R^{n\ell}(r)$ is the same as for the scalar Coulomb problem, so the bound state radial wavefunctions are hydrogenic, labelled by orbital angular momentum $\ell$ and overtone number $n$. The leading-order energy levels are hydrogen-like,

[TABLE]

Figure 1 shows the $A_{i}$ fields at a moment in time, for a set of $j=1$ bound states. These illustrate how the amplitude varies with radius in a hydrogen-like way (controlled by $\ell$), and how the total angular momentum is a combination of spin and orbital components. As we will see in the next section, the fastest-growing bound state around a fast-spinning BH will be the $\ell=0$, $j=m=1$ state with $n=0$. The real time-dependent fields for this state have leading-$\alpha$ form

[TABLE]

[TABLE]

where $a\equiv 1/(\mu\alpha)$ is the characteristic size of the bound state.

II.3 Superradiance rates

If a bound state satisfies the $\Omega_{H}>\omega/m$ superradiance condition, then — assuming that gravitational perturbations to the Kerr background are small — it will grow with time, exponentially increasing in amplitude as it extracts energy and angular momentum from the BH. In particular, this process is possible for any mass $\mu$ (though for large $\mu$, only large-$m$ modes will satisfy the superradiance condition, Eq. 1). However, for superradiance to be important for astrophysical BHs, the growth rate of bound states needs to be fast, at the very least compared to the age of the BH. As in the case of scalar bound states Dolan (2007); Arvanitaki and Dubovsky (2011), it will turn out that the growth rate is a rapidly-varying function of $\alpha$, falling off as a high power of $\alpha$ for $\alpha\ll 1$, and exponentially suppressed for superradiant states with $\alpha\gg 1$.

There are no full analytic solutions for scalar or vector bound states in a general Kerr background, and in most cases, no approximate analytic solution that is close to the true solution at all radii. So, to obtain bound state growth and decay rates, one either needs to solve the wave equation numerically, or match together approximate solutions with different regimes of validity. In Appendix A, we review these different approaches, and perform an analytic matching calculation which gives the leading form of the growth and decay rates in the small-$\alpha$ approximation. Far away from the BH, the bound state wavefunctions are hydrogen-like, as described above. Close to the BH, the mass term in the Proca wave equation becomes a sub-leading correction, allowing us to use the behaviour of massless waves in a Kerr background to determine the energy flux across the horizon. Matching these two regimes together, the leading-$\alpha$ growth rate for a vector bound state scales as

[TABLE]

(where we have also taken $m\Omega_{H}-\omega$ to be a small parameter, to show the behavior near the $\omega=m\Omega_{H}$ crossover point) compared to the scalar form

[TABLE]

In both cases the rate scales as $\alpha^{2\ell+2j+5}$ — the difference in the vector case is that we do not always have $j=\ell$. Roughly speaking, the more localized the bound state wavefunction is close to the BH, and the smaller the total angular momentum barrier it faces to reach the horizon, the larger the growth or decay rate.

The difference between scalar and vector rates has the important consequence that the fastest-superradiating mode for a vector, the $\ell=0,j=m=1$ mode, has leading-$\alpha$ growth rate 333These are the growth rates for the energy density of the bound states (equivalently, for their ‘occupation number’); the growth rate for the field amplitude is a factor $1/2$ smaller, with ${\rm Im}(\omega)=\Gamma/2$.

[TABLE]

where $a_{*}\equiv J_{BH}/(GM_{BH}^{2})\in[0,1)$ is the BH spin, which is related to the horizon angular velocity by $\Omega_{H}=\frac{1}{2}\left(\frac{a_{*}}{1+\sqrt{1-a_{*}^{2}}}\right)\mu\alpha^{-1}$. In contrast, the fastest-growing level for a light scalar, the $\ell=m=1$ level, has $\Gamma\simeq\frac{1}{24}a_{*}\alpha^{8}\mu$. Vector bound states can therefore, in the non-relativistic limit, grow significantly faster than scalar ones, as illustrated in Figure 2.

The strong $\alpha^{6}$ dependence of the superradiance rate implies that, for a given BH, vectors with $\mu\ll r_{g}^{-1}$ will have very small superradiance rates. For $\mu\gg r_{g}^{-1}$, states must have $m\gg 1$ to be superradiant. The ingoing component must therefore tunnel through a large angular momentum barrier, which results in an exponential suppression of the tunnelling rate Zouros and Eardley (1979); Arvanitaki and Dubovsky (2011). Thus, vector masses much larger than $r_{g}^{-1}$ also have very small superradiance rates. Combining these results, we expect astrophysically significant bound state growth rates to occur for vectors whose Compton wavelength is of approximately the same order as the size of the BH.

Figure 2 shows the superradiance rates for the fastest-growing bound states at each $m$ value, around a very-fast-spinning ($a_{*}=0.99$) Kerr BH, illustrating the properties discussed above. For comparison, we also show the superradiance rate for the fastest-growing scalar mode, including both the small-$\alpha$ analytic approximation and numerical results from Dolan (2007). East (2017) performs a full numerical computation of the (time-domain) evolution of vector bound states in a Kerr background. Comparing their results for the fastest-growing vector mode to our analytic calculation, we find good agreement for the superradiant growth rate of the fastest-growing level (to within $\mathcal{O}(1)$ across the entire range).

II.4 Black hole evolution by superradiance

The growth of superradiant bound states can begin spontaneously, without the need for any pre-existing density of the vector field. Though the calculations in this paper are all on the level of classical fields, at the quantum level one can view superradiance as amplifying vacuum fluctuations, with decoherence turning the resulting state into a mixture of coherent, classical-like bound state oscillations. Consequently, while the light, weakly-coupled vector fields we consider could be dark matter candidates, superradiance will occur whether or not there is a pre-existing astrophysical abundance.

If the superradiance time of the fastest-growing bound state — with $j=m=\ell+1$— is significantly smaller than the age of the BH, then its exponential growth will extract enough energy and angular momentum from the BH to bring $\Omega_{H}$ down to the superradiance boundary for that level. For a BH starting with spin $a_{*,0}$, if the cloud extracts a fraction $\Delta a_{*}/a_{*,0}$ of its initial angular momentum, the final occupation number444The bound states grow into coherent classical oscillations of the vector fields, so are not Fock states with a definite occupation number. However, the occupation numbers involved are so large that the fractional variance in occupation number for the coherent states is tiny, so we can sensibly discuss bound state occupation numbers. of the cloud will be

[TABLE]

this behavior is borne out by full numerical time-domain simulations East and Pretorius (2017). Since the angular momentum of fast-spinning stellar-mass BHs is so high, it takes the cloud $\sim\log(N_{m})\sim 180$ e-folds of growth, starting from vacuum fluctuations, to extract an order-1 fraction of the BH’s angular momentum. For small $\alpha$, this corresponds to the extraction of a small proportion of the black hole’s mass, since $\Delta M_{BH}\simeq\mu N_{m}\simeq\alpha M_{BH}\Delta a_{*}/m$, and $\alpha/m<\frac{1}{2}\frac{a_{*}}{1+\sqrt{1-a_{*}^{2}}}$ for a superradiant level. If annihilation rates of vectors to gravitational waves (Section IV) are sufficiently rapid, the cloud may not reach maximum size before annihilations begin to deplete it; however, we do not expect this to be the case in the non-relativistic limit.

After the growth of a given level stops, higher-$m$ modes will still be superradiant. If the fastest-growing of these has a sufficiently small superradiance time, the process will repeat, with the BH being spun down to the next superradiance boundary. Since it takes $\sim 180/m$ e-folds of growth for a bound state level to extract a significant fraction of a roughly-stellar-mass BH’s spin, a BH spends most of its lifetime very near each superradiance boundary, with the transits between them occurring mostly in the last $e$-fold of each level’s growth.

Figure 3 illustrates this evolution for a BH with initial mass $20M_{\odot}$ and initial spin $a_{*}=0.9$, assuming the existence of a weakly-coupled vector of mass $\mu=10^{-12}\text{~{}eV}$ (giving $\alpha=0.15$). The strong dependence of the superradiance rate on $m$ means that, though the $m=1$ and $m=2$ levels grow quickly compared to astrophysical timescales, the $m=3$ level does not have time to grow in the age of the universe, so the BH spends almost its entire lifetime on the $m=2$ superradiance boundary, at $a_{*}\simeq 0.3$ and $M\simeq 18.5M_{\odot}$. The solid curves plotted in Figure 6 are the isocontours of $\Gamma_{m}=\frac{1}{\tau_{bh}}\log N_{m}$, where $\tau_{bh}$ is the lifetime of the BH (here taken to be $5\times 10^{6}\text{~{}yr}$). If a black hole is born above these curves, then it will be spun down by the growth of successive superradiant levels until it reaches a superradiance boundary on or below the curves. If it is born below the curves, it will be unaffected by superradiance. Figure 6 illustrates the results of these processes on a selection of initial BH masses and spins, showing how superradiance depopulates the region above the superradiance rate contours, and leaves a population of BHs scattered along the superradiance boundary curves (with the slight complication that the figure considers BHs in binary systems, where a sufficiently close companion can disrupt superradiance, as discussed in Section II.5).

As we will see in Section IV, bound states are depleted by annihilations to gravitational radiation once they have grown to large occupation number. In most circumstances, the difference between superradiant growth rates for successive $m$ is large enough that, on a timescale short compared to the growth timescale of the next level, the previous level’s occupation number has mostly annihilated away. Once the next level has extracted enough spin from the BH, the previous level goes from being almost on the superradiance boundary, to decaying into the BH. The result is that the small occupation number that did not annihilate away falls back into the BH, delaying the spin-down process very slightly (though not making a significant impact).555This pattern of levels sequentially growing from, and then falling back into, a BH, is explored by Bosch et al. (2016); Sanchis-Gual et al. (2016) for the case of charged field superradiance of a Reissner-Nordstrom BH.

II.5 Perturbations

The presence of a perturbation to the BH gravitational potential can affect superradiance by mixing a growing level with a decaying one. In particular, a binary companion sources a dipole-type gravitational potential (to leading order in the ratio of BH size to companion distance), which can mix levels which have a small energy splitting and different angular momentum numbers.

For vectors, the $j=1,\ell=0$ level — which dominates the phenomenology — is especially robust against perturbations because there are no bound states that decay with a rate faster than its superradiance rate. Furthermore, it is the state with the lowest energy, with a large $\mathcal{O}(\alpha^{2})$ energy difference from other levels, which further suppresses the mixing. For example, for binary systems the largest mixing term is the dipole term with $\Delta j=\Delta\ell=1$. So, at leading order, the fastest-growing $j=1$ level can only mix with the $N=2$, $j=0$, $\ell=1$, $m=0$ state, where $N=n+\ell+1$ is the principal quantum number. Because these states have different energies ($\Delta N\neq 0$), the mixing is further suppressed, and an external perturbation has to be close to order one to significantly affect whether the state is superradiant. For typical binary systems, this translates into the condition that the cloud is sufficiently tightly bound, $\alpha\gtrsim 10^{-3}$, and does not affect the signals or constraints. More details for this estimate, as well as perturbations of $j>1$ levels, are covered in Appendix B.

II.6 Effects of interactions

The BH evolution described above only takes into account the gravitational interactions of the light vector. The vector may also have non-gravitational interactions with SM matter. As long as such couplings always have sub-dominant effects to the gravitational potential of the BH, superradiance proceeds broadly as we have described.

There is also the possibility of the vector having non-gravitational interactions with other non-SM states — for example, if it gets its mass via a Higgs mechanism. Such states will induce self-interactions of the vector, and if the superradiant cloud reaches a high enough density, may be produced on-shell (in analogy to Schwinger pair production). The energy density of a cloud that has extracted $\mathcal{O}(1)$ of a BH’s spin is

[TABLE]

so, for the purely gravitational interactions to dominate, the new states must be sufficiently heavy and/or weakly coupled. If there is some astrophysical abundance of states which interact with the vector, their interactions with the cloud could also change the evolution around BHs.

III Black hole spin measurements

III.1 X-ray Binaries and Supermassive Black Holes

In recent years, an increasing number of stellar BH Miller and Miller (2014) and supermassive BH (SMBH) spins Reynolds (2014) have been measured. The stellar BHs are in binary systems with a companion star. For both stellar BHs and SMBHs, their accretion disks are used to extract information about their spin. The methods used fall into two categories: continuum fitting McClintock et al. (2014) and X-ray relativistic reflection Reynolds (2014), which measure properties of the accretion disk to find the radius of the innermost stable orbit and extract the BH spin. Because a BH that satisfies the superradiance condition can lose its spin quickly on astrophysical time scales, these spin measurements place new limits on vector masses.

Similar to Fig. 3, Fig. 4 shows regions where superradiance of the $j=1,...,5$ levels for a $\mu=10^{-11}$ eV vector will spin down a BH. Each region satisfies the superradiance condition, and also has a superradiance rate fast enough to grow a maximally-filled cloud within the relevant BH time scale, $\tau_{bh}$; $\tau_{bh}$ varies between systems so the regions shown in Figure 4 are approximate. $\tau_{bh}$ is the shortest time scale on which superradiance can be disrupted: for stellar BHs, we use the shorter of the age and the shortest timescale on which the spin can change by accretion, while for SMBHs, the accretion timescale is the relevant one. We conservatively use $\tau_{bh}<\tau_{\text{Salpeter}}/10$, to account for possible periods of super-Eddington accretion, where $\tau_{\text{Salpeter}}\sim 4.5\times 10^{7}$ years Shankar et al. (2009); Shapiro and Teukolsky (1983).

Figure 4 also shows BH spin and mass data with $2\sigma$ errors.666For the more slowly spinning BHs, we do not have a reliable 2$\sigma$ estimate so we assume a Gaussian error on the mass and spin for visualization. A BH excludes the vector mass if it lies in the shaded region, within experimental error. The faster superradiance rates of vectors compared to scalars allows more BHs to be useful in excluding vector masses. The two most rapidly spinning BHs clearly exclude this mass for the vector particle.

As the mass of the vector increases, the affected regions shift to the left, towards lower BH masses. For a given $j$ level and a given BH data point, the superradiance condition is no longer met at some maximum vector mass. Decreasing the mass of the vector causes the affected regions to shift to the right. The minimum BH mass for exclusion is determined by requiring that the superradiance rate remains fast enough to grow a maximally-filled cloud within the relevant BH time scale,

[TABLE]

In the presence of non azimuthally-symmetric perturbations to the cloud, e.g. from the accretion disk or the companion star, a given level can mix with decaying bound states, which can relax the constraint at low vector mass Arvanitaki et al. (2015). However, for the stellar mass BH systems we consider, the perturbations are a sub-leading effect, as discussed in Appendix B.

The resulting constraints are shown in Figure 5. In Figure 5 (right), each stellar BH places a limit on the vector mass. The overlapping rectangles correspond to different $j$ levels which set a constraint. For example, the BH Cyg X-1 places limits on vector masses with the $j=1,...,5$ levels.

For SMBHs shown in Figure 5 (left), we conservatively use the $j=1$ level to set constraints. The $j>1$ levels may place an additional constraint at higher vector masses. However, the density of the accretion disk grows quickly with BH mass, so the effect of accretion disks around SMBHs on the cloud is relatively larger than those around stellar mass BHs. Although the mass in the accretion disk is small compared to the BH mass, estimating the size and form of the perturbation relies on details of the disk Arvanitaki et al. (2015). The $j=1$ level is robust against small perturbations, since the decay rate of the fastest decaying bound state also scales as $\sim\alpha^{7}$, similar to the $j=1$ superradiance rate; we save the extension of limits using higher levels to future work.

The constraints set using spin measurements depend on our computation of the superradiance rates (see Appendix A) and thus have some theoretical uncertainty. In Figure 5, we indicate the uncertainty with gray bands. The right and left edges of these gray bands are set by superradiance rates $\frac{1}{2}\times$ and $2\times$ the calculated value. The right edge of each $j$ level is cut off by the superradiance condition; one source of uncertainty is the energy of the bound state that goes into the superradiance condition. We expect the uncertainty in the energy to be less than the lowest order correction in $\alpha$, $\delta\omega\sim\frac{\alpha^{2}}{2N^{2}}\mu$.

In Table 1, we list details about the stellar BHs used to set limits in Figure 5. In addition to their spins, which are determined by both methods, these BHs have precise measurements for their masses, as well as estimated ages for the binary systems. They set limits on the vector mass in the range $5\times 10^{-14}$ eV $<\mu<2\times 10^{-11}$ eV. The lower limit is about an order of magnitude better than that for scalars ($6\times 10^{-13}$ eV – $2\times 10^{-11}$ eV) Arvanitaki et al. (2015), by virtue of the faster superradiance rates for vectors and weaker level-mixing of the leading $j=1,\ell=0$ bound state.

Details about the SMBH masses and spins that we used to set limits in Figure 5 are listed in Table 2. We assume approximately Gaussian errors to estimate the $90\%$ interval on the masses. We use one-tenth of the Salpeter time to set the constraints, ensuring that superradiance is the fastest process affecting the BH spin evolution. The SMBHs exclude the mass range $6\times 10^{-20}$ eV $<\mu<2\times 10^{-17}$ eV (with the exception of a small gap around $1\times 10^{-19}$ eV).777This mass range is slightly less constraining than found in Pani et al. (2012b); we use superradiance rates for the $j=1,\ell=0$ mode (while Pani et al. (2012b) use the slower $j=1,\ell=1$ mode for the conservative bound), take into account the $\sim 180$ e-folding times needed to grow the cloud, and use an updated set of BH measurements (which includes a new BH Reynolds (2014), but has a smaller spin for Fairall 9 than used in Pani et al. (2012b)). In addition, we assume purely gravitational interactions, while Pani et al. (2012b) consider the SM photon with a mass, the dynamics of which would be strongly affected by the plasma surrounding a SMBH. We consider these limits preliminary since the properties of SMBHs are less well-known, and the spin measurements only come from X-ray reflection Reynolds (2014).

As more BHs are measured with better precision, we expect the limits to extend further.

III.2 Binary black hole mergers

Black hole spins have only recently been measured in another way: via the GW signals from binary BH (BBH) mergers, as observed at aLIGO Abbott et al. (2016a). Compared to the spin measurements in X-ray binary systems, the GW measurements of the two pre-merger spins have very large errors, (roughly $\sigma_{a_{*}}\gtrsim 0.3$ Abbott et al. (2016a)). However, as the sensitivity of gravitational wave observatories improves, they will detect an increasing rate of such mergers. Advanced LIGO may build up a catalogue of hundreds of spin measurements — at design sensitivity, aLIGO is expected to detect 80-1200 BBH merger events per year of data Abbott et al. (2016b, c, a). This raises the possibility that, while individual measurements are not very informative, the statistical properties of the whole set may provide evidence for the altered BH mass vs. spin distribution predicted by superradiance (as proposed for spin-0 superradiance signals in Arvanitaki et al. (2017)).

To give an example, Fig. 6 shows the effects of BH spin-down through a weakly-coupled vector of mass $10^{-12}$ eV, assuming (for illustrative purposes) that BHs are born with a uniform spin distribution. As described in Section II.4, black holes that are born with high spin end up on one of the superradiance boundary lines. For levels with $j>1$, the mixing due to the companion BH may be significant, and can limit the growth of the cloud at small $\alpha$ (see Appendix B); this effect is included in the curves of Fig. 6.888Since bound levels will, in general, annihilate away most of their energy to gravitational radiation (Section IV) before the BHs merge, we do not need to worry about the gravitational effects of the clouds once the BHs get close to each other. This neat picture is spoiled by the large spin measurement errors, which blur the measured distribution (Fig. 6, right panel) — however, there is still an overdensity along the superradiance boundary lines, and an under-density above the superradiance rate isocontours.

In reality, it is unlikely that BHs are born with a uniform spin distribution. We also have no way of knowing the formation history for a given BH binary; any given pair may not have been able to superradiate, due to merging too quickly or forming too close together. A proper estimate of the statistical detectability of superradiance-induced spindown would, in addition to taking into account detector-related measurement errors, need to scan over an astrophysically-plausible ranges of BH formation histories and birth distributions in the mass-spin plane.

We will not attempt such an estimate here. Instead, we show projections for the number of events required to detect that the BH spin distribution varies with mass (as per Fig. 6), without the statistical test assuming a particular form for the spin distribution. Such variation could come from astrophysical mechanisms, as well as from superradiance. Our estimates simply indicate the number of events required to discover mass-dependent structure in the BH mass-spin distribution — if this were to be seen, a detailed study of possible astrophysical explanations would be required. In constrast, if we had strong priors on the astrophysical processes involved, then significantly fewer events may be necessary to point to superradiance; for example, simply seeing a deficit of high-spin BHs may suffice.

Figure 7 shows our estimates for the number of events required to detect mass-dependent structure in the BH spin distribution. We have assumed a uniform initial spin distribution Vitale et al. (2014); *Vitale, and a power-law BH mass distribution $\rho(M)\propto~{}M^{-2.35}$ Abbott et al. (2016b), as were assumed in LIGO analyses. Our measurement error estimates are based on studies of intermediate mass BBHs; at design-sensitivity LIGO/Virgo detectors, one expects to obtain a 90% confidence interval of width $|\Delta a_{*}|<0.8$ for total masses up to   $600M_{\odot}$,999For the most pessimistic case of equal BH masses and misaligned spins; better measurements are possible for dissimilar masses or aligned spins. and a $\sim 10\%$ error in mass determination for an order-one fraction of primary BH masses Vitale et al. (2016). These measurement errors are broadly comparable to those in the spin and mass posterior distributions published by the LIGO collaboration for the three BBH events observed in the O1 aLIGO run Abbott et al. (2016a).

The ‘gap’ at intermediate masses in Figure 7 occurs because superradiance there is too efficient — if such a vector existed, practically all BHs in the observable mass range would be spun down to very small spins. As a result, there would be no mass-dependent structure in the BH spin distribution, and it would not be distinguishable from all of the BHs in such binaries having been born with small spin. Of course, this would be interesting in itself. Figure 7 illustrates that there is a strong possibility of observing structure in the BH mass-spin distribution for vector masses an order of magnitude below the exclusion from the X-ray binary spin measurements (Section III), after $\sim$ hundreds of observed events. Such a statistical test would be an interesting hint of a light vector’s presence, and could not on its own easily be used to exclude or confirm it.

As shown in Fig. 7, reducing the assumed merger time from $10^{10}\text{~{}yr}$ to $10^{7}\text{~{}yr}$ (the range suggested by BBH formation models Abbott et al. (2016c)) reduces the sensitivity of the statistical search at small vector masses, since superradiance becomes too slow to spin down low-mass BHs, and the effect of perturbations for $j>1$ levels is increased. Conversely, third-generation gravitational wave observatories may achieve much higher signal-to-noise BBH merger detections, and consequently much higher precision progenitor spin measurements (e.g. a $90\%$ interval of $|\Delta a_{*}|<0.1$ for a majority of high-SNR events Vitale and Evans (2016)). Along with the increased event rate, the higher precision would improve the reach of the statistical search, and enable detailed investigation of any features found by Advanced LIGO.

IV Gravitational wave signals

So far, we have considered the evolution of a massive vector field on a fixed background. However, since superradiance can extract $\mathcal{O}(\alpha)$ of the BH’s energy into the bound states, the gravitational effects of the bound states may not be negligible. In particular, since superradiance builds up a coherent classical wave in the bound states, it results in an oscillatory stress energy tensor, which sources gravitational radiation. As we will see in this Section, this gravitational radiation provides one of the most important potential observational signatures of a superradiant cloud. Such signatures were analyzed for the case of scalar bound states in Arvanitaki et al. (2010a); Arvanitaki and Dubovsky (2011); Arvanitaki et al. (2015, 2017). Here, we will discuss the corresponding signals for vectors.

Decomposing $A_{\mu}$ into a combination of oscillations in different bound state levels,

[TABLE]

the stress energy tensor is

[TABLE]

where we have only displayed the $A_{\mu}A_{\nu}$ term to demonstrate the schematic form of $T_{\mu\nu}$. The components of $T_{\mu\nu}$ with quadrupole or higher angular dependence, and non-zero time dependence, source gravitational radiation in the usual way. At the particle level, the high-frequency $\omega_{i}+\omega_{j}\simeq 2\mu$ terms correspond to annihilations of two bound-state vectors to a single graviton, while the $\omega_{i}-\omega_{j}\sim\alpha^{2}\mu$ terms correspond to transitions from a higher-energy level to a more deeply bound one.

Both of these processes result in monochromatic gravitational waves, emitted in a characteristic angular pattern from the BH, at a frequency related to the mass of the vector. Gravitational wave observatories, present and future, can search for such signals.

IV.1 Transitions

Transition signals will be sizable only in the case where two levels are populated simultaneously at nearly maximum cloud sizes. Since the occupation number attained by slower-growing levels is exponentially suppressed by the difference between their superradiance rates and that of the fastest-growing level, this only occurs for levels where the superradiance rates of the ‘excited’ and ‘ground’ states are similar. That is, we need $\Gamma_{g}\sim\Gamma_{e}$, where the excited state is the one with a larger real energy component, $\omega_{e}>\omega_{g}$; furthermore, the peak signal strain scales and the superradiance rate Arvanitaki and Dubovsky (2011); Arvanitaki et al. (2015).

The frequency of the emitted radiation is suppressed by the $\alpha^{2}/N^{2}$ splitting between levels, so its wavelength is long compared to the cloud size. As reviewed in Appendix C, this means that it is simple to reliably estimate the rate of GW emission due to transitions. This is approximately given by the quadrupole formula, so is dominantly affected by the size of the cloud, controlled by $\ell$. For the vector case, the fastest-superradiating level at any time is one with $j=\ell+1$, with growth rate $\Gamma\sim\alpha^{4\ell+6}\mu$. In comparison to the scalar case, where the fastest-growing level has $\Gamma\sim\alpha^{4\ell+4}\mu$, a vector bound state with equivalent transition rate (i.e. the same $\ell$) has a smaller superradiance rate than its scalar equivalent. For a given vector mass, we therefore expect to observe fewer transition signals than we would for a scalar of the same mass. Since Advanced LIGO has marginal sensitivity to scalar transition signals, with optimistic assumptions giving $\sim 0.3$ detectable signals Arvanitaki et al. (2015), we do not expect it to be sensitive to transition signals from vector superradiance.

Nevertheless, it is possible that at large $\alpha$ the superradiance rates have a different dependence on $n$ than the form we can extract in the $\alpha\ll 1$ limit. For example, in the scalar case, numerical results Yoshino and Kodama (2015) seem to indicate that the $n=0$ level for $\ell=m=3$ has a superradiance rate that crosses over with those for other $n$ at large $\alpha$ and $a_{*}$. This may result in other levels giving rise to transition signals, which could be more promising for observational purposes. We leave the investigation of such issues to future work.

IV.2 Annihilation rates

Annihilation signals have the advantage that they only require a large occupation number in a single level, so can occur without the coincidence of superradiance rates required for transitions. However, as reviewed in Appendix C, the annihilation rates are parametrically smaller than transition rates, and are rather more subtle to evaluate. The momentum of the emitted gravitational wave is $\sim 2\mu$, which is much larger than the ‘typical momentum’ in a bound state, $\sim\alpha\mu$; the high-momentum behaviour of the bound state is dominated by the discontinuity from the $e^{-r/a}$ wavefunction dependence near the origin. This gives an emission rate that scales as a higher power of $\alpha$ compared to a quadrupole formula estimate, and also means that the small-$r$ distortion of the gravitational wave in the BH background can affect the rates at the same order as the flat-space estimate.

In the flat-space approximation, we obtain an emitted GW power of $P\simeq\frac{32}{5}\alpha^{12}\frac{GN^{2}}{r_{g}^{4}}$ for the fastest-growing $\ell=0$, $j=m=1$ level. As discussed in Appendix C, and similarly to the scalar case, it is necessary to take into account the modified propagation of the gravitational wave due to the $r_{g}/r$ part of the metric in order to obtain the leading-$\alpha$ rate. We estimate that this can increase the power to $\sim 10\times$ larger than the flat-space result,

[TABLE]

This compares to a GW power of $\sim 0.02\alpha^{16}GN^{2}r_{g}^{-4}$ for the fastest $\ell=m=1$ scalar bound state Brito et al. (2015a). For both scalars and vectors, the GW power emitted from a bound state scales as $\alpha^{4\ell+12}GN^{2}r_{g}^{-4}$, with the faster rates possible for vectors reflecting the existence of a $\ell=0$ superradiant level.

We use Eq. (19) as the central value for reach and event rates and take the difference between the flat- and curved-space rates as an estimate of the theoretical uncertainty in the event rates plots. Numerical simulations of the GW emission rate from the fastest-growing level East (2017) approximately agree with our estimate, though they do not extend to small enough $\alpha$ to see the $\alpha^{12}$ scaling in the $\alpha\ll 1$ limit.

The angular distribution of GW emission is set, at leading order in $\alpha$, by the form of the $T_{ij}$ source. For a $j=m=1$ state, this corresponds to an emission profile that is the same as for a standard quadrupole GW source, such as a spinning neutron star with an asymmetry. The emitted power per solid angle is

[TABLE]

where $\theta_{k}$ is the angle of the line-of-sight from the BH’s spin axis. While the power varies by a factor of a few across different angles, this does not make a qualitative difference to the phenomenology, and we ignore this variation in the following sections.

IV.3 Annihilation signals at LIGO

The monochromatic GW signals from bound state annihilations around stellar-mass BHs are in the right frequency range to be detected at Advanced LIGO. A vector of mass $\mu$ results in a GW signal at frequency

[TABLE]

In particular, whatever the mass of the BH, the GW signals from a given vector will be clustered within the angular frequency interval $\sim(1.7\mu,2\mu)$; numerical results indicate that this lower bound is correct even for relativistic bound states East (2017). Writing the emitted GW power as $P\equiv\Gamma_{\rm ann}N^{2}\omega_{\rm GW}$, the occupation number $N$ of the cloud evolves as $\dot{N}=\Gamma_{\rm sr}N-\Gamma_{\rm ann}N^{2}$. When the annihilation term dominates, as it generally will after the cloud has reached maximum size and spun down the BH, $N$ decreases as

[TABLE]

where $N_{0}\equiv N(t=0)$. The observed strain therefore decreases like $1/t$ for $t\gg(N_{0}\Gamma_{\rm ann})^{-1}$, with the timescale of the signal set by

[TABLE]

Hence, signals can last anywhere between seconds and millions of years, depending on the value of $\alpha$.

The mass of the cloud itself contributes to its own gravitational binding energy. As the cloud annihilates away, this contribution to the binding energy decreases, leading to a positive frequency drift of the annihilation signal,

[TABLE]

Monochromatic searches may become difficult with large frequency drifts. On the other hand, a positive frequency drift is quite unusual for an astrophysical source and could be used to distinguish the signal from e.g. monochromatic radiation from rotating neutron stars.

IV.4 Annihilation reach

The rotation and motion of the Earth complicates the search for narrow-band GW signals at terrestrial GW observatories. A monochromatic signal from a distant source is subjected to a time-dependent redshift in the detector frame, with this time dependence set by its sky location. As a result, it is not possible to find such signals by simply Fourier-transforming the strain data; instead, it is necessary to scan over each possible frequency and sky position for the source Wette (2009). This is a computationally demanding process, and if the source is not well-localized, it is generally too expensive to integrate coherently over the entire observation time. Instead, LIGO searches Abbott et al. (2016e) break up the observation time into shorter segments and integrate coherently within each segment, combining the results from each segment incoherently. This results in a strain sensitivity (for given desired signal-to-noise ratio) of

[TABLE]

where $N_{\rm seg}$ is the number of segments, $T_{\rm coh}$ is the length of each segment, $S_{h}$ is the detector’s noise spectral sensitivity Amaro-Seoane et al. (2012), and $C_{\rm tf}$ is the ‘trials factor’, which quantifies the look-elsewhere effect inherent in searching over many bins in frequency and sky position. This should be compared to a GW signal strain at the detector of

[TABLE]

for a gravitational wave source with power $P$ and angular frequency $\omega_{\rm GW}$ at a distance $r$ away from the Earth.

Figure 8 shows the maximum distance at which Advanced LIGO (at design sensitivity) could detect annihilation signals from bound states around a BH with $a_{*}=0.9$. This distance is set by the strength of annihilations from the $\ell=0,j=m=1$ bound state, which is the fastest-annihilating superradiant state. The reach is calculated assuming a large dataset of 121 segments of 250 hour coherent integration times, and trials factor of 20; the most recent search uses 90 segments of 60 hour coherent times Abbott et al. (2016e).

As Figure 8 illustrates, the existence of the $\ell=0$ superradiant state means that vector annihilation signals can be seen from much further away than those from scalar states. The converse of this is that the signal lifetimes for vector annihilations are correspondingly shorter (since they radiate away approximately the same overall energy from the cloud).101010The brightest signals are typically shorter than the coherent integration time $T_{\text{coh}}$. The reach is thus penalized by a shorter integration time. Consequently, increasing the annihilation power increases the signal strain but reduces the detector strain sensitivity by the same amount, so the reach at fixed $\alpha$ stays the same even if we account for uncertainty in the annihilation rate by varying its value from $\times 10^{-1}$ to $\times 10^{1}$. Because the reach for $j=1$ annihilations extends to cosmological distances, we account for the signals red-shifting due to Hubble expansion. The detected frequencies are reduced compared with the source frequency, $f_{\text{observed}}=(1+z)^{-1}f_{\text{source}}$, and signal times are increased by the same factor. Between $10^{3}$ - $10^{4}$ Mpc in luminosity distance, the signal frequencies are $1.2$ - $2.5$ times lower. The effect of redshift changes the peak of the signal reach toward higher values of $\mu$.

IV.5 All-sky searches

Advanced LIGO performs all-sky searches for monochromatic gravitational radiation from unknown sources (e.g. Abbott et al. (2016e)), and vector annihilation signals provide an excellent target for such searches. Figure 9 shows our estimates for the number of events that might be observed by Advanced LIGO, using the assumptions detailed in Appendix D. To account for the large uncertainties regarding the population of astrophysical BHs, we show a band bounded by optimistic and pessimistic estimates for BH mass and spin distributions, as well as their formation rates. We also account for theoretical uncertainties in our estimate of the annihilation rates by including a range of annihilation power a factor of 10 above and below the estimate in Eq. (19). The resulting uncertainty is much smaller than the uncertainty in astrophysical populations. The non-detection of annihilation signals in an all-sky search at LIGO can in principle constrain vector masses; however, most of the parameter space has already been excluded by BH spin measurements.

While, as illustrated in Figure 8, sufficiently strong annihilation signals could be visible from cosmological distances, such signals last for only a short time. At smaller vector masses, for which strong signals only occur around heavy BHs (i.e. at $\alpha$ not too small), such signals are probably rare enough that it is unlikely a suitable BH was born at precisely the right time for its annihilation signals to reach us now (though as we discuss below, BBH mergers may give rise to a suitable population of heavy, fast-spinning BHs). Consequently, the blind search event rates are dominated by longer, quieter signals originating from inside our galaxy. A typical such signal, for a vector mass of $\mu\sim 10^{-13}\text{~{}eV}$, comes from a BH of mass $\sim 20M_{\odot}$, corresponding to a small $\alpha\sim 10^{-2}$, and consequently a very long signal time $\sim 10^{8}\text{~{}yr}$.

At larger vector masses ($\gtrsim$ few $\times 10^{-13}\text{~{}eV}$), lighter and more abundant BHs are able to give strong annihilation signals, so the expected rate for these is large enough that extra-galactic signals dominate the blind search. The typical frequency drift of the short, bright signals outside the Milky Way (MW) is $\sim 10^{-7}$ Hz/s, which falls well outside the current range covered by the aLIGO search; in Fig. 9, we only include signals with frequency drift smaller than $2\times 10^{-11}$ Hz/s (see App. E). With this constraint, virtually all events from beyond the MW are cut off for having a too-large frequency drift, leaving only events from the MW. Vector masses for which such distant events could dominate the expected event rate are disfavored by X-ray spin measurements.

Finally, given a single vector particle, we expect the fractional variation in signal frequencies to be at most $\sim 17\%$ (ignoring redshift), set by the fractional binding energy, $\sim\frac{\alpha^{2}}{2(n+\ell+1)^{2}}+\mathcal{O}(\alpha^{3})$ East (2017). This clustering of multiple monochromatic signals would be another handle pointing toward a particle source of the GW signal. For signals from the MW, the typical signal has $\alpha\ll 0.1$, so the monochromatic lines would appear very close together, within $\delta f/f<10^{-3}$.

In addition to monochromatic signals from individual BHs, there will also be a ‘stochastic background’ from signals too weak to be individually resolved. As for the resolvable signals, at low vector masses we expect this stochastic background to be dominated by signals from BHs within the galaxy. Due to the very good frequency resolution of GW detectors such as Advanced LIGO, individually-resolvable signals will stand out well above such a background, but whether such a background may be observable warrants further investigation.

IV.6 Follow-up Searches

The blind search described above fits well into the existing LIGO search strategy. Another interesting category of search is to look for monochromatic GWs from the product of a BBH or BH – neutron star (BH-NS) merger, as proposed for scalar superradiance in Arvanitaki et al. (2017). At design sensitivity, Advanced LIGO may observe a number of BH-NS mergers, and up to hundreds of BBH mergers a year Abbott et al. (2016c); Belczynski et al. (2015), each of which could result in the perfect source for superradiance signals: a new, isolated, rapidly spinning BH. If its spin is high enough, and there is a light vector in the right mass range, the new BH will superradiate a cloud of vectors which will subsequently annihilate into gravitational radiation. The cloud will begin to grow once the final BH closely approximates the Kerr metric; this happens quickly since the ringdown time is parametrically short compared to the superradiance time. From the properties of the new BH, and the following monochromatic GWs, it would be possible to check that such a signal was compatible with superradiance.

Figure 10 shows the expected number of ‘follow-up search’ events that could be seen at Advanced LIGO and future GW observatories. The signals which dominate the event rates are very bright, visible up to the edge of the reach sensitivity (Fig. 8), corresponding to redshifts of $z\sim 1-2$ with aLIGO design sensitivity. On the other hand, the signals are quite short (see Eq. 24), and fall outside of the current LIGO search strategy Abbott et al. (2016e). The signals also have a large frequency drift of up to $\sim 10^{-5}\text{~{}Hz}/\sec$, which may make a search more computationally challenging Wette (2009) (see App. E). The events visible at vector masses below the region constrained by X-ray binaries rely on the presence of a population of heavy BHs ($M>100M_{\odot}$). We include an estimate of the uncertainty by using a power-law mass distribution with a varying exponential suppression of BH populations at high mass ($M>60M_{\odot}$) Belczynski et al. (2015). Recent aLIGO searches put the first upper bound on merger rates of intermediate-mass BHs, which is $\mathcal{O}(1)$ consistent with the upper end of the power-law mass distribution at $M\sim 100M_{\odot}$ Abbott et al. (2017); the BBH mass function will be better constrained with further measurements, narrowing the range of expected event rates. In the case of BH-BH mergers, Advanced LIGO design reach for the merger signal and the follow-up signal are comparable; for BH-NS mergers, the rates may be limited by the reach to the merger, which we take to be $0.927$ Gpc Abadie et al. (2010).

For both types of mergers, the event will be localized on the sky, and the follow-up search can be focused on a particular region. For BH-BH mergers there is likely no electromagnetic counterpart, so the localization is poor Abbott et al. (2016f). For BH-NS events, there may be an electromagnetic counterpart, which would localize the event well enough to eliminate the need to scan over sky location. We thus assume longer integration times for the BH-NS follow-up search, although the signal length limits the overall coherence time for the brightest signals Wette (2009).

GW150914, the first observed BBH merger, resulted in a final BH mass of $62.2^{+3.7}_{-3.4}~{}M_{\odot}$ and spin of $0.68^{+.05}_{-.06}$, at a luminosity distance of $440^{+160}_{-180}$ Mpc Abbott et al. (2016g). For an optimal vector mass of $\sim 2\times 10^{-13}\text{~{}eV}$ the following annihilation signal would have been visible from as far as $\sim 1$ Gpc; however this mass is excluded by X-ray binary measurements. A lighter vector with mass outside the exclusion would only be visible at distances of less than $\sim 1$ Mpc. As shown in Fig. 8, a rapidly spinning final BH would be visible for a wider range of vector masses.

As discussed above, although there is theoretical uncertainty in the GW annihilation rates, this does not translate into a large uncertainty in the event rates: the larger strain is compensated by shorter signal times, and therefore shorter coherent integration times. If the annihilation rate is more than a factor of a few larger than Eq. (19), the cloud does not reach maximum size, which leads to smaller than expected signal power.

In Fig. 10, we assume the binary formation mechanism does not allow for superradiance, and estimate the spin of the final BH using Buonanno et al. (2008) starting with an initial flat distribution of spins. For the BBH merger, if the two initial BHs have spin close to zero, the final BH will have spin $a_{*}\sim 0.68$ for equal mass BHs Buonanno et al. (2008); an initial distribution with zero spin reduces the number of BBH follow-up events by a factor of $\sim 3$.

V Discussion

Black hole superradiance provides a unique probe of new light, weakly-coupled bosons. In this paper, we have derived constraints and estimated future signals for a new, light, gravitationally-coupled vector particle. Existing X-ray measurements of spins of stellar BHs set stringent constraints on the parameter space, excluding nearly three orders of magnitude in mass below $2\times 10^{-11}$ eV; if confirmed, supermassive BH spin measurements would constrain another three orders of magnitude at smaller masses. Continuous GW signals from annihilating vector clouds could be bright enough to be seen across cosmological distances; these signals can be seen in Advanced LIGO all-sky searches, as well as in follow-up searches after BBH and BH-NS merger events. The strong spin constraints, combined with poor LIGO sensitivity at low frequencies, leave a range of a factor of a few in vector mass that could be observed through GWs.

Superradiance may spin down a significant fraction of the pre-merger BHs, giving a statistical signal. The statistical search has the additional benefit that it only depends on measuring the broadband burst from the BBH merger, extending its sensitivity to masses an order of magnitude below X-ray spin constraints. The lack of aLIGO sensitivity to frequencies below $\sim 10\text{~{}Hz}$ means that GW searches cannot probe vector masses too far below the X-ray binary bound. As discussed in Section III.2, as long as high-mass BBH events are seen, detector simulations indicate that reasonably accurate spin measurements should be possible — the question then becomes how many of the merging BHs have such high masses.

As we have explored in this paper, the observational signatures of light vector superradiance around astrophysical BHs are analogous to those of light scalars. The parametrically faster superradiance and GW emission rates of vector bound states alter the relationships between different signals and constraints. For scalars, the all-sky monochromatic GW search for superradiant clouds is the most promising observational route Arvanitaki et al. (2017). For vectors, the stronger exclusions coming from spin measurements in X-ray binaries mean that there is only a narrow mass window in which such searches could discover a light vector (Figure 9). The strong constraints derived in this paper give one handle to distinguish vector and scalar signals: observing monochromatic radiation pointing to a particle of mass in the range $5\times 10^{-14}-6\times 10^{-13}$ eV would suggest scalar rather than vector superradiance. In addition, the vector $j=1,\ell=0$ mode is more tightly bound to the BH, so while for scalars we expect a spread of $\lesssim 4\%$ in frequency below $2\mu_{a}$ due to BH binding energy, vector monochromatic lines will cluster within a broader range, of up to $\sim 17\%$ (not taking into account the red-shifting of signals). This difference will only be apparent for large-$\alpha$ signals, such as those giving rise to extra-galactic sources. In all-sky searches where weak, low-$\alpha$ signals dominate, the expected frequency differences between different sources are much less than $1\%$ for both scalars and vectors.

Follow-up searches for newly formed BHs are more promising for vectors than scalars at current GW observatories, due to the larger reach. Signals could be seen from distances beyond a Gpc — however, their short duration and relatively large positive frequency derivative requires a dedicated LIGO analysis. Measuring the parameters of the newly formed BH determines, for a given particle mass and spin, the expected superradiance and annihilation rates. It would therefore be possible to check whether any follow-up GW observations were consistent with vector or scalar superradiance.

Statistical searches for structure in the mass/spin distribution of merging BBHs are similarly powerful for scalars and vectors, through the more stringent X-ray spin limits on vectors mean that sensitivity in the non-excluded region relies on the observations of heavier BBH mergers. The existence of the superradiant $\ell=0$ mode for a vector means that BH spin-down can be significantly more efficient. Because of slower superradiance rates and a larger impact of companion perturbations, we do not expect scalar superradiance to reduce BH spins close to zero on astrophysical timescales, while this could be possible for vectors.

Our discussion of GW signals was limited to those from stellar-mass BHs. Lower-frequency observatories, both space- and ground-based Dimopoulos et al. (2008); Graham et al. (2013); Audley et al. (2017), could extend the reach to higher black hole masses, and therefore lower vector masses. Once lower-frequency GW detectors such as LISA Audley et al. (2017) are operational, mergers involving SMBHs will become observable, as well as possible monochromatic signals from low-mass vector clouds. Estimating the sensitivity of LISA to monochromatic signals from vector clouds, we find that, if clouds carrying a significant fraction of a SMBH’s spin can be built up, then their GW signals could be visible from most of the way across the observable universe. The relatively more complicated environments of SMBHs and the importance of processes of accretion and mergers on SMBH spin evolution Arvanitaki and Dubovsky (2011) requires further dedicated study for reliable event rate estimates.

A number of papers have considered whether exponential amplification though superradiance may occur via the SM photon. It was claimed in e.g. Pani et al. (2012b) that, if the SM photon had a very small mass, then the growth of gravitationally bound photon states could lead to significant spin-down of supermassive BHs. However, this analysis neglects the effects of the surrounding plasma on the photon equation of motion, which for the extremely low photon masses considered, $\lesssim 10^{-18}\text{~{}eV}$, will be very important (for example, the plasma frequency corresponding to a typical interstellar medium electron density of $0.01\,\text{cm}^{-3}$ is $\simeq 4\times 10^{-12}\text{~{}eV}$).

Another possibility is that, even for a zero-mass SM photon, the effects of an astrophysical plasma medium on the photon’s dispersion relation could result in superradiant amplification of bound states. The original ‘black hole bomb’ paper Press and Teukolsky (1972) noted that, if there was an ‘evacuated cavity’ in the plasma around a BH, this could act as a ‘mirror’ to confine superradiant waves. However, we are not aware of plausible astrophysical mechanisms to generate suitable ‘cavities’. More recently, papers including Pani and Loeb (2013) and Conlon and Herdeiro (2017) have claimed that for BH within astrophysical plasmas (the hot early universe plasma for Pani and Loeb (2013), and the diffuse interstellar medium for Conlon and Herdeiro (2017)), the ‘effective plasma mass’ could gravitationally confine EM oscillations around the BH, leading to exponential superradiant growth.

Both of these papers make the assumption that EM oscillations within a uniform plasma are governed by the Proca equation, $D_{\mu}F^{\mu\nu}=\omega_{p}^{2}A^{\nu}$. However, only the transverse photon modes have a dispersion relation approximately corresponding to a free particle with mass $\omega_{p}$ — the longitudinal mode in a non-relativistic plasma propagates at a much lower speed $\mathcal{O}(T/m_{e})$ Raffelt (1996). As reviewed in Appendix A, only the ‘axial’ ($j=\ell$) bound states for a massive vector correspond to a superposition of purely transverse waves, so we do not expect the bound state properties for polar states of a Proca field to reflect the behaviour of the SM photon in a plasma. Also, an increased plasma density close to the BH (as expected from e.g. accretion onto the BH) will generally result in a suppressed bound state density close to the BH, leading to suppressed growth rates. It should be kept in mind that there are other astrophysical situations in which energy and angular momentum could be extracted from Kerr BHs via electromagnetic effects, e.g. the Blandford-Znajek process Blandford and Znajek (1977). However, more careful calculations would be required to determine whether setups closer to those considered in this paper could be realised for the SM photon.

In this paper, we have only considered the gravitational interactions of light vectors. Since the vector masses we are interested in are so small, the couplings of the vector to SM matter are generally constrained to be significantly sub-gravitational in strength Adelberger et al. (2003); Wagner et al. (2012). The effects of such couplings should therefore be a small perturbation to the gravitational effects we have worked out. An exception is if long-distance fifth-force or equivalence principle tests do not apply. A clear example is a very light kinetically-mixed dark photon, for which the effective coupling to SM matter is suppressed by at least $(m_{A^{\prime}}/\omega_{p})$, where $m_{A^{\prime}}$ is the dark photon mass, and $\omega_{p}$ is the ‘effective plasma frequency’ of the SM medium. The charge density around Earth-based experiments is high enough that they impose weak constraints on the kinetic mixing at the very light vector masses we consider Essig et al. (2013). In contrast, the BH superradiance signals in this paper will apply for any kinetic mixing below some upper threshold. Calculating this threshold, and whether any additional observational signals arise near to it, would be very interesting, though we leave it to future work.

VI acknowledgements

We thank Joseph Bramante, Savas Dimopoulos, Sergei Dubovsky, Luis Lehner, and Riccardo Penco for many informative and clarifying discussions, Junwu Huang, Anthony Lasenby and Ken Van Tilburg for comments on the manuscript, Sam Dolan for correspondence regarding vector bound states around Schwarzschild BHs, Warren Morningstar for BH spin measurements discussions, and especially Asimina Arvanitaki for conversations throughout the completion of this work and comments on the manuscript. We are grateful to William East and Solomon Endlich for many extensive conversations about their respective spin-1 superradiance rate calculations. This research was supported in part by Perimeter Institute for Theoretical Physics. Research at Perimeter Institute is supported by the Government of Canada through the Department of Innovation, Science and Economic Development and by the Province of Ontario through the Ministry of Research, Innovation and Science. M.T. is supported in part by the Stanford Graduate Fellowship.

Appendix A Superradiance rates for non-relativistic bound states

As mentioned in Section II, there are no analytic approximations to scalar or vector bound states in a Kerr background that are valid at all radii. Finding bound state superradiance rates, which depend on the wavefunction near the horizon, therefore requires either a numerical solution of the wave equation, or some matching between approximate solutions with different regimes of validity.

Both approaches have been used to find scalar superradiance rates. The Klein-Gordon equation is separable in the Kerr background, so the problem can be reduced to solving the one-dimensional radial equation. This has been done numerically in a number of papers (e.g. Dolan (2007)). At small $\alpha$, Detweiler (1980) finds approximate solutions at large and small radii in terms of hypergeometric functions, and matches at intermediate radius to obtain the superradiance rates to leading order in $\alpha$ (as discussed below, there appears to be a factor of 2 error in Detweiler (1980)’s calculation).

We use a similar matching argument to find the growth/decay rates for vector bound states at small $\alpha$. Compared to the scalar case, the problem is complicated further by the fact that the Proca equation is not separable in the Kerr background. Instead, we find separable approximations which apply at large and small radii, which can be matched in an intermediate regime.

In the limit of small mass, the Proca field can be decomposed into two transverse modes, which obey the photon EoM, and a longitudinal mode, which obeys the massless KG equation. The evolution of massless scalars, and the behavior of some of the gauge-invariant functions of the photon field, can be solved exactly on the Kerr background via the Teukolsky equation Teukolsky (1973). This helps us because, for $r\ll\mu^{-1}$, the mass term in the Proca wave equation is sub-dominant, and solutions to the full wave equation are well approximated by solutions to the massless equations. Conversely, for $r\gg r_{g}$, the wavefunction is hydrogenic in the small-$\alpha$ limit, as reviewed in Section II. Consequently, if the hydrogenic wavefunction in the regime $r_{g}\ll r\ll\mu^{-1}$ can be approximated as a superposition of massless transverse and scalar modes, then we can use the behavior of these massless modes to find the energy flux across the BH horizon for the bound state.

The absorption probability for a long-wavelength massless (bosonic) spin-$s$ wave carrying total angular momentum $j$ is Page (1976)

[TABLE]

where $A_{H}$ is the BH horizon area, $\kappa\equiv(4\pi(r_{+}-r_{g}))/A_{H}$, and a negative absorption probability corresponds to superradiant amplification. A massless wave of frequency $\omega\simeq\mu$ has an absorption probability scaling as $\alpha^{2j+1}$, to leading order in $\alpha$. As we will see below, part (in some cases all) of the energy flux into the BH can be obtained by matching the hydrogenic bound state at small $r$ onto a massless wave of this frequency. The relevant scaling of the hydrogenic bound state’s squared amplitude at small $r$ is $\alpha^{2\ell}$, giving an overall growth rate of $\Gamma_{\mathrm{SR}}\propto\alpha^{2j+2\ell+4}\mu$ (where the 4 comes from choosing $\mu$ as our normalization; see below).

For a scalar field, where $j=\ell$, this reproduces the known $\alpha^{4\ell+4}\mu$ scaling. As reviewed in Section A.7, previous works observed the $\alpha^{2j+2\ell+4}\mu$ scaling in numerical solutions of the Proca equation in a BH background, and derived it analytically for the $j=\ell$ mode in the small-$a_{*}$ limit. Here, we demonstrate how the scaling arises for all of the bound states, and analytically determine the leading-$\alpha$ growth/decay rates for some of the low-$j$ modes, including the fastest-growing $\ell=0,j=m=1$ mode. In the main text, when we use higher-$\ell$ superradiance rates for phenomenological purposes, we adopt an estimate that should be $\mathcal{O}(1)$ correct at small $\alpha$ — the full leading-order expression could be derived in a similar fashion to our calculations here.

A.1 Hydrogenic wavefunctions

For $r\gg r_{g}$, the bound state wavefunctions are hydrogen-like,

[TABLE]

where $f_{0}(r)$ is fixed by the Lorentz condition. The vector spherical harmonics are given by

[TABLE]

where $\xi^{0}=\hat{z},\,\xi^{\pm 1}=\mp\frac{1}{\sqrt{2}}(\hat{x}\pm i\hat{y})$ are unit vectors Thorne (1980). The labels run over $\ell=0,1,2,\dots$, with $j=\ell-1,\ell,\ell+1$ ($j=0,1$ for $\ell=0$), and $m=-j,\dots,j$; $\ell$ can be identified as the orbital angular momentum, since $Y_{i}^{\ell,jm}(\theta,\phi)$ are eigenfunctions of the orbital angular momentum operator, $-r^{2}\nabla^{2}Y_{i}^{\ell,jm}=\ell(\ell+1)Y_{i}^{\ell,jm}$.

For $a\gg r\gg r_{g}$,

[TABLE]

where $R_{0}$ is the first coefficient in the Taylor expansion of $R^{n\ell}(r)$ around the origin.

A.2 ‘Massless’ wavefunctions

As noted above, for $r\ll\mu^{-1}$ the mass term is a sub-leading correction to the Proca wave equation. This is particularly easy to see in a Schwarzschild background, where the Proca equation separates into angular and radial parts. The angular equations do not depend on the mass, and the mass only enters into the radial equations via the operator Rosa and Dolan (2012)

[TABLE]

where $r_{*}$ is the standard tortoise coordinate. As a result, for $r\ll\mu^{-1}$, the $\mu^{2}$ term is a sub-leading correction to the EoM. This is also true in the full Kerr background.

For $r\gg r_{g}$, the solutions of the massless wave equation are close to the flat-space solutions. A (flat-space) transverse plane wave has $\partial_{i}A_{i}=0$, so by the Lorentz condition, $A_{0}=0$. A basis for the transverse modes with frequency $\omega$ that are regular at the origin is given by

[TABLE]

[TABLE]

where the $j_{n}$ are spherical Bessel functions of the first kind. Near the origin, for $r\ll\omega^{-1}$,

[TABLE]

where $!!$ denotes the double factorial.

The longitudinal solutions of the flat-space wave equation are, in the massless limit, pure-gauge; that is, $-\partial_{i}A_{0}-\partial_{t}A_{i}=0$ and $\nabla\times\vec{A}=0$. In terms of the pure-orbital vector spherical harmonics, the regular-at-the-origin longitudinal solutions are

[TABLE]

[TABLE]

Together with the transverse modes $\Psi^{E}$ and $\Psi^{B}$, these form a complete basis for regular-at-the-origin solutions to the flat-space wave equation for vanishing mass.

The reason for taking the solutions that are regular at the origin is that these correspond to an ingoing wave which ‘passes through the origin’ undisturbed to become an outgoing wave. At large radii, $r\gg\omega^{-1}$, we can decompose $j_{n}(\omega r)$ into ingoing and outgoing parts,

[TABLE]

A.3 $j=\ell$ bound states

The $j=\ell$ hydrogenic bound states have wavefunction

[TABLE]

In the interval $r_{g}\ll r\ll\mu^{-1}$ we can match this form onto the $\Psi^{B}$ transverse solution, with

[TABLE]

Thus, near the BH, the wavefunction looks like a superposition of ingoing and outgoing massless waves. In that physical situation, the energy flux through the horizon is given by the energy flux in the ingoing wave, multiplied by the BH absorption probability from equation 28. This gives us the energy flux into the BH from the ‘gauge-invariant’ components of the stress-energy tensor, i.e. the first and third terms in

[TABLE]

Since $A_{0}\simeq 0$ in the hydrogenic regime, the $\mu^{2}A_{\mu}A_{\nu}$ contribution to the energy flux into the BH can be shown to be sub-leading in $\alpha$ compared to the $F_{\mu\alpha}F_{\nu}^{\,\,\alpha}$ ‘Poynting’ contribution (the $g_{\mu\nu}$ parts of $T_{\mu\nu}$ do not give an ingoing energy flux).

The magnitude of the Poynting flux can be found via equation 28. To set the normalization of the ‘massless waves’, we can evaluate the energy flux in the ingoing part of $\Psi^{B}$ by calculating the Poynting vector at large radii. There, the ingoing part has wavefunction $\Psi_{i}\simeq C\frac{e^{i\omega r}}{2\omega r}Y^{j,jm}_{i}$, for which the time-independent part of the Poynting vector is

[TABLE]

where $\perp$ indicates components perpendicular to the radial direction $\hat{r}$. This gives an energy flux through a sphere of large radius of $\sim|C|^{2}/(4\omega)$. Thus, the energy flux across the BH horizon for the bound state is, to leading order in $\alpha$, $|C_{B}|^{2}\mathbb{P}_{s=1,j}/(4\mu)$. Performing the analogous calculations for a scalar bound state, we obtain a flux of $|C_{B}|^{2}\mathbb{P}_{s=0,j}/(4\mu)$, so

[TABLE]

This agrees with the expression derived by Pani et al. (2012a), and is accurate to leading order in $\alpha$ for any $a_{*}$. Our matching calculation gives a scalar bound state decay rate that is two times smaller than that found in Detweiler (1980), though it agrees with the expression given in Pani et al. (2012a), and also matches the numerical results of Dolan (2007).

A.4 $\ell=j-1$ bound states

Vector bound states with $j\neq\ell$ have $A_{0}\neq 0$ at the level of the hydrogenic wavefunctions, so the matching procedure is slightly more complicated than for the $j=\ell$ states. For the $\ell=j-1$ hydrogenic bound state, we have

[TABLE]

Since $\Psi_{i}^{H}\sim(r/a)^{j-1}Y^{j-1,jm}_{i}$ near the origin, the leading term in the brackets cancels, and we obtain

[TABLE]

Taking the example of the $\ell=0,j=1,n=0$, mode, we have

[TABLE]

In the matching region, where only the first term of (35) is significant at leading order in $\alpha$, the bound state can be approximated by the purely-transverse even-parity mode $\Psi^{E}$ (35). As for the $j=\ell$ modes, this matching tells us the energy flux into the BH from the $E$ and $B$ fields. However, for the $\ell=0,j=1$ mode, the contribution from $\mu^{2}A_{\mu}A_{\nu}$ is no longer sub-leading. In a Schwarzschild background, we can see this from the fact that the radial dependence of $\Psi_{0}$ is approximately constant near to the BH, taking the form from equation 47 even at small radii.111111Using the Proca equations of motion, which are separable in Schwarzschild coordinates Rosa and Dolan (2012), we have $\hat{\mathcal{D}}_{2}(rA_{0})-\frac{2r_{g}f}{r}\left(E_{r}+(rA_{0})/r^{2}\right)=0$, where $f\equiv 1-2r_{g}/r$, and $E_{r}$ is the radial electric field. Since $A_{0}$ is approximately constant at small radii in the hydrogenic regime, $A_{0}\simeq-r_{g}E_{r}$ there. In a Schwarzschild background, $E_{r}$ obeys its own wave equation, and it can be checked that $A_{0}+r_{g}E_{r}$ remains small even when $r$ becomes comparable to $2r_{g}$. Thus, $A_{0}$ approximately obeys the $\ell=0$ scalar radial equation, so is approximately constant all the way to the horizon. The energy flux through the BH horizon is

[TABLE]

where $k^{\mu}$ is the timelike Killing vector field, and $\xi^{\mu}$ ($=k^{\mu}$ in Schwarzschild) is the normal to the horizon. Using ingoing Eddington-Finkelstein coordinates, we obtain

[TABLE]

(where the angle brackets indicate time averaging) to leading order in $\alpha$, compared to $\frac{32}{3}\alpha^{7}\mu^{2}$ for the Poynting term from the $\Psi^{E}$ matching. Numerical simulations of Proca bound states in a Schwarzschild background show a ratio of Poynting to $\mu^{2}A_{\mu}A_{\nu}$ fluxes that is very close to 2, as expected. The combined decay rate for the bound state is $\Gamma\simeq 16\alpha^{7}\mu$  , (where we take ‘decay rate’ to mean the energy density decay rate, i.e. twice the field amplitude decay rate).

The $\Psi^{E}$ matching computation does not depend on whether the background is Schwarzschild or Kerr, since the $\sim(r_{g}/r)^{2}$ differences between the two metrics give corrections that are sub-leading in $\alpha$. So, to obtain the Poynting energy flux into a Kerr BH, we can multiply the Schwarzschild flux by the pre-factor $\frac{\mathbb{P}_{j=1}|_{a_{*}}}{\mathbb{P}_{j=1}|_{a*=0}}$. For the $\mu^{2}A_{\mu}A_{\nu}$ contribution, the energy flux into the BH scales like $\frac{\mathbb{P}_{j=0,m=1}|_{a_{*}}}{\mathbb{P}_{j=0,m=1}|_{a*=0}}$ (where we have abused notation slightly by evaluating (28) with $j=0$ but $m=1$). Combining these scalings, we can obtain the leading-$\alpha$ growth rate for the bound state around a Kerr BH of any spin,

[TABLE]

where the last equality is in the $\alpha\ll a_{*}$ limit. While the $\omega$ terms in the $\prod_{n=1}^{j}$ part of $\mathbb{P}$ (equation 28) are technically higher order in $\alpha$, we keep them in the expressions we use, since they give the correct behavior when $\omega-m\Omega$ is a small parameter (i.e. near the superradiance boundary for that level).

For $\ell=j+1>0$ bound states, the situation is similar: the Poynting flux can be obtained by matching onto the transverse $\Psi^{E}$ mode, while there is an additional ‘scalar’ flux that can be determined from the behavior of the $A_{0}$ part of the wavefunction. It can be demonstrated that both of these result in a $\Gamma\sim\alpha^{2j+5}\mu$ scaling (for the Schwarzschild decay rate), while their coefficients could be calculated by arguments similar to those above.

A.5 $j=0$, $\ell=1$ bound states

For $j=0$, the longitudinal mode is

[TABLE]

(there being no $Y^{-1,00}$ harmonic) which matches the form of the hydrogenic bound state,

[TABLE]

This is as expected, since massless transverse states must have $j\geq 1$, so cannot match the form of a $j=0$ state. Matching onto the massless scalar absorption rate, we obtain

[TABLE]

as the leading-$\alpha$ Schwarzschild absorption rate. We can derive the Kerr rates via multiplying by the $\frac{\mathbb{P}_{j=0,m=0}|_{a_{*}}}{\mathbb{P}_{j=0,m=0}|_{a*=0}}$ factor as above.

A.6 $\ell=j+1$ bound states

Similarly to the $\ell=j-1$ bound states, the $\ell=j+1$ bound state growth/decay rates receive same-order contributions from the Poynting flux obtained by matching onto the $\Psi^{E}$ mode, and from the ‘scalar’ flux. The $\Psi_{i}^{H}\sim(r/a)^{j+1}$ behavior near the origin results in a Schwarzschild decay rate $\Gamma\sim\alpha^{4j+7}\mu$, which can then by translated in a Kerr growth/decay rate as above.

A.7 Comparison to literature

A number of papers have considered the growth/decay rates of vector bound states around BHs. The first work to notice the $\Gamma\sim\alpha^{2j+2\ell+5}\mu$ scaling, Rosa and Dolan (2012), performed numerical computations of bound state decay rates around a Schwarzschild BH, where the Proca equation is separable (as noted in Section II, their labels for $j$ and $\ell$ are switched relative to ours). They observed that the numerical results were well fit by the form we derived above, and gave an analytic argument for the $j=\ell$ case. In Figure 11, we compare our analytic leading-$\alpha$ decay rates for the Schwarzschild background, to Rosa and Dolan (2012)’s numerical results. The numerical results converge to the analytic expressions at small $\alpha$, and deviate at finite $\alpha$ in a manner similar to the scalar case, which is plotted in the lower panel of Figure 11.

A different analytic matching approach was pursued in Endlich and Penco (2016). While we match at the level of wavefunctions, Endlich and Penco (2016) posit an effective Hamiltonian that describes the interaction of long-wavelength field modes with a small BH, and then use this to compute bound state growth/decay rates in the small-$a_{*}$ approximation. Given the assumed form of their effective Hamiltonian, there appears to be no choice for the operator coefficients that leads to the correct scaling of growth rates / amplification factors as a function of $\alpha$ for all of the bound states / scattering states. For example, the choice of coefficients presented in Endlich and Penco (2016) gives the correct scaling $\alpha^{2\ell+2j+5}\mu$ for $j\geq 1$ but gives a decay rate of $\sim\alpha^{9}\mu$ in Schwarzschild for the $j=0,\ell=1$ mode, whereas we find $\sim\alpha^{7}\mu$, which matches the numerical computations of Rosa and Dolan (2012) and Pani et al. (2012a). Moreover, though this choice gives the correct $\alpha$ dependence for other bound states, it gives different numerical coefficients to ours. Comparing the leading-$a_{*}$, leading-$\alpha$ superradiance rates presented in Endlich and Penco (2016) with those derived in this work, we have

[TABLE]

where $\Gamma_{\rm EP}$ denotes the rate from Endlich and Penco (2016). In particular, Endlich and Penco (2016) obtains a rate for the $j=\ell=1$ level a factor of 2 larger than this work, where the latter agrees with the analytic calculations of Rosa and Dolan (2012) and Pani et al. (2012a).

Another approach is to expand in small $a_{*}$, obtaining a coupled set of one-dimensional equations are then solved numerically Pani et al. (2012a, b) . Whereas our analytic approach is under control at small $\alpha$ but arbitrary $a_{*}$, this semi-analytic method Pani et al. (2012a, b) is theoretically under control at arbitrary $\alpha$ but small $a_{*}$. The numerics are challenging at small $\alpha$, resulting in the $\alpha^{2\ell+2j+4}$ scaling being obscured due to numerical errors for the leading $j=1,\ell=0$ level Pani et al. (2012a). Comparing the two approaches, we find good agreement in the small-$\alpha$, small-$a_{*}$ regime, as expected. Figure 12 compares the results for the $\ell=0$, $j=m=1$ level. At $a_{*}=0$, the small difference between the result from Pani et al. (2012a), which is a numerical calculation to all orders in $\alpha$, and our leading-$\alpha$ result, comes from $\alpha$ needing to be smaller for our expansion to accurately converge, as per Figure 11. As Figure 12 illustrates, the difference is significantly greater at larger $a_{*}$. It would be interesting to compare to a full numerical calculation, to see how small $\alpha$ needs to be before our leading-$\alpha$ form gives the correct $a_{*}$ dependence.

To perform a controlled calculation at large $a_{*}$ and large $\alpha$, it appears to be necessary to do a full 2D (in $r$ and $\theta$) computation in the Kerr background. Witek et al. (2013) performs a numerical time-domain simulation of this kind, and finds that the amplitude of an initial wave packet grows over time around a fast-spinning BH, though they do not appear to cleanly separate the different bound states. East (2017) performs a similar numerical simulation, and extracts the fastest-growing vector bound states in a Kerr background. As described in Section II.3 and illustrated in Figure 2, their results correspond well to our analytic approximation.

Appendix B Perturbations due to companion

The presence of a star or BH as a binary companion to a superradiating BH can hinder superradiance of a growing level because such an asymmetric perturbation can cause the growing level to mix with decaying ones. In this appendix, we show that this perturbation does not affect the bounds set with BH spin measurements, and estimate the size of the effect on spins in BBH systems.

The orbits of companion stars in X-ray binaries have a period on the order of days or more, at distances of $R\sim 10^{7}$ km away from the BH. Since the orbital periods are much smaller than the energies of the bound states, the perturbing potential changes adiabatically, and we can consider a static perturbing potential $\delta V(\vec{r})$. This is also the case for BBH systems at the early stages of their formation, where they are separated by distance

[TABLE]

assuming two equal mass BHs that merge in time $T$ through GW emission.

The condition that a particular superradiant level has a net positive growth is given by

[TABLE]

where $\Gamma_{\text{decay}}$ and $\Gamma_{\text{sr}}$ are the decay and superradiance rates, $A_{\text{decay}}$ and $A_{\text{sr}}$ are the bound states, and $\Delta E$ is the energy difference between the states Arvanitaki et al. (2015).

At small $\alpha$ (large cloud size) where this perturbation is most important, the ratio of the rates is dominated by the $\alpha$ power, and scales as

[TABLE]

To estimate the perturbation by the companion on the bound state, we decompose its gravitational potential into multipoles:

[TABLE]

where $M^{\prime}$ is the companion mass, $R$ the distance between BH and companion, and primed coordinates correspond to the companion’s location. The leading contribution comes from the dipole term. The dipole mixing term allows mixing of the superradiating ($N=\ell+1$, $j$, $\ell=j-1$, $m=j$) state with any other bound state with $\Delta\ell=1$.

The fastest superradiating level has $j=m=1,\ell=0,N=1$ and can undergo dipole mixing with the $N=2$, $j=m=0$, $\ell=1$ level. Since $\ell+j$ is the same for both states, $\Gamma_{\text{decay}}/\Gamma_{\text{sr}}\sim 1$. Because $\Delta N\neq 0$, the energy difference is set by the difference in the bound state energies $\Delta E\sim\frac{3\mu\alpha^{2}}{8}$. For perturbations to be unimportant for this level, the condition is

[TABLE]

For typical X-ray binary and BBH systems, this condition saturates at $\alpha\sim 10^{-3}$. At such low values of $\alpha$, the superradiance rates are too slow to form a maximally filled cloud within a binary lifetime, so mixing does not impose an additional constraint.

The story for the $j>1$ levels is slightly different. For the strongest perturbation, consider mixing with the “lower” ($N^{\prime}=N$, $j^{\prime}=j-1$, $\ell^{\prime}=\ell-1$, $m^{\prime}=m-1$) bound state. Because the energy quantum number $N$ stays the same but $\ell$ changes, the energy difference is given by the “fine” splitting $\sim\mu\alpha^{4}$ in the hydrogenic limit; including higher-order corrections to the bound states and metric may lift degeneracies and increase the splitting, which would reduce the effect of the mixing. The ratio of decay to SR rate is $\sim\alpha^{4}$.

For $j>1$ levels, the resulting mixing does not affect superradiant growth if

[TABLE]

For the X-ray binary systems used to constrain vector masses, the values of $\alpha$ that saturate this condition are small enough such that the “lower” bound state still satisfies the SR condition. (Note that these “lower” bound states we considered are not the same as the fastest superradiating bound states with a lower $j$; they have different energy quantum number $N$.) Thus, for these BHs, perturbations due to a companion star are not the limiting factor. Instead, the limiting factor is that the superradiance time has to be faster than the BH lifetime.

For BBH systems, the typical time for superradiance to occur is longer, $10^{7}$ – $10^{10}$ years, and the BH companion may start at a closer distance. This means that for the $j>1$ levels the mixing is important in limiting the spin-down of BHs in black hole binaries, and we take this into account when estimating the power of statistical searches at aLIGO (Section III.2).

For $j=1$, the values of $\alpha$ affected by companions are also typically too small to produce observable GW signals, so blind searches for GWs may be sensitive to BHs in binary systems as well as isolated BHs; the populations are expected to be comparable Shakura and Sunyaev (1973).

Appendix C Gravitational radiation rates

From Misner et al. (1973), if the metric $g_{\mu\nu}^{(B)}$ corresponds to a background vacuum ($R_{\mu\nu}=0$) spacetime, then small metric perturbations

[TABLE]

are sourced by a small stress-energy tensor $T_{\mu\nu}$ through the linearized wave equation

[TABLE]

where $\bar{h}_{\mu\nu}=h_{\mu\nu}-\frac{1}{2}hg^{B}_{\mu\nu}$ is the trace-reversed metric perturbation (working in ‘Lorentz gauge’ $\bar{h}_{\mu\nu}^{;\nu}=0$). Indices here are raised and lowered by $g^{(B)}_{\mu\nu}$, and $\nabla$s denote covariant derivatives.

In flat space, this reduces to the usual wave equation

[TABLE]

This can be solved using the flat-space Green’s function. For a periodic $T_{\mu\nu}$, we obtain an emitted power Weinberg (1972)

[TABLE]

where $\tilde{T}_{\mu\nu}$ is the spatial Fourier transform of $T_{\mu\nu}$, and $\tilde{T}_{TT}^{ij}$ is the transverse (relative to $\vec{k}$) traceless part of $\tilde{T}^{ij}$.

If the background space is weakly curved,

[TABLE]

where $\epsilon_{\mu\nu}$ is small, then we can expand $h_{\mu\nu}$ to each order in $\epsilon$, writing $h_{\mu\nu}=h^{(0)}_{\mu\nu}+h^{(1)}_{\mu\nu}+\dots$. The homogeneous wave equation at first order in $\epsilon$ is then (schematically)

[TABLE]

where we have suppressed indices.

For gravitational radiation from transitions between levels, the above story works well. It can be checked that, in the small-$\alpha$ limit, the contribution to gravitational radiation from the strongly-curved space near the BH ($r\sim{\rm few\,}\times r_{g}\ll\omega^{-1}$) is small. In the weakly-curved background far away from the BH, the $h^{(1)}$ correction is sub-leading in $\alpha$ to $h^{(0)}$, where $h^{(0)}$ is the result obtained from the non-relativistic bound state $T_{\mu\nu}$ sourcing gravitational waves through their flat-space Green’s function.

However, the small-curvature expansion breaks down for gravitational waves from bound state annihilations. The problem is that we have only a single small parameter, $\alpha$, which controls both the small curvature and the properties of $T_{\mu\nu}$. From the continuity properties of the (non-relativistic) $T_{\mu\nu}$ near the origin, its Fourier transform must fall off at least as fast as some particular power of $k$, at large $k\gg a^{-1}$. However, due to a cancellation that occurs for hydrogenic bound states, it actually falls off faster than this, resulting in $h^{(0)}$ being suppressed by extra powers of $1/(ak)\sim\alpha$. If we try to solve Eq. (70) for $h^{(1)}$, there is a suppression from the small non-flat parts of the metric, but the cancellation from the radial profile of the source term is lifted. The result is that $h^{(1)}$ is of the same order or lower order in $\alpha$ as $h^{(0)}$, and the expansion $h_{\mu\nu}=h^{(0)}_{\mu\nu}+h^{(1)}_{\mu\nu}+\dots$ stops making sense.

In the case of the $2p$ scalar bound state Brito et al. (2015a), the cancellation suppresses the flat space amplitude by two powers of $\alpha$ relative to the ‘naive’ expectation, while (in the formal expansion) the $h^{(1)}$ part is only suppressed by $\alpha$ from the small curvature. The result is that the ‘Schwarzschild background corrections’ actually determine the leading-order-in-$\alpha$ annihilation rate. For higher-$\ell$ scalar bound states, the cancellation only suppresses the flat space amplitude by one power of $\alpha$, resulting in the Schwarzschild correction being of the same order.

For all of the vector bound states, the latter situation applies; the Schwarzschild correction gives a same-order contribution to the annihilation rate. To actually compute the leading-$\alpha$ form of the rate, it is generally simpler to use the wave equation for perturbations to the Riemann tensor, rather than to the metric as above. Using the Newman-Penrose scalar $\Psi_{4}$, which corresponds to the appropriate components of the Weyl tensor describing gravitational radiation, it is possible to write down a one-dimensional radial equation, sourced by appropriate combinations of $T_{\mu\nu}$ and its derivatives Brito et al. (2015a); Poisson (1993). We have calculated an approximation to this solution for the $\ell=0,j=1,m=1$ level (the fastest-growing vector bound state), finding that it leads to a rate $\sim 10$ times the naive flat-space estimate in the small-$\alpha$ limit. This gives an emitted GW power of

[TABLE]

where $N_{M}\equiv GM^{2}$ is of order the maximum occupation number of the level (as per Eq. (14)). Thus, for small $\alpha$, the cloud annihilates on a timescale much larger than its oscillation time. However, as discussed in Section II.4, this timescale can be shorter than other astrophysically-relevant ones, including the growth timescale for the next superradiant level. For a general bound state, the gravitational radiation power scales like $\alpha^{4\ell+12}GN^{2}r_{g}^{-4}$. Table 3 shows the ‘flat-space’ results for the first few levels; as described above, one would expect the full leading-$\alpha$ results to have the same $\alpha$ dependence, but different (generally larger) numerical coefficients.

For a spherically-symmetric background metric, the angular distribution of the emitted gravitational radiation is set by the angular form of the $T_{ij}$ source. This will be the case for the leading-$\alpha$ GW emission from bound states, since as per above, the leading-$\alpha$ annihilation power can be obtained purely from the Schwarzschild part of the metric. Thus, the emitted tensor modes from a $Y^{\ell,jm}$ state will be of the $T^{E2,(2j)(2m)}$ form (in the notation of Thorne (1980)), giving the same polarization and angular emission profiles as a standard quadrupole GW source, such as a spinning asymmetric neutron star. For a $j=m=1$ state, this corresponds to an angular emission power profile of

[TABLE]

In particular, the emitted power in the most intense direction is $5/2$ times the angle-averaged power, and in the least intense direction is $5/16$ times the average power. Since this angular variation is not too strong, we have simply used the averaged value for the calculations in the main text.

Appendix D Calculation for number of detectable sources

For each vector mass $\mu$, we estimate the number of detectable sources by integrating over various astrophysical distributions of stellar mass BHs:

[TABLE]

where $P(M)$, $P(a_{*})$, and $P(r)$ are the normalized probability distributions of finding a BH with mass $M$ and spin $a_{*}$ at a distance $r$ away, reach($M,a_{*}$) is the maximum distance of such a BH whose signal is above the detector noise floor, $\tau_{\text{sig}}(M,a_{*},r)$ is the signal duration, and BHFR is the BH formation rate.

We follow the astrophysical distributions used in Appendix C of Arvanitaki et al. (2015) with a few changes to the BH mass distributions in light of recent LIGO observations confirming the existence of $\sim 30M_{\odot}$ BHs. To account for astrophysical uncertainties, the parameters in the distributions are varied to produce optimistic, realistic, and pessimistic estimates (see Table 4).

For the mass distribution, we use

[TABLE]

with $M_{\text{min}}=4.1M_{\odot}$, and width $M_{0}$ that is varied in our estimates. We impose a maximum cutoff mass with two different values: $M_{\text{max}}=(80,160)M_{\odot}$. Using a distribution that falls off more slowly, e.g. power law, increases the number of detectable sources at low vector mass. We use the spin distributions stated in Table 4. For the distance distribution, we use a $\delta$-function at 8 kpc for BHs in the Milky Way, and follow the blue luminosity distribution (Figure 6 in Kopparapu et al. (2008)) for those beyond our galaxy.

Appendix E Frequency drift

As mentioned in Section IV.3, the gravitational self-energy of the cloud will affect the frequency of GWs emitted through annihilations. As the cloud annihilates away, this contribution to the binding energy decreases, resulting in the emitted frequency increasing with time. The rate at which the frequency changes can affect the detectability of the annihilation signal.

LIGO continuous wave searches currently cover a range of positive to negative frequency derivatives, $3.39\times 10^{-10}$ Hz/s through $-2.67\times 10^{-9}$ Hz/s with $T_{\text{coh}}=60$ hours Abbott et al. (2016e). For $T_{\text{coh}}=250$ hours, and the same amount of computational time, the maximum positive frequency drift searched over is $\dot{f}_{\text{max}}\sim 2\times 10^{-11}$ Hz/s Wette (2009).

To estimate the expected frequency drift $\dot{f}$, we start by calculating the gravitational self-energy (per particle) of a bound state,

[TABLE]

where

[TABLE]

For the $j=1$ bound state, $\psi$ is the $n=0$, $\ell=0$ hydrogen wave function, and the frequency drift is

[TABLE]

where we have used $\dot{N}_{m}=-\Gamma_{\rm ann}N_{m}^{2}$. For annihilation signals that come from beyond the Milky Way, which are stronger and have shorter annihilation times than those from the Milky Way, the frequency drifts are too large ($\gtrsim 10^{-8}$ Hz/s), and are thus beyond the range searched by aLIGO. Without this constraint on the frequency drift, there would be a comparable number of expected events from outside the Milky Way peaking at $\sim 10^{-12}$ eV.