Impact of multiple modes on the black-hole superradiant instability

TL;DR

This paper investigates how multiple modes influence the superradiant instability of black holes caused by ultralight bosonic fields, revealing complex dynamics and conditions that affect the instability's development and observational constraints.

Contribution

It introduces a multimode analysis of black-hole superradiance, showing the impact of mode interactions and seed energy on the instability's evolution and end-state.

Findings

Multimode interactions can suppress or enhance superradiant instability.

Seed energy levels significantly influence the instability's growth.

Nonsuperradiant modes have negligible effects at very low seed energies.

Abstract

Ultralight bosonic fields in the mass range $\sim (10^{-20}-10^{-11})\,{\rm eV}$ can trigger a superradiant instability that extracts energy and angular momentum from an astrophysical black hole with mass $M\sim(5,10^{10})M_\odot$, forming a nonspherical, rotating condensate around it. So far, most studies of the evolution and end-state of the instability have been limited to initial data containing only the fastest growing superradiant mode. By studying the evolution of multimode data in a quasi-adiabatic approximation, we show that the dynamics is much richer and depend strongly on the energy of the seed, on the relative amplitude between modes, and on the gravitational coupling. If the seed energy is a few percent of the black-hole mass, a black hole surrounded by a mixture of superradiant and nonsuperradiant modes with comparable amplitudes might not undergo a superradiant unstable phase, depending on the value of the boson mass. If the seed energy is smaller, as in the case of an instability triggered by quantum fluctuations, the effect of nonsuperradiant modes is negligible. We discuss the implications of these findings for current constraints on ultralight fields with electromagnetic and gravitational-wave observations.

Tables

Table 1. Table 1: Initial mode configuration considered in this work.

Model	$(l,m)$		$\Psi$
	Mode-$1$	Mode-$2$
I	$(1,1)$	$(1,-1)$	Eq. (11)
II	$(1,1)$	$(2,2)$	Eq. (12)
III	$(1,1)$	$(2,1)$	Eq. (13)

Table 2. Table 2: Configurations considered in the adiabatic evolution for Models I-III.

Case	$m_B/\mathrm{eV}$	$M_0/M_{\odot }$	$\mu M_0$	$J_0/M^2$	$M_{S0}/M_0$
A	${10}^{-18}$	${10}^7$	$0.075$	$0.8,0.95$	${10}^{-9}$
B	${10}^{-18}$	${10}^7$	$0.075$	$0.8,0.95$	$0.025$
C	$4\times {10}^{-18}$	${10}^7$	$0.299$	$0.95$	$0.025$

Equations

$$S=\int d^{4}x\sqrt{-g}\left(\frac{R}{16\pi }-\frac{1}{2}{\partial }_{\mu }\Psi {\partial }^{\mu }\Psi -\frac{{\mu }^2}{2}{\Psi }^2\right),$$

$$\Psi (t,r,\vartheta ,\phi )=\Re \left[\int d\omega e^{-i\omega t+im\phi }{}_0{S}_{lm\omega }(\vartheta ){\psi }_{lm\omega }(r)\right],$$

$${\mathcal{D}}_{\vartheta }[{}_{0}S]+$$

$${\mathcal{D}}_r[\psi ]+$$

$$-\Delta ({\mu }^2{r}^2+a^2{\omega }^2+\lambda )]\psi =0,$$

$$\psi (r)\propto \frac{r^{\nu }{e}^{-\sqrt{{\mu }^2-{\omega }^2}r}}{r}\mathrm{as}r\to \infty ,$$

$$\omega \sim$$

$${\gamma }_{lm}=$$

$$\chi >$$

$$\psi \left(r;\mu ,\chi ,M\right)\propto g_l(r),$$

$$g_l(r)={\left(\frac{2rM{\mu }^2}{l+1}\right)}^l\exp \left(-\frac{rM{\mu }^2}{l+1}\right)L_0^{2l+1}\left(\frac{2rM{\mu }^2}{l+1}\right).$$

$$r_{\mathrm{cloud}}\sim \frac{l\left(l+1\right)}{{\left(M\mu \right)}^2}M,$$

$$\Psi =\sum _{lm}{A}_{lm}{g}_l(r)\cos \left(m\varphi -{\omega }_{R}t\right)P_{lm}(\cos \theta ),$$

$$\Psi$$

$$M_S=\int -T_0^0{r}^2\sin \theta \mathrm{d}r\mathrm{d}\theta \mathrm{d}\varphi ,$$

$$A_{11}^2$$

$${\lambda }_1=\frac{A_{1-1}}{A_{11}},{\lambda }_2=\frac{A_{22}}{A_{11}},{\lambda }_3=\frac{A_{21}}{A_{11}}.$$

$${\dot{E}}_{GW}$$

$${\dot{J}}_{GW}$$

$${\dot{E}}_{GW}$$

$${\dot{J}}_{GW}$$

$${\dot{E}}_{GW}$$

$${\dot{J}}_{GW}$$

$${\dot{E}}_S$$

$${\dot{J}}_S$$

$$⟨{\dot{E}}_S⟩$$

$$⟨{\dot{J}}_S⟩$$

$$⟨{\dot{E}}_S⟩$$

$$⟨{\dot{J}}_S⟩$$

$$⟨{\dot{E}}_S⟩$$

$$⟨{\dot{J}}_S⟩$$

$${\lambda }_1(t)$$

$${\lambda }_2(t)$$

Full text

I Introduction

In classical general relativity, where gravity is minimally coupled to massive bosonic fields, Kerr black holes (BHs) can be unstable against the superradiant instability (for an overview, see Superradiance ). This process was discovered almost 50 years ago zeldo1 ; Damour:1976kh ; Teukolsky:1974yv ; Press:1972zz but only recently it has been subject to intense scrutiny, including rigorous mathematical proofs Shlapentokh-Rothman:2013ysa ; Moschidis:2016wew . It was realized that this instability effectively turns astrophysical BHs into detectors of axion-like particles Arvanitaki:2009fg ; Arvanitaki:2010sy and of ultralight, beyond-standard model bosons in general.

For a BH of mass $M$ and an ultralight boson with mass $m_{B}\equiv\hbar\mu$, the instability is efficient only when the gravitational coupling $\frac{GM\mu}{c}\sim{\cal O}(1)$, i.e. when the Compton wavelength of the particle is comparable to the BH radius. Since astrophysical BHs are expected to exist at least in the mass range $\sim(5,10^{10})M_{\odot}$, the superradiant instability is effective for bosons approximately in the mass range $m_{B}\in(10^{-20}-10^{-11})\,{\rm eV}$, i.e. for ultralight bosons. The latter are compelling dark-matter candidates111The collapse and collision of compact objects composed of these dark matter candidates has been studied in Refs. Helfer:2016ljl ; Helfer:2018vtq . and are predicted in a multitude of beyond-standard model scenarios Jaeckel:2010ni ; Essig:2013lka ; Hui:2016ltb ; Irastorza:2018dyq .

The superradiant instability of a Kerr BH has been investigated perturbatively for scalar fields Detweiler:1980gk ; Cardoso:2006wa ; Dolan:2007mj ; Dolan:2012yt , including recent exploration of its phenomenological implications Hannuksela:2018izj ; Isi:2018pzk ; DAntonio:2018sff ; Baumann:2018vus ; Boskovic:2018lkj ; Ikeda:2018nhb , more recently for vector and tensor fields in a small-rotation expansion Pani:2012vp ; Pani:2012bp ; Brito:2013yxa ; Pani:2013hpa ; Endlich:2016jgc ; Cardoso:2017kgn , and for vector fields around BHs with arbitrary spin in an analytical Newtonian approximation valid for small gravitational coupling Baryakhtar:2017ngi , and numerically for generic values of the BH spin and gravitational coupling Witek:2012tr ; Cardoso:2018tly , also using a novel perturbation scheme Frolov:2018ezx ; Dolan:2018dqv . Recently, nonlinear simulations222The superradiant instability affects also Kerr BHs in asymptotically anti de Sitter spacetime; see Refs. Bosch:2016vcp ; Chesler:2018txn ; Sanchis-Gual:2015lje for nonlinear simulations in this context. of the Einstein equations minimally coupled to complex single mode vector fields East:2017ovw ; East:2017mrj ; East:2018glu have confirmed the analysis of previous quasi-adiabatic and perturbative evolution Brito:2014wla (see also Ref. Herdeiro:2017phl ). The latter is justified by the long instability time scale as compared to the dynamical time scale of the BH.

Note that in the case of complex massive bosonic fields or, equivalently, multiple real fields such that the resulting energy-momentum tensor respects the symmetry of the Kerr spacetime, there exist stationary spinning BH solutions surrounded by an oscillating condensate Herdeiro:2014goa ; Herdeiro:2016tmi . These solutions interpolate between boson stars and Kerr BHs and are formed during the evolution of the superradiant instability of Kerr BHs against complex bosons East:2017ovw ; Herdeiro:2017phl . These solutions are unstable against higher-order azimuthal modes Ganchev:2017uuo and, at least in some region of their parameter space Degollado:2018ypf , the instability time scale is comparable to that of Kerr. We deal here with a single real bosonic field, so the only stationary BH configuration is the Kerr metric, as guaranteed by the no-hair theorems Chrusciel:2012jk ; Herdeiro:2015waa .

