Black hole binaries: ergoregions, photon surfaces, wave scattering, and quasinormal modes

Closed photon orbits around isolated black holes are related to important aspects of black hole physics, such as strong lensing, absorption cross section of null particles and the way that black holes relax through quasinormal ringing. When two black holes are present -- such as during the inspiral and merger events of interest for gravitational-wave detectors -- the concept of closed photon orbits still exists, but its properties are basically unknown. With these applications in mind, we study here the closed photon orbits of two different static black hole binaries. The first one is the Majumdar-Papapetrou geometry describing two extremal, charged black holes in equilibrium, while the second one is the double sink solution of fluid dynamics, which describes (in a curved-spacetime language) two"dumb"holes. For the latter solution, we also characterize its dynamical response to external perturbations, and study how it relates to the photon orbits. In addition, we compute the ergoregion of such spacetime and show that it does not coincide with the event horizon.


I. INTRODUCTION
High-frequency electromagnetic or gravitational waves propagating around extremely compact objects -such as black holes (BHs) -can travel in closed orbits. These can be either stable or unstable, according to whether or not small deviations from the orbit grow away from it as time evolves. In the case of a Schwarzschild BH (of mass M ), the union of all such orbits forms a sphere, the so called Schwarzschild photon sphere. Any inward-directed lightray emitted inside the photon sphere will eventually end in the singularity, while any outward-directed lightray emitted outside eventually reaches future null infinity. Hence the photon sphere is not a horizon (for which no lightray would ever be able to exit) but has a similar fundamental role.
Several BH phenomena are directly related to the presence of a photon sphere [1][2][3][4][5][6][7][8][9][10][11][12][13]. First, the characteristic modes of oscillation of a spherically symmetric BH can be related to the parameters of the null closed orbits [1-4, 14, 15]. In the eikonal limit, the oscillation frequency of such modes is a multiple of the orbital frequency while their characteristic damping time is related to the instability timescale of such orbits [2][3][4][5]. Second, the absorption cross section of a highly energetic scalar field scattered by the gravitational field of a spacetime endowed with a photon sphere is simply described by the orbital frequency and instability timescale of the associated null closed orbits [6]. Third, the extreme bending of light near the photon sphere also allows the formation of multiple (possibly infinite) images on the optical axis [7,8] -in addition to the primary image typically formed by any weak gravitational lens. In particular, strong gravitational lenses produce an infinite number of Einstein rings, and the astrophysical observation of a sequence of such rings would be another successful test of General Relativity in the strong field regime [8]. Last but not least, the presence of stable photon spheres implies the existence of a trapping region in spacetime where matter can accumulate, potentially leading to nonlinear instabilities [9][10][11]16].
The concept of a photon sphere, useful for Schwarzschild BHs, can be generalized to the concept of a photon surface in an arbitrary spacetime [17]. The definition of a photon surface, in simple terms, requires that any null geodesic tangent to the surface will remain everywhere tangent to it. More specifically, in this work our main interest is in closed null hypersurfaces, i.e. non-spacelike, chronal hypersurfaces S that admit a null geodesic (i) which is tangent to S and (ii) whose projection onto a given spatial cross-section of S forms a "closed orbit" on the cross-section 1 . The notion of a photon surface is independent of any symmetry that the spacetime may possess, and is applicable to the situation of two very compact objects that inspiral into each other and eventually merge, forming a single compact object. These events, notably the merger phase, being important sources of gravitational-waves, are of interest for detectors like LIGO and LISA. Photon surfaces and closed photon orbits around such binaries, although never stud-ied before, could provide important astrophysical information about BH-BH (and BH-neutron star) collisions, specially now that the gravitational-wave astronomy era has begun [18].
The main objective of this work is to provide the first step towards understanding closed photon orbits around binary systems of astrophysical compact objects. Due to the complexity of simulating and studying collisions of BHs and other compact objects, we focus on simpler geometries of BH binaries that are static and, therefore, never merge. We expect that our results are useful for low-frequency binaries whose velocities are much smaller than the speed of light. This work is also a first step towards understanding how other peculiarities -such as ergoregions -extend to dynamical spacetimes. There is a clear motivation for that, since it could lead to Penroselike phenomena or interesting superradiant effects [19].
The first toy model we consider is that of two electrically charged BHs in equilibrium [20][21][22] (the Majumdar-Papapetrou solution of Einstein's equations), for which we study null geodesics both on the symmetry plane between the BHs and on the plane containing both BHs. Our second model is based on analogue gravity, the fact that sound propagation in moving fluids is formally equivalent to a scalar field evolving in curved spacetime [23,24]. Points where the fluid velocity exceeds the local sound speed behave as ergoregions; if this happens on a normal to a certain surface, such surface is an analogue horizon. By considering a double sink hydrodynamical fluid flow (an exact solution of the fluid-dynamic equations) we are able to study the analogue of a binary BH system. We should note that extending the notion of ergoregions to dynamical curved spacetimes is nontrivial in general, but quite straightforward in acoustic setups. Our results may thus help in understanding whether phenomena such as superradiance exist in BH-binary spacetimes [19] We note that partial results concerning the structure of lightray or scalar-field propagation in the background of BH binaries were reported earlier. The evolution of scalar fields in the background of an inspiralling BH binary was investigated in Ref. [26]. The shadow of a static BH binary (unphysical, since it contains a strut holding the two BHs apart) was also recently reported [27,28]. However, unlike our analysis, these references do not address the question of the structure and timescales associated to the different possible closed photon surfaces. Therefore, through our work, we hope to pave the way for future investigations of closed null orbits in more realistic situations of interest for LIGO and LISA. 2 The existence of ergoregions in static spacetimes describing two Kerr BHs held apart by a strut was established in Ref. [25]. A generalization to dynamical regular spacetimes has not been done. An exact solution in General Relativity which describes two or more static BHs is known as the Majumdar-Papapetrou (MP) solution [20][21][22]. In a cylindrical coordinate system the two-BH version of the MP solution is written as with This solution represents two maximally charged BHs in equilibrium, each with mass M and charge Q = M . Here, and throughout this work, we use geometric units. In these coordinates, their horizons are shrunk to two points at z = ±a (hence, the parameter a measures the distance between them). The spacetime ADM mass is 2M . The Lagrangian for geodesic motion in the Majumdar-Papapetrou spacetime is given by where dots refer to derivatives with respect to an affine parameter. The constant δ on the right handside is zero (δ = 0) for null geodesics and one (δ = 1) for timelike geodesics. The two conserved quantities are E =ṫ/U 2 and L =φU 2 ρ 2 , which we identify as energy and angular momentum per unit rest mass, respectively.
A. Geodesics in the meridian φ = 0 plane We can fix the azimuth angle by setting φ = 0 (this is an arbitrary choice as any fixed angle yields the same behavior) and study geodesic motion confined to the meridian plane. For a generic separation between the BHs, the geodesic equations are not separable, as noted in Ref. [29]. However, when the two BHs coalesce, a = 0 and it is easy to solve the corresponding geodesic equations [2]. We find that there is a null geodesic at r = z 2 + ρ 2 = 2M .
With the aid of standard numerical integration methods, three types of closed null and closed timelike orbits were identified, as shown in Figs. 1 and 2: i. An orbit which encloses both horizons, shown in blue in both figures. For small separations, i.e. a M , the coordinate radius of such an orbit is 2M , as expected from the analysis of the a = 0 case. The corresponding coordinate period is ∼ 16πM . At large separations, we find that the period of null outer orbits is T ∼ 20πM + 4a, for values of M in the range 1 M 10.
ii. Orbits which enclose only one horizon, shown in red in both figures. At large separations, such geodesics only "see" the gravitational field of one BH. Indeed, we find that then the orbit is nearly circular with radius ∼ M , and a coordinate period ∼ 8πM , as expected from the above.
iii. An "8-shaped" orbit which encloses each of the horizons once, shown in black in both figures.
For the two-BH MP geometry these three types of orbits seem to exist for any separation a. These results are in complete agreement with both Chandrasekhar's and Shipley and Dolan's recent results [29,30]. A full description of the equations and normalization conditions is given in Appendix A. We found, numerically, that all these null geodesics are unstable: the numerical integration requires fine tuning to follow them. On the other hand, timelike geodesics are more stable and complete a larger number of full orbits before collapsing to the singularity (some "8-shaped" trajectories are actually longterm stable). For highly relativistic particles, timelike geodesics approach null geodesics. It is impossible to find closed timelike orbits for arbitrary energies. Half way between the two BHs, there is a symmetry plane (z = 0) in which the geodesic equations are separable. By setting a = 1, which simply amounts to rescaling units as ρ → ρ/a and M → M/a, Eq. (2) becomes Substituting E and L and setting δ = 0 in the normalization condition, the equation for null geodesics becomeṡ Provided that M 2 > 27/8, the effective potential V eff has two critical points, which correspond to a stable (in- 0367020. The red curve shows the portion of the motion between two consecutive apastrons, for which the accumulated angle is given by ∆φ = 2πq. When q is a rational number, the motion is periodic, and the geometric properties of the plot can be extracted from q. From left to right, q is given by (1 + 7/10), (1 + 33/47), and (1 + 45/64). Notice that q increases monotonically with E. Geometrically this means that the motion is divided into more "leaves", which creates more vertices in the figures. The stable circular orbit has a radius ρin = 1.167a, and the unstable orbit is at ρout = 9.740a. ner) and an unstable (outer) closed orbit. 3 The analytic expressions for these two radii are [31] ρ out = 4M 2 9 1 + 2 cos When M a we have ρ out → 2M and ρ in → 0, so that the single hole solution is recovered. Besides these circular orbits, more general closed photon orbits also exist in the symmetry plane, as shown in Fig. 4. There is, however, no apparent connection between these closed orbits and the ones found in the meridian plane.
For timelike geodesics, the number of roots of the equation dV eff /dρ = 0 (and consequently the number of circular orbits) depends not only on the mass of the BHs but also on the angular momentum of the particle due to an extra term in the effective potential [32]. The Lyapunov exponent λ associated with the unstable circular orbits can be computed in a straightforward way [2,5]. A detailed calculation for null geodesics is given in Appendix A. In our context, the Lyapunov exponent λ describes the time scale in which a deviation δρ from a circular orbit grows exponentially, i.e. λ is characterized by δρ ∼ e λt . The Lyapunov exponent attains its maximum value 0.03158 at M = 2.197, and approaches 0 for very large M . For instance, at M = 10 6 , the exponent is 8.38 × 10 −8 . Besides the circular orbits, more general closed timelike orbits also exist, as shown in Fig. 5.
Closed orbits in the symmetry plane exhibit a "zoomwhirl" behavior, as seen in Figs. 4 and 5. A brief remark concerning these orbits is in order. Applying the formalism presented in Ref. [33], we associate geometric properties of symmetric orbits in the plane z = 0 with the accumulated angle between two consecutive apastrons ∆φ given by where ρ p and ρ a are the radial coordinates of the periastron and apastron respectively. We may look for orbits for which ∆φ = 2πq, where q is a rational number. We can always write any non-integer rational number uniquely as q = w + v/z, where w = q ∈ Z, and v < z are coprime positive integers. The whirl number w gives the number of times the orbit whirls between two consecutive apastrons and its contribution to the accumulated angle is evidently 2πw. The number z gives the number of "leaves" drawn by the solution before completing a full closed orbit. The path described by the geodesic will not necessarily follow consecutive "leaves" throughout the motion, so the number v tells how many vertices the orbit will skip before drawing the next "leaf".

III. THE DOUBLE SINK SOLUTION
Our second toy model for a BH binary is inspired in Analogue Gravity, a formal equivalence between the propagation of sound waves in fluids and a scalar field in a (analogue) curved spacetime [24,34,35]. In particular, a sound wave Φ associated with a velocity field v = ∇Φ is governed by the equation The curved spacetime geometry associated with the D'Alembertian is fixed by the background fluid flow, and contains horizons if and when the fluid flow exceeds, at some point, the local sound speed [24,34,35]. We consider a generalization of the 2 + 1 "draining bathtub" geometry [35], which here describes two acoustic "dumb" holes. Our background flow is two-dimensional, consists on two sinks of unit strength at x = ±a, and can be written as [36], with r 2 = x 2 + y 2 . The constant A > 0 is arbitrary and fixes the fluid speed at some radius. For simplicity, we assume A = 1 throughout the paper. Some notions of curved spacetime geometry can be used to formally describe the propagation of sound waves, with velocity c s , in irrotational fluids like the one determined by Eq. (9). Among these, we find useful the following definitions. Ergocurves are curves at which v 2 = c 2 s , and play the role of spacetime ergosurfaces. Acoustic horizons, on the other hand, are curves at which n · v = c s , where n is the normal vector at each point in the curve. Acoustic horizons play the role of event horizons when the acoustic geometry is stationary. Finally, soundcurves (which are analog light rings or analog photon surfaces) are closed curves for sound waves; a sound wave of velocity v s ≡ (dx/ds, dy/ds) (with s an affine parameter) therefore satisfies, on a soundcurve, ||v s − v|| = c s . These analog light rings were recently investigated experimentally in a rotating vortex representing the draining bathtub geometry [37].
These properties can all be recovered using the effective metric experienced by sound waves (for simplicity we set the background density to unity), When there is no separation between the sinks (a = 0), the metric above simplifies to Transforming t and r into new variablest andr, defined by r =r/c s and dt = dt/c 2 s + 2rdr/ c 2 s (r 2 − 4) , brings the acoustic metric above to the canonical form [38] When the two holes are sufficiently close to one another, the metric reduces to (12), and the soundcurves are located at r = 2 √ 2 [2] (here and in all numerical calculations we assume c s = 1). The coordinate time that sound takes to travel across this single "merged sink" is T = 8π (corresponding to a frequency Ω = 1/2), and the associated Lyapunov exponent is λ = 1/(2 √ 2). For double sinks, we find that there are different sound curves, shown in Fig. 6. Global outer orbits (shown in blue) surrounding both BHs exist for any value of separation a. Inner orbits around each sink (shown in red) and "8-shaped" orbits (shown in black) exist only for values of separation a > 1. It is worth noticing that these orbits are not symmetric, as opposed to the ones found in the Majumdar-Papapetrou spacetime. Furthermore, the global outer orbits can cross the other two types of orbits for small values of a, in contrast to what happens in MP spacetimes.
For double sinks separated by a large distance a 1, the coordinate period of the global geodesic (encircling both holes) is T ∼ 15.6 ± 0.3 + 4a. This is in rough agreement with expectations: the total crossing time is of order 4a plus two semi-circles, each of which takes 2π to complete. In the same regime, the coordinate period of the inner closed lines is ∼ 12.5 ∼ 4π, as expected from the single-hole analysis.
For separations a 5 the travel time does not strongly depend on the separation itself, and is ∼ 8π. Thus, depending on the dynamical aspect under investigation, the two sinks can be considered to be merged at these separations. Defining φ = y (x), the normal to the curve y(x) can be written as

Acoustic horizon
It is straightforward to show that the equation n · v = c s determining the horizon becomes a first order differential equation for y(x): where the background velocities are given by (9). Since any acoustic horizon engulfing both sinks must necessarily intercept the symmetry line (y axis), we look for solutions of (13) that pass through x = 0. We find that a common horizon engulfs both sinks only if a ≤ 1. When a > 1, the sinks are sufficiently apart and we recover a single horizon around each sink. For a ≤ 1, by searching the parameter space of the initial conditions of Eq. (13), we have found that the acoustic horizon is always a smooth curve with y(0) = ± 1 + √ 1 − a 2 and y (0) = 0. In Fig. 7 we plot the shape of the horizon at the critical separation a = 1.

Ergoregion
Although a notion and definition of ergoregion for dynamical spacetimes does not exist, the acoustic version described above is readily extended to any geometry. Solving the equation v 2 = c 2 s for the double sink flow given by Eq. (9) yields The function y(x) informs on the location of the ergocurve. In general, the acoustic horizon and the ergocurve do not coincide. To show this, assume that they do and that v 2 = c 2 s everywhere along the horizon surface. At any given point x where y(x) is a smooth function, and v x and v y are non-zero, the horizon equation (13) implies that which clearly contradicts the assumption that v x = 0 and v y = 0. Consequently, the ergocurve and the acoustic horizon can only coincide at points for which either v x = 0 or v y = 0 or y(x) is not smooth. This fact is illustrated in Fig. 7, which exhibits four points where they coincide. In other words, for acoustic binaries the horizon lies inside the ergocurve, raising the interesting prospect of phenomena such as superradiance [19,39] to occur. These results suggest that similar phenomena may be present in gravitational BH binaries which may lead to new effects during the inspiral.

B. Wave scattering and quasinormal modes
There is an established connection between photon surfaces and the dynamical response of perturbed, single-BH spacetimes [2,13,40,41]. In a nutshell, the photon surface works as a trapping region where scalar, electromagnetic or gravitational perturbations can linger. The instability timescale of null geodesics is then related directly to the lifetime of massless perturbations.
We wish to understand if such connection carries over to binaries. The analysis of wave dynamics in BH binaries is notoriously difficult. For single BH spacetimes with certain symmetries, a linearized study where the relevant PDEs are transformed into ODEs is possible. Sound waves Φ(t, r, φ) ∼ Φ(t, r)e imφ in single acoustic dumb-hole spacetimes were studied years ago, and it was concluded that [38,42]: i. The early-time response is strongly dependent on the initial conditions, but is followed by quasinormal ringdown [4]. In the eikonal regime, such ringdown is described by sound waves circling the soundcurve and decaying away. For the geometry (12), the fundamental (longest-lived) frequencies ω m for the first few azimuthal numbers m are [38,42] 4 (see [43] for a ringdown analysis in presence of vorticity): ii. The ringdown stage eventually gives way to a latetime power law tail, where the field decays as Appealing to the correspondence between quasinormal modes and geodesic properties [2,44], one expects a relation of the type to hold at large m. This relation indeed describes extremely well the results above for the quasinormal modes of acoustic holes 5 . Nevertheless, as shown in Ref. [44] for the draining bathtub geometry, expansions of ω m can be obtained to arbitrary order in powers of m.

Numerical approach
For the numerical evolutions of (8), we used the Ein-steinToolkit infrastructure [45][46][47], which is based on the Cactus Computational Toolkit [48], a software framework for high-performance computing. The spacetime metric is fixed according to Eq. (10), and we evolve the Klein-Gordon equation on top of this background using the same numerical code as the one employed in [49], herein accordingly adapted to 2+1 evolutions. To reduce the system to a first order form, we introduce the "canonical momentum" of the scalar field Φ where L denotes the Lie derivative, to write our evolution system in the form where D i is the covariant derivative with respect to the 2-metric γ ij and K is the trace of the extrinsic curvature K ij = 1 2α L β γ ij . To numerically evolve this system, 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. A complication arising from metric (10) is the presence of curvature singularities at x = ±a. To deal with these we use a simple excision procedure, whereby our evolved quantities are multiplied by a mask function M (x, y) where u = Φ, K Φ and M (x, y) is a smooth transition function that evaluates to 1 outside the horizon and 0 inside. 6 We evolve the system with initial conditions given by The constants c m characterize the multipole of order m. The function Φ(x, y) represents a Gaussian pulse of width σ centered at (x 0 , y 0 ) and peaked at R = R 0 . We extract the scalar field Φ at x = 30, y = 0 and at x = 0, y = 30. The first point lies along the symmetry line; the second extraction point lies along the line joining the two sinks, and is closer to one than to the other. The general features of our results remain the same at other extraction points. Figure 8 shows the result of evolving a Gaussian, quadrupolar initial data, in the background of a single sink (a = 0). The signal is exponentially damped at intermediate times, and we find excellent agreement (to within ∼ 2%) with prediction (16b), and hence also with the soundcurve analysis. The late-time tail, also clear in Fig. 8, is well described by Φ ∼ t −4.9 , in good agreement with prediction (17). Scattering of a sound wave off a binary dumb-hole with a = 1. The Gaussian, centered at the origin, has a width σ = 1/2 and is localized at R0 = 5. The field Φ is extracted along the x-axis (black solid line) and along the y-axis (red, dashed line) at a distance of 30 from the origin. Top: Monopolar wave, with c0 = 1. For monopolar waves the prompt signal quickly gives way to a power-law tail Φ ∼ t −1.1 at late times, in agreement with prediction (17). Bottom: Quadrupolar wave, with c2 = 1. We find a ringdown ω = 0.472 − 0.170i and a power law ∼ t −1.2 . Figure 9 shows the result of evolving monopolar and quadrupolar waves in the background of a binary with finite but small separation, a = 1. We find that the signal still has all the features of the single-BH spacetime: a single ringdown stage, followed by a late-time power law tail. Both the ringdown and the tail are well-described by single BH geometries. In particular, due to mode-mixing caused by the lack of axisymmetry, an initial pure-m data will evolve to a superposition of different m's. In particular, the m = 0 mode seems to be present at late times in the waveform and causes the latetime power-law tail to be well described by Φ ∼ 1/t (cf. Eq. (17)). In other words, binaries separated by a distance a 1 behave very similarly to a single BH. This is not completely surprising and is most likely connected to such binaries having a common horizon (cf. Fig. 7). This feature is also present in nonlinear evolutions of Einstein equations for the BH-binary problem, where ringdown is seen to be triggered as soon as merger occurs, and in a very linear way [50][51][52][53][54].

Results
When the separation is made to increase, the different scales in the problem become noticeable, as seen in Fig. 10. For a = 5, the coordinate period of the global geodesic is of order ∼ 35, larger by about a factor three than the period of the geodesics encircling each BH.
Focus on the top panel of Fig. 10. The initial data is localized around only one of the BHs (located at x = 5). Therefore the initial data excites this BH first and to a higher extent than the second BH, farther apart. On the x-axis, the initial ringdown is well described by the modes of this BH. After a time ∼ 25 the ringdown mode of the perturbed rightmost BH reaches the extraction point x = 30 on the x-axis (black solid curve). Later, starting at t ∼ 55 the ringdown of the leftmost BH can be seen, followed by the decaying power-law Φ ∼ t −1 . The global mode is encoded in the envelope of this signal, and the only reason it is not more clearly visible is because the powerlaw tail sets in very quickly. Such behavior is similar to the setting in of global modes of anti-de Sitter spacetime, or of ultracompact exotic objects [12,13,55,56].
The two lower panels in Fig. 10 show the result of evolving a monopolar and quadrupolar data, centered at the origin and localized away from the two holes. The signal is peaked at some intermediate time, most likely the result of interference between different modes. We do not have a simple and elegant explanation for such feature. For any initial data we looked at, late-time tail is always dominated by t −1 , that of the circularly symmetric mode.

IV. ON THE STRUCTURE OF CIRCULAR NULL GEODESICS FOR ETERNAL AND FOR COALESCING BINARIES
One may think that even though our geometries of BH binaries described by the MP metric are stationary, they can be interpreted as a "snapshot" of a dynamically coalescing process of two BHs, and therefore the sequences of the circular null geodesics obtained in this paper could mimic a photon surface around a dynamical geometry of colliding BHs in a more astrophysically realistic situation. However, there is a sharp contrast between the stationary MP binary and a dynamically colliding BH binary. One should first note that the BH binaries in the MP solution never merge to form a single BH; they have two mutually disconnected components of the event horizon and therefore allow for closed null/timelike orbits outside the event horizon (see left panel of Fig. 11). In contrast, when two BHs dynamically coalesce, the event horizon as a three-dimensional null hypersurface has only a single connected component from the beginning.
Let us recall that the future event horizon of a BH is formed by null geodesic generators which have no end points in the future but can in general admit past end points. For stationary BHs, the event horizon is a Killing horizon and the horizon null geodesic generators have no past end points. As for dynamical BH system, such as those formed by gravitational collapse or by the merger of two BHs, (some of) the horizon null geodesic generators have past end points, which form a subset in the spacetime. An important feature is that such a pastend-point set-called crease set-must be spacelike and achronal rather than timelike or null, since the event horizon itself must be an achronal null hypersurface. Now when we say, "a single" BH or "two" BHs, we refer to as the number of the cross-section components of the event horizon and a Cauchy surface. It has been shown in [57] that this notion of the "number of cross-sections" is in fact time-slice dependent when the event horizon admits a crease set as a part, since in that case one can choose a Cauchy slice such that it crosses the crease-set as many times as one wishes (middle panel of Fig. 11). Furthermore, depending upon the choice of Cauchy slicing, even the topology of the horizon cross-section changes (for instance, when the crease set is 2-dimensional, it can even be temporarily toroidal [58][59][60]). For the mathematical structure of the crease set and the classification of the possible (temporal) horizon cross-section topologies, see Ref. [57]. For the dynamically coalescing "two" BH binary, the temporarily "two BHs" are connected by a crease set, and the achronal nature of this crease set prevents any causal geodesic curve (either null or timelike) from forming a closed orbit around either "each single" BH or the "8-shaped" closed orbit that would encompass the "two BHs." An exception is a global outer closed orbit. Therefore one has to be careful when attempting to make a physical interpretation of closed null orbits found in static MP solution.
The observation above would be manifest if, for example, one examines closed null geodesics in Kastor-Traschen (KT) solution [61] with two BHs. The KT solution also describes configuration of multiple BHs with maximal charge, i.e., M = |Q|, but it includes a positive cosmological constant. Therefore those multiple BHs are on a de Sitter expanding universe. The time-reverse of the solution can describe coalescing BHs, and their event horizon must have crease sets. Since each BH has its own apparent horizon inside the event horizon, one can presumably find e.g., "8-shaped" closed null orbits that go around two "apparent horizons" but the orbits should entirely be inside the event horizon.
If we consider, instead of a coalescing BH binary, a coalescing compact star binary, each of the stars is compact enough to admit its own photon sphere, then we can find 8-shaped closed null/timelike orbits, besides the global null closed orbits. See right-most panel of Fig. 11.
The above notes of caution refer to the purely formal  11. Left: Null geodesic curves in two MP BHs, whose event horizon consists of two disconnected components. Depicted are two circular orbits around each of the disconnected components (blue lines) and a "8-shaped" orbits (green line) around the event horizon. Middle: Two BHs merge to form a single BH. The event horizon has only a single connected component and admits a "crease set" (depicted by thick red-line) consisting of the past end points of some of the null geodesic generators of the event horizon. Here "two black holes" is meant to be the number of cross-sections of the event horizon on a particularly chosen Cauchy time slice as illustrated. Note however that since the event horizon has, as a part, the crease set, which is spacelike, the notion of the number of cross-sections (as well as the topology of the cross-section) is dependent on the choice of Cauchy time slice. Also because of the spacelike nature of the crease-set, null/timelike orbits cannot be of "8-shaped." Inside the event horizon there are apparent horizons, for which null geodesic curves can possibly form circular orbits (blue lines). Right: Two ultra-compact stars collide to form a single BH. The system can admit a 8-shaped null orbit (green line) as well as circular null orbits (blue lines) around each of the ultra-compact stars if they are sufficiently compact. Further study is needed to see whether these 8-shaped and circular null orbits eventually can go around outside the event horizon (as illustrated here) or must go inside the event horizon, after the BH formation.
description of such geodesics and event horizons. Our detectors are placed some finite distance away from the black hole binary and are active for a short amount of time: they have no access to a teleological quantity such as eternal horizons. Instead, they measure properties related to local quantities, such as apparent horizons. Therefore, the geodesics that we describe should capture basically the physics at play in gravitational-wave detectors.

V. CONCLUSION
We have given the first steps in our understanding of wave dynamics around BH binaries. We found three types of null geodesics, which should play a role in the decay of perturbed binaries. Our results for sound waves in an acoustic geometry indeed indicate that the different timescales can be associated to such geodesics. There are many open issues, concerning the specific features of the wave dynamics, and how it generalizes when the binary itself has dynamics. In particular, defining and exploring the notion of closed null orbits and ergoregions for dynamical spacetimes is challenging, but potentially important, as it raises the possibility of gravitational lensing effects, superradiant effects and Penrose energy extraction.

Integration of geodesic equations in the MP metric
The geodesic equations are found by varying the action with respect to the Lagrangian in Eq. (3). The equations of motion are the well know Euler-Lagrange equations: where x µ = {t, ρ, φ, z}, and λ is an affine parameter. The two conserved quantities are E =ṫ/U 2 and L = ρ 2 U 2φ , which we identify as energy and angular momentum. We can find an energy conservation equation by substituting E and L into Eq. (3): Here and throughout the appendix, a dot means derivative with respect to the affine parameter λ. In terms of E and L, the Euler-Lagrange equations reduce to the following system of coupled differential equations for ρ and z:ρ To find symmetric orbits we set initial conditions as ρ(0) = 0, z(0) = z 0 , andż(0) = 0. By virtue of Eqs. (A2) and (A3) we finḋ With the above initial conditions, we are able to find various orbits by choosing values of E, L, and z 0 . Setting L = 0 means constraining the motion to the plane that contains both holes, as explained in the main text. On the other hand, to look at the symmetry plane z = 0 we need to choose non-zero values for ρ(0) and L.

Oscillating particle in the MP metric
We start with a special case of the three-body problem in which one the bodies is much lighter than the other two. The small body is a test particle following a timelike geodesic motion at z = 0, while the two massive bodies are the two extremely charged BHs in the Majumdar-Papapetrou metric. Setting δ = 1 in Eq. (3) yieldṡ where U = U (ρ, z = 0). The effective potential 1/U 2 tends to 1 as ρ → ±∞ and z = 0. Thus, the solution is bounded and oscillatory around ρ = 0 for 1/(1 + 2M ) 2 < E 2 < 1. The first inequality comes from the fact that the initial velocity must be real. The motion can be parametrized by the coordinate time t by dividing both sides of Eq. (A7) byṫ 2 : For small-amplitude oscillations, Thus, to lowest order E 2 = 1/(1 + 2M ) 2 and the motion is sinusoidal with (coordinate) frequency with V 0 = 1/(1 + 2M ) 2 . For large separations (in our units, when M → 0) we recover the Newtonian result.

Lyapunov exponent
One can compute the Lyapunov exponent associated with the unstable circular orbits by considering perturbations in the radial coordinate ρ and in the canonical momentum p ρ = ∂L/∂ρ [2,5]. Let p ρ = δp ρ and ρ = ρ out + δρ, then the equation for the perturbation reads d dt where The Lyapunov exponent λ is the eigenvalue of the matrix in Eq. (A11), namely √ K 1 K 2 . It is worth noting that the Lyapunov exponent depends on the choice of parametrization. For that reason the perturbed equation is written as an ODE with respect to the coordinate time t, so that associated time scale is the one measured by an observer at infinity. By virtue of the Euler-Lagrange equations, Taking the derivative of Eq. (6) with respect to ρ yieldṡ where V = dV /dρ. From Eqs. (6) and (A14), ρ d dρ (g ρρρ ) = g ρρρ dρ dρ +ρ 2 g ρρ = − 1 2 For circular geodesicsρ = 0, which implies that V = E 2 and that the radius is a critical point of the potential. Eqs. (A13) and (A15) imply d dρ Combining Eqs. (A12) and (A16), the Lyapunov exponent is given by We recall that all quantities are evaluated at the unperturbed (circular) solution. With that in mind, we can arrive at an expression for λ computed for ρ = ρ out , the outermost and unstable circular null geodesic. The conditions V = E 2 and V = 0 yield By substituting these expressions in Eq. (A17) and using the definition of E, we find: It is possible to solve Eq. (A11) by applying the standard procedure to solve coupled differential equations with constant coefficients. The solution is: δρ(t) = δρ 0 cosh( K 1 K 2 t) + K 2 K 1 δp 0 sinh( K 1 K 2 t) , where δp 0 = δp ρ (0) and δρ 0 = δρ(0). We can compute λ numerically from the previous expressions. To simplify the procedure, we set δp 0 = 0 and find that where ρ num (t) is the solution of the radial equation found with a standard numerical integration technique.