What does a binary black hole merger look like?
We present a method of calculating the strong-field gravitational lensing caused by many analytic and numerical spacetimes. We use this procedure to calculate the distortion caused by isolated black holes and by numerically evolved black hole binaries. We produce both demonstrative images illustrating details of the spatial distortion and realistic images of collections of stars taking both lensing amplification and redshift into account. On large scales the lensing from inspiraling binaries resembles that of single black holes, but on small scales the resulting images show complex and in some cases self-similar structure across different angular scales.
pacs:95.30.Sf, 04.25.dg, 42.15.Dp
Black holes are the most compact gravitating objects in the universe, with such strong gravitational fields that not even light can escape them. In the vicinity of a black hole, light rays can be very strongly deflected from a straight-line path, sometimes orbiting around the black hole before continuing on their way. It is now well-known that the bending of light by massive objects like galaxy clusters can create brightness amplification , deformed images, or even multiple images  of background objects such as quasars. These signatures have so far only been directly observed in cases where the deflection of light is very slight, up to approximately arc seconds [3; 4]. However, here we are interested in the lensing effects associated with much more extreme bending of light near single or binary black holes, where the deflection angle is unbounded.
The lensing effects near general-relativistic bodies were first studied in the 1970s, with Cunningham and Bardeen  looking at a star on an orbit in a Kerr spacetime, and Luminet  studying an accretion disk around a Schwarzschild black hole. More recently, open-source codes such as GYOTO  and GeoViS  have produced images of lensing in the neighborhood of various compact objects. While the lensing caused by an isolated black hole has been understood analytically, the case of lensing by a binary black hole (BBH) is much more challenging because of the difficulty of solving for the geometry of the spacetime. With some arguably unrealistic assumptions (e.g., two maximally charged black holes in static equilibrium), analytic solutions can be found and subsequently used for lensing [9; 10; 11; 12; 13; 8].
For astrophysically relevant binaries, however, we must instead rely on numerical solutions. Solving these binary spacetimes numerically to high accuracy has been possible for the last decade (see [15; 16] for a review), motivated by the need to provide gravitational-wave templates used by experiments such as LIGO, VIRGO, and KAGRA to make detections. By using the spacetimes computed in such simulations, we gain the ability to solve for the lensing effects in BBH systems.
In this paper, we focus on the question of what an observer in the vicinity of a BBH would actually see as the black holes orbit, spiral inward, and merge, with an example shown in figure 1. This is in contrast to most BBH visualizations, in which the positions or horizons of the two black holes are simply shown as a function of time in some coordinate system. We instead compute the paths of light rays that enter the observer’s eye or camera to find what would actually be seen. Furthermore, this path must be computed in the fully time-dependent spacetime, as the orbital velocities for a black-hole binary are typically large enough that the system cannot be approximated as time-independent during the time taken by the photons to travel across it.
Because the black holes themselves do not emit light (we ignore Hawking radiation, which is significant only for microscopic black holes), the observer would see nothing unless there is some additional light source. For illustrative purposes, we will take an artificial background “painted on” at infinity (figure 3) as the light source for most of our examples; this will allow us to study in detail where each light ray originates.
We set up the problem with our black hole(s) near the center of our chosen coordinate system. While any physical system representable by a spacetime metric can be used, we specialize in this paper to single and binary black holes. The observer (henceforth taken to be a camera) can be located anywhere in the space and is typically chosen to look towards the origin. A sphere with our light source encloses the black hole(s) and camera, infinitely far away.
To recreate the image taken by the camera in this configuration, we must find the properties of the light that arrives at each point on the camera’s image plane. A naïve approach would be to trace all possible light rays (i.e., null geodesics) emanating from the light source to determine which rays reach the camera and from what directions they arrive, but this is computationally infeasible. A more efficient approach is to reverse the problem by tracing light rays away from the camera and backwards in time (the computer graphics community calls this a ray-casting algorithm). This method identifies the origin of any light ray that illuminates the camera, from which we infer the color and intensity of the corresponding photons as detected by the camera. When black holes are present, some of the null geodesics traced backwards in time from the camera may approach arbitrarily close to an event horizon as ; these geodesics correspond to dark image regions.
In what follows we describe how the light rays are traced from the camera using the geodesic language from general relativity. We show how we initialize these geodesics based on camera parameters such as position and viewing angle. Finally, we show how the origin of each light ray is determined and describe how the simulated image is constructed.
ii.1 Geodesic tracing
Our code can trace geodesics independently through either numerical or analytic metric data. It is common for numerical simulations to use the decomposition , so we express the metric in the form
where is the lapse function, is the shift vector, and is the spatial metric.111 Our convention is that Greek indices, as in , denote temporal or spatial components, while Latin indices, as in , denote only spatial components. We obtain numerical data from simulations performed using the Spectral Einstein Code (SpEC) [18; 19; 20; 21; 22]. The geodesics are traced by evolving a solution to the geodesic equation
where is the four-position of the geodesic, is an affine parameter, and are the Christoffel symbols describing the effective force caused by spacetime curvature.
To facilitate the numerical geodesic evolution, we split this second-order differential equation into two first-order differential equations using an intermediate, momentum-like variable such as . As we have some freedom in the definition of this momentum variable, we look for one that helps to minimize computational time and numerical errors when evolving through spacetimes with black holes.
We initially explored using the variable from Hughes et al. , along with converting the evolution equations from affine parameter to the coordinate time of SpEC evolutions through the use of . Although the resulting evolution equations are concise and have no time derivatives of metric variables, the variables and grow exponentially near black hole horizons in typical coordinate systems used by SpEC simulations. This forces our time-stepper to take prohibitively small steps in order to achieve the desired accuracy.
We therefore choose a momentum variable slightly different than to mitigate this time-stepping problem. Null geodesics satisfy , which can be rewritten as using the metric (1). This expression shows that and scale similarly, so we can eliminate the exponential behavior of these variables by evolving the ratio. Our intermediate variable thus becomes
where is the extrinsic curvature (see, e.g., ) and is defined via the inverse spatial metric as . Note that the geodesic equation consists of four second-order equations, yet we only have three pairs of coupled first-order equations in (4). Because we are evolving a normalized momentum (3), we have lost information about during evolution. Compared to Hughes et al. , we have introduced a time derivative of the three-metric inside , but we have significantly sped up the evolution near black holes by removing the exponential growth of and .
The equations in (4) are similar to those in (28) of Vincent et al. . In fact our intermediate evolution variable is related to their variable by the three-metric, such that . But our (4) has a reduced number of both temporal and spatial derivatives of metric quantities compared to Vincent’s (28).
During the backwards-in-time geodesic evolution, many geodesics are traced until they are far from the strong-field region, but some are traced until they encounter a black hole. These latter geodesics slowly converge towards the black hole’s event horizon, but as they can in principle be evolved indefinitely, we need some way of identifying them in finite time. We do this by monitoring for each geodesic, which (as discussed above) grows large near black hole horizons. Since our evolution equations (4) do not evolve , we must evolve another equation to keep track of it. However, we would still like to avoid the exponential growth of near the horizon. This can be accomplished by evolving the logarithm of . As was done in (3), we multiply by the lapse to reduce the number of terms in the resulting equation, which gives the evolution variable . This leads to the evolution equation
When becomes too large, signaling a large energy, we flag the geodesic as originating from the black hole and we stop evolving it.
The remaining geodesics are those that originate from infinity, so we need to determine the (,) location at infinity where they come from. In section II.3, we will need the gravitational redshift of each photon, which can be calculated from the ratio of the photon’s energy at the two ends of its trajectory via
where is the photon’s energy at infinity, and is the photon’s energy as measured by the camera. Therefore we will need to compute the energy that each photon would have at infinity. In practice, these geodesics are traced backwards in time until they reach a large distance from the black hole(s), chosen so that the metric at is equal to the flat space metric within about a percent error. We use the approximation that the metric is exactly flat for . Under this approximation, the geodesic’s direction and at infinity are the same as at . The direction is used to calculate a (,) location on the sky, while is the photon’s energy at infinity, .
ii.2 Initial data
Here we outline how we initialize our geodesic evolution variables. Because the geodesics are traced away from the camera, backwards in time, we initialize each geodesic’s evolution variables to their values at the camera. We have seven variables to set: three each for the initial position and momentum in (4), and one for the initial redshift in (5).
The initial position for every geodesic is simply the camera’s position. The initial momentum, however, is different for each geodesic and is dependent on the angle at which it enters the camera. We express the momenta in terms of an orthonormal tetrad defined as
The camera’s four-velocity, a timelike vector. For stationary cameras ;
The direction in which the camera is pointing;
The “upward” direction for the camera;
The “rightward” direction for the camera.
The four-vectors , , and are all spacelike, and their orientations in the camera’s reference frame are illustrated in figure 2.
In order to specify this tetrad, we give guesses for the vectors , , and , with the condition that the guessed time components of and must be zero. We then apply the Gram-Schmidt process to the sequence , , and to transform these vectors into an orthonormal set. The final vector, , is found by calculating the generalized cross product of the other three; explicitly,
where is the Levi-Civita tensor (see [25, p. 202] for more details).
Given the four orthonormal unit vectors, we can construct a null vector tangent to the geodesic that enters the camera from a given direction. The vector will be proportional to the four-momentum of a photon following the geodesic; that is, for some positive constant . We define by
where give the ray’s arrival direction in terms of fractions of the image’s horizontal and vertical lengths, respectively, and are the angular sizes of the camera aperture (field of view angles) in the vertical and horizontal directions. For the sign convention chosen in (8), corresponds to a photon seen at the bottom left corner of the image. We find by requiring that is null, i.e., :
We then use the metric to lower the index on , and we compute the initial value of our evolution variable using . Note that is independent of the proportionality constant relating and the actual photon momentum ; physically, this is because the photon trajectory is independent of the photon energy. The only place where enters is in the initial value of in (5). We fix the value of by demanding that the energy of the photon in the frame of the camera be unity when the photon strikes the camera, so in (6).
ii.3 Image generation
We create our image of the physical system by dividing the image plane into rectangular regions corresponding to the pixels of the output image and assigning an appropriate color to each region. Because each region has an extended size, there is no single source point we can look at to obtain its color, so we must adopt some prescription for assigning a single color to each pixel. We use two different prescriptions, based on the nature of the light source illuminating the system.
For extended sources, such as the artificial grid in figure 3, we use a subpixel sampling method. On each pixel we construct an evenly spaced grid of points, and at each of these points we determine where incident light rays originate, either from one of the holes or a location at infinity. We assign a color to each grid point based on that of the corresponding source point; the color of the pixel is then the average of these. We find that a grid of sample points gives sufficiently smooth images without too much computational cost. For these images, we neglect the effects of redshift and focus on the spatial distortions.
To create more astronomically relevant images, we wish to use a collection of point sources (i.e., stars) as our illumination. In this case we cannot determine a pixel’s color using sampling, but must instead sum the contributions from all the point sources contributing light there. For our list of sources, we use about stars from the Two Micron All Sky Survey (2MASS) . To simplify computations, we approximate each star as a thermal source with temperature and brightness determined by fitting to the photometric information in the catalog. When we calculate the contribution of each star to the light arriving at the camera, we must account not only for its properties as a light source, but also for the effects of the spacetime curvature encountered by the photon. These effects come in two forms. First, the observed energies of photons at the camera will be modified by redshift effects, changing sources’ apparent brightnesses and temperatures. Second, the spatial convergence or divergence of nearby geodesics produces an overall adjustment to each source’s apparent brightness without affecting its spectrum. Both of these effects are discussed in detail in Mollerach and Roulet . After we have drawn the entire image in this manner, we convolve it with a blurring function to make the stars more visible. This has the effect of transforming each star into a fuzzy circle with size dependent on its brightness.
The result of this scheme can be seen in figure 1, which shows the BBH image from figure 11 in front of a background of stars. Note that by generating our starfield images from a catalog of point sources, we obtain a substantially more realistic image than would be generated by applying the lensing deformation to a raster image of the unlensed Milky Way stars. In such a raster image, each star is usually represented (whether as a result of camera optics or software rendering) as a blurred circle whose area depends on the star’s brightness. These circles are typically hundreds of arc seconds wide, and therefore lensing distortions applied to the image tend to produce stars that appear as smeared ellipses. In contrast, the angular sizes of real stars are many orders of magnitude smaller, so we expect them to remain as unresolved points under all but the most extreme lensing magnifications. These unresolved points can then be rendered as previously described, giving stars that better portray what an observer would actually see (as in figure 1). The difference between these methods lies in the non-commutativity between the lensing deformations and the blurring of each star. A minor shortcoming of our method arises at Einstein rings (discussed in section III.1), where the magnification diverges. There a star could in principle (though with very low probability) appear as an extended object, but in our treatment it would remain point-like. On the other hand, blurring first and then lensing is almost guaranteed to produce unphysical extended streaks at the Einstein ring.
Before applying our lensing code to binary black hole systems, we generate images of simpler analytic spacetimes. These serve both to provide checks that our images are consistent with earlier work, and also to illustrate general features of lensing around black holes that will appear again in BBH images. We then proceed to show two different configurations of BBH mergers.
To help visualize the lensing, we divide our light source at infinity into colored quadrants with a superimposed grid. An external view of this sphere is shown in figure 3. In addition to the colored sections, our light source has a bright reference spot in the direction towards which we point our camera. This spot will prove useful in illustrating an important feature of black hole lensing called an Einstein ring.
iii.1 Analytic spacetimes
In figure 4, we compare a flat space image with the images obtained by lensing our light source through Schwarzschild and Kerr black hole spacetimes. The top row from left to right shows flat Minkowski space and a Schwarzschild black hole. These spacetimes are spherically symmetric, so viewing them from different angles produces the same lensing effects. The bottom row shows a Kerr black hole, where in the left frame the spin vector is pointing out of the page and in the right frame it is pointing up. Here the spin breaks the spherical symmetry of the spacetime, leading to different lensing effects from different viewing directions.
In Minkowski space in the top left image we expect no deflection of light, which is what we observe. The camera sees an upright image of the portion of the grid near the white dot. The bowing of the grid lines is an expected geometric effect of viewing a latitude-longitude grid.
In the top right image, we see the lensing effects of a non-spinning black hole. The black circle in the center of the image is called the shadow of the black hole, where the hole prevents any light from reaching the camera. Alternatively, a shadow is a region of the image where geodesics are traced backwards in time from the camera to a black hole. Another prominent feature is that the white dot on our grid at infinity has been lensed into a large ring, called an Einstein ring . Light from the point situated directly on the opposite side of the black hole, the antipodal point, will by symmetry be lensed into a ring around the black hole as observed by our camera. Regions inside the Einstein ring correspond to photons that are deflected by larger angles than are the Einstein ring photons; this results in an inverted image of the reference grid inside the Einstein ring. A second Einstein ring can be seen near the shadow, corresponding to light from a source behind the camera wrapping around the hole on its way to the camera. In fact, photons can wind an arbitrarily large number of times around the black hole, resulting in an infinite number of Einstein rings.
The bottom row of figure 4 shows a single black hole with a large dimensionless spin of . As in the Schwarzschild case, there is an Einstein ring around the black hole shadow as well as image inversion inside the Einstein ring. However, for the case of a Kerr spacetime, the light coming from the Einstein ring does not originate from a single point directly behind the black hole, but from a small region (unless the camera is pointing directly along the spin axis). The spin of the black hole causes frame dragging, where space is dragged in the direction of the rotation [29; 30]. In the bottom left image, the spin axis of the black hole is pointing out of the page, so space is dragged in a counterclockwise motion. The effect of the frame dragging on the photon trajectories produces an image in which the grid itself appears to be dragged by the spin, as is evident when compared to the non-spinning black hole in the top right image. The strength of frame dragging increases closer to the black hole, which can also be inferred from the deformation of the background grid.
Frame dragging manifests differently in the bottom right image, where the spin axis is pointing up. The direction of frame dragging is out of the page on the left of the shadow of the black hole and into the page on the right. A photon traveling in the direction of the frame dragging can orbit closer to the black hole without being captured than a photon traveling opposite the frame dragging direction, resulting in an asymmetrical shadow about the spin axis. This causes the shadow to appear offset relative to the shadow of a Schwarzschild hole.
iii.2 Binary black hole spacetimes
Astrophysical black hole binaries are expected to radiate energy via gravitational waves, leading to a long inspiral followed by a merger, and finally a ringdown to a steady-state single black hole. Lensing by a final, steady-state black hole will look like the single black holes already seen in figure 4. However, the situation becomes more interesting when viewing these systems before merger. The first images we will present show an equal-mass BBH with non-spinning black holes—one of the simplest binary inspiral spacetimes to analyze—shortly before merger. The simulation we use is case 1 of Taylor et al. .
Figure 5 shows the image of our reference grid in the presence of this BBH, where the camera is situated such that the orbital angular momentum is pointing out of the page. This image bears a striking resemblance to the bottom left frame of figure 4, excluding the details near the shadows. This shows that, away from the shadows, the spacetime looks similar to a single rotating black hole, where the lensing is dominated by the mass monopole with corrections caused by the angular momentum of the system. In the single-hole case, the spin is responsible for frame dragging, whereas here the orbital angular momentum is responsible.
Focusing on the inner portion of the image, we observe that the binary lensing is markedly different from the Schwarzschild or Kerr cases. Figure 6 shows a cropped version of figure 5, emphasizing the structure of the shadows. As might be expected, there are two prominent shadows visible, each associated with one of the two black holes. We also see a narrow secondary shadow (an “eyebrow” ) close to the outside of each primary shadow. These secondary shadows correspond to one black hole (BH) casting a shadow which is lensed by the other BH on the way to the camera. Equivalently, they are image regions where geodesics are traced backwards from the camera to a BH, but bend around at least one BH on the way there. The first pair of eyebrows is evident in figure 6; however, we can resolve a pair of smaller eyebrows, shown in the inset.
We show another view of the same system in figure 7. Here the camera is looking at the system edge on, such that the orbital angular momentum is pointing up. We see again an overall similarity with the corresponding orientation of Kerr spacetime (the bottom-right frame of figure 4), indicating the dominant effects of the mass and angular momentum in these images. We can see a primary shadow for each black hole, but in this configuration one black hole is located roughly behind the other and as a result its shadow gets lensed into a dark ring. Extending along the right side of this ring we see a long thin eyebrow, which is shown in the inset, along with another, smaller, eyebrow.
To illustrate how photon trajectories behave near shadows, we plot trajectories of a few geodesics on the horizontal line passing through the middle of figure 7 near the eyebrow. Figure 8 shows four snapshots of these trajectories in time, with their current locations in each frame denoted by large dots. It is easiest to consider these trajectories as we evolve them, out of the camera and backwards in time, to see where they came from. In frames A–C, we see the trajectories under consideration start close together then diverge significantly, demonstrating how nearby pixels on the image can correspond to vastly different physical locations. In frame D we see the entire trajectories. A few extend to infinity, but most terminate on the black holes; these are denoted by solid lines and dotted lines, respectively. Only the trajectories extending to infinity result in a photon reaching the camera; those that reach the hole on the right of frame D correspond to the primary ring-like shadow in figure 7, while those that reach the left hole correspond to the larger eyebrow visible on the right side of figure 7. Note that the black holes are orbiting rapidly, so they move significantly while the photons pass through the system.
We can also uniquely identify which black hole casts each shadow, which enables us to show in figure 9 the origin of the photons along the horizontal line across the center of figure 7. We arbitrarily label the large shadow in the middle of figure 7 as BH 2, and the ring-like shadow as BH 1. Regions where photons reach the camera from infinity are labeled . The top plot in figure 9 shows the origin of the photons that reach the camera along the entire middle horizontal line in figure 7. We see that each transition from to either of the BHs includes transitions to the other BH. Even though we cannot resolve them numerically, each vertical line in principle contains infinitely many transitions. To illustrate this idea, the second plot in figure 9 investigates the group of shadows indicated by the zoomed inset of figure 7. Here we find a structure which resembles the first plot. The third plot in figure 9 zooms to a similar group of shadows on the right side of the second plot to again reveal the same structure. This figure clearly shows evidence of self-similarity in the structure of BBH lensing, where the smaller length scales explore more photon orbits through the system. Furthermore, the structure of shadows in BBH lensing is more complex than figures 6 and 7 appear to suggest. The shadows these images focus on are merely some of the largest visible shadows, associated with simpler geodesic orbits around the binary.
If we consider this equal-mass BBH earlier in the inspiral when its separation is large, the black holes are only weakly interacting. Therefore most camera viewpoints of this binary will yield images with two primary shadows, one for each black hole. Each shadow will be similar to an isolated Schwarzschild or Kerr shadow but with the addition of small eyebrows. However, when the binary is viewed edge-on and the black holes are nearly aligned with the camera, we see an interesting image.
Figure 10 shows the equal-mass binary in this configuration, hundreds of orbits before merger. Just as in figure 7, the more distant black hole is lensed into a ring-like shadow; however, the ring is thinner here, primarily because of the large separation of the binary. The angular momentum causes the lensed grid outside the shadows to strongly resemble lensing by a Kerr black hole rather than lensing by a Schwarzschild black hole. In addition to the usual primary Einstein ring, another ring is visible between these shadows. Both of these rings correspond to the same source of light, which is in front of the camera and behind the BBH. The second Einstein ring is caused by photons following an “S”-shaped trajectory through the system.
The second binary system we consider is a fully generic black hole binary with a mass ratio of and black hole spins of and in arbitrary directions. This is case 4 of Taylor et al. . In figure 11 we see a top view of this system, in analogy with what is presented in figure 6. Away from the shadows, the lensing is similar to a single black hole with spin, as was seen with the equal-mass binary images. This appears to be a generic feature of lensing from orbiting BBHs. We can clearly see that the symmetry present in the equal-mass system is gone. The unequal masses evidently change the relative sizes of not only the primary shadows, but all additional shadows as well. The inset in figure 11 zooms to show two successively smaller eyebrows near the small black hole’s primary shadow. However, the effects of the black holes’ spins are not at all clear from this viewpoint.
In figure 12 we see the same binary as in figure 11, viewed with the orbital angular momentum pointing upward, in analogy with figure 7. We again see that, away from the shadows, the system looks like a Kerr black hole. The unequal mass ratio is apparent here, with the smaller black hole lensing the shadow of the larger black hole into a partial ring. If it were not for the black hole spins, the lensing by the binary would be symmetric, giving either a ring-like shadow similar to figure 7 or a shadow and a very thick eyebrow. In this particular BBH, the effect of the individual black hole spins on the image depends strongly on the camera position.
In this paper, we present the first images of gravitational lensing by astrophysically relevant binary black holes, thereby providing a realistic representation of what an observer near such a system would actually see. To accomplish this, we have developed a new set of equations that evolve photons efficiently near black hole horizons. Our images show there is a primary shadow—a region where the black hole prevents light from reaching the camera—for each black hole, as well as multiple secondary shadows (or eyebrows).
We have found that, early in the inspiral, images of a BBH look similar to two separate Kerr black hole shadows, unless viewed when the holes are nearly collinear with the camera. Shortly before the merger, all camera angles yield interesting images of not just one shadow for each black hole, but a handful of smaller visible shadows. We showed for an equal-mass binary viewed edge-on that the lensing structure exhibits self-similarity on smaller scales, corresponding to photons taking an increasing number of orbits through the system. Lensing by a fully generic BBH illustrated that the spin of black holes in a binary can have a clear effect on the lensed shadows.
We chose not to classify eyebrows and shadows into a hierarchy in this paper. In the inset of figure 6, for instance, identifying the largest eyebrow as the primary eyebrow and the next largest as the secondary eyebrow feels very natural, but the exact definition of such a hierarchy is not immediately clear. For example, simply specifying a geodesic winding number around each black hole is likely not to be sufficient. In addition to the trajectories not lying in a plane, the order that a geodesic orbits the black holes does not commute. Furthermore, the black holes are moving at comparable speeds to the geodesics. For these reasons, we leave the task of classifying shadows as future work.
We have also shown in this paper that, away from the shadows, an image of a binary black hole system looks like that of an isolated black hole. Thus it is necessary to resolve individual shadows in order to discern the unique visual characteristics present in such images, which places limits on our ability to observe them.
For systems involving matter, however, the combination of the lensing effects of strong gravity with the disruption and distortion of radiation-emitting matter might yield a unique optical signature. Generating lensed images of black hole-neutron star and neutron star-neutron star mergers is an avenue of future investigation. The techniques presented here would allow us to produce detailed visualizations of these mergers; integrating over such images, we could predict the optical signature of an unresolved system.
Acknowledgements.We would like to thank Curran Muhlberger for providing the temperature fits to the 2MASS photometric data. This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. We would like to thank Daniel Hemberger and Saul Teukolsky for comments on an earlier version of this paper. The authors from Cornell would also like to thank Saul Teukolsky and Lawrence Kidder for general advice while writing this paper. This work was supported in part by NSF Grants PHY-1306125 and AST-1333129 at Cornell University, by NSF Grants PHY-1440083, AST-1333520, PHY-1005655, and DMS-1065438 at the California Institute of Technology, and by a grant from the Sherman Fairchild Foundation. FH acknowledges support by the NSF Graduate Research Fellowship under Grant No. DGE-1144153. DB acknowledges support from the LIGO Laboratory, with funding from the National Science Foundation under cooperative agreement PHY-0757058 and NSF REU award PHY-1062293. The binary black hole simulations were performed using the Zwicky computer system operated by the Caltech Center for Advanced Computing Research and funded by NSF MRI No. PHY-0960291 and the Sherman Fairchild Foundation.
-  Sjur Refsdal. The gravitational lens effect. Monthly Notices of the Royal Astronomical Society, 128(4):295--306, 1964.
-  R. D. Blandford and R. Narayan. Cosmological applications of gravitational lensing. Annual Review of Astronomy and Astrophysics, 30(1):311–358, 1992.
-  Naohisa Inada, Masamune Oguri, Bartosz Pindor, Joseph F Hennawi, Kuenley Chiu, Wei Zheng, Shin-Ichi Ichikawa, Michael D Gregg, Robert H Becker, Yasushi Suto, et al. A gravitationally lensed quasar with quadruple images separated by 14.62 arcseconds. Nature, 426(6968):810–812, 2003.
-  Naohisa Inada, Masamune Oguri, Tomoki Morokuma, Mamoru Doi, Naoki Yasuda, Robert H. Becker, Gordon T. Richards, Christopher S. Kochanek, Issha Kayo, Kohki Konishi, Hiroyuki Utsunomiya, Min-Su Shin, Michael A. Strauss, Erin S. Sheldon, Donald G. York, Joseph F. Hennawi, Donald P. Schneider, Xinyu Dai, and Masataka Fukugita. SDSS J1029+2623: A gravitationally lensed quasar with an image separation of 225. The Astrophysical Journal Letters, 653(2):L97, 2006.
-  C. T. Cunningham and J. M. Bardeen. The Optical Appearance of a Star Orbiting an Extreme Kerr Black Hole. Astrophys. J. Lett., 173:L137, May 1972.
-  J.-P. Luminet. Image of a spherical black hole with thin accretion disk. Astron. Astrophys., 75:228–235, May 1979.
-  F.H. Vincent, T. Paumard, E. Gourgoulhon, and G. Perrin. GYOTO: a new general relativistic ray-tracing code. Class.Quant.Grav., 28:225011, 2011, 1109.4769.
-  Thomas Müller. GeoViS—Relativistic ray tracing in four-dimensional spacetimes. Computer Physics Communications, 185(8):2301 – 2308, 2014.
-  Sudhansu Datta Majumdar. A class of exact solutions of Einstein’s field equations. Phys. Rev., 72:390–398, Sep 1947.
-  Achilles Papapetrou. A static solution of the equations of the gravitational field for an arbitrary charge-distribution. Proc. Roy. Irish Acad., A 51:191–204, 1947.
-  David Kastor and Jennie Traschen. Cosmological multi-black-hole solutions. Phys. Rev. D, 47:5370–5375, Jun 1993.
-  Daisuke Nitta, Takeshi Chiba, and Naoshi Sugiyama. Shadows of colliding black holes. Phys. Rev. D, 84:063008, Sep 2011.
-  Akifumi Yumoto, Daisuke Nitta, Takeshi Chiba, and Naoshi Sugiyama. Shadows of multi-black holes: Analytic exploration. Phys. Rev. D, 86:103001, Nov 2012.
-  http://www.black-holes.org/lensing.
-  Joan Centrella, John G. Baker, Bernard J. Kelly, and James R. van Meter. Black-hole binaries, gravitational waves, and numerical relativity. Rev.Mod.Phys., 82:3069, 2010, 1010.5260.
-  Harald P. Pfeiffer. Numerical simulations of compact object binaries. Class.Quant.Grav., 29:124004, 2012, 1203.5166.
-  R. Arnowitt, S. Deser, and Charles W. Misner. The dynamics of general relativity. In L. Witten, editor, Gravitation: An Introduction to Current Research, pages 227–265. Wiley, New York, 1962, gr-qc/0405109.
-  http://www.black-holes.org/SpEC.html.
-  Simulating eXtreme Spacetimes. http://www.black-holes.org/.
-  Béla Szilágyi, Lee Lindblom, and Mark A. Scheel. Simulations of binary black hole mergers using spectral methods. Phys. Rev. D, 80:124010, 2009, 0909.3557.
-  Daniel A. Hemberger, Mark A. Scheel, Lawrence E. Kidder, Béla Szilágyi, Geoffrey Lovelace, Nicholas W. Taylor, and Saul A. Teukolsky. Dynamical excision boundaries in spectral evolutions of binary black hole spacetimes. Class. Quantum Grav., 30(11):115001, 2013, 1211.6079.
-  http://www.black-holes.org/waveforms.
-  S. A. Hughes, C. R. Keeton, P. Walker, K. T. Walsh, S. L. Shapiro, and S. A. Teukolsky. Finding black holes in numerical spacetimes. Phys. Rev. D, 49:4004, 1994.
-  F.H. Vincent, E. Gourgoulhon, and J. Novak. 3+1 geodesic equation and images in numerical spacetimes. Classical and Quantum Gravity, 29(24):245005, 2012.
-  Charles W. Misner, Kip S. Thorne, and John Archibald Wheeler. Gravitation. Freeman, New York, New York, 1973.
-  M. F. Skrutskie, R. M. Cutri, R. Stiening, M. D. Weinberg, S. Schneider, J. M. Carpenter, C. Beichman, R. Capps, T. Chester, J. Elias, J. Huchra, J. Liebert, C. Lonsdale, D. G. Monet, S. Price, P. Seitzer, T. Jarrett, J. D. Kirkpatrick, J. E. Gizis, E. Howard, T. Evans, J. Fowler, L. Fullmer, R. Hurt, R. Light, E. L. Kopan, K. A. Marsh, H. L. McCallon, R. Tam, S. Van Dyk, and S. Wheelock. The two micron all sky survey (2MASS). The Astronomical Journal, 131(2):1163, 2006.
-  Silvia Mollerach and Esteban Roulet. Gravitational Lensing and Microlensing. World Scientific, 2002.
-  Albert Einstein. Lens-like action of a star by the deviation of light in the gravitational field. Science, 84(2188):506–507, 1936.
-  H. Thirring. Über die Wirkung rotierender ferner Massen in der Einsteinschen Gravitationstheorie. Physikalische Zeitschrift, 19:33, 1918.
-  H. Thirring. Berichtigung zu meiner Arbeit: “Über die Wirkung rotierender Massen in der Einsteinschen Gravitationstheorie”. Physikalische Zeitschrift, 22:29, 1921.
-  Nicholas W. Taylor, Michael Boyle, Christian Reisswig, Mark A. Scheel, Tony Chu, Lawrence E. Kidder, and Béla Szilágyi. Comparing gravitational waveform extrapolation to Cauchy-characteristic extraction in binary black hole simulations. Phys. Rev. D, 88:124010, Dec 2013, 1309.3605.