The general properties of this process do not depend strongly on the nature of the bosonic field: the fundamental unstable mode has a frequency $\omega_{R}\sim\mu$ (henceforth we use $G=c=1$ units) and must satisfy the superradiant condition, $\omega_{R}<m\Omega_{H}$, where $m$ is the azimuthal number of the perturbation and $\Omega_{H}$ is the BH angular velocity. As a result of the instability, a single mode with $m>0$ and arbitrarily small amplitude grows exponentially near the BH, extracting energy and angular momentum on a time scale $\tau\equiv 1/\omega_{I}\gg M$, and forming a nonspherical, rotating condensate of characteristic size $r_{\rm cloud}\gg M$. Thus, if ultralight bosonic fields exist in nature, they would produce two generic signatures Arvanitaki:2009fg ; Arvanitaki:2010sy ; Superradiance : i) they would favor slowly-spinning BHs against highly-spinning ones, since BHs would lose their angular momentum over a time scale $\tau$ which can be much shorter than the typical BH accretion rate; and ii) they would produce a continuous gravitational-wave (GW) signal at a frequency set by the boson mass. The first signature translates into the existence of “gaps” in the BH “Regge plane”, i.e., in its spin–mass plane Arvanitaki:2010sy ; Pani:2012vp ; Pani:2012bp ; Brito:2013yxa ; Brito:2014wla , whereas the second signature can be directly searched for in LIGO/Virgo (and in the future LISA) data, both as isolated resolvable sources Arvanitaki:2014wva ; Baryakhtar:2017ngi ; Brito:2017zvb ; DAntonio:2018sff ; Isi:2018pzk ; Ghosh:2018gaw or through the GW stochastic background of a population of BH-boson condensates Brito:2017wnc .

Most phenomenological studies so far have focused on the idealized case in which the BH is initially surrounded by a single-mode superradiant seed (see Ref. Arvanitaki:2010sy , where the generic setup for the evolution of multiple modes has been laid down, although without discussing the phenomenology). However, more realistic configurations are likely to contain a superposition of modes, both superradiant (i.e., satisfying the $\omega_{R}<m\Omega_{H}$ condition) and nonsuperradiant. This is particularly important if the initial seed is due to quantum fluctuations, since in that case modes with different values of $m$ are expected to be produced with comparable amplitude.

Full-fledged $3+1$ numerical simulations including the backreaction of massive scalar Okawa:2014nda or vector fields Zilhao:2015tya ; Witek2018inprep onto the spacetime employed such multimode initial data either through explicit superposition or mode mixing due to the construction of metric initial data. These simulations, furthermore, assumed the presence of an appreciable bosonic cloud, i.e. a condensate of a few percent of the BH mass. Those are formed naturally via the superradiant evolution with small seeds Brito:2014wla . External effects such as a binary companion or the inspiral and merger of two such BH-condensate systems will cause mode mixing. The merger remnant would form in an environment containing a single cloud with complex multipolar structure; see e.g. Baumann:2018vus for work in this direction.

In those cases, the BH was shifted out of the superradiant regime by absorbing a counter-rotating mode with sufficiently large amplitude. This essentially switches off the superradiant instability, leaving a rotating BH surrounded by a slowly decaying bosonic condensate. These results indicate that the presence of multiple modes might crucially change the dynamics of the system. However, it is unclear whether this conclusion would persist for arbitrarily small initial seeds. Ideally, one wishes to follow the nonlinear evolution of a small initial seed at least for a few instability $e$-folding times, $\tau\sim 10^{6}M$ (resp. $\tau\sim 10^{4}M$), in the most favorable cases for scalars (resp. vectors). These type of simulations are numerically expensive and, hence, only a small number of cases with time scales of $\sim\mathcal{O}\left(10^{3}\right)M$ were analyzed. Instead, a quasi-adiabatic treatment along the lines of Ref. Brito:2014wla can provide crucial new insight into the evolution of (multimode) massive bosonic clouds surrounding BHs.

That is precisely the goal of this paper: study the impact of multiple modes on the evolution of the superradiant instability. As we shall show, the impact of an initial mixture of nonsuperradiant and superradiant modes with comparable amplitude depends strongly on the energy of the initial seed and on the value of the gravitational coupling $M\mu$. If this energy is initially much smaller than the BH mass (as expected in the most natural scenarios) the effect of multiple modes is negligible. On the other hand, if the energy is at least a few percent of the BH mass and $M\mu\sim{\cal O}(0.1)$ the presence of nonsuperradiant modes might affect the evolution and quench the instability completely. As we shall discuss, this latter scenario is more speculative and may comprise only a small fraction of the BH-scalar condensates expected in the universe.

The rest of this paper is organized as follows. In Sec. II we introduce our setup and different multimode models, and calculate their energy and momentum fluxes. In Sec. III we present the quasi-adiabatic evolution of our systems. We discuss their implications for current electromagnetic and GW-based bounds on the mass of axion-like particles in Sec. IV. We conclude in Sec. V.

II Setup

We focus on the action describing a real scalar field $\Psi$ with mass $m_{B}=\mu\hbar$ minimally coupled to gravity,

[TABLE]

where $g$ is the determinant of the spacetime metric $g_{\mu\nu}$, and $R$ is the Ricci curvature scalar. Minimization of this action yields the Klein-Gordon equation $\nabla_{\mu}\nabla^{\mu}\Psi=\mu^{2}\Psi$ and Einstein’s equations coupled to the stress energy tensor $T_{\mu\nu}=\partial_{\mu}\Psi\partial_{\nu}\Psi-\frac{1}{2}g_{\mu\nu}\left(\partial_{\alpha}\Psi\partial^{\alpha}\Psi+\mu^{2}\Psi^{2}\right)$.

Our setup will be the same as that of Ref. Brito:2014wla . In particular, we study the quasi-adiabatic evolution of the instability, i.e., neglecting the backreaction of the scalar field and instead employing energy and angular momentum balance argument. Although the total mass of the condensate can reach a few percent of the black hole mass, the stress-energy tensor (e.g. the energy density) remains small and our approximation valid as shown in Brito:2014wla . Furthermore, the energy and angular momentum extraction occurs over the instability time scale which is much longer than the BH dynamical time scale; this justifies a quasi-adiabatic evolution Brito:2014wla ; Herdeiro:2017phl . In this regime the dynamics is governed by the scalar-field equation on a fixed Kerr geometry, the mass and spin of which evolve adiabatically through energy and angular momentum fluxes.

The linearized dynamics of a Klein-Gordon field on the Kerr background with mass $M$ and spin $J=aM=\chi M^{2}$ is described by the Teukolsky equation for a spin-[math] perturbation, whose general solution can be written as

[TABLE]

where a sum over harmonic indices $({l},\,m)$ is implicit, and ${}_{s}Y_{{l}m\omega}(\vartheta,\varphi)={}_{s}S_{{l}m\omega}(\vartheta)e^{im\varphi}$ are the spin-weighted spheroidal harmonics of spin weight $s$ which, for $s=0$, reduce to the scalar spheroidal harmonics Berti:2005gp . The radial and angular functions satisfy the following coupled system of differential equations

[TABLE]

where for simplicity we omit the $({l},\,m)$ subscripts, $r_{\pm}=M\pm\sqrt{M^{2}-a^{2}}$ denotes the coordinate location of the inner and outer horizons, $\Delta=(r-r_{+})(r-r_{-})$, ${\cal D}_{r}=\Delta\partial_{r}\left(\Delta\partial_{r}\right)$, and ${\cal D}_{\vartheta}=(\sin\vartheta)^{-1}\partial_{\vartheta}\left(\sin\vartheta\partial_{\vartheta}\right)$.

II.1 Unstable modes

Imposing appropriate boundary conditions, namely purely ingoing waves at the horizon and exponential decay of the scalar field at infinity, a quasi-bound solution to the above coupled system can be obtained numerically, e.g. using continued fractions Cardoso:2005vk ; Dolan:2007mj or a shooting method Pani:2013pma . The eigenspectrum contains an infinite, discrete set of complex quasi-bound modes Berti:2009kk , $\omega=\omega_{R}+i\omega_{I}$. We will consider only fundamental modes with overtone number $n=0$, i.e. eigenfunctions with zero nodes. In particular, this system admits unstable ($\omega_{I}>0$) quasi-bound states satisfying the superradiant condition $\omega_{R}<m\Omega_{H}$ Detweiler:1980uk ; Dolan:2007mj , with $\Omega_{H}=a/(2Mr_{+})$ being the angular velocity at the event horizon. For these solutions the eigenfunctions are exponentially suppressed at spatial infinity:

[TABLE]

