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

TL;DR

This paper explores the fundamental physics of black hole binaries, demonstrating the existence of photosurfaces and quasinormal modes, and proposing a new energy extraction mechanism that could operate in systems with spinning components or stars.

Contribution

It establishes the presence of global photosurfaces and quasinormal modes in black hole binaries and introduces a novel energy transfer mechanism similar to Fermi acceleration.

Findings

Existence of binary photosurfaces confirmed.

Binary quasinormal modes exhibit high quality factors.

Energy transfer via a cascading process analogous to Fermi acceleration.

Abstract

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.

Figures

Tables

Table 1. Table 1: 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 0.5 0.5 . T 20 subscript 𝑇 20 T_{20} is the period of oscillations of the l = 2 , m = 0 formulae-sequence 𝑙 2 𝑚 0 l=2,m=0 component at late-times. Our results show that the l = m = 2 𝑙 𝑚 2 l=m=2 has a similar period. The oscillation period T 20 subscript 𝑇 20 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 subscript 𝑇 20 T_{20} is highly modulated and hard to extract any characteristic frequency.

Name	$L$	$\sigma$	$\omega$	$r_0$	$T_{20}$
BHB1	$20$	$40$	$0.1$	$100$	26.6
BHB2	$20$	$40$	$0.05$	$100$	28.0
BHB3	$10$	$40$	$0.1$	$100$	19.0
BHB4	$40$	$40$	$0.1$	$100$	49.6
BHB5	$20$	$40$	$0.01$	$100$	27.2
BHB6	$20$	$40$	$0.02$	$100$	27.0
BHB7	$20$	$80$	$0.2$	$100$	28.6
BHB8	$20$	$80$	$0.5$	$100$	*
BHB9	$10$	$4$	$0.1$	$100$	18.6

Equations

$$\Phi (0,\vec{x})$$

$${\partial }_t\Phi (0,\vec{x})$$

$$T/M=(1.03\pm 0.04)L/M+8\pm 1,$$

$${{\rm Y}}_{\max }=\frac{1+v}{1-v}.$$

$$f=f_h+O\left(h^n\right),$$

$$Q\equiv \frac{f_{h_c}-f_{h_m}}{f_{h_m}-f_{h_f}}=\frac{h_c^n-h_m^n}{h_m^n-h_f^n}.$$

$$d{s}^2=-\frac{d{t}^2}{U^2}+U^2\left(d{\rho }^2+{\rho }^{2}d{\varphi }^2+d{z}^2\right),$$

$$U(\rho ,z)=1+\frac{M}{\sqrt{{\rho }^2+(z-a{)}^2}}+\frac{M}{\sqrt{{\rho }^2+(z+a{)}^2}}.$$

$$r_1^2\equiv {\rho }^2+(a-z{)}^2=a^2(\chi +\eta {)}^2,$$

$$r_2^2\equiv {\rho }^2+(a+z{)}^2=a^2(\chi -\eta {)}^2.$$

$$0$$

$${\partial }_{\chi }\left(({\chi }^2-1){\partial }_{\chi }\Psi \right)-{\partial }_{\eta }\left(({\eta }^2-1){\partial }_{\eta }\Psi \right)$$

$${\partial }_{\eta }\left((1-{\eta }^2){\partial }_{\eta }S\right)+\left(-a^2{\omega }^2{\eta }^2-\frac{m^2}{1-{\eta }^2}+\Lambda \right)S=0,$$

Full text

I. Introduction. Einstein’s theory of General Relativity (GR) is the most accurate known description of gravity Will (2018a, b, 2014); Berti et al. (2015). One of its outstanding predictions is the existence of black holes (BHs), vacuum spacetimes defined by an event horizon, a null one-way surface. Isolated BHs have been studied for decades. They are extremely simple 111This simplicity is best epitomized in Chandrasekhar’s words: “In my entire scientific life, extending over forty-five years, the most shattering experience has been the realization that an exact solution of Einstein’s equations of GR provides the absolutely exact representation of untold numbers of black holes that populate the universe” Cha ., and fully characterized by their mass and angular momentum Kerr (1963); Newman et al. (1965); Bekenstein (1997); Carter (1999); Chrusciel et al. (2012); Robinson (1975, 2009); Cardoso and Gualtieri (2016). These properties are instrumental to build templates of gravitational wave (GW) signals generated by dynamical BHs, which eventually led to the first direct detection of GWs Abbott et al. (2016a). 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) Berti et al. (2006, 2009); Cardoso and Gualtieri (2016); Abbott et al. (2016b); Cardoso et al. (2016); Cardoso and Pani (2017, 2019).

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, is known very well through post-Newtonian expansion techniques at large separations Blanchet (2014). 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 Pretorius (2007); Cardoso et al. (2015); however, such techniques only probe relatively small timescales using finely tuned initial data. In particular, efforts to date focus mostly 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 timescales can be at most of order of a few thousands $GM/c^{3}$, with $M$ the total spacetime mass. For stellar-mass components, these are of 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 timescales. 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.

