Black hole binaries and light fields: Gravitational molecules
BBlack hole binaries and light fields: gravitational molecules
Taishi Ikeda,
1, 2
Laura Bernard, Vitor Cardoso, and Miguel Zilh˜ao Dipartimento di Fisica, “Sapienza” Universit´a di Roma, Piazzale Aldo Moro 5, 00185, Roma, Italy CENTRA, Departamento de F´ısica, Instituto Superior T´ecnico – IST,Universidade de Lisboa – UL, Avenida Rovisco Pais 1, 1049 Lisboa, Portugal LUTH, Observatoire de Paris, PSL Research University, CNRS,Universit´e de Paris, 5 place Jules Janssen, 92195 Meudon, France
We show that light scalars can form quasi-bound states around binaries. In the non-relativisticregime, these states are formally described by the quantum-mechanical Schr¨odinger equation fora one-electron heteronuclear diatomic molecule. We performed extensive numerical simulations ofscalar fields around black hole binaries showing that a scalar structure condensates around the binary– we dub these states “gravitational molecules”. We further show that these are well-described bythe perturbative, non-relativistic description.
I. INTRODUCTION
More than one century after Einstein wrote down thefield equations of General Relativity, black holes (BHs)remain one of its most outstanding and intriguing predic-tions. Among all of its features, the inherent simplicityof BHs is possibly the most remarkable one: just twonumbers (mass and angular momentum) suffice to fullycharacterize these objects in vacuum [1–4].This simplicity and fundamental nature has led toanalogies being drawn between BHs in General Relativ-ity and the hydrogen atom in Quantum Mechanics. Infact, in the context of fundamental fields in curved space-time, BHs do behave as atoms: massive scalar fields canform long-lived states which are, in a certain limit, math-ematically described by the non-relativistic Schr¨odingerequation for the Hydrogen atom [5–8]; such states havebeen dubbed “gravitational atoms” [8].The horizon of non-spinning BHs acts as a dissipa-tive surface, hence these scalar states are in general“quasi-bound”. When the host BH is spinning theseconfigurations may grow via superradiance, extractinga substantial fraction of the BH rotation energy to abosonic “cloud” in the BH exterior. The process slowsthe BH spin down and releases monochromatic gravi-tational waves, giving rise to very particular imprints.These states may form, through a different process, as aconsequence of boson star collisions [9] or collisions be-tween axion stars and BHs [10]. Thus BHs can be usedas efficient particle detectors of ultralight field across awide range of mass [5, 11, 12]. For fine tuned conditions,the states can become truly bound-states, and new BHsolutions become possible [13].Black hole binaries were recently shown to have char-acteristic vibration modes [14]. Together with the abovediscussion on gravitational atoms, one is led to askwhether the program can be taken a step further: Doesit make sense to talk about “gravitational molecules”?We explore this question here, and give a positive answerproviding (under certain mild conditions), an equivalencebetween BH binaries and di-Hydrogen molecules.We will study this issue by looking at the dynamics of massive scalar fields in a BH binary (BHB) spacetime.We thus consider a Klein-Gordon scalar, governed by theequation (cid:3) φ = µ φ , (1)in a nontrivial background describing a binary, in partic-ular a BHB. We always neglect backreaction of the field inthe spacetime geometry, which in all but extreme situa-tions should be a very accurate approximation. Here, themass parameter µ is related to the boson mass m B = (cid:126) µ .We use geometric units G = c = 1 and a ( − + ++)convention for the metric signature throughout. II. NON-RELATIVISTIC SCALARS AROUNDBINARIES
We first explore the problem in a non-relativistic set-ting. This means that the spacetime is taken to be weaklycurved, and could describe for instance two non-compactstars. It also means that the rest-mass of the scalar dom-inants over its kinetic energy.
A. Equivalence with di-Hydrogen molecules
In a non-relativistic setting, the Klein-Gordon equa-tion around a single BH reduces to the Schr¨odinger equa-tion [5, 7, 8]. There is thus a quasi-bound state structurefor scalar fields around a BH which is identical to thespectrum of the Hydrogen atom. Let us now considerthe non-relativistic limit of a scalar field around a BHB.To lowest order in a post-Newtonian expansion, the ge-ometry of a binary (including that of a BHB, if we’re notinterested in near-horizon phenomena) can be written inthe form ds = − (1 + 2Φ N ) dt + (1 − N ) δ ij dx i dx j , (2)where Φ N ( t, x i ) = − M | (cid:126)r − (cid:126)r ( t ) | − M | (cid:126)r − (cid:126)r ( t ) | , (3) a r X i v : . [ g r- q c ] S e p is Newton’s potential. Here, M i ( i = 1 ,
2) are the indi-vidual component masses, and (cid:126)r i ( t ) their position vec-tor. There are higher-order terms which depend on thespecifics of the system, and which become relevant forrelativistic and strongly gravitating systems, but whichdo not affect the physics we are interested in.Using standard non-relativistic limit procedures, wedefine the complex field Ψ( t, (cid:126)r ) as φ = 1 √ µ (cid:0) Ψ e − iµt + Ψ ∗ e iµt (cid:1) . (4)Assuming that the binary components are widely sep-arated and that the angular frequency is so small thattime-dependent terms in the Newtonian potential can beneglected, the Klein-Gordon equation (1) reduces to theSchr¨odinger equation i∂ t Ψ( t, x i ) = (cid:18) − ∇ µ + µ Φ N (cid:19) Ψ( t, x i ) , (5)where we neglect the sub-leading (for weakly gravitating,non-relativistic systems) terms2 µ ∂ t Φ N ∂ t Ψ , i∂ t Φ N Ψ , ∂ t Ψ2 µ , µ Φ N ∇ Ψ , We can also recover Eq. (5) via a Lagrangian ap-proach [15].Equation (5) is written in the lab frame, x µ =( t, r, θ, ϕ ). It will be useful to write it in the binary restframe (co-rotating frame), ¯ x µ = (¯ t, ¯ r, ¯ θ, ¯ ϕ ), which we cando with the usual coordinate transformation ∂ t = ∂ ¯ t − Ω ∂ ¯ ϕ , ∂ ϕ = ∂ ¯ ϕ , (6)where ∂ ¯ ϕ = − ¯ y∂ ¯ x + ¯ x∂ ¯ y spans the orbital plane, and Ω isthe BHB orbital angular velocity. The frames are relatedthrough¯ t = t , ¯ r = r , ¯ θ = θ , ¯ ϕ = ϕ − Ω t . (7)Denoting ¯Ψ(¯ t, ¯ r, ¯ θ, ¯ ϕ ) ≡ Ψ(¯ t, ¯ r, ¯ θ, ¯ ϕ + Ω¯ t ) = Ψ( t, r, θ, ϕ ),and remembering that ∇ = ¯ ∇ , we can then rewriteEq. (5) in the co-rotating frame as i∂ ¯ t ¯Ψ(¯ t, ¯ x i ) = H ¯Ψ(¯ t, ¯ x i ) + i Ω ∂ ¯ ϕ ¯Ψ(¯ t, ¯ x i ) , (8)where H = − µ ¯ ∇ + V , (9)and the (time-independent) potential V is given by V = − µM r − µM r , (10)where r , = (cid:112) (¯ x ∓ a ) + ¯ y + ¯ z is the distance to BH“i” and ¯ x = ± a are the positions of each BH. Equation (8) can be treated perturbatively when Ω issmall. Let us first consider the unperturbed system, i∂ ¯ t ¯Ψ = H ¯Ψ . Since the potential V is not time dependent, we can con-sider the energy eigenstate problem; writing¯Ψ(¯ t, ¯ x i ) = ¯ ψ (¯ x i ) e − i ¯ E ¯ t , (11)we then have ¯ E ¯ ψ = − µ ¯ ∇ ¯ ψ + V ¯ ψ . (12)Introducing prolate spheroidal coordinates¯ x ( ξ, η, χ ) = a η ξ ¯ y ( ξ, η, χ ) = a (cid:112) − η (cid:112) ξ − χ )¯ z ( ξ, η, χ ) = a (cid:112) − η (cid:112) ξ − χ ) , (13)where 2 a is the separation between the two BHs on the¯ x -axis and − ≤ η ≤
1, 1 ≤ ξ , and 0 ≤ χ < π , Eq. (12)becomes separable. Using the ansatz ¯ ψ ( ξ, η, χ ) = e im χ χ √ π R ( ξ ) S ( η ) , (14)with m χ = 0 , ± , ± , . . . , we then find0 = ∂ η (cid:0) (1 − η ) ∂ η S (cid:1) + (cid:32) A − µa ¯ Eη + 2 aµ ∆ αη + m χ η − (cid:33) S , (15a)0 = ∂ ξ (cid:0) ( ξ − ∂ ξ R (cid:1) + (cid:32) − A − m χ ξ − αµaξ + 2 µa ¯ Eξ (cid:33) R , (15b)where we defined α i = M i µ , α = α + α , ∆ α = α − α . Here, A is a separation constant. ¯ E and A are labeled bythree integers m ξ , m η , m χ , characterizing the propertiesof the solutions of the coupled system. We will focus onbound-state solutions, for which ¯ E < Z + Z = µα, Z − Z = µ ∆ α then theequations above are mathematically equivalent to thequantum-mechanical problem (in atomic units), wherethe nuclei have atomic numbers Z i each and are sepa-rated by a fixed distance D ≡ a [16–20]. In particu-lar, this system also describes the ionized di-Hydrogenmolecule [16, 17]. We thus have a formal equivalence be-tween two similar systems, that of a molecule governedby electromagnetism and a simple binary system in a TABLE I. Eigenvalues for equal-mass binaries correspondingto α = 0 .
2, and to two different binary separation 2 a = 10 , E and angular separation A are labeled by threeintegers ( m ξ , m η , m χ ) in a manner analogous to the eigen-values of the Hydrogen atom in Quantum Mechanics. Alldimensionful quantities are normalized by the total mass M . a = 5 a = 30( m ξ , m η , m χ ) A × ¯ E A × ¯ E (0,0,0) − . − . − . − . − . − . − . − . − . − . − . − . . − . . − . . − . . − . . − . . − . Newtonian setting. We will see below that the inclusionof full general-relativistic effects alters this picture onlyslightly.Equations (15) are of spheroidal type. The first is an“angular”-type scalar spheroidal equation [20, 21], and iscoupled to the second (radial) equation through the (un-known) energy ¯ E and separation constant A . We havesolved this system to find the characteristic energies ¯ E ,with two different methods. We used direct integrationof the ordinary differential equations, shooting to the en-ergy; in addition, we used a high accuracy continued frac-tion approach to solve the same problem [20]. Our nu-merical results for selected values of the separation areshown in Table I. B. The single black hole limit
At zero separation, we are effectively dealing with onesingle BH, for which the energy levels are known to highprecision. In spherical polar coordinates (¯ r, ¯ θ, ¯ ϕ ) theeigenfunction is¯ ψ = e − i ¯ E ¯ t γ (cid:96) Y m(cid:96) (¯ θ, ¯ ϕ ) e − γ/ L (cid:96) +1 n ( γ ) , (16) γ ≡ M µ ¯ r(cid:96) + n + 1 , (17)with Y m(cid:96) the scalar spherical harmonics and L (cid:96) +1 n a gen-eralized Laguerre polynomial (it’s a polynomial of order n in its argument). Note that M = M + M is the BHBmass. With these definitions and conventions, the energyeigenvalue is ¯ E = − µα (cid:96) + n + 1) , (18)up to O ( α ). Higher order expansions in α can be cal-culated using well-known techniques and are shown inRef. [8], with a different state label. Notice that we fol-low the state labeling of Ref. [5]. The scalar profile ofmode ( n, (cid:96), m ) decays spatially as r (cid:96) + n e − γ/ and the an-gular profile is dictated by the corresponding spherical harmonic. The spatial extend of the scalar configura-tion is defined by the exponential decay, and is of order S ∼ / ( M µ ).Note that the prolate coordinates (13) are also givenby ξ = r + r a , η = r − r a . (19)Defining spherical coordinates such that¯ x = ¯ r cos ¯ θ , ¯ y = ¯ r sin ¯ θ sin ¯ ϕ , ¯ z = ¯ r sin ¯ θ cos ¯ ϕ , one finds, when a → ξ → ¯ ra , η → cos ¯ θ , χ = ¯ ϕ . The relation between quantum numbers is then, in thislimit [22] m ξ = n , m η = (cid:96) − | m | , m χ = m . (20)For small separations aµ , our results – shown in Table I– for the fundamental ( n = 0) mode are compatible with(up to terms of order O ( a µ ))¯ E = − µα (cid:96) + 1) × (cid:32) a µ (cid:0) (cid:96) ( (cid:96) + 1) − m χ (cid:1) (cid:0) α − ∆ α (cid:1) (cid:96) ( (cid:96) + 1) (2 (cid:96) − (cid:96) + 1)(2 (cid:96) + 3) (cid:33) . (21)This analytic expansion was derived by Bethe, and insubsequent analytical work for the di-Hydrogen moleculeand generalizations thereof [18, 19, 22–25]. Our own nu-merical results in the limit of small binary separation arein perfect agreement with such expression.We note that for our numerical simulations shown inSec. IV, we will use a frame rotated by 90 degrees, whichtherefore superposes different m states. In addition,and more importantly, at finite separation a , we can nolonger use Eq. (20): a given state in the ( m ξ , m η , m χ )basis involves, generically, a superposition of all modesin a spherical-harmonic decomposition, such as that ofSec. IV. C. Relation with classical closed orbits
Particle analogies of wave equations often play an im-portant role in several physical phenomena. In partic-ular, a lot of BH physics described by wave equationshave particle descriptions which helps us to understandthe system from a different point of view. For example,BH quasi-normal modes (QNMs) can be related to closednull orbits around the BH based on the WKB approxi-mation, where the QNM frequencies can be estimatedfrom the orbital periods [26]. Conversely, the existenceof closed null orbits around relativistic objects may hintat the possibility of existence of QNMs [14].We will then here study the particle analog of the sys-tem we have been considering of a massive scalar fieldaround a BHB. As we will now see, this particle de-scription can be derived from the WKB approximationof Schr¨odinger’s equation and can be applied to statesaround a general separable spacetime.In Quantum Mechanics, it is well known that the firstorder of the WKB approximation corresponds to theHamilton-Jacobi equation for a classical particle, and theenergy spectrum can be related to the particle motionthrough the Bohr-Sommerfeld condition [27] (cid:73) p k dq k = 2 πn k , (22)where q k and p k are the canonical coordinates and n k is an integer. The integral is performed over a closedclassical orbit which is described by the corresponding Hamilton-Jacobi equation. Since this is derived from theWKB approximation, we can use this relation in the large n k limit to estimate the frequency of the bound states ofSec. II A.Furthermore, from the classically allowed regions of theHamilton-Jacobi equation, we can have an idea about thespatial profile of the wave function. Since, as shown inthe previous subsection, our system can be described bySchr¨odinger’s equation, we can apply these arguments toget the relation between (truly) bound states and closedorbits of a massive particle. Let us then consider the classical motion under New-ton’s potential (10), which corresponds to the motion of amassive particle (of mass µ ) around a BHB, in the binaryrest frame. The Lagrangian for the particle is L (¯ x i (¯ t ) , ∂ t ¯ x i (¯ t )) = µ (cid:0) ∂ ¯ t ¯ x (¯ t ) + ∂ ¯ t ¯ y (¯ t ) + ∂ ¯ t ¯ z (¯ t ) (cid:1) − V , where (¯ x (¯ t ) , ¯ y (¯ t ) , ¯ z (¯ t )) are the particle’s coordinates. Us-ing the coordinates ( ξ, η, χ ), the Hamilton-Jacobi equa-tion reads2 µa ∂ ¯ t S + (cid:0) ξ − (cid:1) ( ∂ ξ S ) + (cid:0) − η (cid:1) ( ∂ η S ) − µ a [( M + M ) ξ + ( M − M ) η ]( ξ − η )( ξ + η ) + ( ∂ χ S ) (1 − η ) ( ξ −
1) = 0 , (23)where S (¯ t, ξ, η, χ ) is Hamilton’s principal function. Inthese coordinates the Hamilton-Jacobi equation is sepa-rable; we thus write S (¯ t, ξ, η, χ ) = S ξ ( ξ ) + S η ( η ) + m χ χ − ¯ E ¯ t , and after substituting into (23) we obtain( S (cid:48) η ) = 2 a µ − η (cid:32) C − ¯ Eη + ∆ αa η − m χ a µ − η (cid:33) , ( S (cid:48) ξ ) = 2 a µξ − (cid:32) − C + ¯ Eξ + αa ξ − m χ a µ ξ − (cid:33) , (24)where C is a separation constant. Comparing withEq. (15) we see that 2 a µC corresponds to A . For clas-sical motion, the right hand side of Eqs. (24) must bepositive, and the parameter space where this happenscorresponds to the classically allowed region. The topol-ogy of this region can change depending on the parame-ters; for simplicity, we will focus on equal mass binaries(∆ α = 0) and on bound states around the binary ( ¯ E < m χ = 0. Strictly speaking, the Bohr-Sommerfeld condition is (cid:72) p k dq k =2 π (cid:0) n k + (cid:1) ; however, since we are only interested in its large n k limit, we have neglected the term in Eq. (22). These are actual bound states because we are using a Newtonianapproximation, and horizons are not present.
Considering first the C > η is unconstrained ( − ≤ η ≤ ξ is determined from | ¯ E | ξ − αa ξ + C < . (25)When | ¯ E | − αa + C <
0, the allowed region for ξ is ξ ∈ [1 , ξ + ], and when | ¯ E | − αa + C >
0, the allowed regionis ξ ∈ [ ξ − , ξ + ]. Here, ξ ± are solutions of the equation | ¯ E | ξ − αa ξ + C = 0, ξ ± = 1 | ¯ E | (cid:32) α a ± (cid:114) α a − | ¯ E | C (cid:33) . (26)Therefore, when C >
0, the particle motion is an orbitaround the binary (see first and second panel in Fig. 1).We now focus on the C < η is η ∈ (cid:104) − , − (cid:113) | C || ¯ E | (cid:105) ∪ (cid:104)(cid:113) | C || ¯ E | , (cid:105) for | C | < | ¯ E | .Existence of solutions to Eq. (25) then implies that | ¯ E | − αa − | C | <
0, and the corresponding allowed region for ξ is then [1 , ξ + ] (see third panel in Fig. 1). In this case,the motion is an orbit around each individual BH.Since the WKB approximation is valid only in thehigh-frequency limit, we cannot easily apply this classicalpicture to the states from the previous subsection. Wecan, however, see that there are two distinct profiles forthese states: one corresponding to configurations boundto each individual BH, and another to a configurationbound to the binary as a whole. FIG. 1. Typical shape of the classically allowed region for a massive particle with m χ = 0 around an equal-mass BHB. Blackcircles denote the BHs. The first panel shows the allowed region for C > | ¯ E | − α/a + C <
0, and the second panelshows the allowed region for C > | ¯ E | − α/a + C >
0. These describe orbits around the binary. The third panel standsfor C < | ¯ E | − α/a + C <
0, and describes orbits around each individual BH.
In order to compare the particle picture to the statesof the previous subsection, let us discuss the spectrum of these states using the Bohr-Sommerfeld conditionEq. (22). From this condition we immediately see that m χ must be the integer n χ and we further get (cid:90) η + η − (cid:115) a µ − η (cid:18) C − ¯ Eη + ∆ αa η − m χ a µ − η (cid:19) dη = 2 πn η , (27a) (cid:90) ξ + ξ − (cid:115) a µξ − (cid:18) − C + ¯ Eξ + αa ξ − m χ a µ ξ − (cid:19) dξ = 2 πn ξ , (27b)where n ξ and n η are integers, and the integration limits η ± and ξ ± are the roots of the expressions inside eachsquare root. Since Eqs. (27) cannot be integrated ana-lytically let us consider the Hydrogen atom limit ( a → (cid:96) = 2 a µC and r = aξ in the equal mass case(∆ α = 0). We obtain the spectrum (cid:96) − | m χ | = n η , − α (cid:114) µ | ¯ E | − (cid:96) = n ξ . These expressions are in good agreement with the spec-trum of the Hydrogen atom in the large n η and n ξ limit,confirming the validity of this particle picture. D. Corrections induced by orbital motion
Having solved Eq. (8) to zeroth order in Ω, let us nowconsider, perturbatively, the effect of rotation. The firstorder correction in Ω to the unperturbed energy levelsthat we have just computed are given by the expectationvalue of the rotation operator L ¯ z = − i∂ ¯ ϕ for the system in the unperturbed state. In the coordinates ( ξ, η, χ ) thisoperator takes the form L ¯ z = i ¯ ya ( ξ − η ) ( ξ∂ η + η∂ ξ ) + i ¯ z ξηa ( ξ − − η ) ∂ χ . (28)Since ¯ z ∼ cos χ and ¯ y ∼ sin χ , we can easily see that theexpectation value of this operator in an eigenstate (14)is zero, (cid:104) L ¯ z (cid:105) ¯ ψ = 0. This shows that rotation effects onlymanifest themselves at order Ω and will therefore beneglected. III. SETUP FOR TIME EVOLUTIONS
In the previous Section we discussed, in a perturbativesetup, the non-relativistic limit of a massive scalar fieldaround a BHB and found molecule-like states which werelabeled therein with three parameters ( m ξ , m η , m χ ). Wewill now numerically solve the Klein-Gordon equationaround a BHB to construct these (quasi-)bound statesthrough time evolutions. A. Numerical implementation
For our numerical implementation we ignore the back-reaction of the massive scalar field on the BHB spacetimeand follow the approach described in [14]. In this ap-proach one builds an approximate BHB background met-ric using the construction outlined in detail in Mundimet al. [28], to where we refer the interested reader (seealso Ref. [29] for the equivalent construction used in thecontext of generating BHB initial data). It is importantto mention that we turn off radiation effects and thereforealways consider binaries with constant separation.To solve Eq. (1) in this spacetime we write the equationin a first order form by introducing K φ ≡ − N ( ∂ t − L β ) φ , (29)where N and β i are the lapse function and shift vectorrespectively. The resulting system is numerically evolvedusing the Einstein Toolkit infrastructure [30–32] withCarpet [33, 34] for mesh-refinement capabilities and themultipatch infrastructure Llama [35]. The scalar fieldequations are evolved in time by adapting the ScalarE-volve code available in [36], which was first used and de-scribed in [37]. In our evolutions we further excise the BHinterior with the procedure outlined in [38]. The overallinfrastructure and evolution is essentially the same as theone used and tested in [14]. B. Initial data
We will evolve two different types of scalar field initialdata, which will be referred to as “non-spinning” and“spinning” initial data. The first of these consists of amomentarily static Gaussian profile given by φ = Ae − r / (2 σ ) , K φ = 0 , (30)where A , σ denote the amplitude, width of the Gaussianpulse, respectively.Secondly, we will evolve configurations which, un-like the previous construction, have angular momentum.These are the “spinning” initial data, for which φ = R ( r ) A ( t, θ, ϕ ) , (31)with R ( r ) = rσ e − r σ , (32) A ( t, θ, ϕ ) = A , (cid:114) π sin θ cos( ϕ + ωt )+ A , − (cid:114) π sin θ cos( − ϕ + ωt ) , (33)and the initial configuration for the K φ field can be triv-ially obtained from Eq. (29). Here, A l,m is the amplitude of each ( l, m ) mode, which can be freely specified, and σ is a typical width of the clump. ( r, θ, ϕ ) are the stan-dard spherical coordinates around the center of mass ofthe system. The time dependence mϕ + ωt ( m = ±
1) ineach term introduces angular momentum on the z -axis.Initial data with negative m is co-rotating with the BHB,while initial data with positive m is anti-rotating. C. Frequency extraction
To analyze our results, we decompose the evolved fieldat fixed radial distance into spherical harmonics Y ml ( θ, ϕ )as follows φ ( t, r = r ext , θ, ϕ ) = ∞ (cid:88) l =0 l (cid:88) m = − l Y ml ( θ, ϕ ) φ lm ( t ) . (34)Note that, since φ is real, φ lm ( t ) = ( − ) m φ ∗ l, − m ( t ) where ∗ denotes complex conjugation. We will further Fouriertransform φ lm ( t ) to check its frequency spectra and com-pare with the results obtained in Sec. II A. In order to doso, however, we must note that the frequencies computedin Sec. II A were computed in the rest frame of the bi-nary, whereas here we will be extracting the data in thelab frame. To compare the data we then need to changeframes once again, which can be done as follows.The Fourier spectra in the lab frame is computedthrough ( φ = φ ( t, r ext , θ, ϕ )) F [ φ lm ]( ω ) = (cid:90) dt e iωt φ lm ( t )= (cid:90) dt e iωt sin θ dθdϕ Y ml ∗ ( θ, ϕ ) φ . (35)Given that Y ml ( θ, ϕ ) = N e imϕ P ml (cos θ ), where P ml arethe Legendre polynomials and N is the normalizationconstant, using (6) we can write F [ φ lm ]( ω ) = (cid:90) d ¯ t sin ¯ θ d ¯ θ d ¯ ϕ e i ¯ t ( ω − m Ω) × N e − im ¯ ϕ P ml (cos ¯ θ ) φ (¯ t, r ext , ¯ θ, ¯ ϕ + Ω¯ t )= (cid:90) d ¯ t sin ¯ θd ¯ θ d ¯ ϕ e i ¯ t ( ω − m Ω) × Y ml ∗ (¯ θ, ¯ ϕ ) ¯ φ (¯ t, r ext , ¯ θ, ¯ ϕ )= F [ ¯ φ lm ]( ω − m Ω) , (36)relating the frequencies computed in the lab frame to theones computed in the co-moving frame.Since in the following Section we will be focusing onthe real part of the multipolar components let us alsowrite F (cid:60) = (cid:90) dt e iωt (cid:60) ( φ lm ( t ))= 12 (cid:90) dt e iωt ( φ lm ( t ) + φ ∗ lm ( t ))= 12 ( F [ φ lm ]( ω ) + ( − ) m F [ φ l, − m ]( ω ))= 12 (cid:0) F [ ¯ φ lm ]( ω − m Ω)+ ( − ) m F [ ¯ φ l, − m ]( ω + m Ω) (cid:1) , (37)where F (cid:60) ≡ F [ (cid:60) ( φ lm )]( ω ) and Eq. (36) was used in thelast step. To connect with our upcoming results we fur-ther need to take the real part of Eq. (37), (cid:60) ( F (cid:60) ) = 14 (cid:0) F [ ¯ φ lm ]( ω − m Ω) + ( − ) m F [ ¯ φ l, − m ]( ω + m Ω)+ ( − ) m F [ ¯ φ l, − m ]( − ω + m Ω)+ F [ ¯ φ lm ]( − ω − m Ω) (cid:1) where we have used F [ φ ∗ lm ]( ω ) = ( − ) m F [ φ l, − m ]( − ω ) . In conclusion, in the lab frame we expect to see a super-position of ω ± m Ω frequencies for each m mode of theco-rotating frame. IV. GRAVITATIONAL MOLECULES
We are now in a position to discuss the evolution ofscalar fields on a background describing realistic BHBs.We evolved a battery of different configurations with theinitial data of Eqs. (30) and (31) on background space-times described by equal-mass binaries, M = M = M/
2, for different values of binary separation D . In thisSection we will report on the results obtained with a sub-set of these runs, summarized in Tables II and III. Theseresults, as described below, are in good agreement withthe non-relativistic analysis of Sec. II, and provide strongevidence to the existence and formation of gravitationalmolecules with scalar fields. TABLE II. List of simulations analyzed for the momentarilystatic Gaussian initial data of Eq. (30). Note that we haverun a much larger set of simulations – listed here are onlythose that are analyzed later on.Run
D/M µM σ/M nonspin1
60 0.5 12 nonspin2
10 0.2 25 nonspin3
60 0.2 25
It is useful to keep in mind the different length scalesinvolved in this problem. As we saw in Sec. II B, the scale S of a state around an isolated BH of mass M i is of order TABLE III. List of simulation analyzed for the spinning initialdata of Eq. (31).Run
D/M µM σ/M A , A , − ωM spin1
10 0.2 25 0 1 0.2 spin2
60 0.2 25 0 1 0.2 spin3
60 0.2 25 1 0 0.2 S i ∼ / ( M i µ ) in the small M i µ limit. In a binary ofcomponent masses M , M , we thus have scales S i , anda global scale S BHB ∼ µ − /M where the total BHB massis M = M + M . If S i is much smaller than D = 2 a , thequasi-bound state can be formed around each BH, andfeels a tidal force from the companion object. On theother hand, if S i is much larger than D , the companionBH strongly disturbs such state, destroying it. However,as discussed in Sec. II, we can expect that a quasi-boundstate forms around the BHB. Furthermore, if this statearound the binary is stable, it should be formed startingfrom generic initial conditions.Note also that timescales are important. A light-crossing timescale is of order D/M , whereas an orbitaltimescale is of order 2 π (cid:112) D /M ∼ M, M for bi-naries separated by D = 10 M or 60 M , respectively. A. Quasi-bound states around individual BHs
Based on the scales above, we expect that for largeenough couplings
M µ the scale S i obeys S i < D , in whichcase the scalar cloud localizes around each BH but notaround the binary. In other words, we expect that BHssufficiently far apart can support clouds as if in isola-tion. To test this, we evolved momentarily static Gaus-sian initial data, Eq. (30), corresponding to configuration nonspin1 in Table II. For these parameters, the typicalsize of the cloud around each BH is S = S ∼ M ,which is smaller than the BH separation. We thereforeexpect that quasi-bound states form around each BH andour results confirm this, as shown in Fig. 2. As can beseen, localized structures are apparent with the dominantmode being an l = 1 , m = 0 state around each individualBH. B. Global quasi-bound states from evolution ofstatic initial data
We will now discuss molecular-like structures, i.e.,scalar clouds around BHBs. We focus on BH separations D = 10 M, M , and we fix the mass of the scalar field to µM = 0 .
2. This corresponds to configurations nonspin2 and nonspin3 , respectively. The scales are such that nowthe size of the quasi-bound states, if present, would en-compass the binary when D = 10 M . We will see thateven for D = 60 M such global states exist. FIG. 2. Profile of a massive scalar field around an equal-mass BHB of total mass M . Color codes are unimportant to interpretthe results (and will be omitted throughout this work), since the Klein-Gordon equation is linear in the scalar field, and there isno backreaction on the metric. The field has a dipolar, l = 1 , m = 0 profile around each BH, which was obtained from evolvingconfiguration nonspin1 in Table II. Left panel shows contours of constant scalar field. Purple, green, and red lines representlines of constant scalar field which are 0 . − .
05, and − .
16, respectively. Right panel shows a contour plot of (scalar field)energy density. The snapshots were taken at t = 1600 M . The evolution of these configurations is show in Fig. 3.Perhaps the most evident aspect of these simulations isthat there is a persistent structure, a cloud or quasi-bound state of scalar field around the binary for rela-tively long time. In the context of Fig. 3, the evolutiontimescale is large enough that the binary performed overten periods. This timescale is orders of magnitude largerthan the free-fall time, and yet the scalar structure per-sists.The second noteworthy aspect is that the state keeps,roughly, the symmetry of the initial conditions. This isclearly seen in the energy density plots, Fig. 3. How-ever, it is clear also that fine structure arises, clearlyseen for larger binary separations, excited by the pres-ence of the binary, an asymmetric perturber. In particu-lar, we observe the excitation of the quadrupole modeby the BHB. This is clearly shown in Fig. 4, wherewe show the monopole l = m = 0 and quadrupole l = m = 2 component of the field at selected “extrac-tion” r ex /M = 40 , , × M for D = 10 M , and 10 M for D = 60 M .In accordance with our analysis in Sec. II, the scalarfield is oscillating with a frequency ∼ µ . To quan-tify the agreement with the non-relativistic analysis, we TABLE IV. Spectrum content of waveforms and comparisonagainst non-relativistic results. The third column of this tableshows the location of the dominant peak of the waveform,from time evolutions, in Fourier space (for l = m = 1 , E should always be present in each spherical harmonic ( l, m )basis used for the simulations. Other components are alsopresent, but we find those to be subdominant. The agreementbetween both is very good and lends strong support to theinterpretation that these are “molecular” gravitational quasi-bound states.Run ( l, m ) Mω M ( µ + E ± m Ω) nonspin3 (0 ,
0) 0.1976 0.1973 spin2 (1 ,
1) 0.1992 0.19940.1948 0.1951 nonspin3 (2 ,
2) 0.2012 0.20160.1930 0.1930 computed the Fourier spectrum of the monopole andquadrupole modes at different radii, and averaged theresult. We find clear peaks at different frequencies, thedominant ones are shown in Table IV for D = 60 M and compared against the non-relativistic prediction of FIG. 3. Snapshots depicting the evolution of a scalar field around an equal-mass BHB. The first two rows correspond toevolutions of configuration nonspin2 (cf. Table II), whereas the last two depict evolutions of configuration nonspin3 . First andthird rows: scalar field. Second and fourth row: energy density. At late times, the scalar and energy density profile rotatecounter-clockwise at a frequency equal to the binary Keplerian orbital frequency.FIG. 4. Monopolar l = m = 0 and dipolar l = m = 2 components of the scalar, at selected extraction radius r ex are shown forconfiguration nonspin2 (left panel) and configuration nonspin3 (right panel) of Table II. The signal modulation is not due tobeating of higher overtones, but to the binary orbital motion. The modulation frequency is 2 m Ω to a good approximation. The agreement is of order 0.1% or better, pro-viding strong evidence that bound, molecular-like gravi-tational states do form around BHBs. Notice also that wefind two (or more) peaks for l = m modes. For D = 60 M we can read from the Table that their separation in fre-quency ∆ ωm ≡ M ( ω m − ω − m ) is ∆ ωm =2 = 0 . ωm = 2 mM Ω ( M Ω ∼ . x -axis is depicted in thethree panels of this figure. The late time profile aroundthe binary is well described by a density ∼ e − r/ (25 M ) .This profile is in accord with the non-relativistic, singleBH expression (16) for the quasi-bound state. Our re-sults indicate this is also a good expression even for themoderately large separations that we studied. Figure 5shows that the exponential falloff gives rise to a slowlydecaying but small tail of energy density for distances r (cid:38) M . This looks to be a slow, radiative componentpart of the signal.In fact, this gravitational system behaves like a rotat-ing molecule in Quantum Mechanics. The snapshots ofthe scalar field and energy density shown in Fig. 3 showthat the scalar field is not only localized around the BHB,but that it is dragged by it, rotating counter-clockwise,along the orbital motion of the binary. The field showsmodulations at low-frequency, in particular at 2 m Ω, ascan be seen in the figures. Such modulation is the ex-pected for an equal-mass binary. Notice that the signalis almost equally modulated in amplitude at different ex-traction radii, thus the low-frequency envelope is not theresult of beatings caused by overtone excitation. In otherwords the scalar is indeed gravitationally dragged by thebinary. The period of such pattern is the same as theorbital period of the binary.
C. Co-rotating dipolar global states
The previous evolutions referred to spherically sym-metric, momentarily static initial data. We now considertime-asymmetric initial data. Let us focus on dipolar ini-tial data co-rotating with the binary, i.e., configurations spin1 and spin2 of Table III. The evolution of the spa-tial profile of the field and its energy density is shown inFigs. 6-7. A dipolar pattern stands out from these plots.Figures 6-7 show that for compact binaries (when theBH separation is smaller than the typical size of thecloud), the energy density of the state acquires a torus-like shape, centered on the binary, and supported by its One needs to pay attention when comparing against the resultsof Sec. II, since a rotation of axes is involved and, as we men-tioned, at finite separations a spherical-harmonic mode maps intoa superposition of the ( m ξ , m η , m χ ) states used in Table I. angular momentum. This is similar to the topology ofquasi-bound states around single BHs. On the otherhand, for large BH separations as in Fig. 7, the profile isno longer connected; the torus “breaks up”, and leavestwo over-density clumps of scalar field co-rotating withthe BHB.The presence of the binary excites other modes withsimilar symmetries as that in the initial data. In par-ticular, we see a strong octupole l = m = 3 mode,shown in Fig. 8 for D = 60 M . Notice that the octupo-lar mode grows from a negligible value to roughly 10%in amplitude of the dipolar component. Notice also thelarge timescales involved: the amplitude of these compo-nents is approximately constant up to timescales of order ∼ M or larger. The modulation in the signal is, aswe explained before, due to the motion of the binary,and has a frequency 2 m Ω as expected for an equal-massbinary.These results indicate that the evolution drove the sys-tem to a quasi-bound state, a relativistic analogue ofthe molecular solutions discussed in Sec. II A for thenon-relativistic system. Together with the previous re-sults, and as we will insist below, these features indi-cate that the formation of quasi-bound states is a ro-bust result for general initial conditions. The analysisof the Fourier-decomposed signal shows a frequency con-tent which is peaked at
M ω = 0 . , . . , . ωm =1 = 0 . re − r/ M ) . Note that the non-relativisticprediction Eq. (16) for the scalar profile is of the form ∼ e − Mµ r/ ( (cid:96) +1) = e − r/ (50 M ) for the fundamental modewhen M µ = 0 .
2. This is in perfect agreement with ourtime-evolution results around a BHB.
D. Counter-rotating dipolar global states
Similar results hold for the evolution of counter-rotating initial data (with respect to the binary, thusdata rotating clockwise). We have evolved configuration spin3 of Table III and the corresponding profiles of thefield and energy density are shown in Fig. 10.Since the initial profile has angular momentum whichis opposite to that of the binary, the scalar field rotatesin a direction opposite to the BHB. However, the energydensity of the scalar cloud does rotate in the same direc-tion as that of the binary. The asymptotic late-time spa-tial distribution of the energy density is well describedby a radial dependence ( re − r/ M ) , consistent with anon-relativistic analysis of quasi-bound states.1 FIG. 5. Energy density of the scalar field on the x -axis for configuration nonspin2 (left panel) and nonspin3 (right panel).Dashed lines are best fits to the numerical results, and agree well with the analytical, non-relativistic predictions (see text).FIG. 6. Snapshots depicting the evolution of a scalar aroundan equal-mass BHB for configuration spin1 of Table III.The initial conditions are therefore those of a dipolar l =1 , m = − As we discuss in Appendix A, selection rules imply thatother modes must be excited, in particular the l = m = 3mode. Figure 11 shows precisely this. It is also clear thata very long-lived quasi-bound state forms (the lifetime isso large, in fact, that we were unable to estimate it). FIG. 7. Same as Fig. 6, but for configuration spin2 of Ta-ble III.
E. General initial data
We investigated other types of initial conditions, suchas narrower pulses, with smaller width σ . This changesthe amount of field that is dissipated in the early stages,but a quasi-bound state always ends up forming, with aphenomenology similar to what we described.We also studied data with higher initial frequency con-tent: since quasi-bound states have ω ∼ µ , it is conceiv-able that high frequency initial data just dissipates away.Our results show that this does not happen. Instead,again, the system evolves towards a frequency content ω ∼ µ and a spatial distribution described well by anon-relativistic quasi-bound state. It seems that thesequasi-bound states are an attractor in the phase-space.We should also add that quasi-bound states decay ex-2 FIG. 8. Dipolar ( l = m = 1) and octupolar ( l = m = 3)components of the scalar field, measured at various r ex , forconfiguration spin2 of Table III. Notice that a dipolar modeis already present initially. This is the same initial data asthat of Fig. 7. The signal modulation is not due to beatingof higher overtones, but by the binary orbital motion. Themodulation frequency is 2 m Ω to a good approximation.FIG. 9. Energy density of a quasi-bound state around a BHBseparated by D = 60 M . The density is measured along the x -axis. Initial data is that of configuration spin2 in Table III,corresponding to a co-rotating dipole. The density was mea-sured at different instants, and our results show that the late-time profile is in accordance with the non-relativistic predic-tion for a bound state. ponentially in space, far away from the binary. At largedistances, the dominant behavior is controlled by power-law tails [39–42]. Asymptotically, our results are consis-tent with a late-time decay φ lm ∼ t − p sin( µt ) where theexponent p = l + 3 / p = 5 / FIG. 10. Evolution of configuration spin3 in Table III, corre-sponding to a dipolar ( l = m = 1) scalar field counter-rotatingwith the binary. Top panels: evolution of scalar field. Bottompanels: evolution of energy density.FIG. 11. Multipolar l = m = 1 , spin3 . These are thesame initial conditions as those of Fig. 10. The signal mod-ulation is not due to beating of higher overtones, but by thebinary orbital motion. V. CONCLUSIONS
Light scalar fields are interesting solutions to some ofthe most pressing problems in physics. One example isthe dark matter problem. It is tempting to introducefields with a scale (the Compton wavelength) similar tothe size of galactic cores. One would thus be dealingwith fields of mass 10 − eV or similar, in what is knownas fuzzy dark matter models [43, 44]. Understandingthe physics and evolution of compact binaries in suchenvironments is crucial to model their evolution and tosearch for such fields [45, 46].3Our results indicate that, in the presence of a back-ground scalar, the scalar field dynamics close to a BHBparallels very closely that of an electron in a one-electronheteronuclear diatomic molecule. We do note that thereare very important differences between realistic BHBsand molecules. In particular, a BHB is a dissipative sys-tem. Our numerical results for the time evolution of ini-tial data show that non-relativistic bound states turn intoquasi-bound states, via absorption at the horizon.A BHB is dissipative in another way, not included(for simplicity) in our simulations: the system loses en-ergy through gravitational wave emission. We focusedon timescales much shorter than the typical scale for BHcoalescence. Naturally there are systems for which thisassumption is not justified. The typical time scale untilcoalescence is t c ∼ D /M [47] (for simplicity, we as-sume a circular orbit and an equal mass binary). Thetime scale is t c ∼ M for D = 60 M , and t c ∼ M for D = 10 M . Our numerical simulations clearly showthat the quasi-bound state lifetime is at least O (10 ) M .Thus, BHB evolution via gravitational-wave emission isindeed relevant for the evolution of these states, speciallyat small separations, and left for future work.We have not dealt with eccentricity, nor did we con-sider unequal-mass binaries, although it is straightfor-ward to apply our formalism and methods to these situa-tions. Unequal-mass binary evolution might lead, dueasymmetric accretion and drag, to substantial center-of-mass velocities making it specially interesting tostudy [48].In the context of gravitational-wave imprints, dynam-ical friction caused by such fields and its impact on thegravitational-wave phase was recently described [45, 46].However, such description does not include possiblequasi-bound state formation. The clouds have a size1 / ( M µ ) and according to general considerations andspecific calculations [49] they should contribute an ex-tra attracting force which scales linearly with the scalardensity. Since this is a conservative effect, its only conse-quence is a slight renormalization of the binary mass, andwe do not expect changes to the dephasing introduced bydissipative effects [45, 46] (in particular, a dephasing ap-pearing at post-Newtonian order “ −
6” with respect tothe leading vacuum General Relativity prediction).Still in the context of ultralight dark matter, consider aBHB evolving within a “cloud” of coherently oscillatingscalars. This cloud could have primordial origin – andbe a component of ultralight dark matter, or could sim-ply arise as a consequence of superradiant instabilities,and be localized around a supermassive, spinning BH.Now, when this cloud is much larger than any scale in aBHB system, the corresponding boundary conditions aredifferent. It is possible that molecular-like states arise,but their study requires understanding the time evolu-tion of scalar fields with coherent oscillating boundaryconditions [50].Tidal disruption of scalar clouds by orbiting compan-ions was recently discussed [51]. Our results raise the interesting possibility that the final state of such a dis-ruption can be a gravitational molecule.Finally, we note that the response of BHBs to externalfluctuations has been studied recently. In particular, theresponse to high frequency and low frequency scalars wasstudied in toy models [38, 52, 53]. A realistic BHB config-uration revealed already universal ringdown for binaries,and hints of superradiance [5, 14, 54]. Together with theresults we discussed here, these studies show that com-pact binaries are a fertile ground to new phenomenology.
ACKNOWLEDGMENTS
V.C. acknowledges financial support provided un-der the European Union’s H2020 ERC ConsolidatorGrant “Matter and strong-field gravity: New frontiersin Einstein’s theory” grant agreement no. MaGRaTh–646597. M.Z. acknowledges financial support pro-vided by FCT/Portugal through the IF programme,grant IF/00729/2015. T.I. acknowledges financial sup-port provided under the European Union’s H2020ERC, Starting Grant agreement no. DarkGRA–757480.This project has received funding from the EuropeanUnion’s Horizon 2020 research and innovation pro-gramme under the Marie Sklodowska-Curie grant agree-ment No 690904. We acknowledge financial supportprovided by FCT/Portugal through grant PTDC/MAT-APL/30043/2017 and through the bilateral agreementFCT-PHC (Pessoa programme) 2020-2021. The authorswould like to acknowledge networking support by theGWverse COST Action CA16104, “Black holes, grav-itational waves and fundamental physics.” Computa-tions were performed on the “Baltasar Sete-Sois” clus-ter at IST and XC40 at YITP in Kyoto University.This work was granted access to the HPC resources ofMesoPSL financed by the Region Ile de France and theproject Equip@Meso (reference ANR-10-EQPX-29-01) ofthe programme Investissements d’Avenir supervised bythe Agence Nationale pour la Recherche.
Appendix A: Mode excitation by a binary spacetime
To understand which modes are excited by the BHBlet us again consider the metric (2). Expanding the New-tonian potential (3) we obtainΦ (cid:39) − M + M r + M r ( t ) + M r ( t ) r − M ( (cid:126)r · (cid:126)r ( t )) + M ( (cid:126)r · (cid:126)r ( t )) r . (A1)The Klein-Gordon equation on this spacetime can bewritten as (cid:0) − ∂ t + ∇ − µ (cid:1) φ = − ∂ t Φ ∂ t φ − ∇ φ + 2 µ Φ φ . φ − φ hom = (cid:90) dt (cid:48) d (cid:126)x (cid:48) G ( t − t (cid:48) , (cid:126)x − (cid:126)x (cid:48) ) × (cid:0) − ∂ t (cid:48) Φ ∂ t (cid:48) φ − ∇ (cid:48) φ + 2 µ Φ φ (cid:1) ( t (cid:48) , (cid:126)x (cid:48) )where φ hom is the homogeneous solution, fixed by theinitial data. G ( t − t (cid:48) , (cid:126)x − (cid:126)x (cid:48) ) is the Green’s function, G ( t, (cid:126)x ) = (cid:90) d (cid:126)k (2 π ) ω k e − iω k t + i(cid:126)k · (cid:126)x i θ ( t ) + c . c . , (A2)where ω k = (cid:113) µ + (cid:126)k , and c . c . is a complex conjugate.For weak couplings, φ ∼ φ hom , φ − φ hom = (cid:90) dt (cid:48) d (cid:126)x (cid:48) G ( t − t (cid:48) , (cid:126)x − (cid:126)x (cid:48) ) (A3) × (cid:0) − ∂ t (cid:48) Φ ∂ t (cid:48) φ hom − ∇ (cid:48) φ hom + 2 µ Φ φ hom (cid:1) ( t (cid:48) , (cid:126)x (cid:48) ) . Let us now expand the right hand side of this equation in spherical harmonics, φ hom = (cid:88) lm φ hom lm ( t, r ) Y ml ( θ, ϕ ) , (A4)Φ = (cid:88) lm Φ lm ( t, r ) Y ml ( θ, φ ) = Φ ( t, r ) Y ( θ, φ )+ (cid:88) l =2 , , , ··· (cid:88) m = ± Φ lm ( t, r ) Y ml ( θ, φ ) , (A5)whereΦ = M + M r + M R + M R r − M R + M R r Φ l, ± = − M R + M R r c l e ∓ i Ω t and c = (cid:113) π , c = (cid:112) π , c = (cid:113) π , . . . . TheGreen function can also be expanded in spherical har-monics through the plane wave expansion e i(cid:126)k · (cid:126)x = 4 π ∞ (cid:88) l =0 l (cid:88) m = − l i l j l ( kr ) Y ml (ˆ (cid:126)k ) Y ml ∗ (ˆ (cid:126)x ) , where the j l are the spherical Bessel functions and thehatˆdenotes a unit vector. Thus, φ ( t, (cid:126)x ) − φ hom ( t, (cid:126)x ) = (cid:90) dt (cid:48) d (cid:126)x (cid:48) d (cid:126)k (2 π ) ω k (cid:26) e − iω k ( t − t (cid:48) ) i θ ( t − t (cid:48) ) e i(cid:126)k · (cid:126)x (cid:88) l,m (cid:88) l (cid:48) ,m (cid:48) (cid:88) l (cid:48)(cid:48) ,m (cid:48)(cid:48) π ( − i ) l (cid:48)(cid:48) j l (cid:48)(cid:48) ( kr (cid:48) ) Y m (cid:48)(cid:48) l (cid:48)(cid:48) ∗ (ˆ (cid:126)k ) Y m (cid:48)(cid:48) l (cid:48)(cid:48) ( θ (cid:48) , ϕ (cid:48) ) Y ml ( θ (cid:48) , ϕ (cid:48) ) Y m (cid:48) l (cid:48) ( θ (cid:48) , ϕ (cid:48) ) A lml (cid:48) m (cid:48) ( t (cid:48) , r (cid:48) ) + c . c . (cid:27) , where A lml (cid:48) m (cid:48) ( t, r ) = − ∂ t Φ lm ( t, r ) ∂ t φ hom l (cid:48) m (cid:48) ( t, r ) − lm ( t, r ) (cid:18) ∂ rr φ hom l (cid:48) m (cid:48) ( t, r ) + 2 r ∂ r φ hom l (cid:48) m (cid:48) ( t, r ) − l (cid:48) ( l (cid:48) + 1) r φ hom l (cid:48) m (cid:48) ( t, r ) (cid:19) + 2 µ Φ lm ( t, r ) φ hom l (cid:48) m (cid:48) ( t, r ) . Let us consider the modes of the scattered field, (cid:0) φ ( t, (cid:126)x ) − φ hom ( t, (cid:126)x ) (cid:1) lm . Performing the integration in the angulardirections of (cid:126)k we obtain (cid:0) φ ( t, (cid:126)x ) − φ hom ( t, (cid:126)x ) (cid:1) lm = (cid:90) dt (cid:48) dr (cid:48) r (cid:48) dk k πω k e − iω k ( t − t (cid:48) ) i θ ( t − t (cid:48) ) j l ( kr ) j l ( kr (cid:48) ) (cid:88) l (cid:48) m (cid:48) (cid:88) l (cid:48)(cid:48) m (cid:48)(cid:48) (cid:110) Λ m,m (cid:48) ,m (cid:48)(cid:48) l,l (cid:48) ,l (cid:48)(cid:48) A l (cid:48)(cid:48) m (cid:48)(cid:48) l (cid:48) m (cid:48) ( t (cid:48) , r (cid:48) ) + ( − ) m Λ − m,m (cid:48) ,m (cid:48)(cid:48) l,l (cid:48) ,l (cid:48)(cid:48) A l (cid:48)(cid:48) m (cid:48)(cid:48) l (cid:48) m (cid:48) ( t (cid:48) , r (cid:48) ) (cid:111) , (A6)where Λ m,m (cid:48) ,m (cid:48)(cid:48) l,l (cid:48) ,l (cid:48)(cid:48) ≡ (cid:90) d Ω Y ml ( θ, ϕ ) Y m (cid:48) l (cid:48) ( θ, ϕ ) Y m (cid:48)(cid:48) l (cid:48)(cid:48) ( θ, ϕ ) . This integral is nonzero only when m (cid:48) + m = − m (cid:48)(cid:48) . Since, from Eq. (A5), A l (cid:48)(cid:48) m (cid:48)(cid:48) l (cid:48) m (cid:48) ( t, r ) is nonzero only for m (cid:48)(cid:48) = 0 , ± (cid:0) φ ( t, (cid:126)x ) − φ hom ( t, (cid:126)x ) (cid:1) lm is nontrivial only when m = ± m (cid:48) , m = m (cid:48) ± m = − m (cid:48) ±
2, for nontrivialvalues of φ hom l (cid:48) m (cid:48) .5 [1] J. D. Bekenstein, Phys. Rev. D , 1239 (1972).[2] D. Robinson, Phys. Rev. Lett. , 905 (1975).[3] J. D. Bekenstein, in Proceedings of the Second Interna-tional A.D. Sakharov Conference on Physics: Moscow,Russia 20-24 May 1996 , edited by I. M. Dremin andA. M. Semikhatov (World Scientific, 1997) arXiv:gr-qc/9605059 [gr-qc].[4] P. T. Chrusciel, J. L. Costa, and M. Heusler, Living Rev.Relativ. , 7 (2012), arXiv:1205.6112 [gr-qc].[5] R. Brito, V. Cardoso, and P. Pani, Lect. Notes Phys. , pp.1 (2015), arXiv:1501.06570 [gr-qc].[6] S. L. Detweiler, Phys. Rev. D , 2323 (1980).[7] A. Arvanitaki and S. Dubovsky, Phys. Rev. D , 044026(2011), arXiv:1004.3558 [hep-th].[8] D. Baumann, H. S. Chia, J. Stout, and L. ter Haar,JCAP , 006 (2019), arXiv:1908.10370 [gr-qc].[9] N. Sanchis-Gual, M. Zilh˜ao, C. Herdeiro, F. Di Giovanni,J. A. Font, and E. Radu, (2020), arXiv:2007.11584 [gr-qc].[10] K. Clough, T. Dietrich, and J. C. Niemeyer, Phys. Rev.D , 083020 (2018), arXiv:1808.04668 [gr-qc].[11] R. Brito, V. Cardoso, and P. Pani, Class. Quant. Grav. , 134001 (2015), arXiv:1411.0686 [gr-qc].[12] H. Yoshino and H. Kodama, PTEP , 043E02 (2014),arXiv:1312.2326 [gr-qc].[13] C. A. R. Herdeiro and E. Radu, Phys. Rev. Lett. ,221101 (2014), arXiv:1403.2757 [gr-qc].[14] L. Bernard, V. Cardoso, T. Ikeda, and M. Zilh˜ao, Phys.Rev. D , 044002 (2019), arXiv:1905.05204 [gr-qc].[15] D. Baumann, H. S. Chia, and R. A. Porto, Phys. Rev.D , 044001 (2019), arXiv:1804.03208 [gr-qc].[16] Ø. Burrau, Danske Vidensk. Selskab. Math.-fys. Meddel.(in German) M 7 , 1 (1927).[17] A. H. Wilson, Proceedings of the Royal Society of London
A118 , 635 (1928).[18] W. B. Brown and E. Steiner, The Jour-nal of Chemical Physics , 3934 (1966),https://doi.org/10.1063/1.1726555.[19] B. Nickel, Journal of Physics A: Mathematical and The-oretical , 395301 (2011).[20] E. W. Leaver, Journal of Mathematical Physics , 1238(1986).[21] E. Berti, V. Cardoso, and M. Casals, Phys. Rev. D73 ,024013 (2006), [Erratum: Phys. Rev.D73,109902(2006)],arXiv:gr-qc/0511111 [gr-qc].[22] D. Bates and R. Reid (Academic Press, 1968) pp. 13 –35.[23] H. Bethe (Springer, 1933) p. 527.[24] D. Abramov and S. Slavyanov, Journal of Physics B:Atomic and Molecular Physics , 2229 (2001).[25] R. J. Damburg, R. K. Propin, S. Graffi, V. Grecchi, E. M.Harrell, J. c. v. ˇC´ıˇzek, J. Paldus, and H. J. Silverstone,Phys. Rev. Lett. , 1112 (1984).[26] V. Cardoso, A. S. Miranda, E. Berti, H. Witek,and V. T. Zanchin, Phys. Rev. D79 , 064016 (2009),arXiv:0812.1806 [hep-th].[27] S. Weinberg,
Lectures on Quantum Mechanics , 2nd ed.(Cambridge University Press, 2015).[28] B. C. Mundim, H. Nakano, N. Yunes, M. Campanelli,S. C. Noble, and Y. Zlochower, Phys. Rev.
D89 , 084008(2014), arXiv:1312.6731 [gr-qc]. [29] N. K. Johnson-McDaniel, N. Yunes, W. Tichy, andB. J. Owen, Physical Review D (2009), 10.1103/phys-revd.80.124039.[30] F. L¨offler, J. Faber, E. Bentivegna, T. Bode, P. Diener,R. Haas, I. Hinder, B. C. Mundim, C. D. Ott, E. Schnet-ter, G. Allen, M. Campanelli, and P. Laguna, Class.Quantum Grav. , 115001 (2012), arXiv:1111.3344 [gr-qc].[31] M. Zilh˜ao and F. L¨offler, Proceedings, Spring Schoolon Numerical Relativity and High Energy Physics(NR/HEP2): Lisbon, Portugal, March 11-14, 2013 , Int.J. Mod. Phys.
A28 , 1340014 (2013), arXiv:1305.5299 [gr-qc].[32] M. Babiuc-Hamilton, S. R. Brandt, P. Diener, M. Elley,Z. Etienne, G. Ficarra, R. Haas, H. Witek, M. Alcubierre,D. Alic, G. Allen, M. Ansorg, L. Baiotti, W. Benger,E. Bentivegna, S. Bernuzzi, T. Bode, B. Bruegmann,G. Corvino, R. De Pietri, H. Dimmelmeier, R. Doo-ley, N. Dorband, Y. El Khamra, J. Faber, T. Font,J. Frieben, B. Giacomazzo, T. Goodale, C. Gund-lach, I. Hawke, S. Hawley, I. Hinder, S. Husa, S. Iyer,T. Kellermann, A. Knapp, M. Koppitz, G. Lanferman,F. L¨offler, J. Masso, L. Menger, A. Merzky, M. Miller,P. Moesta, P. Montero, B. Mundim, A. Nerozzi, C. Ott,R. Paruchuri, D. Pollney, D. Radice, T. Radke, C. Reis-swig, L. Rezzolla, D. Rideout, M. Ripeanu, E. Schnet-ter, B. Schutz, E. Seidel, E. Seidel, J. Shalf, U. Sper-hake, N. Stergioulas, W.-M. Suen, B. Szilagyi, R. Taka-hashi, M. Thomas, J. Thornburg, M. Tobias, A. Tonita,P. Walker, M.-B. Wan, B. Wardell, M. Zilh˜ao, B. Zink,and Y. Zlochower, “The Einstein Toolkit,” (2019), tofind out more, visit http://einsteintoolkit.org.[33] E. Schnetter, S. H. Hawley, and I. Hawke, Class. Quant.Grav. , 1465 (2004), arXiv:gr-qc/0310042 [gr-qc].[34] Carpet, Carpet: Adaptive Mesh Refinement for the Cac-tus Framework.[35] D. Pollney, C. Reisswig, E. Schnetter, N. Dorband, andP. Diener, Physical Review D (2011), 10.1103/phys-revd.83.044045.[36] H. Witek, M. Zilhao, G. Ficarra, and M. Elley, “Canuda:a public numerical relativity library to probe fundamen-tal physics,” (2020).[37] P. V. P. Cunha, J. A. Font, C. Herdeiro, E. Radu,N. Sanchis-Gual, and M. Zilh˜ao, Phys. Rev. D96 , 104040(2017), arXiv:1709.06118 [gr-qc].[38] T. Assumpcao, V. Cardoso, A. Ishibashi, M. Richartz,and M. Zilhao, Phys. Rev.
D98 , 064036 (2018),arXiv:1806.07909 [gr-qc].[39] S. Hod and T. Piran, Phys. Rev. D , 044018 (1998),arXiv:gr-qc/9801059.[40] H. Koyama and A. Tomimatsu, Phys. Rev. D , 044014(2001), arXiv:gr-qc/0103086.[41] H. Koyama and A. Tomimatsu, Phys. Rev. D , 084031(2002), arXiv:gr-qc/0112075.[42] H. Witek, V. Cardoso, A. Ishibashi, and U. Sperhake,Phys. Rev. D87 , 043513 (2013), arXiv:1212.0551 [gr-qc].[43] V. H. Robles and T. Matos, Mon. Not. Roy. Astron. Soc. , 282 (2012), arXiv:1201.3032 [astro-ph.CO].[44] L. Hui, J. P. Ostriker, S. Tremaine, and E. Witten,Phys. Rev.
D95 , 043541 (2017), arXiv:1610.08297 [astro-ph.CO]. [45] L. Annulli, V. Cardoso, and R. Vicente, (2020),arXiv:2007.03700 [astro-ph.HE].[46] L. Annulli, V. Cardoso, and R. Vicente, (2020),arXiv:2009.00012 [gr-qc].[47] P. Peters, Phys. Rev. , B1224 (1964).[48] V. Cardoso and C. F. Macedo, (2020), arXiv:2008.01091[astro-ph.HE].[49] R. Feynman, Phys. Rev. , 340 (1939).[50] K. Clough, P. G. Ferreira, and M. Lagos, Phys. Rev. D , 063014 (2019), arXiv:1904.12783 [gr-qc]. [51] V. Cardoso, F. Duque, and T. Ikeda, Phys. Rev. D ,064054 (2020), arXiv:2001.01729 [gr-qc].[52] K. Nakashi and T. Igata, Phys. Rev. D , 124033 (2019),arXiv:1903.10121 [gr-qc].[53] K. Nakashi and T. Igata, Phys. Rev. D , 104006(2019), arXiv:1908.10075 [gr-qc].[54] L. K. Wong, Phys. Rev. D100