where $\nu=M(2\omega^{2}-\mu^{2})/\sqrt{\mu^{2}-\omega^{2}}$. In the small-coupling limit, $M\mu\ll 1$, these solutions are well approximated by a hydrogenic spectrum Detweiler:1980uk ; Dolan:2007mj with angular dependence governed by the spherical harmonics $Y_{lm}(\theta,\phi)$, angular separation constant $\lambda\simeq{l}({l}+1)$, and frequency

[TABLE]

where we introduced the dimensionless spin parameter $\chi=a/M=J/M^{2}$ and the coefficient $\gamma_{lm}$ is defined by the following relation:

[TABLE]

with $C_{l}=\frac{2^{4l+1}\left(2l+1\right)!}{\left(l+1\right)^{2l+4}}\left[\frac{l!}{2l!\left(2l+1\right)!}\right]^{2}$, $C_{1}=1/48$ for the dominant unstable ${l}=1$ mode. From Eq. (4) it is clear that these modes become unstable ($\omega_{I}>0$) whenever $\omega_{R}<m\Omega_{H}$, with an instability time scale roughly given by the $e$-folding time, $\tau_{lm}=1/\omega_{I}$, which strongly depends on the gravitational coupling $M\mu$, dimensionless spin $\chi$, and quantum numbers $(l,m)$.

The critical value of the spin that saturates the superradiance condition reads

[TABLE]

In particular, for positive frequencies and spin, $m>0$ (i.e., a mode corotating with the BH) is a necessary but not sufficient condition for the instability.

In the small-$M\mu$ limit, the radial eigenfunctions read Brito:2014wla ; Detweiler:1980gk ; Yoshino:2015nsa :

[TABLE]

where $g_{l}(r)$ can be written in terms of Laguerre polynomials:

[TABLE]

The eigenfunction peaks at

[TABLE]

and thus extends well beyond the horizon, where rotation effects can be neglected.

II.2 Multiple modes

Using an ansatz of the form (2), we consider a generic superposition of monochromatic modes as

[TABLE]