II. Setup: a black hole binary in post-Newtonian theory. Consider a BHB, whose dynamics are governed by vacuum Einstein equations. For the reasons outlined (computational expenses, dependence on initial data, etc), instead of numerically-generated spacetimes, we use the approximate BHB spacetime discussed in Ref. Mundim et al. (2014). The construction relies on the theory of matched asymptotic expansions and proceeds as follows. The spacetime is divided in three different regions, see Fig. 2 in Ref. Mundim et al. (2014). There are two inner zones sufficiently close to each BH, characterized with 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 to 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 Blanchet (2014). Such spacetime was implemented and used to investigate the physics of accretions disks in the presence of BHBs Noble et al. (2012); Zilhão et al. (2015); Bowen et al. (2018, 2017).

Here onwards we set $G=c=1$, and focus exclusively on equal mass binaries, of total 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\pi$ to an accuracy which increases with separation $L$ and which self-intersect at least once) null geodesics in such a construction for a vacuum BHB are shown in Fig. 1. We found three types of null geodesics: the “standard” null circular geodesic surrounding each BH (which in Schwarzschild coordinates sits at a coordinate distance $\sim r_{H}/2$ away from the horizon, with $r_{H}$ the Schwarzschild radius of a single BH), a global non-intersecting geodesic and an eight-shaped trajectory. Such null geodesics were found in toy models in the past, in extremal BH or analog spacetimes Assumpcao et al. (2018); Shipley and Dolan (2016); Chandrasekhar (1989). The distance of closest approach of the global non-intersecting geodesic to each BH horizon is $(0.6\pm 0.03)r_{H}$ and its period $T=2L+31\,(\pm 1)$ for the separations we studied ($L=[20,40]M$). Its distance of closest approach to the center of mass is not strongly sensitive on $L$ and is $y/M=3.5\pm 0.07$ for these separations. It is natural to speculate that such global geodesic can be identified with the null circular geodesic of the final BH, when such binary merges. 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 semi-trapping of massless waves; accordingly, they provide useful information on the characteristic modes of vibration (QNMs) of BHs Ferrari and Mashhoon (1984); Cardoso et al. (2009, 2016); Cardoso and Pani (2017). 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 $\square\Phi(t,\vec{x})=0$, in a known (albeit time-dependent) background. We use purely ingoing initial data of the form

[TABLE]

where $r$ is the radial coordinate and $W(r)$ is a window function that smooths the $r=0$ behavior. Here, $r_{0},\,\sigma$ characterize the typical radius and width of the initial ingoing scalar field, while $\omega$ 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. Cunha et al. (2017); Assumpcao et al. (2018), which makes use of the EinsteinToolkit infrastructure Löffler et al. (2012); EinsteinToolkit ; Zilhão and Löffler (2013) with the Carpet package Schnetter et al. (2004); Carpet .

We project the scalar field onto scalar harmonics, $\Phi=\sum_{lm}\Phi_{l,m}Y_{lm}$. 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 excited in the current setup (see also Ref. Bentivegna et al. (2008)). 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). The binary is separated by $L=10M$ 222Such separation is at the limits of applicability of the post-Newtonian approximation Mundim et al. (2014); results for larger separations show the same qualitative features.. Our results show fourth-order convergence (see Appendix).

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\sim 100M$ the observer starts receiving signals which interacted with the BHB. The first clear signal is a strongly damped sinusoid, associated to 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 Berti et al. (2009); GRI . After $t\sim 100M$, the leading monopolar component is now outgoing and drives also higher multipoles.

After the initial driving monopolar mode dies away at $t\sim 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 on 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 Ferrari and Mashhoon (1984); Cardoso et al. (2009, 2016); Cardoso and Pani (2017). A reassuring check on the correspondence is that our results are well described by the linear fit

[TABLE]

in the range of separations we studied, where we made use of the simulations listed in the Supplement. The geodesic calculation would have predicted, for $l=m=2$ modes $T\sim 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 timescale during this global ringdown phase are less accurate, but suggest that the relaxation timescale $\tau$ 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 timescale of order $10M$. In other words, the binary only traveled 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 Kokkotas and Schmidt (1999); Cardoso and Pani (2019); although all 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, $\Phi\sim t^{-\gamma}$. For isolated, non-spinning BH spacetimes and for scalar fluctuations $\gamma=2l+3$ Ching et al. (1995). Our results indicate $\gamma\sim 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 Zenginoglu et al. (2014); Burko and Khanna (2014).

