Gravitational lensing by spinning black holes in astrophysics, and in the movie Interstellar

At a summer school in Les Houches, France, in summer 1972, Bardeen [1], building on earlier work of Carter [2], initiated research on gravitational lensing by spinning black holes. Bardeen gave a thorough analytical analysis of null geodesics (light-ray propagation) around a spinning black hole; and, as part of his analysis, he computed how a black hole's spin affects the shape of the shadow that the hole casts on light from a distant star field. The shadow bulges out on the side of the hole moving away from the observer, and squeezes inward and flattens on the side moving toward the observer. The result, for a maximally spinning hole viewed from afar, is a D-shaped shadow; cf figure 4. (When viewed up close, the shadow's flat edge has a shallow notch cut out of it, as hinted by figure 8.)

Despite this early work, gravitational lensing by black holes remained a backwater of physics research until decades later, when the prospect for actual observations brought it to the fore.

There were, we think, two especially memorable accomplishments in the backwater era. The first was a 1978 simulation of what a camera sees as it orbits a non-spinning black hole, with a star field in the background. This simulation was carried out by Palmer et al [3] on an Evans and Sutherland vector graphics display at Simon Fraser University. Palmer et al did not publish their simulation, but they showed a film clip from it in a number of lectures in that era. The nicest modern-era film clip of this same sort that we know of is by Riazuelo (contained in his DVD [4] and available on the web at [5]); see figure 3 and associated discussion. And see [6] for an online application by Müller and Weiskopf for generating similar film clips. Also of much interest in our modern era are film clips by Hamilton [7] of what a camera sees when falling into a nonspinning black hole; these have been shown at many planetariums, and elsewhere.

The other most memorable backwater-era accomplishment was a black and white simulation by Luminet [8] of what a thin accretion disk, gravitationally lensed by a nonspinning black hole, would look like as seen from far away but close enough to resolve the image. In figure 15(c), we show a modern-era colour version of this, with the camera close to a fast-spinning black hole.

