Ultralight boson cloud depletion in binary systems

TL;DR

This paper investigates how tidal interactions in binary systems can deplete ultralight boson clouds around black holes, especially for larger values of the gravitational fine structure constant, using perturbation theory.

Contribution

It extends previous models by computing instability rates for a wider range of parameters and includes eccentric orbits, revealing larger depletion effects and broader observational prospects.

Findings

Depletion estimates increase for α ≥ 0.1.

Resonances can occur at all harmonics in eccentric binaries.

Extended the analysis to generic α and eccentric orbits.

Abstract

Ultralight scalars can extract rotational energy from astrophysical black holes through superradiant instabilities, forming macroscopic boson clouds. This process is most efficient when the Compton wavelength of the boson is comparable to the size of the black hole horizon, i.e. when the "gravitational fine structure constant" $\alpha\equiv G \mu M/\hbar c\sim 1$. If the black hole/cloud system is in a binary, tidal perturbations from the companion can produce resonant transitions between the energy levels of the cloud, depleting it by an amount that depends on the nature of the transition and on the parameters of the binary. Previous cloud depletion estimates considered binaries in circular orbit and made the approximation $\alpha\ll 1$. Here we use black hole perturbation theory to compute instability rates and decay widths for generic values of $\alpha$, and we show that this leads to much larger cloud depletion estimates when $\alpha \gtrsim 0.1$. We also study eccentric binary orbits. We show that in this case resonances can occur at all harmonics of the orbital frequency, significantly extending the range of frequencies where cloud depletion may be observable with gravitational wave interferometers.

Equations

$$\left(\square -{\mu }^2\right)\Psi (t,\text{r})=0,$$

$$\Psi =\sum _{\ell ,m}{e}^{im\varphi }{e}^{-i\omega t}{S}_{\ell m}(\theta ){\psi }_{n\ell m}(r),$$

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

$$\Psi (t,\text{r})=\frac{1}{\sqrt{2\mu }}\left[\psi (t,\text{r})e^{-i\mu t}+{\psi }^{\ast }(t,\text{r})e^{i\mu t}\right],$$

$$i\frac{\partial }{\partial t}\psi (t,\text{r})=\left[-\frac{1}{2\mu }{\nabla }^2-\frac{\alpha }{r}\right]\psi (t,\text{r}),$$

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

$$\delta {\omega }_{n\ell m}\simeq \left(-\frac{{\alpha }^4}{8n^4}+\frac{(2\ell -3n+1){\alpha }^4}{n^4(\ell +1/2)}+\frac{2\tilde{a}m{\alpha }^5}{n^3\ell (\ell +1/2)(\ell +1)}\right).$$

$$r_{\mathrm{Bohr}}\simeq \left(\frac{n^2}{{\alpha }^2}\right)M.$$

$${\Gamma }_{n\ell m}=\frac{2r_{+}}{M}{C}_{n\ell m}(a,\alpha )\left(m{\Omega }_H-\omega \right){\alpha }^{4\ell +5},$$

$$C_{n\ell m}(a,\alpha )$$

$$\prod _{j=1}^{\ell }\left[j^2\left(1-{\tilde{a}}^2\right)+{\left(\tilde{a}m-2{\tilde{r}}_{+}\alpha \right)}^2\right].$$

$$|{\Gamma }_d^{(i)}|=\frac{B_{(i)}}{24}\frac{{\alpha }^{10}}{M}{\left(\frac{1-4{\alpha }^2}{1+4{\alpha }^2}\right)}^2\left(\frac{2}{1+4{\alpha }^2}+{\tilde{r}}_{+}\right).$$

$${\Gamma }_{n\ell -1}=-B_n\frac{{\alpha }^{10}\left[16{\alpha }^4{\tilde{r}}_{+}^2+4{\alpha }^2\left({\tilde{r}}_{+}^2+4{\tilde{r}}_{+}+1\right)+1\right]}{{\left(4{\alpha }^2+1\right)}^2},$$

$$\frac{a}{M}=\frac{4mM\omega }{m^2+4(M\omega {)}^2}.$$

$$\Omega =\sqrt{\frac{M+M_{\ast }}{R_{\ast }^3}},$$

$$R_{\ast ,\mathrm{cr}}={\left[\frac{0.49q^{-2/3}}{0.6q^{-2/3}+\ln \left(1+q^{-1/3}\right)}\right]}^{-1}{r}_{\mathrm{Bohr}},$$

$${\tau }_0=\frac{5}{256}\frac{M_{\mathrm{tot}}}{\nu }{\left(\frac{1}{M_{\mathrm{tot}}{\Omega }_0}\right)}^{8/3},$$

$${\Omega }_0>0.042\frac{{\left(1+q\right)}^{1/8}}{(M/M_{\odot })q^{3/8}}{\alpha }^{45/8}\mathrm{Hz}\equiv {\Omega }_c.$$

$${ϵ}_h=\frac{\mu }{12}\tilde{a}{\alpha }^5,{ϵ}_b^{(n)}=-\frac{n^2-4}{16n^2}\mu {\alpha }^2,$$

$$|\psi (t)⟩=c_g(t)|{\psi }_g⟩+c_d^{(h)}(t)|{\psi }_d^{(h)}⟩+c_d^{(b)}(t)|{\psi }_d^{(b)}⟩,$$

$$|c_g(t){|}^2+|c_d^{(h)}(t){|}^2=1.$$

$$|c_d^{(h)}(t){|}^2=\left[1-{\left(\frac{{ϵ}_h\mp \Omega }{{\Delta }_R^{(h)}}\right)}^2\right]{\sin }^2\left[{\int }_{t_0}^{t}d{t}^{\prime }{\Delta }_R^{(h)}(t^{\prime })\right],$$

$${\Delta }_R^{(h)}=\sqrt{{\left(9\eta \right)}^2+{\left({ϵ}_h\mp \Omega \right)}^2},$$

$$\eta ={\alpha }^{-3}\left(\frac{q}{R_{\ast }}\right){\left(\frac{M}{R_{\ast }}\right)}^2.$$

$$|c_d^{(b)}(t){|}^2=\left[1-{\left(\frac{{ϵ}_b^{(3)}\mp \Omega }{{\Delta }_R^{(b)}}\right)}^2\right]{\sin }^2\left[{\int }_{t_0}^{t}d{t}^{\prime }{\Delta }_R^{(b)}(t^{\prime })\right],$$

