The physics of black hole binaries: geodesics, relaxation modes and energy extraction

Black holes are the simplest macroscopic objects, and provide unique tests of General Relativity. They have been compared to the Hydrogen atom in quantum mechanics. Here, we establish a few facts about the simplest systems bound by gravity: black hole binaries. We provide strong evidence for the existence of `global' photosurfaces surrounding the binary, and of binary quasinormal modes leading to exponential decay of massless fields when the binary spacetime is slightly perturbed. These two properties go hand in hand, as they did for isolated black holes. The binary quasinormal modes have high quality factor and may be prone to resonant excitations. Finally, we show that energy extraction from binaries is generic and we find evidence of a new mechanism -- akin to the Fermi acceleration process -- whereby the binary transfers energy to its surroundings in a cascading process. The mechanism is conjectured to work when the individual components spin, or are made of compact stars.


I. INTRODUCTION
Einstein's theory of general relativity (GR) is the most accurate known description of gravity [1][2][3][4]. One of its outstanding predictions is the existence of black holes (BHs), vacuum spacetimes defined by an event horizon, i.e., a null one-way surface. Isolated BHs have been studied for decades. They are extremely simple [5] and fully characterized by their mass and angular momentum [7][8][9][10][11][12][13][14]. These properties are instrumental to building templates of gravitational-wave (GW) signals generated by dynamical BHs, which eventually led to the first direct detection of GWs [15]. The lack of complex multipolar structure of BH geometries is crucial to perform strong-field tests of the theory, for example, through the late-time relaxation of BHs, as a superposition of quasinormal modes (QNMs) [14,[16][17][18][19][20][21].
By contrast, and due to their inherent complexity, BH binaries (BHBs) are less well studied and understood: their GW output and dynamics, when in isolation, are known very well through post-Newtonian expansion techniques at large separations [22]. In this approach the individual binary components are stationary vacuum BHs, slightly deformed in response to the companion's field. The dynamical behavior of the BHB itself is poorly understood. Such knowledge can in principle be obtained using numerical methods [23,24]; however, such techniques only probe relatively small time scales using finely tuned initial data. In particular, efforts to date have mostly focused on purely vacuum spacetimes describing isolated BHBs which have been evolving solely through GW emission, leading to an inspiral and merger, possibly observable by current or future GW detectors. The simulation time scales can be at most of the order of a few thousand GM=c 3 , where M is the total spacetime mass. For stellar-mass components, these are of the order of one second or less, but BHBs can live for millions of years on tight orbits. Thus, new effects may be triggered and relevant on large time scales. Do perturbed BHBs also have characteristic ringdown modes, and can they be resonantly excited? Do BHBs amplify incoming, low-frequency radiation? Here, we provide a framework for studying these open questions and answer some of them.
BHB spacetime discussed in Ref. [25]. The construction relies on the theory of matched asymptotic expansions and proceeds as follows. The spacetime is divided into three different regions; see Fig. 2 in Ref. [25]. There are two inner zones sufficiently close to each BH, characterized by BH perturbation techniques (the metric perturbation is described by tidal fields generated by the companion; the tidal moments include the quadrupole and octupole deformation and their time derivatives). The inner zone is "stitched" to a near zone described by a post-Newtonian expansion (to second post-Newtonian order), itself matched to a far zone which is described using a multipolar, post-Minkowskian formalism. In addition, there are two buffer zones defined as the region of spacetime where the three main zones overlap. The existence of such overlapping regions is crucial for constructing the matched metric. In particular, the multipolar post-Newtonian formalism used to build the near-and far-zone metrics ensures that they are matched by construction [22]. Such a spacetime was implemented and used to investigate the physics of accretions disks in the presence of BHBs [26][27][28][29].
Henceforth we set G ¼ c ¼ 1, and focus exclusively on equal-mass binaries of total Arnowitt-Deser-Misner (ADM) mass M separated by a coordinate distance L (proper separations are very close to the coordinate distances we discuss). "Nearly closed" (i.e., curves of period 2π to an accuracy which increases with separation L and which selfintersect at least once) null geodesics in such a construction for a vacuum BHB are shown in Fig. 1. We find three types of null geodesics: the "standard" null circular geodesic surrounding each BH (which in Schwarzschild coordinates sits at a coordinate distance ∼r H =2 away from the horizon, where r H is the Schwarzschild radius of a single BH), a global nonintersecting geodesic, and a figure-of-eightshaped trajectory. Such null geodesics were found in toy models in the past, in extremal BH or analog spacetimes [30][31][32]. The distance of closest approach of the global nonintersecting geodesic to each BH horizon is ð0.6 AE 0.03Þr H and its period T ¼ 2L þ 31ðAE1Þ for the separations we study (L ¼ ½20; 40M). Its distance of closest approach to the center of mass is not strongly sensitive to L and is y=M ¼ 3.5 AE 0.07 for these separations. It is natural to speculate that such global geodesics can be identified with the null circular geodesic of the final BH, when such a binary merges. The instability of such orbits is clear from their numerical search (fine-tuning is necessary). A complete picture of geodesic motion in BHB spacetimes is outside the scope of this work.

III. SCATTERING AND BINARY RELAXATION: INDIVIDUAL AND GLOBAL QNMs, AND POWER-LAW TAILS
The closed null geodesics of isolated BHs correspond to a semitrapping of massless waves; accordingly, they provide useful information on the characteristic modes of vibration (QNMs) of BHs [19,20,33,34]. It is thus natural to associate the previous global null geodesics with binary relaxation, i.e., with the hitherto unknown QNMs of BHBs. An analytic understanding of these issues for BHBs is challenging, and we turn instead to the numerical simulation of massless scalar fields in the BHB background spacetime described above.
The dynamics of the BHB is governed by vacuum GR, as described previously. We ignore the backreaction of the test scalar field on the BHB spacetime, an approximation which is valid for all realistic setups known to us (except when the scalar field mimics GWs and the BHB is in the last stages of inspiral). Thus, the scalar field is governed by the Klein-Gordon equation □Φðt; ⃗ xÞ ¼ 0, in a known (albeit timedependent) background. We use purely ingoing initial data of the form where r is the radial coordinate and WðrÞ is a window function that smooths the r ¼ 0 behavior. Here, r 0 and σ characterize the typical radius and width of the initial ingoing scalar field, while ω characterizes its frequency. The initial field amplitude is irrelevant, since the Klein-Gordon equation is linear. To numerically evolve the scalar field, we employ the code presented in Refs. [30,35], which makes use of the EINSTEIN TOOLKIT infrastructure [36][37][38] with the CARPET package [39,40]. We project the scalar field onto scalar harmonics, The initial data is spherically symmetric around the center of coordinates, but the presence of the BHB guarantees that upon evolution other components will exist. Note that for symmetry reasons only even modes are FIG. 1. Illustration of closed null geodesics in the background of nearly static, nonspinning BHBs. Each BH, depicted by a full black circle, is surrounded by its own light ring, lying at r ¼ 3M in standard Schwarzschild coordinates, when their separation is large. Two other null closed trajectories are possible: a global nonintersecting null geodesic and a figure-of-eight-shaped trajectory. Such null geodesics were found in extremal BH or analog spacetimes [30][31][32]; we also find them in the matched spacetime describing a vacuum BHB. A similar illustration also describes some timelike geodesics. excited in the current setup (see also Ref. [41]). We studied a variety of different initial parameters and BHB separations. A typical outcome of the evolution is shown in Fig. 2 (further details are provided in the Appendix A). The binary is separated by L ¼ 10M [42]. Our results show fourthorder convergence.
As seen in Fig. 2, the dominant mode is the monopolar one and it drives the dynamics. Initially, the observer sitting at r ¼ 100M sees the field heading towards the BHB. At t ∼ 100M the observer starts receiving signals that interacted with the BHB. The first clear signal is a strongly damped sinusoid, associated with the ringdown of each individual BH in the binary. These individual modes are well studied and our results are consistent with theoretical expectations based on linearized calculations [17,43]. After t ∼ 100M, the leading monopolar component is now outgoing and also drives higher multipoles.
After the initial driving monopolar mode dies away at t ∼ 320M, another exponentially damped sinusoid is apparent in the waveform. This is one of our main results: BHBs possess global QNMs which describe the ringdown of the binary as a whole. Our results indicate that the ringdown parameters depend only on the mass and separation of the binary and are (to the extent probed by our simulations) independent of the initial data. It is compelling to associate such BHB modes with the closed null geodesics around the BHB, as can be done formally for isolated BHs [19,20,33,34]. A reassuring check on the correspondence is that our results are well described by the linear fit in the range of separations we studied, where we made use of the simulations listed in the Appendix A. The geodesic calculation would have predicted, for l ¼ m ¼ 2 modes, T ∼ L þ 16. In other words, our numerical results are in suggestive agreement with geodesic expectations, lending further support to the association of the global QNMs with global geodesics. Our results for the decay time scale during this global ringdown phase are less accurate, but suggest that the relaxation time scale τ also increases with the separation and is of order 10M for L ¼ 10M.
Notice that the insight from geodesic motion hinges on having two static BHs. This crude approximation is however robust for BHB for which the post-Newtonian expansion is trustworthy: the orbital period scales like L 3=2 , whereas the light travel time scales like L. Already for L ¼ 10 we find an orbital period of 200M, and a global ringdown time scale of order 10M. In other words, the binary only travels a few degrees during the global ringdown stage. The approximation is even better for larger separations. Characteristic vibration modes are generic properties of dissipative systems: a binary system of two stars may also have global QNMs without allowing closed null geodesics [21,44]. Although all of our results are compatible with such an association, further work is necessary to understand the details and origin of the BHB QNMs.
At very late times, our simulations indicate that the signal dies as a power-law tail in time, Φ ∼ t −γ . For isolated, nonspinning BH spacetimes and for scalar fluctuations, γ ¼ 2l þ 3 [45]. Our results indicate that γ ∼ 7 for l ¼ m ¼ 0, which is presumably an indication of mode mixing during the evolution, of the kind seen in isolated but spinning BH geometries [46,47].
When the initial wave packet has a very large width (corresponding to a constant force on the BHB), the outgoing pulse is modulated at frequencies ω AE kΩ, where k is an integer. Such effects were previously seen in the scattering of electromagnetic waves by periodically moving obstacles [48,49].

IV. ENERGY EXTRACTION AND INSTABILITIES
Compact binaries are astrophysical blenders and potential energy sources, either when surrounded by accretion disks or in the context of fundamental massive fields. Different mechanisms may be associated with energy extraction in the presence of a compact binary: (i) If the individual objects spin, there are ergoregions in the spacetime and each binary component can transfer rotational energy to bosonic fields through superradiance. Such transfer can be turned into an instability by placing the system inside a cavity [50][51][52]. It is unknown whether binary-intrinsic ergoregions exist (but there are arguments suggesting that superradiance does exist for binaries made of nonspinning objects [53]). (ii) Awell-known Newtonian energy extraction processthe gravitational slingshot-transfers kinetic energy from moving planets or stars to scattered probe objects; it is straightforward to show that, within GR, such a mechanism also occurs with light. Light scattering off a BH moving with velocity v can extract (kinetic) energy with an efficiency ϒ ≡ E final =E initial up to The maximum efficiency occurs when the photon scatters with a 180°angle off of an oppositely moving BH, and is identical to the energy gain of a photon scattering off of a moving mirror. For velocities associated with orbital motion in a compact binary, the efficiency can be 1.2 or higher. We verified such a result via explicit scatterings of photons off of moving BH and BHB geometries [54]. It is conceivable that one photon suffers multiple scatterings with the binary (specially if confined). (iii) A binary provides a periodic force on external fields; a similar lower-dimensional toy model is known to give rise to instabilities in trapped radiation [56], akin to the Fermi acceleration process of cosmic rays [57]. These effects or others may all be part of the astrophysics of compact binaries (and hence are all part of realistic simulations, although perhaps not easily identifiable [58,59]). We will be interested in confined binaries (which are of interest to some dark matter scenarios) where the above mechanisms are expected to trigger instabilities. Unfortunately, a numerical investigation of these issues using the previous (3 þ 1) setting is challenging: time scales for energy extraction are expected to be very large, and for nonspinning BHBs (the ones we are currently able to simulate in our setup) absorption at the horizon will likely quench or strongly suppress any possible energy extraction mechanism (but see Ref. [53]); a BH of mass M BH has an absorption cross section (for scalars) of 20kπM 2 BH [60-62] with k ¼ 27=20, 16=20 ¼ Oð1Þ for highand low-frequency radiation, respectively. Because of this, a naive expectation is that a BHB in a cavity of size R ext will contribute to a decrease in the energy inside that cavity at a rate dE=dt ∼ −λE, λ ∼ 10kðM=LÞ 2 =R ext .
To emulate compact binaries of spinning BHs or neutron stars, we take a binary of two reflecting objects in flat (2 þ 1) dimensions, and we evolve a massless scalar with a simple Gaussian initial profile. This setup allows for the evolution of a very large number of orbits and reflections in the confining cavity. We evolve the system numerically for different values of the orbital frequency Ω, separation L, and cavity size R ext with the EINSTEIN TOOLKIT infrastructure using the code described in Ref. [30]. The total integrated energy inside the cavity is shown in Fig. 3 for one such binary. The total energy increases with time, and at late times this increase is exponential. Our results indicate that this growth happens only when the orbital frequency is of the order of the light travel time inside the cavity. Although this is a toy model, this is probably the first example of instabilities triggered by binaries. There is nothing intrinsic to lower-dimensional spacetimes: the arguments above suggest that it may have a counterpart in (3 þ 1) setups as well. Enclosed BHBs are currently being studied [63].

V. DISCUSSION
The role and simplicity of BHs in GR has often been compared to the hydrogen atom in the development of quantum mechanics. It is compelling to draw a parallel between BH binaries and the simplest possible moleculethat of the hydrogen molecule ion [64][65][66][67]. The null geodesics around isolated BHs have a counterpart in wave dynamics as QNMs of the geometry. In a quantummechanical picture they would correspond to the bound states of the hydrogen atom. Likewise, the global geodesics for BHBs (see Fig. 1) may be tightly connected to global QNMs (Fig. 2), the analog of molecular bound states. It is amusing to note that such a correspondence can be taken a step further for an exact BHB solution: the Majumdar Papapetrou geometry describing a pair of extremal BHs. The Klein-Gordon equation in this background assumes the same form as the Schrödinger equation describing a positron in the ionized hydrogen molecule (see Appendix B). The successful, effective-one-body treatment of post-Newtonian theory for the BHB problem uses very explicit hydrogen-atom analogies to construct an inspiraling waveform model [68]. It is tempting to expect that such methods can also yield insight into the problem of BHB relaxation. A quantitative correspondence can be established between geodesics and QNMs of isolated BHs; a FIG. 3. Total energy in a scalar field inside a cavity composed of a circular reflecting boundary and two reflecting circles in circular orbit around each other. The boundary radius is R ext ¼ 30, and each object has a radius 0.5 and is on a circular orbit of radius 5 with angular velocity Ω ¼ 0.14. similar analysis for BHBs is missing, but would be fundamental for the program of understanding the highfrequency QNM regime and its possible connection to null geodesics. BH spectroscopy makes use of the QNMs of isolated BHs to infer their mass and spin [16,69,70]. Perhaps BHB spectroscopy is within reach, where their separation can also be estimated from the global modes. The excitation mechanism of BHB modes could happen, for example, via three-body interactions.
It is a fact that, since the QNMs of isolated BHs are associated with photospheres (and hence to the largest frequency in such a spacetime), the resonant excitation of such modes is impossible in astrophysical setups. However, the new global modes associated with the entire BHB geometry are associated with a new scale (the binary separation), and have a lower frequency than each of the individual BH QNMs. It is therefore conceivable that a particle orbiting one of the BHs can resonantly excite the global modes. The condition for this to happen is that the orbital period of the small particle equals the period of the QNMs. We find that the particle orbiting at radius r p around one of the BHs should satisfy r p =M BH ¼ 0.466ðL=M þ 3π ffiffi ffi 3 p =2Þ 2=3 (M BH is the mass of each individual BH). For L ¼ 38M, for example, resonant excitation is possible when the particle reaches the innermost stable circular orbit of one individual BH. The consequence could be enhanced emission of GWs from the binary.
There are compelling arguments suggesting that compact binaries may be prone to energy transfer to other degrees of freedom (most notably to fundamental scalars, for example) when inside a cavity. This artificial setup accurately mimics physically motivated scenarios consisting of massive degrees of freedom, such as massive scalars or vectors, which are proxies or serious candidates for dark matter (for instance, when dealing with axionlike particles) [52]. Our results strongly suggest that a new type of instability may be active, which could potentially improve constraints on dark matter models. The mechanism resembles a parametric instability [71], but we find complex growth patterns in the field. The instability growth rate is larger for larger angular velocities, and thus a natural question remains open, which our toy model is unable to answer: is the instability relevant for an astrophysical binary driven by GW emission? In other words, is there any regime during which the growth rate is important during a binary's lifetime? [72].  Table I summarizes the parameters of our simulations for the scattering of Gaussian wave packets off BHBs. TABLE I. Summary of our simulations. The parameters specify the initial conditions, as in Eq. (1) in the main text. Here the mass of each BH is fixed at 0.5. T 20 is the period of oscillations of the l ¼ 2, m ¼ 0 component at late times. Our results show that the l ¼ m ¼ 2 case has a similar period. The oscillation period T 20 is extracted using two periods of the late-time oscillation, when available. For the BHB8 run, the late-time behavior of T 20 is highly modulated and it is hard to extract any characteristic frequency.

Numerical procedure and convergence analysis
To numerically evolve the scalar field equations in our prescribed metric background we employ the code presented in Refs. [30,35], which makes use of the EINSTEIN TOOLKIT infrastructure [36][37][38] with the CARPET package [39,40] for mesh-refinement capabilities. We employ the method of lines, where spatial derivatives are approximated by fourth-order finite-difference stencils, and we use the fourth-order Runge-Kutta scheme for the time integration. Kreiss-Oliger dissipation is applied to evolved quantities in order to damp high-frequency noise.
Our simulations use finite-difference techniques, which approximate the continuum solution of the problem with an error that depends polynomially on the grid spacing h, where n is the convergence order. Since the code we employ uses both second-and fourth-order techniques we expect this to be reflected in the convergence properties of our results. Consistency can be checked by evolving the same configuration with coarse, medium, and fine resolutions (h c , h m , and h f ). One can then compute the convergence factor given by To check the convergence of the extracted waveforms we have evolve the configuration of Fig. 2 in the main text, with resolutions h c ¼ 1.6M, h m ¼ 1.28M, and h f ¼ 1.0M (where this refers to the resolution of the outermost refinement level); the corresponding results are shown in Fig. 4 for the l ¼ 2, m ¼ 2 multipole of Φ. We have amplified the differences between the medium-and fineresolution runs by the factor 2.97 expected for fourth-order convergence. Figure 5 shows the waveform extracted at different radii, aligned by their maxima (in other words, aligned by the light-travel-time propagation delay). The consistent overlap between the three signals indicates that the signal is indeed being measured in the wave zone, and that there is no finiteextraction radius artifact. Notice that the early-time response for the monopole does not align because the initial pulse is ingoing. On the other hand, the late-time  Fig. 2 of the main text (BHB3) for the l ¼ m ¼ 0 (top) and l ¼ 2, m ¼ 0 (bottom) modes. Here we align the waveform extracted at three different radii by properly subtracting the light travel time. As one can see, there is a very good overlap of the waveform for the l ¼ 2 mode, indicating that it is purely outgoing and that one is indeed in the wave zone already and capturing the leading-order dependence of the waveform. The driving l ¼ 0 mode, on the other hand, is not aligned at early times, showing that indeed it is initially ingoing. behavior of the monopolar component is perfectly aligned, showing that the pulse is outgoing at these late stages. Figure 6 shows the waveforms for different types of initial conditions, following the definitions in Table I. The waveforms are aligned in time to show a clear universal ringdown, which is the most important result of our work. The ringdown frequency and damping time scale are the same for different initial conditions, and depend only on the binary parameters (mass and separation). The results also clearly indicate that, although the period and damping time scale with separation L, the quality factor seems to be scale independent.

Dependence on initial data
The initial pulse is spherically symmetric. However, on very short time scales the signal develops a quadrupolar (l ¼ 2) component as well. This behavior is expected, since one is specifying spherically symmetric initial conditions on a nonsymmetric background. Thus, the field will very quickly sense the nonsymmetric background metric. It is possible to show that such a nonsymmetric component is weaker at larger distances. There is an exact solution in general relativity describing two or more static BHs. Such a solution is known as the Majumdar-Papapetrou (MP) solution [73][74][75]. In a cylindrical coordinate system the BHB version of the MP solution is written as

Power-law tails
with This solution represents two maximally charged BHs in equilibrium, each with mass M and charge Q ¼ M. In these coordinates, their horizons are shrunk to two points at z ¼ AEa (and hence the parameter a measures the distance between them). The spacetime ADM mass is 2M. Some geodesic properties of these solutions were studied before [30][31][32]. We will now show the rather remarkable result that the Klein-Gordon equation separates. We change FIG. 6. Waveforms for different initial data, following the notation of Table I  to planar "prolate confocal elliptical" coordinates χ, η (keeping the time and azimuthal coordinates) defined as r 2 1 ≡ ρ 2 þ ða − zÞ 2 ¼ a 2 ðχ þ ηÞ 2 ; ðB3Þ r 2 2 ≡ ρ 2 þ ða þ zÞ 2 ¼ a 2 ðχ − ηÞ 2 : The variable χ plays a role similar to r in standard spherical coordinates, while η plays the role of cos θ. The domains of these variables are −1 ≤ η ≤ 1 and 1 ≤ χ ≤ ∞. In these coordinates the Klein-Gordon equation for Ψe imϕ−iωt can be written as When a=M ≫ 1, we find that the above is separable and reduces to With the ansatz Ψ ¼ SðηÞRðχÞ, we finally find where Λ is a separation constant. This same system describes the Schrödinger equation (for a positron) in the ionized hydrogen molecule [64,65]. We thus have a formal equivalence between two similar systems: that of a molecule governed by electromagnetism and a simple binary system in full general relativity. The effectiveone-body treatment of post-Newtonian theory for the BHB problem uses very explicit hydrogen-atom analogies to construct an inspiraling waveform model [68]. Thus, it is interesting that the converse (i.e., recovering the dynamics of a molecule in quantum mechanics) is borne out of an exact solution in general relativity.