When the initial wavepacket has a very large width (corresponding to a constant forcing on the BHB), the outgoing pulse is modulated at frequencies $\omega\pm k\Omega$, with $k$ an integer. Such effects were previously seen in the scattering of electromagnetic waves by periodically moving obstacles van Bladel (1976); Cooper (1980).

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 Zel’dovich (1971, 1972); Brito et al. (2015). It is unknown whether binary-intrinsic ergoregions exist (but there are arguments suggesting that superradiance does exist for binaries made of non-spinning objects Wong (2019)).

ii. A well-known Newtonian energy extraction process, the 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 $\Upsilon\equiv E_{\rm final}/E_{\rm initial}$ up to

[TABLE]

The maximum efficiency occurs when the photon scatters with a $180^{\circ}$ angle off an oppositely-moving BH, and is identical to the energy gain of a photon scattering off a moving mirror. For velocities associated to orbital motion in a compact binary, the efficiency can be $1.2$ or higher. We verified such result via explicit scatters of photons off moving BH and BHB geometries Tho . It is conceivable that one photon suffers multiple scatters with the binary (specially if confined).

iii. A binary provides a periodic forcing to external fields; a similar lower-dimensional toy model is known to give rise to instabilities in trapped radiation Cooper (1993), akin to the Fermi acceleration process of cosmic rays Fermi (1949).

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 van Meter et al. (2010); D’Ascoli et al. (2018)). We will be interested in confined binaries, 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; timescales for energy extraction are expected to be very large, and for non-spinning 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. Wong (2019)): a BH of mass $M_{\rm BH}$ has an absorption cross section (for scalars) of $20k\pi M_{\rm BH}^{2}$ Chandrasekhar (1983); Macedo et al. (2013); Das et al. (1997) with $k=27/20,16/20={\cal O}(1)$ for high and low frequency radiation respectively. Because of this, a naive expectation is that a BHB in a cavity of size $R_{\rm ext}$ will contribute to a decrease in the energy inside that cavity at a rate $dE/dt\sim-\lambda E$, $\lambda\sim 10k(M/L)^{2}/R_{\rm 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 simple gaussian initial profile. This setup allows for a very large number of orbits and reflections in the confining cavity to be evolved. We evolved the system numerically, for different values of orbital frequency $\Omega$, separation $L$ and cavity size $R_{\rm ext}$, with the EinsteinToolkit infrastructure using the code described in Ref. Assumpcao et al. (2018). 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 such 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 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 Ber .

V. Discussion. The role and simplicity of BHs in GR has often been compared to the Hydrogen atom in the development on quantum mechanics. It is compelling to draw a parallel between BH binaries and the simplest possible molecule, that of the Hydrogen molecule ion Burrau (1927); Wilson (1928); Hylleraas (1931); Leaver (1986). The null geodesics around isolated BHs have a counterpart in wave dynamics, as QNMs of the geometry. In a quantum-mechanical 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 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). 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 Buonanno and Damour (1999). It is tempting to expect that such methods can also yield insight on the problem of BHB relaxation. A quantitative correspondence can be established between geodesics and QNMs of isolated BHs; a similar analysis for BHBs is missing, but would be fundamental for the program of understanding the high frequency 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 Berti et al. (2006, 2007, 2016). 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 to photospheres – and hence to the largest frequency in such 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 to a new scale (the binary separation), and have 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_{\rm BH}=0.466(L/M+3\pi\sqrt{3}/2)^{2/3}$ ($M_{\rm 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 mimics well physically motivated scenarios consisting on massive degrees of freedom, such as massive scalars or vectors, which are proxies or serious candidates for dark matter (for instance when dealing with axion-like particles) Brito et al. (2015). Our results strongly suggest that a new type of instability may be active, potentially interesting to improve constraints on dark matter models. The mechanism resembles a parametric instability Yang et al. (2015), but we find complex growth patterns in the field. The instability growth rate is larger for larger angular velocities, 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 lifetime? Hua

Acknowledgments. We thank Ana Carvalho, for useful input and for producing some of the figures for us. We thank David Hilditch for useful advice regarding the numerical implementation. We are indebted to Lorenzo Annulli, Emanuele Berti, Kfir Blum, Miguel Correia, Will East, Pedro Ferreira, Masashi Kimura, Luis Lehner, José Natário, Rodrigo Vicente, Leong Khim Wong and Huan Yang for useful comments on a version of this manuscript, and for discussions on several issues related to this work. T.I. and M.Z. are indebted to the Theory Division of CERN for the warm hospitality. M.Z. also thanks the hospitality of Perimeter Institute for Theoretical Physics, where part of this work was done. V.C. acknowledges financial support provided under the European Union’s H2020 ERC Consolidator Grant “Matter and strong-field gravity: New frontiers in Einstein’s theory” grant agreement no. MaGRaTh–646597. M.Z. acknowledges financial support provided by FCT/Portugal through the IF programme, grant IF/00729/2015. This research was supported in part by Perimeter Institute for Theoretical Physics. Research at Perimeter Institute is supported by the Government of Canada through the Department of Innovation, Science and Economic Development Canada and by the Province of Ontario through the Ministry of Economic Development, Job Creation and Trade. This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 690904. We acknowledge financial support provided by FCT/Portugal through grant PTDC/MAT-APL/30043/2017. The authors would like to acknowledge networking support by the GWverse COST Action CA16104, “Black holes, gravitational waves and fundamental physics.” Computations were performed on the “Baltasar Sete-Sois” cluster at IST and XC40 at YITP in Kyoto University.

Appendix A Simulations, convergence and radial dependence

Table 1 summarizes the parameters of our simulations, for the scattering of gaussian wavepackets off BHBs.

A.1 Numerical procedure and convergence analysis

To numerically evolve the scalar field equations in our prescribed metric background we employ the code presented in Refs. Cunha et al. (2017); Assumpcao et al. (2018), which makes use of the EinsteinToolkit infrastructure Löffler et al. (2012); EinsteinToolkit ; Zilhão and Löffler (2013) with the Carpet package Schnetter et al. (2004); Carpet 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 differencing techniques, which approximate the continuum solution of the problem with an error that depends polynomially on the grid spacing $h$,

[TABLE]

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 resolution $h_{c}$, $h_{m}$ and $h_{f}$. One can then compute the convergence factor given by

[TABLE]

To check the convergence of the extracted waveforms we have evolved the configuration of Fig. 2 in the main text, with resolutions $h_{c}=1.6~{}M$, $h_{m}=1.28~{}M$ and $h_{f}=1.0~{}M$ (where this refers to resolution of the outermost refinement level); the corresponding results are shown in Fig 4 for the $l=2$, $m=2$ multipole of $\Phi$. We have amplified the differences between the medium and fine resolution runs by the factor 2.97 expected for fourth-order convergence.

A.2 Dependence on extraction radii

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 finite-extraction radius artifact. Notice that the early-time response for the monopole does not align: the reason is that the initial pulse is ingoing. The late time behavior of the monopolar component on the other hand is perfectly aligned, showing that the pulse is outgoing at these late stages.

A.3 Dependence on initial data

Figure 6 shows the waveforms for different types of initial conditions, following the definitions in Table 1. The waveforms were aligned in time, to show a clear universal ringdown, the most important result of our work. The ringdown frequency and damping timescale is the same for different initial conditions, and depends only on the binary parameters (mass and separation). The results also indicate clearly that, although the period and damping timescale scale with separation $L$, the quality factor seems to be scale independent.

The initial pulse is spherically symmetric. However, on very short timescales the signal develops a quadrupolar ($l=2$) component as well. This behavior is expected, since one is specifying spherically symmetric initial conditions on a non-symmetric background. Thus, the field will very quickly sense the non-symmetric background metric. It is possible to show that such non-symmetric component is weaker at larger distances.

A.4 Power-law tails

Figure 7 zooms-in on the late-time behavior of the dominant, driving $l=m=0$ simulations for the scattering of gaussian wavepackets. Our results are consistent with a power-law.

Appendix B Majumdar-Papapetrou and the di-Hydrogen ionized molecule

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 Majumdar (1947); Papapetrou (1947); Hartle and Hawking (1972). In a cylindrical coordinate system the BHB version of the MP solution is written as

[TABLE]

with

[TABLE]

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=\pm a$ (hence, the parameter $a$ measures the distance between them). The spacetime ADM mass is $2M$.

Some geodesic properties of these solutions were studied before Chandrasekhar (1989); Assumpcao et al. (2018); Shipley and Dolan (2016). We will now show the rather remarkable result that the Klein-Gordon equation separates. Change to planar “prolate confocal elliptical” coordinates $\chi,\eta$ (keeping the time and azimuthal coordinates) defined as

[TABLE]

The variable $\chi$ plays a role similar to $r$ in standard spherical coordinates, while $\eta$ plays the role of $\cos\theta$. The domains of these variables are $-1\leq\eta\leq 1,\,1\leq\chi\leq\infty$. In these coordinates the Klein-Gordon equation for $\Psi e^{im\phi-i\omega t}$ can be written as

[TABLE]

When $a/M\gg 1$, we find that the above is separable and reduces to

[TABLE]

With the ansatz $\Psi=S(\eta)R(\chi)$, we find finally

[TABLE]

where $\Lambda$ is a separation constant. This same system describes the Schrodinger equation (for a positron) in the ionized Hydrogen molecule Burrau (1927); Wilson (1928). 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 effective-one-body treatment of post-Newtonian theory for the BHB problem, uses very explicit hydrogen-atom analogies to construct an inspiraling waveform model Buonanno and Damour (1999). 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.