$${\Delta }_R^{(b)}=\sqrt{{\left(7.6\eta \right)}^2+{\left({ϵ}_b^{(3)}\mp \Omega \right)}^2}.$$

$$v=l+2\sum _{j=1}^{\infty }\frac{1}{j}\left\{J_j(je)+\sum _{k=1}^{\infty }{\beta }^p\left[J_{j-k}(je)+J_{j+k}(je)\right]\right\}\sin jl,$$

$$R_{\ast }=\frac{a_{\mathrm{SM}}(1-e^2)}{1+e\cos {\Phi }_{\ast }},$$

$$a_{\mathrm{SM}}={\left(\frac{M+M_{\ast }}{{\Omega }^2}\right)}^{1/3}.$$

$$\cos {\Phi }_{\ast }=-e+\frac{2}{e}(1-e^2)\sum _{j=1}^{\infty }{J}_j(je)\cos jl.$$

$$|\psi (t)⟩=c_g(t)|{\psi }_g⟩+c_d(t)|{\psi }_d⟩,$$

Full text

I Introduction

The observation of gravitational waves (GWs) by the LIGO and Virgo collaborations Abbott et al. (2018a) marked the beginning of a new era in astrophysics and fundamental physics Barack et al. (2018). These observations have already provided crucial information on the formation of binary compact objects Abbott et al. (2018b), tested general relativity in the strong, highly dynamical regime Abbott et al. (2016); Yunes et al. (2016); Berti et al. (2015, 2018a, 2018b), and led to new measurements of the cosmological expansion of the universe Abbott et al. (2017); Soares-Santos et al. (2019).

Quite remarkably, GW observations also have the potential to transform our understanding of particle physics. One example is the possibility to discover ultralight bosonic particles, such as axions, through GWs Arvanitaki et al. (2015); Yoshino and Kodama (2015a); Arvanitaki et al. (2017); Baryakhtar et al. (2017); Brito et al. (2017a, b); D’Antonio et al. (2018); Isi et al. (2018); Ghosh et al. (2018); Tsukada et al. (2018). Ultralight bosons can efficiently extract rotational energy from spinning black holes (BHs) through superradiant instabilities and form macroscopic condensates when the Compton wavelength of the boson is comparable to the characteristic size of the BH horizon, i.e. when the “gravitational fine structure constant” $\alpha\equiv G\mu M/\hbar c\sim 1$ Brito et al. (2015a). This possibility can shed light on bosons with masses in the range $\sim[10^{-19},10^{-11}]$ eV, which have Compton wavelengths comparable to the size of astrophysical BHs.

The existence and formation of bosonic clouds can be inferred through several observational channels. A first possibility is to look for gaps in the “Regge plane” of astrophysical BHs: superradiant instabilities could lead to a lack of highly spinning BHs in a BH mass range that depends on the boson mass. Measurements of the spin and mass of astrophysical BHs can then be used to infer or constrain the existence of ultralight bosons Arvanitaki et al. (2015); Brito et al. (2015b); Baryakhtar et al. (2017); Cardoso et al. (2018); Stott and Marsh (2018). Even more exciting is the prospect of direct detection: once formed, boson clouds would slowly decay through the emission of long-lived, nearly monochromatic GWs. This radiation is potentially observable, either as a continuous, nearly monochromatic signal from individual sources or as a stochastic background Arvanitaki et al. (2015, 2017); Baryakhtar et al. (2017); Brito et al. (2017a, b); Ghosh et al. (2018); Tsukada et al. (2018).

Here we study how bosonic clouds around astrophysical BHs can affect the dynamics of a binary system, revisiting and extending the recent, comprehensive analysis of Baumann et al. (2019). The cloud can affect the motion of small compact objects in its vicinity Ferreira et al. (2017); Zhang and Yang (2019). The special nature of the axionic cloud would leave characteristic signatures in the gravitational waveforms from extreme mass-ratio inspirals, that are potentially detectable by LISA Hannuksela et al. (2018). In this work we consider the effect of a binary companion on the bosonic cloud itself Arvanitaki et al. (2015). Under certain conditions, the perturbations induced by the companion can lead to resonant transitions between superradiant and non-superradiant modes, which can, in some cases, deplete the cloud Baumann et al. (2019). Ref. Baumann et al. (2019) assumed that $\alpha\ll 1$ and that the binary is on a circular, equatorial orbit that can be treated in the Newtonian limit. We relax the approximation $\alpha\ll 1$ (which is expected to fail precisely when the superradiant instability is strongest) by using numerical calculations in BH perturbation theory to estimate the decay rate of the non-superradiant modes, and we point out a sign error that affects the decay rates of Ref. Baumann et al. (2019). We also extend their analysis to eccentric binaries, showing that multiple resonant depletion episodes can occur for eccentric BH binaries of interest for LISA.

The plan of the paper is as follows. In Section II we review the hydrogenic structure of the energy levels of boson clouds around rotating BHs. In Section III we discuss resonances in circular and eccentric binary systems, and in Section IV we present our estimates for cloud depletion. In Section V we highlight some limitations of our study and point out directions for future work. To improve readability and to keep this paper self-contained, we relegate some necessary technicalities to the Appendices. Appendix A shows (following Dolan (2007)) how to compute instability rates and decay widths for generic values of $\alpha$ using continued fraction methods. Appendix B deals with level mixing induced by tidal perturbations and with the resulting selection rules, summarizing some important results from Baumann et al. (2019). Finally, Appendix C presents an approximate analytical calculation of the occupation numbers of decaying levels valid for small-eccentricity orbits.

II Hydrogenic structure of the boson cloud

Superradiant instabilities can lead to the formation of ultralight boson clouds around rotating (Kerr) BHs. Consider a scalar field $\Psi$ of mass $\mu$, described by the Klein-Gordon equation on a Kerr background:

[TABLE]

where $\Box=g^{ab}\nabla_{a}\nabla_{b}$ is the d’Alembert operator, $g^{ab}$ is the contravariant Kerr metric, and $\nabla_{a}$ denotes a covariant derivative. The angular dependence of the scalar field can be separated with the ansatz

[TABLE]

