Physics of

Otto Seiskari

Ray paths and gravitational lensing

As usual in raytracing, each sample (pixel) in the image corresponds to a path starting from the position p0 = (x0,y0,z0) of the camera, initially travelling to a direction d0 3 specified by the camera coordinate system. The path represents photons arriving to the camera from that direction. Instead of traveling in straight lines until hitting an object, the paths solve the photonic geodesic equation in the Schwarzschild space-time, which in certain units and spherical coordinates (t,r,𝜃,ϕ) takes the form (cf. [1, §4.2.1])

ds2 = 1 1 rdt2 + 1 1 r1dr2 r2(d𝜃2 sin 2𝜃 dϕ2) = 0.

Without loss of generality, one can fix 𝜃 = π 2 and, even though not very obviously (see, e.g., [3] and [2]), this gives the ordinary differential equation

u(ϕ) = u(ϕ) 1 3 2u2(ϕ) , (1)

for u := 1r. To solve it in the raytracer, I use the initial conditions

u(0) = 1 p0,u(0) = u(0)d0 n d0 t ,


n = p0 p0 and t = (n ×d0) ×n (n ×d0) ×n. (2)

span the plane in which the ray travels.

Equation (1) is integrated using the Leapfrog method until either a maximum number of steps N is reached or uj + Δuj < 0, which signals that the ray escapes at an angle ϕ (ϕj,ϕj + Δϕj). Given ϕ and u, the corresponding spatial coordinates p = (x,y,z) can be computed from (2) as

p = n cos ϕ + t sin ϕu.

The Schwarzschild time coordinate t is resolved by amending the numerical integration with

t(ϕ) = (u )2 + u2 (1 u) u2(1 u) ,

which is needed to compute light travel times for simulating apparent positions of moving objects. In order to do the integration in a visually accurate way, the step size Δϕj is varied using an ad hoc formula, in addition to using the classical approximation dt2 dx2 + dy2 + dz2 at large distances r from the black hole.

When the observer is moving with velocity v w.r.t. the Schwarzschild coordinate system, the ray directions d in the coordinate system of the observer (the camera coordinate system) are transformed by relativistic aberration, which is corresponds to the Lorentz velocity transformation

d0 = 1 1 + vd (d + v)v0 + 1 γd

where v = v0v, d = v0 d, d = d d and

γ = 1 1 v2. (3)

Circular orbits

The observer (unless stationary) and the planet both move on circular orbits with angular velocities given by [2]

dϕ dt = 1 r2(r 1). (4)

On such orbits, the time dilation of the observer relative to the planet is computed by updating the t-coordinate as

Δt = (Δt0 )2 (1 v2 ) 1 1 r ,

where Δt0 is the elapsed wall clock time between consecutive animation frames and v = rdϕdt.

When computing intersection between the planet and light paths, the motion of the planet has to be taken into account. Within the ODE integration steps, a ray-sphere intersection model can be used after mapping the ray to the moving coordinate system of the planet using a Lorentz velocity transformation.

Doppler shift and beaming

The Doppler factor for light arriving to the observer from direction d, emitted by a source moving at velocity v is [1, §1.7]

δ = γ(1 + d v), (5)

where γ is the Lorentz factor (3). The Doppler factor for a ray path connecting a moving light source and a moving observer can be computed as

δ = δoδs = γoγs(1 do vo)(1 + ds vs),

where do and ds are the initial and final (at intersection with the light source) photon directions on the ray path. The velocities vo and vs of the observer and the source can be determined by (4). All of these are given in the ‘‘stationary’’ Schwarzschild coordinate system. For the accretion disk, the velocity vs is computed assuming the disk consists of glowing matter moving in circular orbits given by (4).

The wavelengths of the emitted and observed light are related by λo = δλs and intensities by I(λo) = δ3I(λ s) [1]. The change in intensity is called relativistic beaming and the change in wavelength is Doppler shift, which is also commonly referred to as redshift (when δ > 1) or blueshift (δ < 1).

The perceived (RGB) colors corresponding to any fixed spectrum I can be precomputed in a one-dimensional lookup table (texture) as a function of the Doppler factor. This is done by integrating the shifted spectra with respect to the standard CIE color sensitivity measures [4]. The spectra used for the colors in this simulation are the black-body spectrum

I(λ) 2hc2 λ5 1 e hc λkBT 1

and hydrogen-α: I(λ) δ656.281nm(λ).


[1]   Jeremy Goodman, Topics in High-Energy Astrophysics

[2]   Riccardo Antonelli, Visualizing a Black Hole

[3]   Wikipedia, Schwarzschild Geodesics

[4]   Munsell Color Science Laboratory, CIE data

All web pages retrieved in December 2015.