where $P_{lm}$ are the Legendre polynomials. For concreteness, we focus on the lowest-lying modes with the shortest absorption or instability timescales beyond the $l=m=1$ mode, and consider three different cases ((cf. Table 1):

[TABLE]

where $\Psi_{11}=A_{11}\,g_{1}(r)\cos\left(\phi-\omega_{R}t\right)\sin\theta$. The first case above (dubbed Model I) corresponds to the superposition of two modes with $l=1$ and $m=\pm 1$, respectively. The second one (dubbed Model II) corresponds to two modes with $l=m=1,2$, whereas the third case (dubbed Model III) corresponds to the superposition of two modes with the same $m=1$ but $l=1,2$, respectively. Note that $\omega_{R}\sim\mu$ for all modes when $M\mu\ll 1$.

It is convenient to express the initial amplitudes $A_{lm}$ in terms of the total mass of the condensate and the relative amplitude between modes. By introducing the scalar-condensate mass computed in the flat spacetime approximation (justified in the $M\mu\ll 1$ limit Brito:2014wla ; Herdeiro:2017phl ),

[TABLE]

we obtain

[TABLE]

for the three above cases, respectively, where we have introduced the relative amplitudes

[TABLE]

Thus, each initial state is defined by $M_{S}$ and by one of the $\lambda_{i}$’s, with $\lambda_{i}\to 0$ being the single $l=m=1$ mode limit.

II.3 Energy and angular momentum fluxes

II.3.1 GW emission from the scalar condensate

The scalar condensate is a source of GWs. Even though the cloud is nonrelativistic, the quadrupole approximation does not apply because the emission is incoherent Brito:2014wla ; Yoshino:2015nsa . Indeed, a dipolar scalar condensate would emit quadrupolar GWs at a frequency $\omega=2\omega_{R}\sim 2\mu$, whose wavelength $\sim 1/\omega$ is generally smaller than the size of the source, $r_{\mathrm{cloud}}$. Thus, computing the GW emission requires a fully relativistic computation using the Teukolsky formalism Brito:2014wla ; Yoshino:2015nsa . This is performed in Appendix A, we report here the final result.

For Model I ($l=1$, $m=\pm 1$), we get

[TABLE]

for the GW energy and angular-momentum fluxes, respectively. Clearly, if $\lambda_{1}\to 0$, both expressions reduce to those of the single-mode case with $l=m=1$ Brito:2014wla . In the opposite limit, $\lambda_{1}\gg 1$, the $m=-1$ mode dominates. This corresponds to the same energy flux but, from Eq. (20), the angular-momentum flux has the opposite sign, as expected for a counter-rotating mode. Thus, the angular momentum variation can be negative when $\lambda_{1}>1$, i.e. when the initial amplitude of the counter-rotating mode is bigger than that of the corotating one.

For Model II ($l=m=1,2$), we get

[TABLE]

In this case the angular momentum variation is always positive because both modes corotate with the BH. When $\lambda_{2}\to 0$ we retrieve the $l=m=1$ single-mode case, whilst for $\lambda_{2}\to\infty$ both expressions are suppressed to leading order in $M\mu\ll 1$.

Finally, for Model III ($m=1$, $l=1,2$), we get

[TABLE]

Again, as $\lambda_{3}\to 0$ we obtain the single mode case, whereas if $\lambda_{3}\to\infty$ both expressions are suppressed to leading order in $M\mu\ll 1$. We always neglect the GW energy flux at the horizon, which is typically subdominant Poisson:1994yf .

II.3.2 Superradiant evolution of the scalar condensate

In the quasi-adiabatic approximation, we assume that the energy and angular-momentum fluxes of the condensate at the BH horizon ($\dot{E}_{S},\dot{J}_{S}$) are entirely converted into the growth of the total scalar-cloud mass and angular momentum Brito:2014wla

[TABLE]

where ${L}_{S}$ is the $z$ component of the angular momentum of the cloud. The computation of $\dot{M}_{S}$ and $\dot{L}_{S}$ is performed in Appendix B. The procedure can be summarized as follows:

(i) include an adiabatic time dependence $A_{lm}\to A_{lm}e^{t/\tau_{lm}}$ in the expression for the eigenfunction $\Psi$; Eq. (10). Note, that this yields a time dependence of the relative amplitudes $\lambda_{i}$; cf. Eq. (18) and Eqs. (33)-(35) below;

(ii) compute $M_{S}(t)$ and ${L}_{S}(t)$ and their corresponding time derivative in the $\omega_{I}\ll\omega_{R}$ limit;

(iii) average the final result over several orbital periods, $T=2\pi/\omega_{R}$, of the scalar cloud.

We report here the final result.

For Model I ($l=1$, $m=\pm 1$), we obtain

[TABLE]

where $\left\langle...\right\rangle$ is the time average over several orbital periods and we defined $\omega_{lm}\equiv\omega_{I}$ for a given $(l,m)$. As discussed in detail below, a crucial point is that $\left\langle\dot{E}_{S}\right\rangle<0$ when $\lambda_{1}$ is sufficiently large, because $\omega_{1-1}<0$.

For Model II ($l=m=1,2$), we obtain

[TABLE]

Finally, for Model III ($m=1$, $l=1,2$), we obtain

[TABLE]

As a consistency check, we note that all above expressions reduce to the expected limits when $\lambda_{i}\to 0$ or $\lambda_{i}\to\infty$.

In the adiabatic approximation, the time dependence of $\lambda_{i}$ can be obtained from Eq. (18) and reads

[TABLE]

III Quasi-adiabatic evolution

We are now in the position to study the quasi-adiabatic evolution of the BH-scalar condensate in the presence of multiple modes. Using conservation of the total energy and angular momentum, the evolution of the system is described by Brito:2014wla :

[TABLE]

We integrated this simple set of ODEs using Mathematica with its built-in NDSolve function. Using the first relation of Eq. (36), we define the residuals of this procedure as ${\rm Residuals}\equiv\dot{M}+\dot{M}_{S}+\dot{E}_{{\rm GW}}$. We varied the working precision of the NDSolve function, which corresponds to changing the resolution of the numerical scheme, and verified that the residuals remain $\lesssim\mathcal{O}(10^{-17})$. At variance with Ref. Brito:2014wla , in Eqs. (36) we neglected mass and angular momentum accretion by ordinary matter (e.g., from an accretion disk), since the latter play a marginal role in the evolution of the system and are not crucial for our purposes. Indeed, accretion simply introduces an extra time scale in the problem, associated with the Salpeter time, $\tau_{\rm accretion}\sim\sigma_{T}/(4\pi m_{p})\sim 4.5\times 10^{7}\,{\rm yr}$, where $\sigma_{T}$ and $m_{p}$ are the Thompson cross section and the proton mass, respectively. The BH mass and spin grow through accretion approximately over this time scale. Thus, a slowly-spinning BH that does not satisfy the superradiant condition can be brought into a superradiant phase through accretion Brito:2014wla . Here, for simplicity we consider only systems which at $t=0$ satisfy the superradiant condition; including accretion is a straightforward extension that should not affect our overall conclusion. Within our framework, the evolution depends on the dimensionless parameters $M_{0}\mu\equiv M(t=0)\mu$, $\chi_{0}\equiv\chi(t=0)$, where $\chi(t)\equiv J(t)/M(t)^{2}$ is the dimensionless BH spin parameter, and $\lambda_{i,0}\equiv\lambda_{i}(t=0)$ (with $i=1,2,3$ depending on the model), as well as on the initial scalar cloud mass $M_{S0}\equiv M_{S}(t=0)$. The initial angular momentum of the cloud, $L_{S0}\equiv L_{S}(0)$ is determined in terms of $M_{S0}$, $\lambda_{i,0}$ and $\chi_{0}$, as discussed in Appendix B.

For concreteness333The mass simply sets the scale of the problem and we could have chosen units such that $M=1$. It is straightforward to consider different values of the BH mass by rescaling all dimensionful quantities accordingly. For example, the evolution for $m_{B}=10^{-18}\,{\rm eV}$ and $M=10^{7}M_{\odot}$ is equivalent to the evolution for $m_{B}=10^{-12}\,{\rm eV}$ and $M_{0}=10M_{\odot}$ after rescaling the time coordinate by a factor $10^{-6}$., we choose to present the numerical results in this section for a BH with mass $M_{0}=10^{7}M_{\odot}$ and consider different cases (see also Table 2):

• •

Case A: The scalar field has a mass $m_{B}=10^{-18}\,{\rm eV}$, corresponding to an initial gravitational coupling $M_{0}\mu\sim 0.075$. The initial BH spin is either $\chi_{0}=0.8$ or $\chi_{0}=0.95$. In both cases the BH is initially in a superradiant state, $\Omega_{H}>\mu/m$. The initial seed has $M_{S0}=10^{-9}M_{0}$. This is representative for the case $M_{S0}\ll M_{0}$, which includes seeds due to quantum fluctuations.

• •

Case B: This is same as Case A above but for a seed with larger initial mass $M_{S0}=0.025M_{0}$, which should model a seed of astrophysical origin, since the energy of the perturbation is a sizeable fraction of the initial BH mass. This case is expected when a scalar cloud is present in the environment in which the BH forms, or during the coalescence of a binary BH each one endowed with its own scalar cloud.

• •

Case C: This is same as Case B above but for a scalar field with slightly larger mass, $m_{B}=4\times 10^{-18}\,{\rm eV}$, corresponding to an initial gravitational coupling $M_{0}\mu\sim 0.3$. In this case we set the initial spin to $\chi_{0}=0.95$ in order to satisfy the superradiant condition initially.

Case A was chosen to agree with the case considered in Ref. Brito:2014wla , whereas Case C is representative for the initial data evolved numerically in Ref. Okawa:2014nda . Note that this case is only marginally consistent with our small-coupling approximation, $M\mu\ll 1$. In all cases, the BH and boson-field masses correspond to a range that will be accessible by future LISA observations Arvanitaki:2014wva ; Baryakhtar:2017ngi ; Brito:2017zvb ; Brito:2017wnc .

III.1 Model I: $l=1$ with $m=\pm 1$

III.1.1 Case A: small initial seed

A representative example of the evolution for Model I in Case A is presented in Fig. 1. In this case the evolution is insensitive to the presence of a nonsuperradiant unstable mode, even if the latter has initially a much larger amplitude than the superradiant mode (e.g., $\lambda_{1}(0)=10$). This is due to Eq. (33): since $\omega_{1-1}<0$ for a nonsuperradiant mode and $\omega_{11}>0$ for the superradiant one, the relative amplitude $\lambda_{1}(t)$ decreases exponentially over a time scale $\sim 1/(\omega_{1-1}-\omega_{11})$. The evolution is then only affected by the superradiant mode and proceeds as in the $l=m=1$ single-mode case Brito:2014wla . For the chosen parameters we find $\tau_{11}\sim 7\times 10^{6}\,{\rm yr}$ which is consistent with the exponential growth of the condensate at $t>10^{7}\,{\rm yr}$ shown in Fig. 1. In this particular case, the condensate extracts $\approx 4\%$ of the initial BH mass. However, its energy density, and hence its backreaction, is negligible Brito:2014wla so that it dissipates on a longer time scale through GW emission. An estimate for the latter time scale is

[TABLE]

in agreement with the late-time behavior shown in Fig. 1.

III.1.2 Case B: large initial seed, small coupling

The impact of a nonsuperradiant mode is stronger when the initial seed has a larger amplitude as illustrated in Fig. 2. There, we show the evolution of Model I for fixed $\lambda_{1,0}=1$ and different initial scalar cloud masses $M_{S0}$, including Case A and Case B.

The case $M_{S0}=10^{-9}M_{0}$ corresponds to the case $\lambda_{1}(0)=1$ shown in Fig. 1 whereas, as $M_{S0}$ increases, we observe further features. Parts of the scalar cloud are initially absorbed by the BH, whose mass, in turn, grows in time. This can be understood as follows. Neglecting for the moment GW emission, system (36) reduces to

[TABLE]

and therefore, when $\langle{\dot{E}_{S}}\rangle<0$ and $\langle{\dot{J}_{S}}\rangle<0$, the mass and angular momentum of the condensate decrease while the BH mass and spin increase, respectively. From Eqs. (27) and (28), this can never happen when $\lambda_{1}=0$, because in that case $\langle{\dot{E}_{S}}\rangle,\langle{\dot{J}_{S}}\rangle\propto\omega_{11}$ and the scalar fluxes are positive in the superradiant phase (i.e., when $\omega_{11}>0$). In this case the instability halts as the superradiant condition is saturated (i.e., as $\Omega_{H}\to\mu/m$ or, equivalently, as $\omega_{11}\to 0$). However, the situation is different when $\lambda_{1}\neq 0$. When $\lambda_{1}\ll 1$, Eqs. (27) and (28) reduce to

[TABLE]

and therefore the energy and angular-momentum fluxes are smaller than in the single-mode case as long as $\omega_{11}>|\omega_{1-1}|$.

On the other hand, if $\lambda_{1,0}$ is sufficiently large, it might happen that $\omega_{11}+\lambda^{2}_{1}\omega_{1-1}<0$ and, therefore, the scalar energy flux is negative; see Eq. (27). Even when this happens, Eq. (33) shows that $\lambda_{1}(t)$ decreases exponentially. As shown in Fig. 3 (solid black curve) when the initial seed mass is negligible the scalar flux can be negative at $t=0$, but then it turns positive (on a time scale time scale $1/|\omega_{1-1}-\omega_{11}|$) as $\lambda_{1}(t)\to 0$. Therefore, the usual superradiant evolution is not affected by the presence of a nonsuperradiant mode as long as the initial seed mass is small.

Larger scalar cloud masses, however, enhance the negative energy flux. This dependence is illustrated in Fig. 3, where we show the energy flux (rescaled by $M_{\rm S}\mu$) for $\lambda_{1,0}=1$ and different initial cloud masses. In particular, clouds with intermediate values of $M_{S0}$ such as Case B (green dashed line in Figs. 2 and 3) may be partly absorbed, but despite a small increase in the BH mass and decrease of the BH spin the energy flux becomes positive, i.e., the system reaches the superradiant regime. This picture changes dramatically for larger cloud masses $M_{S0}$; see, e.g., blue dotted line in Figs. 2 and 3. In such case, the negative scalar flux is significantly enhanced, a sizeable fraction of the scalar cloud – including its counterrotating modes – is absorbed by the BH and the superradiant phase can be highly suppressed or entirely absent. The details of the evolution sensitively depend on both the initial relative amplitude $\lambda_{1,0}$ and scalar cloud mass $M_{S0}$.

Meanwhile, the angular momentum of the scalar cloud grows because $\langle\dot{L}_{S}\rangle$ remains positive (see Eq. (28)). As a result, the BH angular momentum decreases irrespective of being in the superradiant phase. Interestingly, the final BH spin seems to be similar to the single-mode case. However, a careful investigation of the spin evolution presented in Fig. 5 shows that this is a coincidence. For large enough $\lambda_{1,0}$ the final spin depends crucially on the initial parameters.

III.1.3 Case C: large initial seed, large coupling

When the gravitational coupling $M\mu$ and spin increase, the impact of the seed mass becomes even more relevant. As shown in Fig. 4, the presence of a counterrotating mode with $m=-1$ reduces the superradiant energy extraction. For sufficiently large values of $\lambda_{1,0}$ – in the present case $\lambda_{1,0}\geq 1$ – the scalar cloud never grows, since the BH absorbs it before superradiance could kick in. At the same time, the BH angular momentum decreases because the absorbed energy is mostly contained in a counter-rotating mode. Comparing to the critical value $\chi_{\rm crit}\sim 0.882$ (see Eq. (6)), we observe that the system can be driven out of the superradiant regime early in the evolution and therefore does not undergo a superradiant phase.

This behaviour agrees with the expectation raised by fully nonlinear simulations Okawa:2014nda . We show their setup KGl_m30_a3 (cf. Table III of Okawa:2014nda ) as blue line in Fig. 4. While the immediate response differs due to nonlinear (backreaction) effects, we find excellent agreement within $\lesssim 0.5\%$ in the BH mass and spin with the adiabatic evolution at late times.

In these cases the final BH spin is not only driven by superradiance but also the absorption of counterrotating ($m=-1$) modes. Hence, it depends on the initial parameters as illustrated in Fig. 5. Here we show the final dimensionless BH spin (relative to its initial value) as a function of $\lambda_{1,0}$ for different initial cloud masses $M_{S0}$, for Cases A and B (Fig. 5a) and Case C (Fig. 5b). Small perturbations (red solid lines) always yield the superradiant evolution, i.e., BHs whose final spin is smaller than its initial one due to the superradiant instability independently of the presence of a counterrotating mode. The dependence is more complex for large initial scalar clouds: in an intermediate regime, around $\lambda_{1,0}\lesssim 1$ accretion of counterrotating (i.e.,nonsuperradiant) modes and superradiant scattering compete, potentially leading to a larger final spin. Instead, if the initial condensate is dominated by the $m=-1$ mode, the evolution is dominated by accretion of the counterrotating component and can yield considerably smaller final spins.

Summary: To summarize, our quasi-adiabatic evolution for Model I reveals that the BH superradiant instability proceeds as in the case of a single superradiant mode whenever the seed’s energy is negligible (as in the case of quantum fluctuations), whereas the dynamics and the final BH spin are strongly affected by the addition of a nonsuperradiant mode if the latter has a large amplitude relative to the superradiant one and if the initial scalar cloud has a nonnegligible energy. In some extreme cases, the absorption of the counter-rotating superradiant mode is sufficient to reduce the BH angular momentum past the superradiant condition, so that the instability is completely quenched.

III.2 Model II: $l=m=2$ and $l=m=1$

The phenomenology of Model II is different from that of Model I. In particular, both modes can trigger the superradiant instability, albeit on vastly different time scales since $\tau\sim(M\mu)^{4l+5}$ depends strongly on $l$.

We first focus on Case A, i.e. small initial fluctuations, whose evolution is presented in Fig. 6 for different relative amplitudes $\lambda_{2,0}$. We observe two unstable phases: the first one occurring on a time scale $1/\omega_{11}$ and the other occurring on longer scales $1/\omega_{22}$. Because of this separation in time scales, the evolution starts with the first superradiant phase in which the scalar cloud grows and the BH spins down. Then, the cloud is dissipated through GW emission. Finally, the $l=m=2$ mode becomes unstable and the scalar cloud grows again, with the BH spin further decreasing since the superradiant threshold $\mu\approx m\Omega$ implies a smaller final spin; cf. Eq. (6).

Looking at Fig. 6 it is easy to notice that, regardless of the value of $\lambda_{2,0}\neq 0$, the end-state of the system remains unchanged, i.e. the values of the final BH spin and mass are an attractor of the dynamics. In order to better understand this effect, we study the time evolution of $\lambda_{2}(t)$, which is shown in Fig. 7. After a first depletion occurring at $t\sim\tau_{11}$ due to the superradiant instability induced by the $l=m=1$ mode, $\lambda_{2}$ undergoes an exponential divergence at $t\sim\tau_{22}$. This can be seen from the definition of $\lambda_{2}(t)$, Eq. (34): At $t\sim\tau_{11}$ the system reaches the superradiant threshold of the $l=m=1$ mode, for which $\omega_{11}=0$. Afterwards, at $t\sim\tau_{22}$, the secondary mode $l=m=2$ kicks in and takes the system to the superradiant threshold for the $m=2$ mode, which is saturated when $\omega_{22}=0$. In this situation, however, $\omega_{11}$ becomes negative and $\lambda_{2}(t)$ diverges.

Note that the time scale associated with GW dissipation is much longer for $l=m=2$ Yoshino:2015nsa , which explains the long time before the condensate disappears (not shown in Fig. 6). For the system under consideration, these long time scales force us to consider an evolution that lasts much longer than the age of the universe, see Fig. 6. However, if we would consider a stellar-mass BH with $M(0)=10M_{\odot}$ the time scales would be $10^{6}$ times smaller, since all dimensionful quantities scale with the initial BH mass. Then also the secondary superradiant phase might occur within the age of the universe. Model II–Case A is therefore a straightforward interpolation between the case of a single mode with $l=m=1$ and that of a single mode with $l=m=2$.

We now focus our attention to the influence of a larger initial scalar cloud. Its evolution is illustrated in Fig. 8 for various $M_{S0}$ and equal initial amplitude of the $l=m=1$ and $l=m=2$ modes, i.e., $\lambda_{2,0}=1$.

As before, we observe the growth of the scalar cloud at the expense of the BH mass and angular momentum on timescales $1/\omega_{11}$, i.e., due to the $l=m=1$ instability. In contrast to Model I, this process is essentially independent of the cloud’s initial mass since the influence of the secondary mode kicks in on significantly longer time scales $1/\omega_{22}$. Once the $m=1$ superradiant threshold is reached, the scalar condensate dissipates via GW emission. Towards the end of this process, after about $t\sim 10^{12}$yr in our setup, the BH mass increases; see Fig. 8b. This indicates that the scalar cloud is accreted onto the BH.

To better understand this process, let us inspect the energy flux $\sim\omega_{11}+81\lambda_{2}^{2}\omega_{22}$; see Eq. (29) and Fig. 9. Let us also recall the difference in time scales $\omega_{11}\sim(M\mu)^{-9}\gg\omega_{22}\sim(M\mu)^{-13}$. That is, during the early evolution $\omega_{11}>0$ dominates, thus triggering the $m=1$ superradiant instability where the scalar flux is positive and peaks around $\tau\sim 10^{6}$yr as shown in Fig. 9. The scalar’s growth stops as the superradiant threshold is reached where $\omega_{11}=0$, and the clouds starts dissipating. Now, although the secondary mode is still growing at a rate $\sim 1/\omega_{22}$ (recall that $\omega_{22}>0$ is positive), the primary mode starts decaying with a rate $\omega_{11}<0$. The latter can dominate and lead to a negative scalar energy flux as shown in the inset of Fig. 9. This is consistent with the observation of increasing BH mass in Fig. 8b.

In the meantime the relative amplitude between the two modes, $\lambda_{2}(t)\sim\exp[(\omega_{22}-\omega_{11})t]$, is growing exponentially, see Eq. (34). So, eventually the second term in Eq. (29), which is positive, will cancel and then dominate over the $m=1$ contribution. At this point the scalar energy flux is positive, see inset of Fig. 9, and the BH-scalar cloud configuration undergoes its second (i.e. $m=2$) superradiant phase. It sets in after about $\tau\sim 10^{13}-10^{15}$yr, with the specific onset depending on the scalar clous mass; see Fig. 8a. As expected, the scalar cloud grows by dipping into the BH mass and angular momentum that further decreases the final BH spin. Eventually the scalar cloud will dissipate via GW emission on time scales much longer than shown in Fig. 8.

III.3 Model III: $l=1,2$ with $m=1$

Model III is qualitatively similar to Model II. In particular, this model interpolates between the single-mode case with $l=m=1$ (when $\lambda_{3}\approx 0$) and the single-mode case with $l=2$, $m=1$ (when $\lambda_{3}\to\infty$). Also in this case the transition occurs for large values of $\lambda_{3}$, because $|\omega_{21}|\ll|\omega_{11}|$. The only qualitative difference is related to the critical value of the spin, which is the same in both regimes, since the modes have the same azimuthal number $m$ and Eq. (6) does not depend on $l$.

IV Implications

Due to BH no-hair theorems for real bosonic fields Chrusciel:2012jk ; Herdeiro:2015waa , the end-state of the evolution must be a Kerr BH, the condensate being eventually dissipated in GWs. However, an interesting question concerns the final value of the BH spin in the new stationary configuration and, more generically, the phase-space (Regge-plane) of the final BH.

Another relevant question concerns the viability of our Cases B and C, where the energy of the initial seed is a sizeable fraction (roughly $2\%$) of the BH mass. This scenario could occur if the BH is formed in a scalar-rich environment, for example if it is formed out of the merger of two previously scalarized BHs. The time scale for GW dissipation of the condensate strongly depends on the coupling Arvanitaki:2014wva ; Brito:2017zvb ; Brito:2017wnc . For the fundamental $l=m=1$ mode,

[TABLE]

Thus, depending on the mass and spin of the BH and on the mass of the bosonic field, $\tau_{\rm GW}$ can easily exceed the age of the universe, in agreement with Eq. (37) above. In that case the condensate will not have enough time to dissipate during the coalescence, and the merger remnant will form in an environment where the energy of the scalar field is not negligible.

Although this scenario might be relevant only for a fraction of sources, the majority of supermassive BHs are believed to form via hierarchical mergers. Thus, they constitute a sizeable fraction of the GW signal from bosonic condensates which may potentially be detected by LISA Brito:2017zvb ; Brito:2017wnc . We leave a more quantitative analysis of possible formation scenarios and event rates for future work.

IV.1 Spin evolution

Single mode: For reference, let us recall the final spin resulting from the evolution of the single, $l=m=1$ mode. Focusing on the initial parameters used in our previous sections, the final spin is $\chi(t\to\infty)\approx 0.28$ and $\chi(t\to\infty)\approx 0.86$ for Case A and Case C, respectively.

Model I: The presence of a counter-rotating mode can significantly change the value of the final spin, and details depend on the initial parameters. To better understand those dependencies we present the ratio between the final and initial spin as function of the initial relative amplitude in Fig. 5. They are shown for different (initial) scalar cloud masses and fixing ($\chi_{0}=0.8$, $M_{0}\mu=0.075$), i.e., Model I–Case A and B (see Fig. 5a) or ($\chi_{0}=0.95$, $M_{0}\mu=0.3$), i.e., Model I–Case C (see Fig. 5b).

If we consider only small fluctuations this ratio remains constant, i.e., is insensitive to the presence of a counterrotating mode and yields the same final spins as in the single ($l=m=1$) mode case reviewed above.

Similarly, the final BH spin appears independent of the initial scalar cloud mass $M_{S0}$ as long as the relative amplitude between counter- and co-rotating modes is sufficiently small, namely $\lambda_{1,0}\lesssim 0.1$.

However, if the initial scalar cloud contains comparable excitations of the $m=\pm 1$ modes, i.e. if $\lambda_{1,0}\sim\mathcal{O}(1)$, and its mass is a few $\%$ of the BH’s initial mass, the dependency of the final spin is more complex. Now, the angular momentum flux (28) is determined by both the superradiant $m=+1$ mode and the counter-rotating $m=-1$ mode. Since $\omega_{11}>0$ and $\omega_{1-1}<0$, both modes will increase the flux and, hence, reduce the final BH spin, although the latter may be slightly larger than in the single $l=m=1$ mode case, depending on the parameters (see Fig. 5).

As we further increase $\lambda_{1,0}\gtrsim 10$ the ratio between the final and initial BH spin approaches a constant value (i.e., independent of $\lambda_{1,0}$) that can, in fact, be smaller than the $l=m=1$ case but depends strongly on the the initial scalar cloud mass $M_{S0}$. Note, that the spin can actually flip sign due to absorption of counter-rotating modes as shown in Fig. 5a.

Model II: Due to the presence of a secondary superradiant mode, the extraction of the BH spin proceeds in two stages: the first one is caused by the dipole mode and yields the same final BH spin as the single ($l=m=1$) mode case discussed above. After $t\sim 1/\omega_{22}$ the system undergoes the $l=m=2$ superradiant phase that yield further extraction of the BH spin. As indicated in Figs. 6 and 8, its final value appears independent of the initial relative amplitude $\lambda_{2,0}$ or scalar cloud mass $M_{S0}$. This is because $\lambda_{2}(t)$ will have acquired the same value (independent of its initial one) by the time the $l=m=2$ mode becomes active; see Fig. 7. For example, in the case studied here (with $\chi_{0}=0.8$), the final spin is $\chi(t\to\infty)\sim 0.14$.

IV.2 Regge planes

Let us now focus on the mass-spin phase-space of the final BH encapsulated in its Regge-plot. To identify it we performed a set of quasi-adiabatic evolutions whose results are shown in Figs. 10, 11 and 12 for Models I and II, respectively. In particular, we considered $m_{B}=10^{-18}\,{\rm eV}$ and $1000$ configurations starting at $t=0$ with a random distribution of the initial BH spin in the range $\chi_{0}\in(0,0.998)$ and masses in the range $\log_{10}M_{0}\in(6,7.5)$, so that the gravitational coupling $M_{0}\mu\in(0.0075,0.24)$.

For comparison, we show the $\lambda_{i,0}=0$ case in the top-left panel of each plot in Figs. 10– 12. Then, the final BH configuration avoids a specific region of the Regge plane, customarily dubbed Regge “gap” Arvanitaki:2010sy . For a single $(l,m)$ mode, the shape of this gap is approximately given by Brito:2014wla

[TABLE]

where the critical spin $\chi_{\rm crit}$ is given in Eq. (6), $M_{c}$ is the value of $M$ that minimizes the spin when $\tau_{lm}=t_{F}$. An approximate formula is $M_{c}=\left({\frac{l+1}{C_{l}\mu^{4l+6}t_{F}}}\right)^{\frac{1}{4l+5}}$ Brito:2014wla .

Model I: We first focus on Model I whose Regge planes are shown in Fig. 10 for initial scalar cloud masses $M_{S0}=10^{-9},0.025,0.05M_{0}$ and different relative amplitudes $\lambda_{1,0}$. For each initial configuration, we followed the evolution of the system up to $t=t_{F}=10^{8}\,{\rm yr}\gg\tau_{11}$.

The evolution of a system containing only small scalar fluctuations $M_{S0}=10^{-9}M_{0}$, depicted in Fig. 10a, is largely independent of the presence of a counter-rotating mode. In particular, it exhibits the same exclusion regions in the Regge plane as those induced by the $l=m=1$ superradiant evolution Brito:2014wla .

If, instead, the scalar cloud already stores a significant fraction of the BH’s mass – of the order of a few percent – the spin–mass phase-space of the final BH exhibits more structure; see Figs. 10b and 10c. We identify three main features:

(i) we still find gaps in the Regge plane consistent with those of the standard superradiant evolution. However, their onset occurs for smaller masses as the scalar cloud mass increases, even in the single-mode case. This can be explained considering that a bigger initial value of the scalar cloud mass implies a larger energy flux rate via Eq. (27) and, consequently, a shorter instability time scale. Furthermore, if $\lambda_{1,0}\gtrsim 10$, we start populating the low-mass end of the Regge gap. This is not surprising as the absorption of (counter-rotating) modes decreases the BH spin while increasing its mass; see e.g. Fig. 2;

(ii) we find additional gaps in the Regge plane, just below the superradiant threshold, if the scalar cloud is dominated by the $m=+1$ mode; see top panels of Figs. 10b and 10c. BHs that are below the threshold will absorb part of the predominantly co-rotating cloud whose mass is $\mathcal{O}(1)\%$ of the BH mass. Hence their mass and spin will increase towards the superradiant threshold. Should they supercede it the superradiant instability will become active and drive the system towards the threshold from above. That is to say, that the superradiant threshold appears to be an attractor if the initial scalar cloud mass is sufficiently large and dominated by potentially superradiant (i.e. $m=+1$) modes;

(iii) if the initially large scalar cloud is instead dominated by counter-rotating modes, i.e. $\lambda_{1,0}\gtrsim 10$, as shown in the right-bottom panels of Figs. 10b and 10c, these additional holes disappear. Instead we find BHs with a negative final spin relative to its initial one. This can be understood as follows: BHs that have almost vanishing initial spin will absorb the counter-rotating modes that further decrease the BHs’ spin. Interestingly, one could now view the system containing a BH with negative spin and cloud with $m=-1$ modes as co-rotating (but in the opposite direction as before), i.e., the situation is the same as that of a BH with positive spin surrounded by a $m=+1$ cloud. This, again, should suffer from the superradiant instability if the BH is spun up sufficiently. Indeed, the bottom-right panel of Fig. 10c seems to exhibit such a new attractor line.

Model II: The Regge planes for Model II are shown in Figs. 11 and 12 for an evolution time of $t_{F}=10^{8}\,{\rm yr}\gg\tau_{11}$ and $t_{F}=10^{15}\,{\rm yr}\gg\tau_{22}$, respectively. Although the latter time scale is larger than the age of the universe it allows us to explore features in the Regge plane due to both the $l=m=1$ and $l=m=2$ instability for the supermassive BHs under consideration. Note, furthermore, that this time scales with the BH mass. So for a much lighter, stellar-mass BH of $\mathcal{O}(10)M_{\odot}$ (and scalar of $m_{\rm B}\sim 10^{-12}{\rm eV}$) we would observe features of Fig. 12 after $t_{F}\sim 10^{9}$yr.

The Regge plots of Fig. 11 only exhibit the $l=m=1$ Regge gap since we evolved the systems for times that are significantly shorter than the $l=m=2$ instability timescale. If the scalar starts off as only a small fluctuation, i.e. Case A, the Regge gaps are identical to those of the single, $l=m=1$ mode. That is, they are independent of the presence of a secondary mode as shown in Fig. 11a. The Regge gap itself is consistent with the estimate (42).

In Fig. 11b we consider a larger cloud with $M_{S0}=0.025M_{0}$. As we already saw in Model I, the onset of the superradiant instability is shifted towards smaller BH masses as we increase the scalar cloud mass. In particular, the threshold $M_{c}$ is a complicate function of the parameters and differs from the expression below Eq. (42). As before, we find an additional gap below the superradiant threshold for sufficiently massive scalar clouds; see top left panel of Fig. 11b.

Let us now turn our attention to the long-time evolution depicted in Fig. 12 where we capture both the $l=m=1$ and $l=m=2$ phases. For small initial seeds, depicted in Fig. 12a, we observe the appearance of two Regge gaps consistent, respectively, with the $l=m=1$ and $l=m=2$ superradiant evolution. The details are independent of the (initial) relative amplitude $\lambda_{2,0}$ and appear as soon as the secondary mode is switched on.

The Regge plane becomes more complex as we increase the scalar cloud’s mass to a few percent of the BH mass; see Fig. 12b. In particular, we see the formation of a second gap below the superradiant threshold for $l=m=1$ if $\lambda_{2,0}=0$ and below the $l=m=2$ threshold as soon as $\lambda_{2,0}\neq 0$. Again, this can be understood as BHs that start outside the superradiant regime, but absorb mass (and angular momentum) from the scalar cloud until they reach the threshold. Finally, the critical mass parametrizing the onset of the superradiant instability decreases, now for the $l=m=2$ case.

To summarize, even if a scalar condensate surrounding a BH contains counterrotating or higher multipole modes, in all cases studied here the holes in the Regge plane persist and yield a larger and more complex exclusion region.

V Discussion

We have investigated the evolution of the BH superradiant instability against ultralight scalar fields with an initial configuration described by a superposition of modes. We focused on the case $M\mu\ll 1$ which allows for a Newtonian description of the condensate and for a quasi-adiabatic approximation due to the separation of scales between the instability time scale and the dynamical time scale of the BH.

Our analysis shows that the evolution of the superradiant instability in the presence of an initial superposition of modes is very rich and diverse. The evolution of the system depends strongly on the energy of the scalar seed and on the gravitational coupling $M\mu$. If the seed energy is a few percent of the BH mass, a BH surrounded by a mixture of superradiant and nonsuperradiant modes with comparable amplitudes might not even undergo a superradiant unstable phase, depending on the value of the boson mass.

Our analysis adds to the numerical results of Refs. Okawa:2014nda ; Witek2018inprep , where the authors explore the interplay between a highly spinning BH and massive scalars or vectors composed of multimode data and with $M\mu\sim\mathcal{O}(0.5)$. Indeed, our simple adiabatic approximation in the small-$M\mu$ limit is in remarkably good agreement with the evolution presented in Ref. Okawa:2014nda . On the other hand, if the seed energy is much smaller than a few percent of the BH mass – as in the most natural and likely scenario in which the instability is triggered by quantum fluctuations – the effect of nonsuperradiant modes is negligible.

This implies that the only case in which the evolution of the superradiant instability is affected by multiple modes is when the BH is initially surrounded by a nonnegligible scalar environment, or if it is formed out of the coalescence of two BHs merging with their own scalar clouds. This latter scenario might be relevant only for a fraction of sources, in particular for massive BHs formed out of the merger of two BHs surrounded by their own condensates. In certain cases the time scale for GW dissipation of the condensate can exceed the age of the universe so a BH might form in a scalar-rich environment. This might be relevant for searches of ultralight fields with LISA Brito:2017zvb ; Brito:2017wnc , since supermassive BHs are expected to form hierarchically. In these cases the initial configuration of ultralight fields around BHs is generically a superposition of (superradiant and nonsuperradiant) modes, and the initial mass of the scalar configuration might be large enough to suppress the instability. We leave a more detailed analysis of such a binary and event rate estimates for future work.

Likewise, the BH Regge plane is also affected by the presence of nonsuperradiant modes when the initial scalar mass is a sizeable fraction of the BH mass. The pattern of the Regge holes is more involved and additional forbidden regions can appear, depending on the parameters. Interestingly, the region forbidden in the single-mode case is also forbidden in the presence of nonsuperradiant modes, i.e. the original Regge holes are not populated even when the superradiant instability is absent. This is due the absorption of large counter-rotating modes which decrease the BH spin.

Our analysis can be extended in several directions. We have neglected mode mixing and possible transfer of energy between modes (e.g., turbulence) which might significantly change the overall picture. We have also neglected scalar self-interactions which – if sufficiently strong – are known to quench the instability and give rise to interesting nonlinear effects such as “bosenovas” Yoshino:2012kn ; Yoshino:2015nsa . Likewise, we have neglected axion-like couplings to the electromagnetic field, which might also quench the instability through a different channel Boskovic:2018lkj ; Ikeda:2018nhb . We have also neglected accretion of ordinary matter; in light of the analysis of Ref. Brito:2014wla , we expect that including accretion should be a straightforward extension that would not give a substantial contribution to the understanding of the problem. Furthermore, although we focused on scalar fields, it is likely that the qualitative features of the evolution will be the same also for massive vector (Proca) and massive tensor fields, as indicated by nonlinear simulations that will appear soon Witek2018inprep .

Finally, a natural extension of our work is to investigate whether the presence of multiple modes can also suppress the ergoregion instability of BH mimickers Cardoso:2007az ; Cardoso:2008kj ; Pani:2010jz ; Maggio:2017ivp ; Maggio:2018ivz , since the latter shares Superradiance many features with the superradiant instability discussed here.

Acknowledgements.

We are indebted to Emanuele Berti, Roberto Emparan, and Vitor Cardoso for important comments on a first draft version of this work and to Richard Brito for relevant comments on the the current draft. G.F. wishes to thank King’s College London for financial support and the Royal Society for the PhD studentship provided under Research Grant RGF\R1\180073. P.P. acknowledges financial support provided under the European Union’s H2020 ERC, Starting Grant agreement no. DarkGRA–757480 and support from the Amaldi Research Center funded by the MIUR program “Dipartimento di Eccellenza” (CUP: B81I18001170001). H.W. was supported by the European Union’s H2020 research and innovation program under Marie Sklodowska-Curie grant agreement No. BHstabNL-655360 and acknowledges financial support provided by the Royal Society University Research Fellowship UF160547 and Royal Society Research Grant RGF\R1\180073. The authors would like to acknowledge networking support by the COST Action CA16104. H.W. thanks the Yukawa Institute for Theoretical Physics at Kyoto University for their hospitality during the workshop YITP-T-17-02 on “Gravity and Cosmology 2018” and the YKIS2018a symposium on “General Relativity – The Next Generation”. We thankfully acknowledge the computer resources at Marenostrum IV, Finis Terrae II and LaPalma and the technical support provided by the Barcelona Supercomputing Center via the PRACE grant Tier-0 PPFPWG, and via the BSC/RES grants AECT-2017-2-0011, AECT-2017-3-0009 and AECT-2018-1-0014. The authors thankfully acknowledge the computer resources and the technical support provided by the PRACE Grant No. 2018194669 “FunPhysGW: Fundamental Physics in the era of gravitational waves”, STFC DiRAC Grant No. ACTP186 “Extreme Gravity and Gravitational Waves” and STFC DiRAC Grant No. ACSP191 “Exploring fundamental fields with strong gravity”.

Appendix A GW emission from the scalar condensate

Owing to the separation of scales between the size of the cloud and the BH size for $M\mu\ll 1$, the GW emission can be approximately analyzed taking the source to lie in a flat 444We note that the flat-spacetime approximation yields a different prefactor for the GW fluxes emitted from the cloud relative to the case in which the background spacetime is described by a Schwarzschild metric Yoshino:2013ofa ; Brito:2013yxa . The difference between the two cases is small and we adopt here a flat-spacetime approximation for simplicity. background Yoshino:2013ofa . Because the source is incoherent $1/\omega\ll r_{\rm cloud}$, the quadrupolar approximation fails. In the fully relativistic regime, the gravitational radiation generated is best described by the Teukolsky formalism for gravitational perturbations Teukolsky:1973ha .

A.1 General two modes case: $(l,m)$ and $(l^{\prime},m^{\prime})$

The gravitational radiation is described by the Newman-Penrose scalar $\psi_{4}$, which, in the flat spacetime approximation, can be decomposed as

[TABLE]

where the radial function $R_{jk\omega}(r)$ satisfies the inhomogeneous Teukolsky equation,

[TABLE]

The source term $T_{jk\omega}$ is given by Poisson:1993vp

[TABLE]

where we have defined $\mathcal{L}\equiv\partial_{r}+i\omega$ and

[TABLE]

where $\Theta_{S}={T}_{nn},\,{T}_{n\bar{m}},{T}_{\bar{m}\bar{m}}$ for $S=0,\,-1,\,-2$, respectively.

The source term $T_{jk\omega}$ is related to the scalar field stress-energy tensor ${T}_{\mu\nu}$ through the tetrad projections

[TABLE]

where

[TABLE]

For a scalar configuration with two modes with $(l,m)$ and $(l^{\prime},m^{\prime})$, the contributions to the source term are given by a sum over several active modes $(j,k)$, defined by the non vanishing contributions of the integrals (46) that are strictly dependent on the values of $(l,m)$ and $(l^{\prime},m^{\prime})$. The contributions will feature two frequencies $\omega=\pm 2\omega_{R}$, due to the fact that $T_{jk\omega}$, computed through Eq. (45), contains only terms $\propto\delta\left(\omega\pm 2\omega_{R}\right)$.

Once the source term is known, the radial equation (44) can be solved using the Green’s function. The latter can be found by considering two linearly independent solutions of the homogeneous equation associated with Eq. (44), with the following asymptotic behavior Sasaki:2003xr ,

[TABLE]

[TABLE]

where $k=\omega-m\Omega_{H}$, $\left\{A,B\right\}_{in,out}$ are constants. Owing to the flat spacetime approximation, the tortoise coordinate usually defined to deal with these kind of problems coincides with the standard radial coordinate.

Imposing ingoing boundary conditions at the horizon and outgoing boundary conditions at infinity, one finds that the solution of Eq. (44) is given by Sasaki:2003xr

[TABLE]

where $W=\left(R^{\infty}\partial_{r}R^{H}-R^{H}\partial_{r}R^{\infty}\right)/r=2i\omega B_{in}$ is the Wronskian, which is a constant by virtue of the homogeneous Teukolsky equation. From the asymptotic solution of Eq. (44) we find

[TABLE]

where $C_{1}$ is an arbitrary constant that we set to unity without loss of generality. The solution $R^{H}$ can be found through

[TABLE]

where $\psi^{H}$ is the Regge-Wheeler function that at small frequencies reads

[TABLE]

where $j_{j}$ are the spherical Bessel functions of the first kind. At radial infinity the solutions reads

[TABLE]

Since the frequency spectrum of the source $T_{jk\omega}$ is discrete with frequencies $\omega=\pm 2\omega_{R}$, $\tilde{Z}_{jk\omega}^{\infty}$ can be written as

[TABLE]

where $\omega_{1}=2\omega_{R}$ and $\omega_{2}=-2\omega_{R}$. Replacing the above equation in Eq. (43), we obtain $\psi_{4}$ at radial infinity,

[TABLE]

which can be written as

[TABLE]

where $h_{+}$ and $h_{\times}$ are the two independent GW polarizations. Then, using Eq. (60) in the previous relation and integrating twice with respect to the time, we obtain the gravitational waveform,

[TABLE]

The energy flux carried by these waves at infinity is given by Teukolsky:1974yv

[TABLE]

Finally, combining the last two equations, we get the energy and angular momentum fluxes at radial infinity Hughes:2001jr

[TABLE]

A.2 Particular cases

We shall now apply the above results to the three models presented in the main text.

A.2.0.1 Model I

. For the scalar configuration (11) the contributions to the source term are given by $(j=2,k=\pm 2)$ with two frequencies $\omega=\pm 2\omega_{R}$, since $T_{jk\omega}$ contains only terms $\propto\delta\left(\omega\pm 2\omega_{R}\right)$. Then, the right-hand sides of Eqs. (64) and (65) take the form

[TABLE]

In this particular case, using Eq. (11) and considering the small $M\mu$-limit, we obtain

[TABLE]

Finally, using Eq. (15) we get Eqs. (19) and (20).

A.2.0.2 Model II

. For the scalar configuration (12) in Eq. (45), we have different contributions relative to $(j=2,k=\pm 2)$, $(j=3,k=\pm 3)$ and $(j=4,k=\pm 4)$. Furthermore, the contributions with $k>0$ are $\propto\delta\left(\omega-2\omega_{R}\right)$, while those with $k<0$ are $\propto\delta\left(\omega+2\omega_{R}\right)$. In this case the right-hand sides of Eqs. (64) and (65) become

[TABLE]

Using Eqs. (12), (68) and (69) for $M\mu\ll 1$, we get

[TABLE]

Finally, using Eq. (16), the above equations reduce to

[TABLE]

with

[TABLE]

The contributions $\propto\lambda_{2}$ in the numerator are subleading when $M\mu\ll 1$. This implies that, in the considered limit, $\dot{E}_{GW},\dot{J}_{GW}\propto(M_{S}/M)^{2}(M\mu)^{14}$. In this limit, the GW energy and angular momentum fluxes are given by Eqs. (21) and (22).

A.2.0.3 Model III

. The stress-energy tensor corresponding to the scalar configuration (12) in Eq. (45) yields contributions corresponding to $(j=2,k=\pm 2)$, $(j=3,k=\pm 3)$ and $(j=4,k=\pm 4)$. Again, the ones with $k>0$ are $\propto\delta\left(\omega-2\omega_{R}\right)$, while for $k<0$ they are $\propto\delta\left(\omega+2\omega_{R}\right)$. The analysis is the same as for Model II above; the right-hand sides of Eqs. (64) and (65) are given by Eqs. (68) and (69).

Finally, considering Eqs. (13), (68) and (69) for $M\mu\ll 1$, we get:

[TABLE]

Using Eq. (17), the equations above can be written as

[TABLE]

with

[TABLE]

Once again, the contributions $\propto\lambda_{3}$ in the numerator are subleading when $M\mu\ll 1$ and can be neglected, finally obtaining Eqs. (23) and (24).

Appendix B Scalar energy and angular momentum fluxes at the horizon

In this section we compute the adiabatic time variation of the mass $M_{S}$ and angular momentum $L_{S}$ of the condensate due to the superradiant instability.

In the Newtonian approximation, the condensate mass is given by Eq. (14), whereas the $z$-component of the angular momentum of the condensate reads

[TABLE]

where the quantity $xT^{0y}-yT^{0x}$ has to be expressed in spherical coordinates.

In order to include an adiabatic time dependence, we make the substitution $A_{lm}\to A_{lm}e^{t/\tau_{lm}}$ in the expression (10) of the mode $\Psi$. Clearly, $\tau_{lm}>0$ in the superradiant phase, whereas $\tau_{lm}<0$ otherwise.

B.1 Model I: $l=1$ and $m=\pm 1$

We start analyzing the case of a scalar cloud described by Eq. (11). Including the time dependence, the expression of the scalar cloud reads

[TABLE]

where we recall that $\omega_{lm}\equiv\omega_{I}$ with a given value of ($l$, $m$). Using Eq. (14), we obtain

[TABLE]

Note that $\omega_{I}\ll\omega_{R}$ in the small-$M\mu$ limit; an important consequence of the latter is the absence of terms proportional to $A_{11}A_{1-1}$ in the above formula. Using Eq. (15) we get

[TABLE]

where $M_{S}(0)$ is the value of $M_{S}$ at $t=0$. Then, from Eq. (25) we obtain

[TABLE]

which can be expressed as a function of $M_{S}(t)$, $\lambda_{1}$ and $\omega_{lm}$ by isolating $M_{S}(0)$ in Eq. (79) and replacing it in the last expression:

[TABLE]

Likewise, using Eq. (76) for the configuration Eq. (77), we obtain the angular momentum of the condensate,

[TABLE]

or, equivalently,

[TABLE]

Note that by using Eqs. (82) and (79) the initial angular momentum of the cloud, $L_{S0}\equiv L_{S}(0)$, is fixed in terms of the the initial mass $M_{S0}\equiv M_{S}(0)$, the initial relative amplitude of the modes, $\lambda_{10}$, the initial BH parameters $M_{0}$ and $\chi_{0}$, and the gravitational coupling $M_{0}\mu$.

In the quasi-adiabatic approximation there should not be any explicit time dependence because all the quantities of interest implicitly vary over the time evolution of the system. To remove the explicit time dependence in Eqs. (81) and (83), we can consider an average over a time $2\pi/\omega_{R}\sim 2\pi/\mu$, which is the characteristic orbital time period of the scalar condensate,

[TABLE]

where $M_{S}$, $\omega_{11}$ are $\omega_{1-1}$ are treated as constants. Performing the integrals we obtain

[TABLE]

Finally, by expanding the above relations in the $M\mu\ll 1$ limit, we obtain Eqs. (27)–(28). Note that the final result would be the same if the averages were performed over several orbital periods.

B.2 Model II: $l=m=2$ and $l=m=1$

We now consider the case of a scalar cloud described by Eq. (12). After including the time dependence, the expression of the scalar condensate reads

[TABLE]

Using Eq. (76) for this configuration, we get

[TABLE]

and the same computation described above for Model I yields

[TABLE]

In this case, the time average gives

[TABLE]

Finally, by expanding the above relations in the $M\mu\ll 1$ limit, we obtain Eqs. (29)–(30).

B.3 Model III: $l=1,2$ with $m=1$

At last let us consider the case of a scalar cloud described by Eq. (13). In this case the scalar condensate reads

[TABLE]

and we obtain

[TABLE]

The last relation is expected because both modes have $m=1$. In this case it is sufficient to average Eq. (91),

[TABLE]

Again, by expanding the above relations in the $M\mu\ll 1$ limit, we finally obtain Eqs. (31)–(32).