which leads to ordinary differential equations for the radial eigenfunctions $\psi_{n\ell m}(r)$, where $n$ is an integer labeling the discrete eigenfrequencies $\omega$. Modes with angular frequency $\omega$ and azimuthal number $m$ will be superradiantly amplified if the BH rotates faster than the field’s phase velocity, i.e.

[TABLE]

where $\Omega_{H}=\frac{a}{2Mr_{+}}$ is a function of the Kerr angular momentum parameter $a=J/M$ (where $M$ and $J$ are the BH mass and angular momentum), and $r_{+}\equiv M+\sqrt{M^{2}-a^{2}}$ denotes the Boyer-Lindquist horizon radius (here and below we will use geometrical units, $G=c=1$). Let us also define $\tilde{a}=a/M$ and $\tilde{r}_{+}\equiv r_{+}/M$ for future use. The mass of the scalar field works as a potential barrier that confines the superradiant modes, leading to a continuous extraction of angular momentum from the BH until the inequality (3) is saturated.

By plugging into Eq. (1) the ansatz

[TABLE]

where $\psi(t,\textbf{r})$ is a complex scalar field that varies on timescales much longer than $\mu^{-1}$ and ∗ denotes complex conjugation, and keeping only terms up to first order in $r^{-1}$ and linear in $\alpha$, we obtain

[TABLE]

where $\nabla^{2}$ is the Laplacian operator and $\alpha\equiv M\mu$ is the equivalent of the fine-structure constant for the hydrogen atom. In fact, Eq. (5) is formally equivalent to the Schrödinger equation for the hydrogen atom, and thus the eigenstates $\psi_{n\ell m}(r)$ re hydrogenic eigenfunctions with principal and orbital quantum numbers $n$ and $\ell$, respectively. Let us remark that we follow the conventions of Baumann et al. (2019), which are more convenient to highlight similarities with the spectrum of the hydrogen atom. In particular, our principal quantum number $n$ is the same as $\tilde{n}=n+\ell+1$ in Dolan’s notation Dolan (2007), and the dominant superradiant mode – a nodeless ($n=0$) solution with $\ell=m=1$ in Dolan’s notation – corresponds to $n=2$ in our conventions. The eigenfrequencies of these states are Detweiler (1980)

[TABLE]

where $\delta\omega_{n\ell m}$ denote higher-order corrections, that (up to fifth-order in $\alpha$) are given by Baumann et al. (2019)

[TABLE]

Finally, the characteristic Bohr radius – i.e., the radius at which the radial profile of the scalar field achieves its maximum value – is well approximated by

[TABLE]

The eigenstates are not stationary because of dissipation, therefore the eigenfrequencies have an imaginary part of the form $i\Gamma_{n\ell m}$, where the coefficients $\Gamma_{n\ell m}$ are the instability rates (decay widths) for unstable (stable) modes, respectively Detweiler (1980); Dolan (2007).

For generic values of the constant $\alpha$, the decay width $\Gamma_{n\ell m}$ must be computed numerically, as explained in Appendix A. However, as first shown by Detweiler Detweiler (1980), in the limit $\alpha\ll 1$ these quantities can be computed analytically, with the result:

[TABLE]

where

[TABLE]

One of the main purposes of this paper is to improve over this small-$\alpha$ approximation, which is significantly inaccurate when $\alpha\gtrsim 0.1$. This is shown in Fig. 1, where – for illustration – we compare the decay width of the (stable) $\Gamma_{21-1}$ mode computed using continued fractions (CF, solid black line) and the Detweiler small-$\alpha$ approximation (dashed red line), for a BH spin that saturates the superradiant amplification of the $\ell=m=1$ mode [cf. Eq.(13)].

When comparing our results with those of Ref. (Baumann et al., 2019) we found large deviations for $\alpha\gtrsim 0.05$. This discrepancy is caused only in part by the use of the Detweiler approximation, and we think that it is partly due to the use of an incorrect equation111Eq. (3.40) of (Baumann et al., 2019) reads

$|\Gamma_{d}^{(i)}|=\frac{B_{(i)}}{24}\frac{\alpha^{10}}{M}\left(\frac{1-4\alpha^{2}}{1+4\alpha^{2}}\right)^{2}\left(\frac{2}{1+4\alpha^{2}}+\tilde{r}_{+}\right)\,.$

(11)

where $B_{(i)}$ is a numerical constant dependent on the transition. From Eqs. (9) and (10), the correct result for the dominant transitions should read instead

$\Gamma_{n\ell-1}=-B_{n}\frac{\alpha^{10}\left[16\alpha^{4}\tilde{r}_{+}^{2}+4\alpha^{2}\left(\tilde{r}_{+}^{2}+4\tilde{r}_{+}+1\right)+1\right]}{\left(4\alpha^{2}+1\right)^{2}}\,,$

(12)

where $B_{2}=1/6$ and $B_{3}=128/2187$ for the hyperfine and the Bohr mixing, respectively. Their equation can be recovered by using $m=1$ in Eq. (10), instead of the correct value $m=-1$. for the decay width $\Gamma_{21-1}$. As shown in Fig. 1, where the dotted blue line shows the prediction from Eq. (3.40) of (Baumann et al., 2019), the different predictions are in good agreement for $\alpha\ll 1$, but there are large discrepancies when $\alpha\gtrsim 0.1$.

III Hyperfine resonance and Bohr resonance

III.1 Circular orbits

At the end of the superradiant process, i.e. when the superradiant amplification saturates, the mass and spin of the final BH are related to the frequency of the dominant superradiant mode of the cloud by Brito et al. (2015b)

[TABLE]

Throughout this paper we will assume that the spin of the BH that carries the cloud is given by Eq. (13).

If the BH that carries the cloud is part of a binary, new cloud instabilities arise due to the existence of resonant orbits. The tidal field of the companion will induce perturbations in the potential of Eq. (5), which in turn induce overlaps between different states $\psi_{n\ell m}$, also known as level mixings. Ref. Baumann et al. (2019) studied level mixings for binaries in quasicircular orbits with orbital frequency

[TABLE]

where $M_{*}$ and $R_{*}$ denote the mass of the companion and the orbital separation, respectively. Here we briefly summarize their main results. Defining $\Phi_{*}$ as the azimuthal angle of $M_{*}$ relative to $M$ and setting $\Omega>0$ without loss of generality, configurations with $\Phi_{*}=\Omega t$ ($\Phi_{*}=-\Omega t$) correspond to orbits corotating (counterrotating) with the cloud.