Gravitational lensing by black holes began to be observationally important in the 1990s. Rauch and Blandford [9] recognized that, when a hot spot, in a black hole's accretion disk or jet, passes through caustics of the Earth's past light cone (caustics produced by the hole's spacetime curvature), the brightness of the hot spot's x-rays will undergo sharp oscillations with informative shapes. This has motivated a number of quantitative studies of the Kerr metric's caustics; see, especially [9–11] and references therein.

These papers' caustics are relevant for a source near the black hole and an observer far away, on Earth—in effect, on the black hole's 'celestial sphere' at radius $r=\infty$. In our paper, by contrast, we are interested in light sources that are usually on the celestial sphere and an observer or camera near the black hole. For this reversed case, we shall discuss the relevant caustics in sections 3.3 and 3.4. This case has been little studied, as it is of primarily cultural interest ('everyone' wants to know what it would look like to live near a black hole, but nobody expects to make such observations in his or her lifetime), and of science-fiction interest. Our paper initiates the detailed study of this culturally interesting case; but we leave a full, systematic study to future research. Most importantly, we keep our camera outside the ergosphere—by contrast with Riazuelo's [12] and Hamilton's [7] recent simulations with cameras deep inside the ergosphere and even plunging through the horizon³ .

In the 1990s astrophysicists began to envision an era in which very long baseline interferometry would make possible the imaging of black holes—specifically, their shadows and their accretion disks. This motivated visualizations, with ever increasing sophistication, of accretion disks around black holes: modern variants of Luminet's [8] pioneering work. See, especially, Fukue and Yokoyama [13], who added colours to the disk; Viergutz [14], who made his black hole spin, treated thick disks, and produced particularly nice and interesting coloured images and included the disk's secondary image which wraps under the black hole; Marck [15], who laid the foundations for a lovely movie now available on the web [16] with the camera moving around close to the disk, and who also included higher-order images, as did Fanton et al [17] and Beckwith and Done [18]. See also papers cited in these articles.

In the 2000s astrophysicists have focused on perfecting the mm-interferometer imaging of black-hole shadows and disks, particularly the black hole at the centre of our own Milky Way Galaxy (Sgr A*). See, e.g., the 2000 feasibility study by Falcke et al [19]. See also references on the development and exploitation of general relativistic magnetohydrodynamical (GRMHD) simulation codes for modelling accretion disks like that in Sgr A* [20–22]; and references on detailed GRMHD models of Sgr A* and the models' comparison with observations [23–26]. This is culminating in a mm interferometric system called the Event Horizon Telescope [27], which is beginning to yield interesting observational results though not yet full images of the shadow and disk in Sgr A*.

All the astrophysical visualizations of gravitational lensing and accretion disks described above, and all others that we are aware of, are based on tracing huge numbers of light rays through curved spacetime. A primary goal of today's state-of-the-art, astrophysical ray-tracing codes (e.g., the Chan et al massively parallel, GPU-based code GRay [28]) is very fast throughput, measured, e.g., in integration steps per second; the spatial smoothness of images has been only a secondary concern. For our Interstellar work, by contrast, a primary goal is smoothness of the images, so flickering is minimized when objects move rapidly across an IMAX screen; fast throughput has been only a secondary concern.

With these different primary goals, in our own code, called Double Negative Gravitational Renderer (DNGR), we have been driven to employ a different set of visualization techniques from those of the astrophysics community—techniques based on propagation of ray bundles (light beams) instead of discrete light rays, and on carefully designed spatial filtering to smooth the overlaps of neighbouring beams; see section 2 and appendix A. Although, at Double Negative, we have many GPU-based workstations, the bulk of our computational work is done on a large compute cluster (the Double Negative render-farm) that does not employ GPUs.

In appendix A.7 we shall give additional comparisons of our DNGR code with astrophysical codes and with other film-industry CGI codes.

Our work on gravitational lensing by black holes began in May 2013, when Christopher Nolan asked us to collaborate on building realistic images of a spinning black hole and its disk, with IMAX resolution, for his science-fiction movie Interstellar. We saw this not only as an opportunity to bring realistic black holes into the Hollywood arena, but also an opportunity to create a simulation code capable of exploring a black hole's lensing with a level of image smoothness and dynamics not previously available.

To achieve IMAX quality (with 23 million pixels per image and adequately smooth transitions between pixels), our code needed to integrate not only rays (photon trajectories) from the light source to the simulated camera, but also bundles of rays (light beams) with filtering to smooth the beams' overlap; see section 2, appendices A.2, and A.3. And because the camera would sometimes be moving with speeds that are a substantial fraction of the speed of light, our code needed to incorporate relativistic aberration as well as Doppler shifts and gravitational redshifts.

Thorne, having had a bit of experience with this kind of stuff, put together a step-by-step prescription for how to map a light ray and ray bundle from the light source (the celestial sphere or an accretion disk) to the camera's local sky; see appendices A.1 and A.2. He implemented his prescription in Mathematica to be sure it produced images in accord with others' prior simulations and his own intuition. He then turned his prescription over to our Double Negative team, who created the fast, high-resolution code DNGR that we describe in section 2 and appendix A, and created the images to be lensed: fields of stars and in some cases also dust clouds, nebulae, and the accretion disk around Interstellar's black hole, Gargantua.

In section 3, we use this code to explore a few interesting aspects of the lensing of a star field as viewed by a camera that moves around a fast-spinning black hole (0.999 of maximal—a value beyond which natural spin-down processes become strong [31]). We compute the intersection of the celestial sphere $(r=\infty )$ with the primary, secondary, and tertiary caustics of the camera's past light cone and explore how the sizes and shapes of the resulting caustic curves change as the camera moves closer to the hole and as we go from the primary caustic to higher-order caustics. We compute the images of the first three caustics on the camera's local sky (the first three critical curves), and we explore in detail the creation and annihilation of stellar-image pairs, on the secondary critical curve, when a star passes twice through the secondary caustic. We examine the tracks of stellar images on the camera's local sky as the camera orbits the black hole in its equatorial plane. And we discover and explain why, as seen by a camera orbiting a fast spinning black hole in our Galaxy, there is just one image of the galactic plane between adjacent critical curves near the poles of the black hole's shadow, but there are multiple images between critical curves near the flat, equatorial shadow edge. The key to this from one viewpoint is the fact that higher-order caustics wrap around the celestial sphere multiple times, and from another viewpoint the key is light rays temporarily trapped in prograde, almost circular orbits around the black hole. By placing a checkerboard of paint swatches on the celestial sphere, we explore in detail the overall gravitational lensing patterns seen by a camera near a fast-spinning black hole and the influence of aberration due to the camera's motion.

Whereas most of section 3 on lensing of stellar images is new, in the context of a camera near the hole and stars far away, our section 4, on lensing of an accretion disk, retreads largely known ground. But it does so in the context of the movie Interstellar, to help readers understand the movie's black-hole images, and does so in a manner that may be pedagogically interesting. We begin with a picture of a gravitationally lensed disk made of equally spaced paint swatches. This picture is useful for understanding the multiple images of the disk that wrap over and under and in front of the hole's shadow. We then replace the paint-swatch disk by a fairly realistic and rather thin disk (though one constructed by Double Negative artists instead of by solving astrophysicists' equations for thin accretion disks [32]). We first compute the lensing of our semi-realistic disk but ignore Doppler shifts and gravitational redshifts, which we then turn on pedagogically in two steps: the colour shifts and then the intensity shifts. We discuss why Christopher Nolan and Paul Franklin chose to omit the Doppler shifts in the movie, and chose to slow the black hole's spin from that, $a/M\simeq 1$, required to explain the huge time losses in Interstellar, to $a/M=0.6$. And we then discuss and add simulated lens flare (light scattering and diffraction in the camera's lenses) for an IMAX camera that observes the disk—something an astrophysicist would not want to do because it hides the physics of the disk and the lensed galaxy beyond it, but that is standard in movies, so computer generated images will have continuity with images shot by real cameras.

Finally, in section 5 we summarize and we point to our use of DNGR to produce images of gravitational lensing by wormholes.

Throughout we use geometrized units in which G (Newton's gravitation constant) and c (the speed of light) are set to unity, and we use the MTW sign conventions [33].

The ray tracing part of DNGR produces a map from the celestial sphere (or the surface of an accretion disk) to the camera's local sky. More specifically (see appendix A.1 for details and figure 1 for the ray-tracing geometry):

• (i)  
  In DNGR, we adopt Boyer–Lindquist coordinates $(t,r,\theta ,\phi )$ for the black hole's Kerr spacetime. At each event we introduce the locally non-rotating observer, also called the fiducial observer (FIDO) in the membrane paradigm [34]: the observer whose four-velocity is orthogonal to the surfaces of constant t, the Kerr metric's space slices. We regard the FIDO as at rest in space, and give the FIDO orthonormal basis vectors, $\{{{{\bf e}}_{{\hat{r}}}},{{{\bf e}}_{{\hat{\theta }}}},{{{\bf e}}_{{\hat{\phi }}}}\}$, that point along the spatial coordinate lines.
• (ii)  
  We specify the camera's coordinate location $({{r}_{c}},{{\theta }_{c}},{{\phi }_{c}})$; its direction of motion relative to the FIDO there, a unit vector ${\bf B}$ in the camera's reference frame; and the camera's speed β relative to the FIDO.
• (iii)  
  In the camera's reference frame, we set up a right-handed set of three orthonormal basis vectors $\{{{{\bf e}}_{x}},{{{\bf e}}_{y}},{{{\bf e}}_{z}}\}$, with ${{{\bf e}}_{y}}\equiv {\bf B}$ along the direction of the camera's motion, ${{{\bf e}}_{x}}$ perpendicular to ${{{\bf e}}_{y}}$ and in the plane spanned by ${{{\bf e}}_{{\hat{r}}}}$ and ${{{\bf e}}_{{\hat{\phi }}}}$, and ${{{\bf e}}_{z}}$ orthogonal to ${{{\bf e}}_{x}}$ and ${{{\bf e}}_{y}}$. See figure 1. And we then set up a spherical polar coordinate system $\{{{\theta }_{{\rm cs}}},{{\phi }_{{\rm cs}}}\}$ for the camera's local sky (i.e. for the directions of incoming light rays) in the usual manner dictated by the camera's Cartesian basis vectors $\{{{{\bf e}}_{x}},{{{\bf e}}_{y}},{{{\bf e}}_{z}}\}$.
• (iv)  
  For a ray that originates on the celestial sphere (at $r=\infty$), we denote the Boyer–Lindquist angular location at which it originates by $\{\theta ^{\prime} ,\phi ^{\prime} \}$.
• (v)  
  We integrate the null geodesic equation to propagate the ray from the camera to the celestial sphere, thereby obtaining the map $\{\theta ^{\prime} ({{\theta }_{{\rm cs}}},{{\phi }_{{\rm cs}}}),\phi ^{\prime} ({{\theta }_{{\rm cs}}},{{\phi }_{{\rm cs}}})\}$ of points on the camera's local sky to points on the celestial sphere.
• (vi)  
  If the ray originates on the surface of an accretion disk, we integrate the null geodesic equation backward from the camera until it hits the disk's surface, and thereby deduce the map from a point on the disk's surface to one on the camera's sky. For more details on this case, see appendix A.6.
• (vii)  
  We also compute, using the relevant Doppler shift and gravitational redshift, the net frequency shift from the ray's source to the camera, and the corresponding net change in light intensity.

Zoom In Zoom Out Reset image size

DNGR achieves its IMAX-quality images by integrating a bundle of light rays (a light beam) backward along the null geodesic from the camera to the celestial sphere using a slightly modified variant of a procedure formulated in the 1970s by Pineault and Roeder [35, 36]. This procedure is based on the equation of geodesic deviation and is equivalent to the optical scalar equations [37] that have been widely used by astrophysicists in analytical (but not numerical) studies of gravitational lensing; see references in section 2.3 of [38]. Our procedure, in brief outline, is this (see figure 1); for full details, see appendix A.2.

• (i)  
  In DNGR, we begin with an initially circular (or sometimes initially elliptical) bundle of rays, with very small opening angle, centred on a pixel on the camera's sky.
• (ii)  
  We integrate the equation of geodesic deviation backward in time along the bundle's central ray to deduce the ellipse on the celestial sphere from which the ray bundle comes. More specifically, we compute the angle μ that the ellipse's major axis makes with the celestial sphere's ${{{\bf e}}_{\theta ^{\prime} }}$ direction, and the ellipse's major and minor angular diameters ${{\delta }_{+}}$ and ${{\delta }_{-}}$ on the celestial sphere.
• (iii)  
  We then add up the spectrum and intensity of all the light emitted from within that ellipse; and thence, using the frequency and intensity shifts that were computed by ray tracing, we deduce the spectrum and intensity of the light arriving in the chosen camera pixel.

Novel types of filtering are key to generating our IMAX-quality images for movies. In DNGR we use spatial filtering to smooth the interfaces between beams (ray bundles), and temporal filtering to make dynamical images look like they were filmed with a movie camera. For details, see appendix A.3.

In appendix A.4 we describe some details of our DNGR implementation of the ray-tracing, ray-bundle, and filtering equations; in appendix A.5 we describe some characteristics of our code and of Double Negative's Linux-based render-farm on which we do our computations; in appendix A.6 we describe our DNGR modelling of accretion disks; and in appendix A.7 we briefly compare DNGR with other film-industry CGI codes and state-of-the-art astrophysical simulation codes.

In this subsection we review well known features of gravitational lensing by a nonspinning (Schwarzschild) black hole, in preparation for discussing the same things for a fast-spinning hole.

We begin, pedagogically, with a still from a film clip by Riazuelo [5], figure 2. The camera, at radius ${{r}_{c}}=30M$ (where M is the black hole's mass) is moving in a circular geodesic orbit around the black hole, with a star field on the celestial sphere. We focus on two stars, which each produce two images on the camera's sky. We place red circles around the images of the brighter star and yellow diamonds around those of the dimmer star. As the camera orbits the hole, the images move around the camera's sky along the red and yellow curves.

Zoom In Zoom Out Reset image size

Images outside the Einstein ring (the violet circle) move rightward and deflect away from the ring. These are called primary images. Images inside the Einstein ring (secondary images) appear, in the film clip, to emerge from the edge of the black hole's shadow, loop leftward around the hole, and descend back into the shadow. However, closer inspection with higher resolution reveals that their tracks actually close up along the shadow's edge as shown in the figure; the close-up is not seen in the film clip because the images are so very dim along the inner leg of their tracks. At all times, each star's two images are on opposite sides of the shadow's centre.

This behaviour is generic. Every star (if idealized as a point source of light), except a set of measure zero, has two images that behave in the same manner as the red and yellow ones. Outside the Einstein ring, the entire primary star field flows rightward, deflecting around the ring; inside the ring, the entire secondary star field loops leftward, confined by the ring then back rightward along the shadow's edge. (There actually are more, unseen, images of the star field, even closer to the shadow's edge, that we shall discuss in section 3.2.)

As is well known, this behaviour is easily understood by tracing light rays from the camera to the celestial sphere; see figure 3.

The Einstein ring is the image, on the camera's sky, of a point source that is on the celestial sphere, diametrically opposite the camera; i.e., at the location indicated by the red dot and labeled 'Caustic' in figure 3. Light rays from that caustic point generate the purple ray surface that converges on the camera, and the Einstein ring is the intersection of that ray surface with the camera's local sky.

(The caustic point (red dot) is actually the intersection of the celestial sphere with a caustic line (a one-dimensional sharp edge) on the camera's past light cone. This caustic line extends radially from the black hole's horizon to the caustic point.)

The figure shows a single star (black dot) on the celestial sphere and two light rays that travel from that star to the camera, gravitationally deflecting around opposite sides of the black hole. One of these rays, the primary one, arrives at the camera outside the Einstein ring; the other, secondary ray, arrives inside the Einstein ring.

Because the caustic point and the star on the celestial sphere both have dimension zero, as the camera moves, causing the caustic point to move relative to the star, there is zero probability for it to pass through the star. Therefore, the star's two images will never cross the Einstein ring; one will remain forever outside it and the other inside—and similarly for all other stars in the star field.

However, if a star with finite size passes close to the ring, the gravitational lensing will momentarily stretch its two images into lenticular shapes that hug the Einstein ring and will produce a great, temporary increase in each image's energy flux at the camera due to the temporary increase in the total solid angle subtended by each lenticular image. This increase in flux still occurs when the star's actual size is too small for its images to be resolved, and also in the limit of a point star. For examples, see Riazuelo's film clip [5].

(Large amplifications of extended images are actually seen in nature, for example in the gravitational lensing of distant galaxies by more nearby galaxies or galaxy clusters; see, e.g. [39].)

For a camera orbiting a spinning black hole and a star field (plus sometimes dust clouds and nebulae) on the celestial sphere, we have carried out a number of simulations with our code DNGR. We show a few film clips from these simulations at stacks.iop.org/cqg/32/065001/mmedia. Figure 4 is a still from one of those film clips, in which the hole has spin $a/M=0.999$ (where a is the hole's spin angular momentum per unit mass and M is its mass), and the camera moves along a circular, equatorial, geodesic orbit at radius ${{r}_{c}}=6.03M$.

In the figure we show in violet two critical curves—analogs of the Einstein ring for a nonspinning black hole. These are images, on the camera sky, of two caustic curves that reside on the celestial sphere; see discussion below.

We shall discuss in turn the region outside the secondary (inner) critical curve, and then the region inside.

As the camera moves through one full orbit around the hole, the stellar images in the outer region make one full circuit along the red and yellow curves and other curves like them, largely avoiding the two critical curves of figure 4—particularly the outer (primary) one.

We denote by five-pointed symbols four images of two very special stars: stars that reside where the hole's spin axis intersects the celestial sphere, $\theta ^{\prime} =0$ and π (${{0}^{{\rm o}}}$ and ${{180}^{{\rm o}}}$). These are analogs of the Earth's star Polaris. By symmetry, these pole-star images must remain fixed on the camera's sky as the camera moves along its circular equatorial orbit. Outside the primary (outer) critical curve, all northern-hemisphere stellar images (images with ${{\theta }_{{\rm cs}}}\lt {{90}^{{\rm o}}}$) circulate clockwise around the lower red pole-star image, and southern-hemisphere stellar images, counterclockwise around the upper red pole-star image. Between the primary and secondary (inner) critical curves, the circulations are reversed, so at the primary critical curve there is a divergent shear in the image flow.

(For a nonspinning black hole (figure 2) there are also two critical curves, with the stellar-image motions confined by them: the Einstein ring, and a circular inner critical curve very close to the black hole's shadow, that prevents the inner star tracks from plunging into the shadow and deflects them around the shadow so they close up.)

3.3.1. Primary and secondary critical curves and their caustics

After seeing these stellar-image motions in our simulations, we explored the nature of the critical curves and caustics for a camera near a fast-spinning black hole, and their influence. Our exploration, conceptually, is a rather straightforward generalization of ideas laid out by Rauch and Blandford [9] and by Bozza [10]. They studied a camera or observer on the celestial sphere and light sources orbiting a black hole; our case is the inverse: a camera orbiting the hole and light sources on the celestial sphere.

Just as the Einstein ring, for a nonspinning black hole, is the image of a caustic point on the celestial sphere—the intersection of the celestial sphere with a caustic line on the camera's past light cone—so the critical curves for our spinning black hole are also images of the intersection of the celestial sphere with light-cone caustics. But the spinning hole's light-cone caustics generically are two-dimensional (2D) surfaces (folds) in the three-dimensional light cone, so their intersections with the celestial sphere are one-dimensional: they are closed caustic curves in the celestial sphere, rather than caustic points. The hole's rotation breaks spherical symmetry and converts non-generic caustic points into generic caustic curves. (For this reason, theorems about caustics in the Schwarzschild spacetime, which are rather easy to prove, are of minor importance compared to results about generic caustics in the Kerr spacetime.)

We have computed the caustic curves, for a camera near a spinning black hole, by propagating ray bundles backward from a fine grid of points on the camera's local sky, and searching for points on the celestial sphere where the arriving ray bundle's minor angular diameter ${{\delta }_{-}}$ passes through zero. All such points lie on a caustic, and the locations on the camera sky where their ray bundles originate lie on critical curves.

For a camera in the equatorial plane at radius ${{r}_{c}}=6.03M$, figure 5 shows the primary and secondary caustic curves. These are images, in the celestial sphere, of the primary and secondary critical curves shown in figure 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The primary caustic is a very small astroid (a four-sided figure whose sides are fold caustics and meet in cusps). It spans just ${{0.45}^{{\rm o}}}$ in both the $\phi ^{\prime}$ and $\theta ^{\prime}$ directions. The secondary caustic, by contrast, is a large astroid: it extends over ${{263}^{{\rm o}}}$ in $\phi ^{\prime}$ and ${{72}^{{\rm o}}}$ in $\theta ^{\prime}$. All stars within ${{36}^{{\rm o}}}$ of the equator encounter it as the camera, at ${{r}_{c}}=6.03M$, orbits the black hole. This is similar to the case of a source near the black hole and a camera on the celestial sphere (far from the hole, e.g. on Earth) [10, 11]. There, also, the primary caustic is small and the secondary large. In both cases the dragging of inertial frames stretches the secondary caustic out in the ϕ direction.

3.3.2. Image creations and annihilations on critical curves

Because the spinning hole's caustics have finite cross sections on the celestial sphere, by contrast with the point caustics of a nonspinning black hole, stars, generically, can cross through them; see, e.g., the dashed stellar path in figure 5. As is well known from the elementary theory of fold caustics (see, e.g., section 7.5 of [41]), at each crossing two stellar images, on opposite sides of the caustic's critical curve, merge and annihilate; or two are created. And at the moment of creation or annihilation, the images are very bright.

Figure 6 (two stills from a film clip available at stacks.iop.org/cqg/32/065001/mmedia) is an example of this. As the star in figure 5(a), at polar angle $\theta ^{\prime} =0.608\pi$, travels around the celestial sphere relative to the camera, a subset of its stellar images travels around the red track of figure 6, just below the black hole's shadow. (These are called the star's 'secondary images' because the light rays that bring them to the camera have the same number of poloidal turning points, one—and equatorial crossings, one—as the light rays that map the secondary caustic onto the secondary critical curve; similarly these images' red track is called the star's 'secondary track'.) At the moment of the right still, the star has just barely crossed the secondary caustic at point $\mathcal{A}$ of figure 5(a), and its two secondary stellar images, #2 (inside the secondary critical curve) and #3 (outside it) have just been created at the point half way between #2 and #3, where their red secondary track crosses the secondary critical curve (figure 6(a)). In the meantime, stellar image #1 is traveling slowly, alone, clockwise, around the track, outside the critical curve. Between the left and right stills, image #2 travels the track counter clockwise and images 1 and 3, clockwise. Immediately after the right still, the star crosses the secondary caustic at point $\mathcal{B}$ in figure 5(a), and the two images #1 (outside the critical curve) and #2 (inside it) annihilate at the intersection of their track with the critical curve (figure 6(b)). As the star in figure 5(a) continues on around the celestial sphere from point $\mathcal{B}$ to $\mathcal{A}$, the lone remaining image on the track, image #3, continues onward, clockwise, until it reaches the location #1 of the first still, and two new images are created at the location between #2 and #3 of the first still; and so forth.

Zoom In Zoom Out Reset image size

Figure 7(a) puts these in a broader context. It shows the tracks of all of the images of the star in figure 5(a). Each image is labeled by its order n: the number of poloidal θ turning points on the ray that travels to it from its star on the celestial sphere; or, equally well (for our choice of a camera on the black hole's equator), the number of times that ray crosses the equator $\theta =\pi /2$. The order-0 track is called the primary track, and (with no ray-equator crossings) it lies on the same side of the equator as its star; order-1 is the secondary track, and it lies on the opposite side of the equator from its star; order-2 is the tertiary track, on the same side of the equator as its star; etc.

Zoom In Zoom Out Reset image size

The primary track (order 0) does not intersect the primary critical curve, so a single primary image travels around it as the camera orbits the black hole. The secondary track (order 1) is the one depicted red in figure 6 and discussed above. It crosses the secondary critical curve twice, so there is a single pair creation event and a single annihilation event; at some times there is a single secondary image on the track, and at others there are three. It is not clear to us whether the red secondary track crosses the tertiary critical curve (not shown); but if it does, there will be no pair creations or annihilations at the crossing points, because the secondary track and the tertiary critical curve are generated by rays with different numbers of poloidal turning points, and so the critical curve is incapable of influencing images on the track. The extension to higher-order tracks and critical curves, all closer to the hole's shadow, should be clear. This pattern is qualitatively the same as when the light source is near the black hole and the camera far away, but in the hole's equatorial plane [11].

And for stars at other latitudes the story is also the same; only the shapes of the tracks are changed. Figure 7(b) is an example. It shows the tracks for a star just ${{1.8}^{{\rm o}}}$ above the black hole's equatorial plane, at $\theta ^{\prime} =0.51\pi$.

The film clips at stacks.iop.org/cqg/32/065001/mmedia exhibit these tracks and images all together, and show a plethora of image creations and annihilations. Exploring these clips can be fun and informative.

The version of DNGR that we used for Interstellar showed a surprisingly complex, fingerprint-like structure of gravitationally lensed stars inside the secondary critical curve, along the left side of the shadow.

We searched for errors that might be responsible for it, and finding none, we thought it real. But Riazuelo [42] saw nothing like it in his computed images. Making detailed comparisons with Riazuelo, we found a bug in DNGR. When we corrected the bug, the complex pattern went away, and we got excellent agreement with Riazuelo (when using the same coordinate system), and with images produced by Bohn et al using their Cornell/Caltech SXS imaging code [43]. Since the SXS code is so very different from ours (it is designed to visualize colliding black holes), that agreement gives us high confidence in the results reported below.

Fortunately, the bug we found had no noticeable impact on the images in Interstellar.

With our debugged code, the inner region, inside the secondary critical curve, appears to be a continuation of the pattern seen in the exterior region. There is a third critical curve within the second, and there are signs of higher-order critical curves, all nested inside each other. These are most visible near the flattened edge of the black hole's shadow on the side where the horizon's rotation is toward the camera (the left side in this paper's figures). The dragging of inertial frames moves the critical curves outward from the shadow's flattened edge, enabling us to see things that otherwise could only be seen with a strong zoom-in.

3.4.1. Critical curves and caustics for ${{r}_{c}}=2.60M$

The higher-order critical curves and the regions between them can also be made more visible by moving the camera closer to the black hole's horizon. In figure 8 we have moved the camera in to ${{r}_{c}}=2.6M$ with $a/M=0.999$, and we show three nested critical curves. The primary caustic (the celestial-sphere image of the outer, primary, critical curve) is a tiny astroid, as at ${{r}_{c}}=6.03M$ (figure 5). The secondary and tertiary caustics are shown in figure 9.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

By comparing these three figures with each other, we see that (i) as we move the camera closer to the horizon, the secondary caustics wrap further around the celestial sphere (this is reasonable since frame dragging is stronger nearer the hole), and (ii) at fixed camera radius, each order caustic wraps further around the celestial sphere than the lower order ones. More specifically, for ${{r}_{c}}=2.60M$, the secondary caustic (figure 9(a)) is stretched out to one full circuit around the celestial sphere compared to $2/3$ of a circuit when the camera is at ${{r}_{c}}=6.03M$ (figure 5), and the tertiary caustic (figure 9(b)) is stretched out to more than six circuits around the celestial sphere!

The mapping of points between each caustic and its critical curve is displayed in film clips at stacks.iop.org/cqg/32/065001/mmedia. For the secondary caustic at ${{r}_{c}}=2.60M$, we show a few points in that mapping in figure 9(a). The left side of the critical curve (the side where the horizon is moving toward the camera) maps into the long, stretched-out leftward sweep of the caustic. The right side maps into the caustic's two unstretched right scallops. The same is true of other caustics and their critical curves.

3.4.2. Multiple images for ${{r}_{c}}=2.60M$

3.4.3. Checkerboard to elucidate the multiple-image pattern

Figure 10 is designed to help readers explore this multiple-image phenomenon in greater detail. There we have placed, on the celestial sphere, a checkerboard of paint swatches (figure 10(a)), with dashed lines running along the constant-latitude spaces between paint swatches, i.e., along the celestial-sphere tracks of stars. In figure 10(b) we show the gravitationally lensed checkerboard on the camera's entire sky; and in figure 10(c) we show a blowup of the camera-sky region near the left edge of the black hole's shadow. We have labeled the critical curves 1CC, 2CC and 3CC for primary, secondary, and tertiary.

Zoom In Zoom Out Reset image size

The multiple images of lines of constant celestial-sphere longitude show up clearly in the blow-up, between pairs of critical curves; and the figure shows those lines being stretched vertically, enormously, in the vicinity of each critical curve. The dashed lines (star-image tracks) on the camera's sky show the same kind of pattern as we saw in figure 7.

3.4.4. Multiple images explained by light-ray trapping

The multiple images near the left edge of the shadow can also be understood in terms of the light rays that bring the stellar images to the camera. Those light rays travel from the celestial sphere inward to near the black hole, where they get temporarily trapped, for a few round-trips, on near circular orbits (orbits with nearly constant Boyer–Lindquist radius r), and then escape to the camera. Each such nearly trapped ray is very close to a truly (but unstably) trapped, constant-r ray such as that shown in figure 11. These trapped rays (discussed in [45] and in chapters 6 and 8 of [40]) wind up and down spherical strips with very shallow pitch angles.

Zoom In Zoom Out Reset image size

As the camera makes each additional prograde trip around the black hole, the image carried by each temporarily trapped mapping ray gets wound around the constant-r sphere one more time (i.e., gets stored there for one more circuit), and it comes out to the camera's sky slightly closer to the shadow's edge and slightly higher or lower in latitude. Correspondingly, as the camera moves, the star's image gradually sinks closer to the hole's shadow and gradually changes its latitude—actually moving away from the equator when approaching a critical curve and toward the equator when receding from a critical curve. This behaviour is seen clearly, near the shadow's left edge, in the film clips at stacks.iop.org/cqg/32/065001/mmedia.

The gravitational lensing pattern is strongly influenced not only by the black hole's spin and the camera's location, but also by the camera's orbital speed. We explore this in figure 12, where we show the gravitationally lensed paint-swatch checkerboard of figure 10(a) for a black hole with spin $a/M=0.999$, a camera in the equatorial plane at radius ${{r}_{c}}=2.60M$, and three different camera velocities, all in the azimuthal ${{{\bf e}}_{{\hat{\phi }}}}$ direction: (a) camera moving along a prograde circular geodesic orbit (coordinate angular velocity ${\rm d}\phi /{\rm d}t=1/(a+r_{c}^{3/2})$ in the notation of appendix A.1); (b) camera moving along a zero-angular-momentum orbit (${\rm d}\phi /{\rm d}t=\omega$, equation (A.2), which is speed $-0.546c$ as measured by a circular, geodesic observer); and (c) a static camera, i.e. at rest in the Boyer–Lindquist coordinate system (${\rm d}\phi /{\rm d}t=0$, which is speed $-0.813c$ as measured by the circular, geodesic observer).

Zoom In Zoom Out Reset image size

The huge differences in lensing pattern, for these three different camera velocities, are due, of course, to special relativistic aberration. (We thank Riazuelo for pointing out to us that aberration effects should be huge.) For prograde geodesic motion (top picture), the hole's shadow is relatively small and the sky around it, large. As seen in the geodesic reference frame, the zero-angular-momentum camera and the static camera are moving in the direction of the red/black dot—i.e., toward the right part of the external Universe and away from the right part of the black-hole shadow—at about half and 4/5 the speed of light respectively. So the zero-angular-momentum camera (middle picture) sees the hole's shadow much enlarged due to aberration, and the external Universe shrunken; and the static camera sees the shadow enlarged so much that it encompasses somewhat more than half the sky (more than $2\pi$ steradians), and sees the external Universe correspondingly shrunk.

Despite these huge differences in lensing patterns, the multiplicity of images between critical curves is unchanged: still three images of some near-equator swatches between the primary and secondary critical curves, and eight between the secondary and tertiary critical curves. This is because the caustics in the camera's past light cone depend only on the camera's location and not on its velocity, so a point source's caustic crossings are independent of camera velocity, and the image pair creations and annihilations along critical curves are independent of camera velocity.

We have used our code, DNGR, to construct images of what a thin accretion disk in the equatorial plane of a fast-spinning black hole would look like, seen up close. For our own edification, we explored successively the influence of the bending of light rays (gravitational lensing), the influence of Doppler frequency shifts and gravitational frequency shifts on the disk's colours, the influence of the frequency shifts on the brightness of the disk's light, and the influence of lens flare due to light scattering and diffraction in the lenses of a simulated 65 mm IMAX camera. Although all these issues except lens flare have been explored previously, e.g. in [8, 13, 14, 17, 18] and references therein, our images may be of pedagogical interest, so we show them here. We also show them as a foundation for discussing the choices that were made for Interstellar's accretion disk.

4.1.1. Gravitational lensing

Figure 13 illustrates the influence of gravitational lensing (light-ray bending). To construct this image, we placed the infinitely thin disk shown in the upper left in the equatorial plane around a fast-spinning black hole, and we used DNGR to compute what the disk would look like to a camera near the hole and slightly above the disk plane. The disk consisted of paint swatches arranged in a simple pattern that facilitates identifying, visually, the mapping from the disk to its lensed images. We omitted frequency shifts and their associated colour and brightness changes, and also omitted camera lens flare; i.e., we (incorrectly) transported the light's specific intensity ${{I}_{\nu }}={\rm d}E/{\rm d}t\;{\rm d}A\;{\rm d}\Omega \;{\rm d}\nu$ along each ray, unchanged, as would be appropriate in flat spacetime. Here E, t, A, Ω, and ν are energy, time, area, solid angle, and frequency measured by an observer just above the disk or at the camera.

Zoom In Zoom Out Reset image size

In the figure we see three images of the disk. The upper image swings around the front of the black hole's shadow and then, instead of passing behind the shadow, it swings up over the shadow and back down to close on itself. This wrapping over the shadow has a simple physical origin: light rays from the top face of the disk, which is actually behind the hole, pass up over the top of the hole and down to the camera due to gravitational light deflection; see figure 9.8 of [40]. This entire image comes from light rays emitted by the disk's top face. By looking at the colours, lengths, and widths of the disk's swatches and comparing with those in the inset, one can deduce, in each region of the disk, the details of the gravitational lensing.

In figure 13, the lower disk image wraps under the black hole's shadow and then swings inward, becoming very thin, then up over the shadow and back down and outward to close on itself. This entire image comes from light rays emitted by the disk's bottom face: the wide bottom portion of the image, from rays that originate behind the hole, and travel under the hole and back upward to the camera; the narrow top portion, from rays that originate on the disk's front underside and travel under the hole, upward on its back side, over its top, and down to the camera—making one full loop around the hole.

There is a third disk image whose bottom portion is barely visible near the shadow's edge. That third image consists of light emitted from the disk's top face, that travels around the hole once for the visible bottom part of the image, and one and a half times for the unresolved top part of the image.

In the remainder of this section 4 we deal with a moderately realistic accretion disk—but a disk created for Interstellar by Double Negative artists rather than created by solving astrophysical equations such as [32]. In appendix A.6 we give some details of how this and other Double Negative accretion disk images were created. This artists' Interstellar disk was chosen to be very anemic compared to the disks that astronomers see around black holes and that astrophysicists model—so the humans who travel near it will not get fried by x-rays and gamma-rays. It is physically thin and marginally optically thick and lies in the black hole's equatorial plane. It is not currently accreting onto the black hole, and it has cooled to a position-independent temperature T = 4500 K, at which it emits a black-body spectrum.

Figure 14 shows an image of this artists' disk, generated with a gravitational lensing geometry and computational procedure identical to those for our paint-swatch disk, figure 13 (no frequency shifts or associated colour and brightness changes; no lens flare). Christopher Nolan and Paul Franklin decided that the flattened left edge of the black-hole shadow, and the multiple disk images alongside that left edge, and the off-centred disk would be too confusing for a mass audience. So—although Interstellar's black hole had to spin very fast to produce the huge time dilations seen in the movie—for visual purposes Nolan and Franklin slowed the spin to $a/M=0.6$, resulting in the disk of figure 15(a).

Zoom In Zoom Out Reset image size

4.1.2. Colour and brightness changes due to frequency shifts

Christopher Nolan, the director and co-writer of Interstellar, and Paul Franklin, the visual effects supervisor, were committed to make the film as scientifically accurate as possible—within constraints of not confusing his mass audience unduly and using images that are exciting and fresh. A fully realistic accretion disk, figure 15(c), that is exceedingly lopsided, with the hole's shadow barely discernible, was obviously unacceptable.

The first image in figure 15, the one without frequency shifts and associated colour and brightness changes, was particularly appealing, but it lacked one element of realism that few astrophysicists would ever think of (though astronomers take it into account when modelling their own optical instruments). Movie audiences are accustomed to seeing scenes filmed through a real camera—a camera whose optics scatter and diffract the incoming light, producing what is called lens flare. As is conventional for movies (so that computer generated images will have visual continuity with images shot by real cameras), Nolan and Franklin asked that simulated lens flare be imposed on the accretion-disk image. The result, for the first image in figure 15, is figure 16.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

This, with some embellishments, is the accretion disk seen around the black hole Gargantua in Interstellar.

All of the black-hole and accretion-disk images in Interstellar were generated using DNGR, with a single exception: when Cooper (Matthew McConaughey), riding in the Ranger spacecraft, has plunged into the black hole Gargantua, the camera, looking back upward from inside the event horizon, sees the gravitationally distorted external Universe within the accretion disk and the black-hole shadow outside it—as general relativity predicts. Because DNGR uses Boyer–Lindquist coordinates, which do not extend smoothly through the horizon, this exceptional image had to be constructed by Double Negative artists manipulating DNGR images by hand.

4.3.1. Simulating lens flare

4.3.2. Modelling the accretion disk for Interstellar

In this section we give a step-by-step prescription for computing (i) the ray-tracing map from a point $\{{{\theta }_{{\rm cs}}},{{\phi }_{{\rm cs}}}\}$ on the camera's local sky to a point $\{\theta ^{\prime} ,\phi ^{\prime} \}$ on the celestial sphere, and also (ii) the blue shift ${{f}_{c}}/f^{\prime}$ of light emitted from the celestial sphere with frequency f ' and received at the camera with frequency f_c. Throughout we set the black hole's mass M to unity, so distances are measured in terms of M.

The foundations for our prescription are: (i) the Kerr metric written in Boyer–Lindquist coordinates⁴

$$\begin{eqnarray}&&{\rm d}{{s}^{2}}=-{{\alpha }^{2}}{\rm d}{{t}^{2}}+({{\rho }^{2}}/\Delta ){\rm d}{{r}^{2}}+{{\rho }^{2}}{\rm d}{{\theta }^{2}}+{{\varpi }^{2}}{{({\rm d}\phi -\omega {\rm d}t)}^{2}},\end{eqnarray} \tag{ A.1 }$$

where

$$\begin{eqnarray}\begin{array}{rcl} \rho & = & \sqrt{{{r}^{2}}+{{a}^{2}}{{{\rm cos} }^{2}}\theta },\ \ \Delta ={{r}^{2}}-2r+{{a}^{2}},\ \ \Sigma =\sqrt{{{({{r}^{2}}+{{a}^{2}})}^{2}}-{{a}^{2}}\Delta {{{\rm sin} }^{2}}\theta }, \\ \alpha & = & \frac{\rho \sqrt{\Delta }}{\Sigma },\quad \omega =\frac{2\;ar}{{{\Sigma }^{2}}},\quad \varpi =\frac{\Sigma {\rm sin} \theta }{\rho }; \\ \end{array}\end{eqnarray} \tag{ A.2 }$$

(ii) the 3 + 1 split of spacetime into space plus time embodied in this form of the metric; and (iii) the family of FIDOs whose world lines are orthogonal to the three-spaces of constant t, and their orthonormal basis vectors that lie in those three-spaces and are depicted in figure 1,

$$\begin{eqnarray}&&{{{\bf e}}_{{\hat{r}}}}=\frac{\sqrt{\Delta }}{\rho }\frac{\partial }{\partial r},\quad {{{\bf e}}_{{\hat{\theta }}}}=\frac{1}{\rho }\frac{\partial }{\partial \theta },\quad {{{\bf e}}_{{\hat{\phi }}}}=\frac{1}{\varpi }\frac{\partial }{\partial \phi }.\end{eqnarray} \tag{ A.3 }$$

We shall also need three functions of r and of a ray's constants of motion b (axial angular momentum) and q (Carter constant), which appear in the ray's evolution equations (A.15) below—i.e., in the equations for a null geodesic:

$$\begin{eqnarray}&&P={{r}^{2}}+{{a}^{2}}-ab,\ \ R={{P}^{2}}-\Delta [{{(b-a)}^{2}}+q],\ \ \Theta =q-{{{\rm cos} }^{2}}\theta \left( \frac{{{b}^{2}}}{{{{\rm sin} }^{2}}\theta }-{{a}^{2}} \right).\end{eqnarray} \tag{ A.4 }$$

And we shall need the function ${{q}_{o}}({{b}_{o}})$, which identifies the constants of motion for rays (photons) that are unstably trapped in constant-r orbits around the black hole. This function is defined parametrically by [45]

$$\begin{eqnarray}&&{{b}_{o}}\equiv -\frac{r_{o}^{3}-3r_{o}^{2}+{{a}^{2}}{{r}_{o}}+{{a}^{2}}}{a({{r}_{o}}-1)},\quad \qquad {{q}_{o}}\equiv -\frac{r_{o}^{3}(r_{o}^{3}-6r_{o}^{2}+9{{r}_{o}}-4{{a}^{2}})}{{{a}^{2}}{{({{r}_{o}}-1)}^{2}}},\end{eqnarray} \tag{ A.5 }$$

where r_o is the radius of the trapped orbit, which runs over the interval ${{r}_{1}}\leqslant {{r}_{o}}\leqslant {{r}_{2}}$, with

$$\begin{eqnarray}&&{{r}_{1}}\equiv 2\left\{ 1+{\rm cos} \left[ \frac{2}{3}{\rm arccos} (-a) \right] \right\},\quad \qquad {{r}_{2}}\equiv 2\left\{ 1+{\rm cos} \left[ \frac{2}{3}{\rm arccos} (+a) \right] \right\}.\end{eqnarray} \tag{ A.6 }$$

The prescription that we use, in DNGR, for computing the ray-tracing map and the blue shift is the following concrete embodiment of the discussion in section 2.1.

• (i)  
  Specify the camera's location $({{r}_{c}},{{\theta }_{c}},{{\phi }_{c}})$, and its speed β and the components ${{B}_{{\hat{r}}}}$, ${{B}_{{\hat{\theta }}}}$, and ${{B}_{{\hat{\phi }}}}$ of its direction of motion relative to the FIDO at its location; and specify the ray's incoming direction $({{\theta }_{{\rm cs}}},{{\phi }_{{\rm cs}}})$ on the camera's local sky. (Note: if the camera is in a circular, equatorial geodesic orbit around the black hole, then

  $$\begin{eqnarray}&&\quad \ \ {{B}_{{\hat{r}}}}=0,\quad {{B}_{{\hat{\theta }}}}=0,\quad {{B}_{{\hat{\phi }}}}=1,\quad \beta =\frac{\varpi }{\alpha }(\Omega -\omega ),\quad {\rm where}\;\ \Omega =\frac{1}{(a+r_{c}^{3/2})}\end{eqnarray} \tag{ A.7 }$$

  is the geodesic angular velocity at the camera's radius r_c, and the other quantities are defined in equations (A.2).)
• (ii)  
  Compute, in the camera's proper reference frame, the Cartesian components (figure 1) of the unit vector $N$ that points in the direction of the incoming ray

  $$\begin{eqnarray}&&{{N}_{x}}={\rm sin} {{\theta }_{{\rm cs}}}{\rm cos} {{\phi }_{{\rm cs}}},\ {{N}_{y}}={\rm sin} {{\theta }_{{\rm cs}}}{\rm sin} {{\phi }_{{\rm cs}}},\ {{N}_{z}}={\rm cos} {{\theta }_{{\rm cs}}}.\end{eqnarray} \tag{ A.8 }$$
• (iii)  
  Using the equations for relativistic aberration, compute the direction of motion of the incoming ray, ${{{\bf n}}_{F}}$, as measured by the FIDO in Cartesian coordinates aligned with those of the camera:

  $$\begin{eqnarray}&&{{n}_{F\,y}}=\frac{-{{N}_{y}}+\beta }{1-\beta {{N}_{y}}},\ \ {{n}_{F\,x}}=\frac{-\sqrt{1-{{\beta }^{2}}}{{N}_{x}}}{1-\beta {{N}_{y}}},\ \ {{n}_{F\,z}}=\frac{-\sqrt{1-{{\beta }^{2}}}{{N}_{z}}}{1-\beta {{N}_{y}}};\end{eqnarray} \tag{ A.9 }$$

  and from these, compute the components of ${{{\bf n}}_{F}}$ on the FIDO's spherical orthonormal basis:

  $$\begin{eqnarray}\begin{array}{lcl} \qquad {\mkern 1mu} {{n}_{F\,\hat{r}}} & = & \frac{{{B}_{{\hat{\phi }}}}}{\kappa }\;{{n}_{F\,x}}+{{B}_{{\hat{r}}}}\;{{n}_{F\,y}}+\frac{{{B}_{{\hat{r}}}}{{B}_{{\hat{\theta }}}}}{\kappa }{{n}_{F\,z}},\quad {{n}_{F\,\hat{\theta }}}={{B}_{{\hat{\theta }}}}\;{{n}_{F\,y}}-\kappa \;{{n}_{F\,z}}, \\ {{n}_{F\,\hat{\phi }}} & = & -\frac{{{B}_{{\hat{r}}}}}{\kappa }\;{{n}_{F\,x}}+{{B}_{{\hat{\phi }}}}\;{{n}_{F\,y}}+\frac{{{B}_{{\hat{\theta }}}}{{B}_{{\hat{\phi }}}}}{\kappa }\;{{n}_{F\,z}},\quad {\rm where}\;\kappa \equiv \sqrt{1-B_{{\hat{\theta }}}^{2}}=\sqrt{B_{{\hat{r}}}^{2}+B_{{\hat{\phi }}}^{2}}. \\ \end{array}\end{eqnarray} \tag{ A.10 }$$
• (iv)  
  Compute the ray's canonical momenta (covariant coordinate components of its four-momentum) with its conserved energy $-{{p}_{t}}$ set to unity as a convenient convention:

  $$\begin{eqnarray*}&&{{p}_{t}}=-1,\quad {{p}_{r}}={{E}_{F}}\frac{\rho }{\sqrt{\Delta }}\;{{n}_{F\,\hat{r}}},\quad {{p}_{\theta }}={{E}_{F}}\rho \;{{n}_{F\,\hat{\theta }}},\quad {{p}_{\phi }}={{E}_{F}}\varpi \;{{n}_{F\,\hat{\phi }}},\quad \end{eqnarray*}$$

  $$\begin{eqnarray}&&{\rm where}\;{{E}_{F}}=\frac{1}{\alpha +\omega \varpi {{n}_{F\,\hat{\phi }}}}\end{eqnarray} \tag{ A.11 }$$

  is the energy measured by the FIDO. (Note: ${{p}_{\alpha }}$ can also be regarded as the ray's wave vector or simply as a tangent vector to the ray.) Then compute the ray's other two conserved quantities: its axial angular momentum b and its Carter constant q:

  $$\begin{eqnarray}&&b={{p}_{\phi }},\quad q=p_{\theta }^{2}+{{{\rm cos} }^{2}}\theta \left( \frac{{{b}^{2}}}{{{{\rm sin} }^{2}}\theta }-{{a}^{2}} \right).\end{eqnarray} \tag{ A.12 }$$

  (Note: if we had not set ${{p}_{t}}=-1$, then b would be $-{{p}_{\phi }}/{{p}_{t}}$ and q would be (Carter constant)$/p_{t}^{2}$.)
• (v)  
  Determine, from the ray's constants $\{b,q\}$, whether it comes from the horizon or the celestial sphere by the following algorithm: (a) Are both of the following conditions satisfied?: (v.1)

  $$\begin{eqnarray}&&{{b}_{1}}\lt b\lt {{b}_{2}},\ \ {\rm where}\;{{b}_{1}}={{b}_{o}}({{r}_{2}}),\ \ {{b}_{2}}={{b}_{o}}({{r}_{1}}),\end{eqnarray} \tag{ A.13 }$$

  with r₁ and r₂ given by equations (A.6) and the function ${{b}_{o}}({{r}_{o}})$ given by equations (A.5); and (v.2) the ray's value of q lies in the range

  $$\begin{eqnarray}&&q\lt {{q}_{o}}(b),\ \ {\rm where}\;{{q}_{o}}(b)\;{\rm is}\;{\rm defined}\;{\rm by}\;{\rm equations}\;{\rm (A}{\rm .5)}.\end{eqnarray} \tag{ A.14 }$$

  (b) If the answer to (a) is yes, then there are no radial turning points for that $\{b,q\}$, whence if ${{p}_{\;r}}\gt 0$ at the camera's location, the ray comes from the horizon; and if ${{p}_{\;r}}\lt 0$ there, it comes from the celestial sphere. (c) If the answer to (a) is no, then there are two radial turning points for that $\{b,q\}$, and if the camera radius r_c is greater than or equal to the radius ${{r}_{{\rm up}}}$ of the upper turning point, then the ray comes from the celestial sphere; otherwise, it comes from the horizon. Here ${{r}_{{\rm up}}}$ is the largest real root of $R(r)=0$, with R(r) defined by equations (A.4).
• (vi)  
  If the ray comes from the celestial sphere, then compute its point of origin there as follows: beginning with the computed values of the constants of motion and canonical momenta, and beginning in space at the camera's location, integrate the ray equations⁵

  $$\begin{eqnarray}\begin{array}{rcl} \frac{{\rm d}r}{{\rm d}\zeta } & = & \frac{\Delta }{{{\rho }^{2}}}\;{{p}_{r}},\quad \frac{{\rm d}\theta }{{\rm d}\zeta }=\frac{1}{{{\rho }^{2}}}\;{{p}_{\theta }},\quad \frac{{\rm d}\phi }{{\rm d}\zeta }=-\frac{\partial }{\partial b}\left( \frac{R+\Delta \Theta }{2\Delta {{\rho }^{2}}} \right), \\ \frac{{\rm d}{{p}_{r}}}{{\rm d}\zeta } & = & \frac{\partial }{\partial r}\left[ -\frac{\Delta }{2{{\rho }^{2}}}\;p_{r}^{2}-\frac{1}{2{{\rho }^{2}}}\;p_{\theta }^{2}+\left( \frac{R+\Delta \Theta }{2\Delta {{\rho }^{2}}} \right) \right], \\ \frac{{\rm d}{{p}_{\theta }}}{{\rm d}\zeta } & = & \frac{\partial }{\partial \theta }\left[ -\frac{\Delta }{2{{\rho }^{2}}}\;p_{r}^{2}-\frac{1}{2{{\rho }^{2}}}\;p_{\theta }^{2}+\left( \frac{R+\Delta \Theta }{2\Delta {{\rho }^{2}}} \right) \right] \\ \end{array}\end{eqnarray} \tag{ A.15 }$$

  numerically, backward in time from $\zeta =0$ to $\zeta ={{\zeta }_{f}}$, where ${{\zeta }_{f}}$ is either $-\infty$ or some very negative value. (This super-Hamiltonian version of the ray equations is well behaved and robust at turning points, by contrast, for example, with a commonly used version of the form ${{\rho }^{2}}{\rm d}r/{\rm d}\zeta =\pm \sqrt{R}$, ${{\rho }^{2}}{\rm d}\theta /{\rm d}\zeta =\pm \sqrt{\Theta }$,...; equations (33.32) of MTW [33].) The ray's source point $(\theta ^{\prime} ,\phi ^{\prime} )$ on the celestial sphere is $\theta ^{\prime} =\theta ({{\zeta }_{f}})$, $\phi ^{\prime} =\phi ({{\zeta }_{f}})$.
• (vii)  
  Compute the light's blue shift, i.e. the ratio of its frequency f_c as seen by the camera to its frequency f' as seen by the source at rest on the celestial sphere, from

  $$\begin{eqnarray}&&\frac{{{f}_{c}}}{f^{\prime} }=\left( \frac{\sqrt{1-{{\beta }^{2}}}}{1-\beta {{N}_{y}}} \right)\left( \frac{1-b\omega }{\alpha } \right).\end{eqnarray} \tag{ A.16 }$$

In this section we give a step-by-step prescription for evolving a tiny bundle of rays (a light beam) along a reference ray, from the camera to the celestial sphere, with the goal of learning the beam's major and minor angular diameters ${{\delta }_{+}}$ and ${{\delta }_{-}}$ on the celestial sphere, and the angle μ from ${{{\bf e}}_{\theta ^{\prime} }}$ on the celestial sphere to the beam's major axis (cf figure 1). As in the previous subsection, we set the black hole's mass to unity, M = 1.

Our prescription is a concrete embodiment of the discussion in section 2.2; it is a variant of the Sachs optical scalar equations [37]; and it is a slight modification of a prescription developed in the 1970s by Pineault and Roeder [35, 36]⁶ .

Our prescription relies on the following functions, which are defined along the reference ray. (a) The components of the ray bundle's four-momentum (wave vector) on the FIDO's orthonormal spherical basis

$$\begin{eqnarray}\begin{array}{rcl} {{p}^{{\hat{r}}}} & = & {{p}_{{\hat{r}}}}=\frac{\sqrt{\Delta }}{\rho }{{p}_{r}},\quad {{p}^{{\hat{\theta }}}}={{p}_{{\hat{\theta }}}}=\frac{{{p}_{\theta }}}{\rho },\quad \\ {{p}^{{\hat{\phi }}}} & = & {{p}_{{\hat{\phi }}}}=\frac{{{p}_{\phi }}}{\varpi }=\frac{b}{\varpi },\quad {{p}^{{\hat{t}}}}=\sqrt{p_{{\hat{r}}}^{2}+p_{{\hat{\theta }}}^{2}+p_{{\hat{\phi }}}^{2}}. \\ \end{array}\end{eqnarray} \tag{ A.17 }$$

(b) A function $\mathcal{M}$ defined on the phase space of the reference ray by

$$\begin{eqnarray}\begin{array}{rcl} \mathcal{M} & = & \frac{1}{4{{\rho }^{3}}{{\Sigma }^{2}}({{p}^{{\hat{r}}}}-{{p}^{{\hat{t}}}})}\left\{ 8{{a}^{3}}\sqrt{\Delta }r{{{\rm sin} }^{2}}\theta {\rm cos} \theta \left( p_{{\hat{\theta }}}^{2}-p_{{\hat{\phi }}}^{2} \right) \right. \\ {} & {} & +2({{p}^{{\hat{r}}}}-{{p}^{{\hat{t}}}})(a{\rm sin} \theta \left[ {{a}^{2}}\left( {{p}_{{\hat{\theta }}}}\;(a-r)(a+r){\rm cos} 2\theta +2\sqrt{\Delta }\;{{p}^{{\hat{t}}}}\;r{\rm sin} 2\theta \right) \right. \\ {} & {} & +\left. {{p}_{{\hat{\theta }}}}\;\left( {{a}^{4}}-3{{a}^{2}}{{r}^{2}}-6{{r}^{4}} \right) \right]-2{{p}_{{\hat{\phi }}}}\;{{\left( {{a}^{2}}+{{r}^{2}} \right)}^{3}}{\rm cot} \theta ) \\ {} & {} & -{{a}^{2}}{{p}_{{\hat{\phi }}}}\;{\rm sin} \theta (2\sqrt{\Delta }\;{{p}_{{\hat{\theta }}}}\;{\rm sin} \theta \left[ {{a}^{2}}(r-1){\rm cos} 2\theta +{{a}^{2}}(r-1)+2{{r}^{2}}(r+3) \right] \\ {} & {} & +{{a}^{2}}\Delta {{p}^{{\hat{t}}}}\;{\rm cos} 3\theta +{\rm cos} \theta \left[ {{a}^{4}}(7{{p}^{{\hat{t}}}}-4{{p}^{{\hat{r}}}})+{{a}^{2}}r({{p}^{{\hat{t}}}}\;(15r-22)-8{{p}^{{\hat{r}}}}\;(r-2)) \right. \\ {} & {} & -\left. \left. 4{{r}^{3}}({{p}^{{\hat{r}}}}\;(r-4)-2{{p}^{{\hat{t}}}}\;(r-3)) \right]) \right\}. \\ \end{array}\end{eqnarray} \tag{ A.18 }$$

(c) A complex component ${{\Psi }_{0*}}$ of the Weyl tensor, evaluated on vectors constructed from the ray's four-momentum; its real and imaginary parts work out to be:

$$\begin{eqnarray*}\begin{array}{rcl} \Re ({{\Psi }_{0*}}) & = & \frac{1}{2(w-1){{({{p}^{{\hat{r}}}}-{{p}^{{\hat{t}}}})}^{2}}}\left\{ {{Q}_{1}}\left( -3w{{\left( p_{{\hat{\theta }}}^{2}-p_{{\hat{\phi }}}^{2} \right)}^{2}} \right. \right. \\ {} & {} & -3p_{{\hat{r}}}^{4}w+2p_{{\hat{r}}}^{3}[3{{p}^{{\hat{t}}}}\;w+{{p}_{{\hat{\phi }}}}\;S(w-1)] \\ {} & {} & +6p_{{\hat{r}}}^{2}\left[ p_{{\hat{\theta }}}^{2}-{{({{p}^{{\hat{t}}}})}^{2}}\;w-{{p}^{{\hat{t}}}}\;{{p}_{{\hat{\phi }}}}\;S(w-1)-p_{{\hat{\phi }}}^{2}(w+1) \right] \\ {} & {} & +2{{p}^{{\hat{r}}}}\left[ {{p}_{{\hat{\phi }}}}\;S(w-1)\left( p_{{\hat{\phi }}}^{2}-3p_{{\hat{\theta }}}^{2} \right) \right. \\ \end{array}\end{eqnarray*}$$

$$\begin{eqnarray}\begin{array}{l} {} & {} & +\left. 3{{({{p}^{{\hat{t}}}})}^{3}}w+3{{({{p}^{{\hat{t}}}})}^{2}}{{p}_{{\hat{\phi }}}}\;S(w-1)-3{{p}^{{\hat{t}}}}\;(w+2)\left( p_{{\hat{\theta }}}^{2}-p_{{\hat{\phi }}}^{2} \right) \right] \\ {} & {} & -3{{({{p}^{{\hat{t}}}})}^{4}}w-2{{({{p}^{{\hat{t}}}})}^{3}}{{p}_{{\hat{\phi }}}}\;S(w-1) \\ {} & {} & +\left. 6{{({{p}^{{\hat{t}}}})}^{2}}\left[ p_{{\hat{\theta }}}^{2}(w+1)-p_{{\hat{\phi }}}^{2} \right]-2{{p}^{{\hat{t}}}}\;{{p}_{{\hat{\phi }}}}\;S(w-1)\left( p_{{\hat{\phi }}}^{2}-3p_{{\hat{\theta }}}^{2} \right) \right) \\ {} & {} & +2\ {{Q}_{2}}\;{{p}_{{\hat{\theta }}}}\left( 3{{p}_{{\hat{\phi }}}}\;w\left( p_{{\hat{\phi }}}^{2}-p_{{\hat{\theta }}}^{2} \right) \right. \\ {} & {} & +p_{{\hat{r}}}^{3}S(1-w)+3p_{{\hat{r}}}^{2}[{{p}^{{\hat{t}}}}\;S(w-1)+{{p}_{{\hat{\phi }}}}\;(w+2)] \\ {} & {} & -{{p}_{{\hat{r}}}}\left[ -S(w-1)\left( p_{{\hat{\theta }}}^{2}-3p_{{\hat{\phi }}}^{2} \right)+3{{({{p}^{{\hat{t}}}})}^{2}}\;S(w-1) \right. \\ {} & {} & +\left. 6{{p}^{{\hat{t}}}}\;{{p}_{{\hat{\phi }}}}\;(w+2) \right]+{{({{p}^{{\hat{t}}}})}^{3}}\;S(w-1)+3{{({{p}^{{\hat{t}}}})}^{2}}{{p}_{{\hat{\phi }}}}(w+2) \\ {} & {} & -\left. \left. {{p}^{{\hat{t}}}}\;S(w-1)\left( p_{{\hat{\theta }}}^{2}-3p_{{\hat{\phi }}}^{2} \right) \right) \right\}, \\ \end{array}\end{eqnarray} \tag{ A.19 }$$

$$\begin{eqnarray}\begin{array}{rcl} \Im ({{\Psi }_{0*}}) & = & \frac{1}{(w-1){{({{p}^{{\hat{r}}}}-{{p}^{{\hat{t}}}})}^{2}}}\left\{ -{{Q}_{2}}\left( 6p_{{\hat{\theta }}}^{2}p_{{\hat{\phi }}}^{2}w+p_{{\hat{r}}}^{3}[3{{p}^{{\hat{t}}}}\;w+{{p}_{{\hat{\phi }}}}\;S(w-1)] \right. \right. \\ {} & {} & +3p_{{\hat{r}}}^{2}\left[ p_{{\hat{\theta }}}^{2}(w+1)-2{{({{p}^{{\hat{t}}}})}^{2}}w-{{p}^{{\hat{t}}}}\;{{p}_{{\hat{\phi }}}}\;S(w-1)-p_{{\hat{\phi }}}^{2} \right] \\ {} & {} & +{{p}_{{\hat{r}}}}\left[ {{p}_{{\hat{\phi }}}}\;S(w-1)\left( p_{{\hat{\phi }}}^{2}-3p_{{\hat{\theta }}}^{2} \right) \right. \\ {} & {} & +\left. 3{{({{p}^{{\hat{t}}}})}^{3}}w+3{{({{p}^{{\hat{t}}}})}^{2}}{{p}_{{\hat{\phi }}}}\;S(w-1)-3{{p}^{{\hat{t}}}}(w+2)\left( p_{{\hat{\theta }}}^{2}-p_{{\hat{\phi }}}^{2} \right) \right] \\ {} & {} & -{{({{p}^{{\hat{t}}}})}^{3}}{{p}_{{\hat{\phi }}}}\;S(w-1) \\ {} & {} & +3{{({{p}^{{\hat{t}}}})}^{2}}\left[ p_{{\hat{\theta }}}^{2}-p_{{\hat{\phi }}}^{2}(w+1) \right]-{{p}^{{\hat{t}}}}\;{{p}_{{\hat{\phi }}}}\;S(w-1)\left( p_{{\hat{\phi }}}^{2}-3p_{{\hat{\theta }}}^{2} \right)) \\ {} & {} & +{{Q}_{1}}\;{{p}_{{\hat{\theta }}}}\left( 3{{p}_{{\hat{\phi }}}}\;w\left( p_{{\hat{\phi }}}^{2}-p_{{\hat{\theta }}}^{2} \right) \right. \\ {} & {} & +p_{{\hat{r}}}^{3}S(1-w)+3p_{{\hat{r}}}^{2}[{{p}^{{\hat{t}}}}\;S(w-1)+{{p}_{{\hat{\phi }}}}\;(w+2)] \\ {} & {} & -{{p}_{{\hat{r}}}}\left[ -S(w-1)\left( p_{{\hat{\theta }}}^{2}-3p_{{\hat{\phi }}}^{2} \right)+3{{({{p}^{{\hat{t}}}})}^{2}}\;S(w-1) \right. \\ {} & {} & +6{{p}^{{\hat{t}}}}\;{{p}_{{\hat{\phi }}}}(w+2)]+{{({{p}^{{\hat{t}}}})}^{3}}\;S(w-1)+3{{({{p}^{{\hat{t}}}})}^{2}}{{p}_{{\hat{\phi }}}}(w+2) \\ {} & {} & -\left. \left. {{p}^{{\hat{t}}}}\;S(w-1)\left( p_{{\hat{\theta }}}^{2}-3p_{{\hat{\phi }}}^{2} \right) \right) \right\}. \\ \end{array}\end{eqnarray} \tag{ A.20 }$$

Here Q₁, Q₂, w, and S are given by

$$\begin{eqnarray}\begin{array}{lll} {{Q}_{1}} & = & \frac{r({{r}^{2}}-3{{a}^{2}}{{{\rm cos} }^{2}}\theta )}{{{\rho }^{6}}},\quad \\ {{Q}_{2}} & = & \frac{a{\rm cos} \theta (3{{r}^{2}}-{{a}^{2}}{{{\rm cos} }^{2}}\theta )}{{{\rho }^{6}}}, \\ w & = & \frac{\Delta {{a}^{2}}{{{\rm sin} }^{2}}\theta }{{{({{r}^{2}}+{{a}^{2}})}^{2}}},\quad \\ S & = & \frac{3a\sqrt{\Delta }({{r}^{2}}+{{a}^{2}}){\rm sin} \theta }{{{\Sigma }^{2}}}. \\ \end{array}\end{eqnarray} \tag{ A.21 }$$

We shall state our prescription for computing the shape and orientation of the ray bundle on the celestial sphere separately, for two cases: a ray bundle that begins circular at the camera; and one that begins elliptical. Then, we shall briefly sketch the Pineault–Roeder [35, 36] foundations for these prescriptions.

A.2.1. Circular ray bundle at camera

A.2.2. Elliptical ray bundle at camera

A.2.3. Foundations for these ray-bundle prescriptions

For completeness we briefly describe the Pineault–Roeder [35, 36] foundations for these ray-bundle-evolution prescriptions.

At an arbitrary location ζ along the reference ray, in the reference frame of a FIDO there, denote by $\{{{n}_{F\,\hat{r}}},{{n}_{F\,\hat{\theta }}},{{n}_{F\,\hat{\phi }}}\}$ the components of the fiducial ray's unit tangent vector. (These are the same n_F's as in the previous two subsections, extended to the fiducial ray's entire world line.) And at this arbitrary ζ, introduce the unit basis vectors ${\bf a}$ and ${\bf b}$ (equations (A.30)) that span the plane orthogonal to the fiducial ray.

Notice that far from the black hole (at the celestial sphere), where ${{n}_{F\,\hat{r}}}=-1$ (negative because the reference ray is incoming), ${{n}_{F\,\hat{\theta }}}={{n}_{F\,\hat{\phi }}}=0$, the only nonzero components of these basis vectors are ${{a}_{{\hat{\theta }}}}=1$ and ${{b}_{{\hat{\phi }}}}=1$; so ${\bf a}={{{\bf e}}_{{\hat{\theta }}}}$ and $b={{{\bf e}}_{{\hat{\phi }}}}$.

We can think of the transverse plane at ζ as the complex plane, and think of the ${\bf a}$ direction as the real axis for complex numbers, and ${\bf b}$ as the imaginary axis. Then any transverse vector ${\bf Y}$ can be represented by a complex number Y, with

$$\begin{eqnarray}&&{\bf Y}={{Y}_{a}}\;{\bf a}+{{Y}_{b}}\;{\bf b},\quad \qquad Y={{Y}_{a}}+{\rm i}{{Y}_{b}}.\end{eqnarray} \tag{ A.31 }$$

To describe the ray bundle and its evolution, (following Pineault and Roeder) we introduce two complex numbers (transverse vectors) ξ and η, whose real and imaginary parts are functions that are evolved by the ray-bundle equations (A.23):

$$\begin{eqnarray}&&\xi =u+{\rm i}v,\quad \qquad \eta =g+{\rm i}h.\end{eqnarray} \tag{ A.32 }$$

The outer edge of the ray bundle, at location ζ along the reference ray, is described by the complex number

$$\begin{eqnarray}&&Y=(\xi \;{{{\rm e}}^{{\rm i}\sigma }}+\eta \;{{{\rm e}}^{-{\rm i}\sigma }}){{{\rm e}}^{{\rm i}\chi }},\end{eqnarray} \tag{ A.33 }$$

where χ is also evolved by the ray-bundle equations. Here σ is a parameter that ranges from 0 to $2\pi$, and as it changes, it carries the head of the vector ${\bf Y}$ around the bundle's edge. That vector, of course, is given by

$$\begin{eqnarray}&&{\bf Y}=\Re [(\xi \;{{{\rm e}}^{{\rm i}\sigma }}+\eta \;{{{\rm e}}^{-{\rm i}\sigma }}){{{\rm e}}^{{\rm i}\chi }}]{\bf a}+\Im [(\xi \;{{{\rm e}}^{{\rm i}\sigma }}+\eta \;{{{\rm e}}^{-{\rm i}\sigma }}){{{\rm e}}^{{\rm i}\chi }}]{\bf b}\ \end{eqnarray} \tag{ A.34 }$$

(cf equation (A.31)), and its tail sits on the reference ray while its head sits on the edge of the bundle. One can verify that the shape of this ${\bf Y}(\sigma )$, with σ varying from 0 to $2\pi$, is an ellipse.

Let us denote by ν and λ the arguments of ξ and η, so $\xi =\mid \xi \mid \;{{{\rm e}}^{{\rm i}\nu }}$, $\eta =\mid \eta \mid {{{\rm e}}^{{\rm i}\lambda }}$. Then equation (A.33) reads $Y=\mid \xi \mid {{{\rm e}}^{{\rm i}(\nu +\sigma +\chi )}}+\mid \eta \mid {{{\rm e}}^{{\rm i}(\lambda -\sigma +\chi )}}$. As σ varies from 0 to $2\pi$, the modulus of Y (equal to the length of the vector ${\bf Y}$) reaches its largest value when the phase of the two terms is the same, i.e. when $\nu +\sigma =\lambda -\sigma$, or equivalently when $\sigma =\frac{1}{2}(\lambda -\nu )$. This maximal value of $\mid Y\mid$ is

$$\begin{eqnarray}&&\mid {\bf Y}{{\mid }_{{\rm max} }}=\mid Y{{\mid }_{{\rm max} }}=\mid \xi \mid +\mid \eta \mid ,\end{eqnarray} \tag{ A.35 }$$

and the argument of Y at the max (the angle of ${\bf Y}$ to ${\bf a}$) is

$$\begin{eqnarray}&&\mu =\nu +\sigma +\chi =\frac{1}{2}(\nu +\lambda )+\chi =\frac{1}{2}{\rm arg} (\xi \;\eta )+\chi .\end{eqnarray} \tag{ A.36 }$$

The major diameter (not angular diameter) of the elliptical ray bundle, as measured by a FIDO at the location ζ, is this $\mid {\bf Y}{{\mid }_{{\rm max} }}$ multiplied by ${{\delta }_{F}}={{\delta }_{{\rm cs}}}({{f}_{c}}/{{f}_{F}})$:

$$\begin{eqnarray}&&{\rm (major}\;{\rm diameter}\;{\rm of}\;{\rm ellipse)}={{\delta }_{F}}\mid {\bf Y}{{\mid }_{{\rm max} }}={{\delta }_{F}}(\mid \xi \mid +\mid \eta \mid ).\end{eqnarray} \tag{ A.37 }$$

When the bundle reaches the celestial sphere, its measured angular diameter is the rate of increase of this major diameter with distance traveled. But at the celestial sphere, the distance traveled is equal to the affine parameter ζ, so the measured angular diameter is ${{\delta }_{F}}({\rm d}\mid {\bf Y}{{\mid }_{{\rm max} }}/{\rm d}\zeta )$. Now, the FIDO-measured light frequency f_F is equal to f' at the celestial sphere, and the real and imaginary parts of $\xi =u+{\rm i}v$ and $\eta =g+{\rm i}h$ are increasing linearly with ζ, so it turns out that the formula ${{\delta }_{F}}({\rm d}\mid {\bf Y}{{\mid }_{{\rm max} }}/{\rm d}\zeta )$ becomes the first of equations (A.25); and the angle μ of the major axis to the real axis ${\bf a}$, equation (A.36), becomes equation (A.27).

By an argument similar to the first part of the paragraph before last, one can deduce that anywhere along the evolving ray bundle

$$\begin{eqnarray}&&{\rm (minor}\;{\rm diameter}\;{\rm of}\;{\rm ellipse)}={{\delta }_{F}}\mid {\bf Y}{{\mid }_{{\rm min} }}={{\delta }_{F}}(\mid \xi \mid -\mid \eta \mid ),\end{eqnarray} \tag{ A.38 }$$

and correspondingly that the minor angular diameter at the celestial sphere is the second of equations (A.25).

The evolution equations (A.23) are the equation of geodesic deviation for the separation vector whose spatial part in a FIDO reference frame is the vector ${\bf Y}$ (equation (A.34)) at fixed σ; see section II of Pineault and Roeder [36], and see [55].

A.3.1. Spatial filtering and how we handle point stars

A.3.2. Temporal filtering: motion blur

For moving images, we need to filter over time as well as space. A traditional film camera typically exposes the film for half the time between individual frames of film, so for a typical 24 fps (frame per second) film, the exposure time will be 1/48 s. (This fraction is directly related to the camera's 'shutter angle': the shutters in film cameras are rotating discs with a slice cut out. The film is exposed when the slice is over the emulsion, and the film moves onto the next frame when covered. The most common shape is a semicircular disk which is a '180° shutter'.) Any movement of the object or camera during the exposure leads to motion blur. In our case, if the camera is orbiting a black hole, stellar images close to a critical curve appear to zip around in frantic arcs and would appear as streaks of light on a 1/48 s photograph. We aim to reproduce that effect when creating synthetic images for films.

This motion blur in computer graphics is typically simulated with Monte Carlo methods, computing multiple rays per pixel over the shutter duration. In order to cleanly depict these streaks with Monte Carlo methods, we would need to compute the paths of many additional rays and the computational cost of calculating each ray is very high. Instead, in DNGR we take an analytic approach to motion blur (cf figure A1 ) by calculating how the motion of the camera influences the motion of the deflected beam:

Zoom In Zoom Out Reset image size

The camera's instantaneous position influences the beam's momentary state, and likewise the camera's motion affects the time derivatives of the beam's state. We augment our ray and ray bundle equations (A.15) and (A.23) to track these derivatives throughout the beam's trajectory and end up with a description of an elliptical beam being swept through space during the exposure. The beam takes the form of a swept ellipse when it reaches the celestial sphere, and we integrate the contributions of all stars within this to create the motion-blurred image⁸ . These additional calculations approximately double the computation time, but this is considerably faster than a naive Monte Carlo implementation of comparable quality.

A.3.3. Formal description of our analytic approach to motion blur

We can put our analytic method into a more formal setting, using the mathematical description of motion blurred rendering of Sung et al [59], as follows. We introduce the quantity

$$\begin{eqnarray}&&{{I}_{xy}}=\mathop{\sum }\limits_{l}{{\int }_{\Omega }}{{\int }_{\Delta T}}r(\omega ,t){{g}_{l}}(\omega ,t){{L}_{l}}(\omega ,t)\;{\rm d}t{\rm d}\omega ,\end{eqnarray} \tag{ A.39 }$$

which represents the resulting intensity of a sample located at coordinates (x, y) on the virtual image plane of the camera, corresponding to a ray direction ω, and subtending a solid angle, Ω. The sum is over each object, l, in the scene, $\Delta T$ is the time the shutter is open, and ${{L}_{l}}(\omega ,t)$ represents the radiance of the object directed back toward the camera along ω from object l. The shutter and optics may filter the incident radiance (specific intensity) and this is accounted for by the term $r(\omega ,t)$. The term ${{g}_{l}}(\omega ,t)$ is a geometric term that takes the value 1 if object l is visible or 0 if occluded.

Considering just the distortion of the celestial sphere to start with, we start with a ray in the direction ${{\omega }_{0}}$ at time t₀, in the middle of the exposure. We determine whether this ray gets swallowed up by the black hole or continues to the celestial sphere using equations (A.13) and (A.14). This corresponds to $g({{\omega }_{0}},{{t}_{0}})$. We evolve the initial direction ω and solid angle Ω subtended by the sample into a new direction at the celestial sphere, ${{\omega }^{\prime }}$, and ellipse ${{\Omega }^{\prime }}$, using the method described in appendices A.1 and A.2.

When the ray strikes the celestial sphere, we integrate over the area ${{\Omega }^{\prime }}$, weighting the integration with our truncated Gaussian filter, corresponding to $r(\omega )$. (We assume there is no time-dependent aspect to this filter and ignore the t term. We also assume the celestial sphere is static.)

Our motion blur method approximates $L(\omega )$ with the following expression:

$$\begin{eqnarray}&&L\left( {{\omega }^{\prime }}+t{{\left( \frac{\partial {{\omega }^{\prime }}}{\partial t} \right)}_{{{\omega }_{0}},{{t}_{0}}}} \right),\end{eqnarray} \tag{ A.40 }$$

where the partial derivative is evaluated at the centre of the beam ${{\omega }_{0}}$ and at t₀. This represents a ray that can be swept over an arc in the celestial sphere by varying the value of t. This is convolved with the beam shape ${{\Omega }^{\prime }}$ and filter and occlusion terms to give:

$$\begin{eqnarray}&&{{I}_{xy}}=g({{\omega }_{0}},{{t}_{0}}){{\int }_{{{\Omega }^{\prime }}}}{{\int }_{\Delta T}}r(\omega )L\left( {{\omega }^{\prime }}+t{{\left( \frac{\partial {{\omega }^{\prime }}}{\partial t} \right)}_{{{\omega }_{0}},{{t}_{0}}}} \right)\;{\rm d}t{\rm d}{{\omega }^{\prime }}.\end{eqnarray} \tag{ A.41 }$$

In this form, the consequences of our approximations become clearer: by evaluating g only at the beam centre ${{\omega }_{0}}$ and shutter-centre t₀, a sample that crosses the black hole's shadow during the exposure would be given a constant value of g (either 0 or 1) instead of the correct expression that transitions during the exposure. This would lead to a hard edge at the shadow if the camera performed a fast pan.

Luckily, in Interstellar, the design of the shots meant this rarely happened: the main source of motion blur was the movement of the spaceship around the black hole, and not local camera motion.

To handle any cases where this approximation might cause a problem, we also implemented a hybrid method where we launch a small number of Monte Carlo rays over the shutter duration and each ray sweeps a correspondingly shorter analytic path.

We used a very similar technique for the accretion disk.

DNGR was written in C++ as a command-line application. It takes as input the camera's position, velocity, and field of view, as well as the black hole's location, mass and spin, plus details of any accretion disk, star maps and nebulae.

Each of the 23 million pixels in an IMAX image defines a beam that is evolved as described in appendix A.2.

The beam's evolution is described by the set of coupled first and second order differential equations (A.15) and (A.23) that we put into fully first order form and then numerically integrate backwards in time using a custom implementation of the Runge–Kutta–Fehlberg method (see, e.g., chapter 7 of [57]). This method gives an estimate of the truncation error at every integration step so we can adapt the step size during integration: we take small steps when the beam is bending sharply relative to our coordinates, and large steps when it is bending least. We use empirically determined tolerances to control this behaviour. Evolving the beam along with its central ray triples the time per integration step, on average, compared to evolving only the central ray.

The beam either travels near the black hole and then heads off to the celestial sphere; or it goes into the black hole and is never seen again; or it strikes the accretion disk, in which case it gets attenuated by the disk's optical thickness to extinction or continues through and beyond the disk to pick up additional data from other light sources, but with attenuated amplitude. We use automatic differentiation [58] to track the derivatives of the camera motion through the ray equations.

Each pixel can be calculated independently, so we run the calculations in parallel over multiple CPU cores and over multiple computers on our render-farm.

We use the OpenVDB library [60] to store and navigate volumetric data and Autodesk's Maya [61] to design the motion of the camera. (The motion is chosen to fit the film's narrative.) A custom plug-in running within Maya creates the command line parameters for each frame. These commands are queued up on our render-farm for off-line processing.

A typical IMAX image has 23 million pixels, and for Interstellar we had to generate many thousand images, so DNGR had to be very efficient. It has 40 000 lines of C++ code and runs across Double Negative's Linux-based render-farm. Depending on the degree of gravitational lensing in an image⁹ , it typically takes from 30 min to several hours running on ten CPU cores to create a single IMAX image. The longest renders were those of the close-up accretion disk when we shoe-horned DNGR into Mantra. For Interstellar, render times were never a serious enough issue to influence shot composition or action.

Our London render-farm comprises 1633 Dell-M620 blade servers; each blade has two ten-core E5-2680 Intel Xeon CPUs with 156 GB RAM. During production of Interstellar, several hundred of these were typically being used by our DNGR code.

Our code DNGR includes three different implementations of an accretion disk.

Thin disk: for this, we adapted our DNGR ray-bundle code so it detects intersections of the beam with the disk. At each intersection, we sample the artist's image in the manner described in section 2.2 and appendix A.3 above, to determine the colour and intensity of the bundle's light and we attenuate the beam by the optical thickness (cf appendix A.4).

Volumetric disk: the volumetric accretion disk was built by an artist using SideFX Houdini software and stored in a volumetric data structure containing roughly 17 million voxels (a fairly typical number). Each voxel describes the optical density and colour of the disk in that region.

We used extinction-based sampling [62] to build a mipmap volume representation of this data. The length of the ray bundle was split into short piecewise-linear segments which, in turn, were traced through the voxels. We used the length of the major axis of the beam's cross section to select the closest two levels in the mipmap volume; and in each level we sampled the volume data at steps the length of a voxel, picking up contributions from the colour at each voxel and attenuating the beam by the optical thickness. The results from the two mipmap levels were interpolated before moving on to the next line segment.

Close-up disk with procedural textures: Side Effects Software's renderer, Mantra, has a plug-in architecture that lets you modify its operation. We embedded the DNGR ray-tracing code into a plug-in and used it to generate piecewise-linear ray segments which were evaluated through Mantra. This let us take advantage of Mantra's procedural textures¹⁰ and shading language to create a model of the accretion disk with much more detail than was possible with the limited resolution of a voxelized representation. However this method was much slower so was only used when absolutely necessary.

Disk layers close to the camera: the accretion-disk images in Interstellar were generated in layers that were blended together to form the final images. Occasionally a layer only occupied space in the immediate vicinity of the camera, so close that the infuences of spacetime curvature and gravitational redshifts were negligible. These nearby layers were rendered as if they were in flat spacetime, to reduce computation time.

Frequency shifts. in our visualizations that included Doppler and gravitational frequency shifts, we modelled the accretion disk as a constant-temperature blackbody emitter and we estimated the star temperatures from published NASA data [64]. The Doppler and gravitational frequency shifts correspond to a temperature shift in the blackbody spectra on arrival at the camera. We convolve those temperature-shifted black body spectra with the published spectral sensitivity curves of a typical motion picture film [63] to generate separate red, green and blue radiance values, R, G and B. For the volumetric disk, these are calculated at each step and used in place of artist-created radiance values. In figure 15(b) we removed the doppler-induced intensity change by dividing figure 15(d)'s RGB triplet of radiance values by a weighted mean that, empirically, leaves the eye-perceived colours unchanged: $\{{{R}^{\prime }},{{G}^{\prime }},{{B}^{\prime }}\}=\{R,G,B\}/(0.30\;R+0.59\;G+0.11B)$. This dimmed the blue side of the disk and brightened the red side.

We set the white balance of our virtual camera to render a 6500 K blackbody spectrum with equal red, green and blue pixel values by applying a simple gain to each colour channel. We did not model the complex, nonlinear interaction between the colour-sensitive layers that occurs in real film.

Near the end of section 1.1, we compared our code DNGR with the state-of-the-art astrophysical visualization code GRay. The most important differences—DNGR's use of light-beam mappings versus GRay's use of individual-ray mappings, and DNGR's use of ordinary processors versus GRay's use of GPUs—were motivated by our different goals: smoothness of images in movies versus fast throughput in astrophysics.

So far as we are aware, our code is unique in using light-beam-based mappings. No other code, astrophysical or film CGI, uses them. However, some film CGI codes use flat-spacetime mappings from the camera's image plane to the source plane that are mathematically equivalent to our light beams. Specifically, they use ray differential techniques described by Igehy [30], which track, along a ray, the transverse derivatives of the locations of nearby rays, thereby producing the matrix $\mid \mid \partial {{\xi }^{i}}/\partial {{x}^{j}}\mid \mid$, where ${{\xi }^{i}}({{x}^{j}})$ is a ray's source-plane location as a function of its image-plane location. Tracking these differentials in flat spacetime is moderately easy because they are constant along straight paths between reflections, refractions, etc; and they change via simple matrix transformation laws at interaction locations. Tracking them through curved spacetime would entail the same geodesic-deviation technique as underlies our ray-bundle mappings.

In appendix A.3.1 we have described our methods of imaging star fields, using stars that are point sources of light: we feed each star's light into every light beam that intersects the star's location, with an appropriate weighting.

Other gravitational-lensing codes deal with star fields differently. For example, Müller and Frauendiener [29] calculate where each point-star ends up in the camera's image plane, effectively producing the inverse of the traditional ray-tracing algorithm. They do this for lensing by a non-spinning black hole (Schwarzchild metric) which has a high degree of symmetry, making this calculation tractable. Doing so for a spinning black hole would be much more challenging.

A common method for rendering a star-field is to create a 2D source picture (environment map) that records the stars' positions as finite-sized dots in the source plane (see e.g. [64]). This source picture is then sampled using the evolved rays from the camera. This has the disadvantage that stars can get stretched in an unrealistic way in areas of extreme magnification, such as near the critical curves described in section 3.3.1 of this paper. As we discussed in appendix A.3.1, our DNGR light-beam technique with point stars circumvents this problem—as does the Müller-Frauendiener technique with point stars.