Let us limit attention to binary separations greater than the Bohr radius ($R_{*}>r_{\rm Bohr}$) to guarantee that the gravitational influence between the two bodies can be analyzed using a multipole expansion, and that corrections to the Kerr metric can be treated perturbatively Baumann et al. (2019). If the Bohr radius $r_{\rm Bohr}$ is greater than the Roche radius the gravitational attraction between the two objects will induce mass transfer from the cloud to the companion. This critical orbital separation can be estimated using Eggleton’s fitting formula Eggleton (1983):

[TABLE]

where $q=M_{*}/M$ is the ratio between the mass $M_{*}$ of the companion and the mass $M$ of the BH-cloud system. We restrict our study to the region $R_{*}>\max(r_{\rm Bohr},R_{*,{\rm cr}})$, so we can neglect mass transfer.

Let us define $\tau_{c}\approx 10^{9}(M/10^{5}M_{\odot})(0.1/\alpha)^{15}$ yr to be the boson cloud lifetime, estimated assuming that GW emission is the only dissipative process and that other effect (such as accretion) are negligible Brito et al. (2015b). The merger time for a binary with orbital frequency $\Omega_{0}$, as estimated from the quadrupole formula for GW emission for quasicircular orbits Peters (1964), is

[TABLE]

where $M_{\rm tot}=M+M_{*}$ is the total mass and $\nu=MM_{*}/M_{\rm tot}^{2}$ is the symmetric mass ratio. If $\tau_{c}>\tau_{0}$, the merger occurs before the cloud is radiated away. This relation can be translated into a bound on the initial orbital frequency:

[TABLE]

For circular orbits, we will only consider initial orbital frequencies greater than the critical frequency $\Omega_{c}$.

The selection rules for transitions induced by the tidal potential of the companion are discussed in Appendix B. Considering only the dominant growing mode $\psi_{211}$, two main resonances are of interest during the orbital evolution of the binary:

• (i)

the hyperfine resonance is caused by an overlap between $\psi_{211}$ and the decaying states $\psi_{210}$ and $\psi_{21-1}$;

• (ii)

the Bohr resonance is caused by an overlap between $\psi_{211}$ and the states $\psi_{n10}$ and $\psi_{n1-1}$, for $n\geq 3$.

These resonances occur when the orbital frequency $\Omega$ matches the energy split between two states $\psi_{n\ell m}$ and $\psi_{\tilde{n}\tilde{\ell}\tilde{m}}$, i.e. when $\Omega\sim\Delta\omega/\Delta m=\epsilon$, where $\Delta\omega=\omega_{n\ell m}-\omega_{\tilde{n}\tilde{\ell}\tilde{m}}$ and $\Delta m=m-\tilde{m}$ Baumann et al. (2019); Zhang and Yang (2019). More precisely, the hyperfine and Bohr resonances will occur when $\Omega$ matches the hyperfine splitting $\epsilon_{h}$ or the Bohr splitting $|\epsilon_{b}|$, which – at leading order in $\alpha$, using Eqs. (6) and (7) – are given by

[TABLE]

where $n\geq 3$. The hyperfine resonance will only occur for corotating orbits (because $\epsilon_{h}>0$), whereas the Bohr resonance will only occur for counterrotating orbits (because $\epsilon_{b}^{(n)}<0$). For orbital frequencies outside these resonances the mixing between the modes is perturbatively small and can, in general, be neglected. For the Bohr mixings, we will only consider the $n=3$ resonance: this is usually dominant because it occurs earlier in the inspiral, and because the decay width decreases with $n$ [cf. Eq. (12)].

Let us now define the occupation densities of the modes as $c_{g}(t)$ for the growing mode $\psi_{211}$, $c_{d}^{(h)}(t)$ for the decaying modes of the hyperfine resonance, and $c_{d}^{(b)}(t)$ for the decaying modes of the Bohr resonance. In general we have a three-level system

[TABLE]

where the occupation densities satisfy the normalization condition $|c_{g}(t)|^{2}+|c_{d}^{(h)}(t)|^{2}+|c_{d}^{(b)}(t)|^{2}=1$.

Consider first the hyperfine mixing. For quasicircular, co-rotating equatorial orbits the growing mode $\psi_{211}$ dominantly couples to the decaying mode $\psi_{21-1}$, while the perturbative coupling to the $\psi_{31-1}$ mode can be neglected. Therefore the occupation densities must satisfy the normalization condition

[TABLE]

Solving the perturbed Schrödinger equation for the coupled states (see Appendix B) with the initial conditions $c_{g}(0)=1$ and $c_{d}^{(h)}(0)=0$ yields the following proportionality relation for the occupation density of the decaying mode:

[TABLE]

where we defined the modified Rabi frequency for the hyperfine splitting as

[TABLE]

where upper sign in $\mp$ stands for the co-rotating orbits and the bottom for the counter-rotating ones, and

[TABLE]

The phase of the oscillations in Eq. (21) has been written as an integral over time to take into account the fact that $\Delta_{R}^{(h)}$ changes as the orbits shrinks due to radiation reaction Baumann et al. (2019).

The Bohr resonance is important only for counter-rotating orbits. For equatorial orbits the mode $\psi_{310}$ decouples, and we only have to consider the decaying mode $\psi_{31-1}$ (neglecting all modes with $n>3$). Near $\Omega\simeq|\epsilon_{b}|$, the phase of the hyperfine mixing oscillates rapidly with a period of the order $\eta^{-1}$. In this region $|c_{d}^{(h)}(t)|^{2}\sim(\eta/\epsilon_{b})\ll 1$ Baumann et al. (2019), and the problem reduces again to a two-level system. Solving the perturbed Schrödinger equation for the remaining states with the initial conditions $c_{g}(0)=1$ and $c_{d}^{(b)}(0)=0$ leads to the following occupation density for the decaying state:

[TABLE]

where $\eta$ was given in Eq. (23), and in this case the modified Rabi frequency reads

[TABLE]

The addition of modes with $n>3$ would require, in general, the solution of an infinite-dimensional system. We can still approximate the system as a two-level system, as these resonances occur at different orbital frequencies for each $n$, and in the vicinity of each resonance only one mode dominates. However, as stated above, higher modes are subdominant in the estimation of the depletion of the cloud, and therefore we neglect resonances with $n>3$.

III.2 Eccentric orbits

The generalization to eccentric orbits can be done by promoting the orbital phase $\Phi_{*}$ to $\Phi_{*}=F(\Omega t,e)$, where $\Omega$ now describes the mean orbital frequency of the orbit and $0\leq e<1$ is the orbital eccentricity. At the Newtonian level, the mean orbital frequency $\Omega$ and the orbital eccentricity $e$ are constants of motion. To find $F(\Omega t,e)$ we use the fact that the true anomaly $v\equiv\Phi_{*}-\Phi_{0}$ and the mean anomaly $l\equiv\Omega(t-t_{0})$, where $\Phi_{0}$ and $t_{0}$ are some initial time and initial orbital phase, are related through the following Fourier series (see e.g. Brouwer and Clemence (1961)):

[TABLE]

where $J_{j}(x)$ denotes the Bessel functions of the first kind and $\beta=(1-\sqrt{1-e^{2}})/e$. Without loss of generality we will set $\Phi_{0}=0$ and $t_{0}=0$, so that $v=\Phi_{*}$ and $l\equiv\Omega t$. For eccentric orbits, the binary separation will depend on the orbital phase through the elliptical orbit equation, which at Newtonian order is given by

[TABLE]

where $a_{\rm SM}$ is the semi-major axis, related to the mean orbital frequency via Kepler’s third law

[TABLE]

The cosine of the orbital phase can also be expanded in a Fourier series (see e.g. Yunes et al. (2009)):

[TABLE]

As in the quasi-circular case, we only consider binary separations greater than the critical (Roche) radius at which mass transfer from the cloud to the companion becomes important. Ref. Dosopoulou et al. (2017) found that in the presence of nonzero eccentricity the Roche radius in Eq. (15) increases by a factor $(1-e)^{-1}$. However gravitational radiation reaction tends to circularize the orbit, and for the orbits considered here the eccentricity is small enough that Eq. (15) is still a very good estimate of the Roche radius.

We now have all the necessary ingredients to generalize the calculation of the occupation densities to eccentric orbits. As in the circular case, let us focus on equatorial orbits and on two-state systems, which (as argued above) describe very well the resonances of interest. In the interaction picture, the wavefunction of the cloud is therefore a linear combination

[TABLE]

where $g$ and $d$ denote the growing and the decaying mode, respectively, and again $|c_{g}(t)|^{2}+|c_{d}(t)|^{2}=1$. In the nonrelativistic limit, and for generic equatorial orbits, the evolution of the coefficients $\mathbf{c}\equiv(c_{g},c_{d})^{T}$ is described by the following Schrödinger equation [cf. Eq. (3.21) of Baumann et al. (2019) and Appendix B]:

[TABLE]

where $A=9$ for the hyperfine resonance and $A=-7.6$ for the $n=3$ Bohr resonance described above. Below we restrict to those resonances and define a parameter $\epsilon=\Delta\omega/\Delta m$, which is given by Eq. (18) for the hyperfine and Bohr transitions.

This system can be written as a single second-order differential equation for $c_{d}(t)$ by taking the derivative of Eq. (31), and by eliminating $c_{g}(t)$ and its derivative from the system. After some algebra we find

[TABLE]

where a dot denotes a derivative with respect to $t$. We can further simplify this equation by writing

[TABLE]

where $B(t)=B_{0}-\Phi_{*}(t)+i\log[\eta(t)]/2$ and $B_{0}$ is an arbitrary integration constant. Without loss of generality we require $B(0)=0$ and therefore $B_{0}=\Phi_{*}(0)-i\log[\eta(0)]/2$. The evolution of the system can be schematically written as

[TABLE]

where the explicit functional form of $V(t)$ is given in Appendix C. For circular orbits $V(t)$ simplifies to $V(t)=\Delta_{R}^{2}$, where $\Delta_{R}$ is the Rabi frequency defined by

[TABLE]

and one can easily recover the solutions first derived in Ref. Baumann et al. (2019) and discussed above, after imposing the initial condition $C(0)=0$ and an initial condition for $\dot{C}(0)$ that can be easily found by imposing $c_{g}(0)=1$ in Eq. (31).

For generic eccentricities, solutions must be found numerically for given values of the orbital frequency $\Omega$ and of the eccentricity $e$. Approximate analytical solutions can be found for small eccentricity: in Appendix C, for illustration, we show a solution valid up to first order in $e$.

So far we neglected gravitational radiation reaction, i.e. we assumed that $\Omega$ and $e$ are constant. The time evolution of the mean orbital frequency and of the eccentricity were first derived in Ref Peters (1964) under the adiabatic approximation (which is valid when the radiation reaction time scale is much longer than the orbital time scale). They are given by

[TABLE]

To compute the equivalent of Eqs. (21) and (24) for eccentric orbits, we must solve Eq. (34) numerically while evolving the orbit adiabatically using Eqs. (36) and (37). The calculation of the occupation numbers can be carried out as follows:

• (i)

For a given set of fixed initial conditions $\Omega(t=0)=\Omega(0)$ and $e(t=0)=e(0)$, solve (34) numerically, and average $2|c_{d}|^{2}$ over several orbits, but over time scales much smaller than the radiation reaction time scale;

• (ii)

Evolve the orbit using Eqs. (36) and (37) with the initial conditions $\Omega(t=0)=\Omega(0)$ and $e(t=0)=e(0)$, and produce a grid of values for $e(\Omega)$.

• (iii)

For each value in the grid, repeat step (i) and construct $|c_{d}|^{2}$ as a function of $\Omega$.

For $e(0)=0$, this procedure reproduces Eqs. (21) and (24) without the oscillatory terms. For $e(0)\neq 0$, our results are shown in Fig. 2. The most important conclusion of this calculation is that for eccentric orbits with $e(0)\neq 0$, resonances can occur whenever $k\Omega=2|\epsilon|$ (where $k\geq 1$ is an integer), in contrast with the quasi-circular case for which resonances only occur at $\Omega=|\epsilon|$. The existence of resonances with $k\geq 2$ implies that significant depletion can occur earlier than in the quasi-circular case.

IV Cloud depletion: numerical results

Following Baumann et al. (2019), we introduce the depletion estimator $\mathcal{A}(t,\,t_{0})$, where we can take $t_{0}$ to be the time for which superradiance has saturated:

[TABLE]

This quantity represents the ratio between the integrated time that the system spends in the decaying modes and the decay timescale $|\Gamma_{n\ell m}|^{-1}$, weighted by the occupation density of each state. The mass of the cloud decays proportionally to $\exp\left(-2\mathcal{A}\right)$, where the factor of $2$ arises from the quadratic dependence of the stress-energy tensor on the scalar field.

The integral in Eq. (38) is more easily performed in the frequency domain. We can make use of Eqs. (36) and (37) to write:

[TABLE]

To perform the integral we also need $e(\Omega)$. To this end we use fits for $e(\Omega)$ from Ref. Yunes et al. (2009), which are valid for any initial eccentricity. In our numerical evaluation of the integral (38) we ignore the oscillatory terms in Eqs. (21) and (24), following Ref. Baumann et al. (2019).

Our Figs. 3 and 4 update the results presented in Figs. 7 and 8 of Ref. Baumann et al. (2019). To facilitate comparisons, in Fig. 3 we select the same examples shown in the lower panels of those figures. For quasi-circular orbits, we choose the initial frequency to be $\Omega_{0}=\Omega_{c}$ [cf. Eq. (17)], and we truncate the integral at a final frequency $\Omega$ such that our approximations break down, i.e. $\Omega$ corresponds to an orbital radius $R={\rm max}(r_{\rm Bohr},R_{\rm*,cr})$.

There are three possibilites, depending on the binary’s mass ratio $q$ and on the gravitational fine structure constant $\alpha$ Baumann et al. (2019):

• (i)

The cloud depletes dramatically during the resonance.

• (ii)

The cloud undergoes a long period of perturbative depletion.

• (iii)

The cloud mostly survives during the entire inspiral.

Strong depletion, i.e. cases (i) and (ii), corresponds to regions where ${\cal A}>0.5$, or $\exp(-2{\cal A})<1/e$: these regions are marked in red in Fig. 4. Most of the qualitative features of these plots can be understood in terms of two competing effects: (1) the decay width increases with $\alpha$; (2) for fixed $M$ and $\alpha$, the orbital evolution due to GW emission is faster (and hence the binary transits through resonances more rapidly) when $q$ increases. The thick white lines in the co-rotating (counter-rotating) cases corresponds to setting $\Omega_{c}=\epsilon_{h}$ ($\Omega_{c}=|\epsilon_{b}|$). The binary experiences resonant depletion only when the initial frequency $\Omega_{0}$ is small enough. Depletion due to co-rotating (hyperfine) resonances occurs only in the top-right region above the white line in the left panels of Fig. 4: regions marked in red in that part correspond to resonant depletion, i.e. case (i) above. The cloud (partly) survives in the top-right region of these panels only because, for very large $q$, the time that the binary spends within the resonance is short compared to the decay time of the decaying state. In the bottom-left region below the thick white lines, the cloud does not experience any resonance: regions marked in red correspond to long periods of perturbative depletion, i.e. case (ii) above.

In contrast, depletion due to counter-rotating (Bohr) resonances occurs almost in the entire parameter range, with the exception of a small region in the bottom-right of the right panels of Fig. 4. However the Bohr resonance occurs at much smaller orbital separations than the hyperfine resonance: for this reason the binary moves through resonance much faster, and the amount of depletion is (typically) significantly smaller than in the hyperfine case.

One of the most important results of this work is that the small-$\alpha$ approximation dramatically underestimates the decay width $\Gamma_{21-1}$ and $\Gamma_{31-1}$ (cf. Fig. 1), and therefore it is inadequate to estimate depletion when $\alpha\gtrsim 0.1$. This is evident in both Figs. 3 and 4, with the difference being particularly striking for the $(\alpha=0.4,\,q=10^{4})$ case in the left panel of Fig. 3.

Fig. 5 shows how cloud depletion proceeds for the same eccentric binary systems that were shown in Fig. 2: a co-rotating orbit with $(\alpha=0.2,\,q=10^{5})$ and a counter-rotating orbit with $(\alpha=0.1,\,q=10^{-2})$. For illustration we start the orbit at $\Omega(0)=0.9\epsilon_{h}/2$ ($\Omega(0)=0.9|\epsilon_{b}|/2$) for the co-rotating (counter-rotating) case. The gradual depletion at the different resonances is apparent. The resonances with $k\geq 2$ imply that significant depletion can occur earlier than in the quasi-circular case, especially for large values of $e(0)$.

V Discussion and conclusions

We reanalyzed the issue of bosonic cloud depletion in binary systems, dropping two crucial assumptions made in Ref. Baumann et al. (2019): the small-$\alpha$ approximation and the assumption of circular orbits. Our study has two important observational implications:

the small-$\alpha$ approximation leads to a significant underestimation of cloud depletion in the region $\alpha\gtrsim 0.1$, where superradiant effects are expected to be stronger.

• 2)

the inclusion of eccentricity significantly extends the frequency band in which depletion effects are potentially observable, because resonances can occur at all integer multiples $k\Omega$ of the binary’s orbital frequency.

To quantify this remark, in Fig. 6 we plot resonant frequencies – defined as $f_{\rm res}=\epsilon/(k\pi)$ with $\epsilon=\epsilon_{h}$ or $\epsilon=|\epsilon_{b}|$ for the hyperfine or Bohr resonance, respectively – for representative values of $\alpha$: $\alpha=0.1$ roughly corresponds to the edge of the region where the small-$\alpha$ approximation is valid, while $\alpha=0.4$ is close to the value $\alpha=0.42$ at which the growth rate of the dominant ($\ell=m=1$) superradiant mode is maximized for near-extremal Kerr BHs Dolan (2007). For each value of $\alpha$, the different lines correspond to the first six harmonics of the orbital frequency ($k=1,\,\dots,\,6$ from bottom to top).

For any fixed $\alpha$, the plot shows that the frequency band in which cloud depletion could be observable is much larger if we allow for eccentric orbits. Second-generation ground-based detectors are severely limited by their seismic cutoff $f\sim 10$ Hz, as shown by the dark gray band: there is a very small chance of observing hyperfine resonances with $\alpha\simeq 0.4$, and the prospects are only slightly better for Bohr resonances (which however, as discussed earlier, would produce much smaller cloud depletion). Third-generation detectors have better chances of observing hyperfine transitions because of their lower seismic cutoff $f\sim 1$ Hz. Observational prospects are much brighter for LISA, where interactions with gas and stellar environments are expected to lead to significant orbital eccentricities Berti (2006); Roedig and Sesana (2012). When superradiance is strongest ($\alpha\simeq 0.4$), hyperfine and Bohr transitions would occur precisely in the BH mass range $M\sim 10^{2}$–$10^{7}~{}M_{\odot}$ targeted by LISA Klein et al. (2016).

Our analysis made several simplifying assumptions. We described the orbit at Newtonian order, but going beyond this approximation is certainly necessary for GW observations. Already at first post-Newtonian order, orbital eccentricity adds an important physical effect that was neglected in our work: periastron precession, which could potentially lead to more orbital resonances. Following Ref. Baumann et al. (2019), we also made several assumptions (e.g. we neglected backreaction, we treated the tidal field in a weak-field approximation, and we assumed that there is no mass transfer) that will inevitably break down at small orbital separations. Full numerical simulations will be necessary to fully understand the dynamics of the system in this regime. Other effects that should be accounted for include self-interactions, which may lead to the suppression of superradiant instabilities and/or bosenovas Arvanitaki and Dubovsky (2011); Yoshino and Kodama (2012, 2015b); the coupling of ultralight bosons with other fields, which may lead to electromagnetic counterparts Rosa and Kephart (2018); Ikeda et al. (2019); Boskovic et al. (2019); dynamical friction and accretion, which could sensibly alter the orbital dynamics of the companion Macedo et al. (2013); Barausse et al. (2014); Hui et al. (2017). All of these effects should be taken into account when modeling gravitational waveforms.

Last but not least, from the point of view of GW data analysis it is particularly important to understand the behavior of the energy flux when binaries are affected by orbital resonances of the kind discussed in this paper. Some techniques to address this problem were developed in the context of the excitation of neutron star quasinormal modes by orbiting companions Pons et al. (2002) and in studies of orbital resonances for extreme mass-ratio inspirals Gair et al. (2011). We plan to address these issues in future work.

Acknowledgments

We wish to thank Vitor Cardoso, Bei-Lok Hu, Lam Hui, Ted Jacobson and Paolo Pani for discussions. This work has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 690904, and it used computational resources at the Maryland Advanced Research Computing Center (MARCC). The authors would like to acknowledge networking support by the GWverse COST Action CA16104, “Black holes, gravitational waves and fundamental physics.” E.B. is supported by NSF Grant No. PHY-1841464, NSF Grant No. AST-1841358, NSF-XSEDE Grant No. PHY-090003, and NASA ATP Grant No. 17-ATP17-0225. C.F.B.M. would like to thank Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), from Brazil, the Johns Hopkins University (JHU) for kind hospitality during the early stages of preparation of this work and the American Physical Society which partially funded the visit to through the International Research Travel Award Program. R.B. acknowledges financial support from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 792862. J.L.R. acknowledges financial support from Fundação para a Ciência e Tecnologia (FCT) - Portugal for an FCT-IDPASC Grant No. PD/BD/114072/2015 and Fulbright Comission Portugal. GR acknowledges financial support provided under the European Union’s H2020 ERC, Starting Grant agreement no. DarkGRA–757480 and by the H2020-MSCA-RISE-2015 Grant No. StronGrHEP-690904.

Appendix A Continued-fraction method

The decay and growth rate of the scalar field eigenmodes can be exactly computed by employing a continued-fraction method Dolan (2007). For completeness, here we summarize this method. The Klein-Gordon equation (1) describing a massive scalar field $\Phi$ with mass $\mu$ can be separated using the ansatz (2) into ordinary differential equations for the radial function

[TABLE]

and the angular function

[TABLE]

where $\Delta=r^{2}+a^{2}-2Mr$. When supplemented by appropriate boundary conditions, the radial and angular equations yield an eigenvalue problem for the angular separation constant $\Lambda_{\ell m}$ and the eigenfrequency $\omega$. The angular solutions are spheroidal harmonics with an angular separation constant that can be accurately computed through a series expansion Berti et al. (2006)

[TABLE]

where $c^{2}=a^{2}(\omega^{2}-\mu^{2})$. Explicit expressions for $f_{k}$ can be found in Berti et al. (2006). Exact values for $\Lambda_{\ell m}$ can also be computed through a continued-fraction method Berti et al. (2006); however, for the modes of interest $c\ll 1$, and therefore (42) provides a very accurate value of the angular eigenvalue. Let us then focus on the calculation of the radial eigenfrequency $\omega$.

At the event horizon the radial function goes as

[TABLE]

where $\sigma=2r_{+}(\omega-m\Omega_{\rm H})/(r_{+}-r_{-})$. For ingoing waves at the horizon only the solution with a negative sign in the exponent is allowed. At spatial infinity the radial function behaves as

[TABLE]

where $q=\pm\sqrt{\mu^{2}-\omega^{2}}$. Here we will be interested in the solutions for which $\psi_{\ell m}(r)$ is regular at infinity, and therefore we are interested in the solutions for which $\text{Re}(q)<0$. These solutions describe quasibound states.

We therefore look for solutions of the form

[TABLE]

where $\chi=(\mu^{2}-2\omega^{2})/q$ with $\text{Re}(q)<0$. We note that choosing $\text{Re}(q)>0$ one would instead find the quasinormal modes of the system: modes described by ingoing waves at the horizon and outgoing waves at infinity.

After inserting this ansatz into the radial ordinary differential equation (40) one obtains a three-term recurrence relation for the coefficients $a_{n}$ given by

[TABLE]

where

[TABLE]

and $c_{0}$, $c_{1}$, $c_{2}$, $c_{3}$ and $c_{4}$ are functions of $\omega$ and $\mu$ that can be found in Eqs. (39)–(43) of Dolan (2007).

The ratio of the coefficients $a_{n}$ satisfy an infinite continued fraction

[TABLE]

This can be further simplified after substituting the $n=0$ term in this expression and noting that $a_{1}/a_{0}=-\beta_{0}/\alpha_{0}$. We then get

[TABLE]

The discrete family of complex values of $\omega$ for which this condition is satisfied correspond to the bound state frequencies. In particular for each choice of $\ell$ and $m$ there is an infinite tower of solutions corresponding to the different overtones. In practice, to find the solutions numerically we truncate the series expansion in Eq. (45) at some large $n$ (say, $n=10^{3}$) and we empirically check that adding higher-order terms (e.g. up to $n=10^{4}$) does not affect the eigenfrequencies within the desired accuracy.

Appendix B Level-mixing due to gravitational perturbations

As discussed in Ref. Baumann et al. (2019), the gravitational perturbations induced by a companion sufficiently far away from the BH-cloud system can be translated into a shift in the potential of the Schrödinger equation (5), causing level mixings. At lowest order in $\alpha$, the tidal perturbation can be written as Baumann et al. (2019)

[TABLE]

where the coordinates ${\bf R}_{*}(t)\equiv\{R_{*}(t),\Theta_{*}(t),\Phi_{*}(t)\}$ describe the position of the companion relative to the isolated BH-cloud system, and $\bar{\bf r}\equiv\{\bar{r},\bar{\theta},\bar{\phi}\}$ are comoving Fermi coordinates with origin at the center of mass of the BH-cloud system.

This gravitational perturbation will induce an overlap between different modes $\Psi_{n\ell m}$ given by

[TABLE]

where

[TABLE]

The angular integral vanishes unless the following selection rules are satisfied:

• (i)

$-m_{j}+m_{i}+m_{*}=0$;

• (i)

$|\ell_{j}-\ell_{i}|\leq\ell_{*}\leq\ell_{i}+\ell_{j}$;

• (iii)

$\ell_{i}+\ell_{j}+\ell_{*}=2p$, for $p\in\mathbb{Z}$.

If one considers, for example, a BH-cloud system that in isolation is only composed of the fastest growing mode $|n\ell m\rangle=|211\rangle$, then the dominant mixings induced by the leading-order quadrupolar perturbations $\ell_{*}=2$ are with the modes $|210\rangle$, $|21-1\rangle$, $|31-1\rangle$ and $|310\rangle$. If one further restricts to equatorial orbits $\Theta_{*}=\pi/2$, the mixing with the modes $|210\rangle$ and $|310\rangle$ is forbidden, because $\langle 211|V_{*}|n10\rangle=0$ when $\Theta_{*}=\pi/2$.

Consider, for simplicity, the mixing between a growing mode $\Psi_{g}$ and a decaying mode $\Psi_{d}$ (the addition of more modes is straightforward Baumann et al. (2019)). Using perturbation theory one finds that the expectation value for the Hamiltonian of the field $\Psi$ is given by

[TABLE]

where $E_{g}$ and $E_{d}$ are the energy eigenvalues of growing and decaying modes for the unperturbed BH-cloud system. The Hamiltonian can be separated into diagonal and non-diagonal parts:

[TABLE]

Since $H_{0}$ is diagonal, the eigenstates are the same as for the isolated BH-cloud system, but with shifted energy states due to the nonzero expectation values $\langle\Psi_{g}|V_{*}|\Psi_{g}\rangle$. Calculations are then more easily performed in the interaction picture, where the evolution of the state $|\Psi_{I}\rangle$ is defined by Baumann et al. (2019)

[TABLE]

In this picture, one can define the operator $H_{1,I}$ as

[TABLE]

such that the Schrödinger equation can now be written as

[TABLE]

The advantage of working in the interaction picture is that $H_{1,I}$ will, in general, be a purely nondiagonal matrix [cf. Eq. (31)].

In the equations above we have assumed that $H_{0}$ is time-independent. When the eccentricity is nonzero, the expectation values $\langle\Psi_{g,d}|V_{*}|\Psi_{g,d}\rangle\propto R_{*}(t)^{-1-\ell_{*}}$ vary on a timescale given by the orbital period, and therefore $H_{0}$ is time-dependent. We can still work in the interaction picture by replacing the unitary propagator $e^{iH_{0}t}$ with $e^{i\int_{0}^{t}dt^{\prime}H_{0}(t^{\prime})}$ (see e.g. Brinkmann (2016)). The integral can be simplified by noting that $H_{0}(t)$ varies over the orbital period, which is much smaller than the typical timescale associated with Rabi oscillations, especially at the resonant frequencies [cf. Eq. (35)]. Therefore the integral can be approximated by $\int_{0}^{t}dt^{\prime}H_{0}(t^{\prime})\approx\langle H_{0}\rangle t$, where $\langle H_{0}\rangle$ denotes a time-average over the orbital period. In addition we assume that radiation reaction can be treated adiabatically, because $\dot{a}_{\rm SM}/a_{\rm SM}\ll\Omega$. When computing the Hamiltonian we are also neglecting the slow decay of the decaying modes. This is usually appropriate since the decay widths are, in general, much smaller than the frequency eigenvalues, i.e. $\Gamma_{n\ell m}\ll\omega_{n\ell m}$.

Appendix C Occupation numbers to first-order in the eccentricity

The potential $V(t)$ in Eq. (34) is given by

[TABLE]

Remarkably, for small eccentricities we find that Eq. (34) can be solved exactly by expanding $V(t)$ up to first order in $e$. In particular, at first order in $e$ one finds

[TABLE]

where we used Eqs. (26),  (27) and (29). At zeroth order in $e$, the solution satisfying the initial condition $C(0)=0$ is $C(t)=A_{0}\sin(\Delta_{R}t)+\mathcal{O}(e)$. For solutions valid up to first order in $e$ we therefore attempt to find a solution of the form

[TABLE]

where we require that $A_{0}$ does not depend on $e$. Requiring the initial condition $c_{g}(0)=1$ and using Eq. (31) to get $\dot{c}_{h}(0)$ we find

[TABLE]

Using these initial conditions for $X(t)$ and inserting (64) in Eq. (34) we find the rather complicated expression for $X(t)$

[TABLE]

where

[TABLE]

and we introduced the quantity

[TABLE]

We find that this solution is in good agreement with numerical solutions for $e\lesssim 10^{-2}$ but the approximation breaks down for larger eccentricities, so we do not use it for the results presented in the main text.