Visualizing Geodesics on Surfaces

Geodesics generalize the idea of straight lines in a plane to curved surfaces. They are the paths taken by objects that move freely, without acceleration. They are determined by the geodesic equation

$\frac{d^{2} x^{i}}{d s^{2}} + Γ_{k m}^{i} \frac{d x^{k}}{d s} \frac{d x^{m}}{d s} = 0$

where s is the invariant interval measured along the geodesics and repeated indices are summed. The Christoffel symbols in the equation are defined in terms of a metric tensor by

$Γ_{k l}^{i} = \frac{1}{2} g^{i m} [\frac{\partial g_{m l}}{\partial x^{k}} + \frac{\partial g_{m k}}{\partial x^{l}} - \frac{\partial g_{k l}}{\partial x^{m}}]$

As a step towards meaningful visualizations of geodesics for general relativistic metrics, this presentation will consider geodesics for two-dimensional surfaces embedded in our three-dimensional space. This means that given a parametrization of a surface in terms of two independent variables

$x = x (u, v) y = y (u, v) z = z (u, v)$

the expression for the square of the invariant interval formed by summing squared differentials of the parametrized dependent variables reduces to two dimensions:

$d s^{2} = d x^{2} + d y^{2} + d z^{2} = E d u^{2} + 2 F d u d v + G d v^{2}$

This expression in general has three terms with coefficient functions depending on both independent variables. The two-dimensional metric tensor and its inverse are

$\begin{array}{l} g = [\begin{array}{c} E & F \\ F & G \end{array}] g^{- 1} = \frac{1}{Δ} [\begin{array}{c} G & - F \\ - F & E \end{array}], Δ = det g = E G - F^{2} \end{array}$

In the general case there are six independent Christoffel symbols, where derivatives with respect to the two variables are denoted by numerical indices:

$\begin{array}{l} Γ_{11}^{1} = \frac{E_{1} G - 2 F_{1} F + E_{2} F}{2 Δ} \\ Γ_{12}^{1} = Γ_{21}^{1} = \frac{- G_{1} F + E_{2} G}{2 Δ} \\ Γ_{22}^{1} = \frac{- G_{1} G + 2 F_{2} G - G_{2} F}{2 Δ} \end{array} \begin{array}{l} Γ_{11}^{2} = \frac{- E_{1} F + 2 F_{1} E - E_{2} E}{2 Δ} \\ Γ_{12}^{2} = Γ_{21}^{2} = \frac{G_{1} E - E_{2} F}{2 Δ} \\ Γ_{22}^{2} = \frac{G_{1} F - 2 F_{2} F + G_{2} E}{2 Δ} \end{array}$

One can now easily construct the two geodesic equations using these symbols, remembering to include a factor of two due to symmetry whenever symbols with differing lower indices occur.

Since many people consider evaluating Christoffel symbols to be tedious, it is worth pointing out that there is an alternative method for deriving geodesic equations. Define a Lagrangian for the metric as

$L = \frac{1}{2} [E {\overset{\cdot}{u}}^{2} + 2 F \overset{\cdot}{u} \overset{\cdot}{v} + G {\overset{\cdot}{v}}^{2}]$

where the dot is differentiation with respect to invariant interval, and evaluate the Lagrange equations

$\frac{d}{d s} \frac{\partial L}{\partial {\overset{\cdot}{q}}_{i}} = \frac{\partial L}{\partial q_{i}}$

This process is entirely equivalent to the one requiring explicit statement of Christoffel symbols and perhaps a bit less error prone. The explicit Lagrange equations are

$\begin{array}{l} \frac{d}{d s} [E \overset{\cdot}{u} + F \overset{\cdot}{v}] = \frac{1}{2} [E_{1} {\overset{\cdot}{u}}^{2} + 2 F_{1} \overset{\cdot}{u} \overset{\cdot}{v} + G_{1} {\overset{\cdot}{v}}^{2}] \\ \frac{d}{d s} [F \overset{\cdot}{u} + G \overset{\cdot}{v}] = \frac{1}{2} [E_{2} {\overset{\cdot}{u}}^{2} + 2 F_{2} \overset{\cdot}{u} \overset{\cdot}{v} + G_{2} {\overset{\cdot}{v}}^{2}] \end{array}$

Expanding the left-hand sides and rearranging

$\begin{array}{l} E \overset{\cdot \cdot}{u} + F \overset{\cdot \cdot}{v} = \frac{1}{2} [- E_{1} {\overset{\cdot}{u}}^{2} - 2 E_{2} \overset{\cdot}{u} \overset{\cdot}{v} + (G_{1} - 2 F_{2}) {\overset{\cdot}{v}}^{2}] \\ F \overset{\cdot \cdot}{u} + G \overset{\cdot \cdot}{v} = \frac{1}{2} [(E_{2} - 2 F_{1}) {\overset{\cdot}{u}}^{2} - 2 G_{1} \overset{\cdot}{u} \overset{\cdot}{v} - G_{2} {\overset{\cdot}{v}}^{2}] \end{array}$

so that the second derivatives of the independent variables are

$\begin{array}{l} \overset{\cdot \cdot}{u} = \frac{- E_{1} G + 2 F_{1} F - E_{2} F}{2 Δ} {\overset{\cdot}{u}}^{2} + \frac{G_{1} F - E_{2} G}{Δ} \overset{\cdot}{u} \overset{\cdot}{v} + \frac{G_{1} G - 2 F_{2} G + G_{2} F}{2 Δ} {\overset{\cdot}{v}}^{2} \\ \overset{\cdot \cdot}{v} = \frac{E_{1} F - 2 F_{1} E + E_{2} E}{2 Δ} {\overset{\cdot}{u}}^{2} + \frac{- G_{1} E + E_{2} F}{Δ} \overset{\cdot}{u} \overset{\cdot}{v} + \frac{- G_{1} F + 2 F_{2} F - G_{2} E}{2 Δ} {\overset{\cdot}{v}}^{2} \end{array}$

One can easily verify that the quantities on the right-hand side are the negatives of the Christoffel symbols previously calculated, allowing for the extra factor of two mentioned above for those with differing lower indices.

For the purposes of visualization, there is no need to simplify any further: all the functions involved can be written in code and evaluated as is. Since the set of resulting equations is nonlinear for nontrivial metrics, numerical integration is in any case generally necessary for a solution.

A good starting point is a spherical surface of radius a that can be parametrized with

$\begin{array}{l} x = a sin u cos v \\ y = a sin u sin v \\ z = a cos u \end{array} 0 \leq u \leq π, 0 \leq v \leq 2 π$

The corresponding two-dimensional metric is

$d s^{2} = d x^{2} + d y^{2} + d z^{2} = a^{2} (d u^{2} + {sin}^{2} u d v^{2})$

Not only is the off-diagonal component zero, but the first metric coefficient is constant and the other depends on only one independent variable. The geodesic equations are relatively simple:

$\begin{array}{l} \overset{\cdot \cdot}{u} = \frac{G_{1}}{2 E} {\overset{\cdot}{v}}^{2} = sin u cos u {\overset{\cdot}{v}}^{2} \\ \overset{\cdot \cdot}{v} = - \frac{G_{1}}{G} \overset{\cdot}{u} \overset{\cdot}{v} = - 2 cot u \overset{\cdot}{u} \overset{\cdot}{v} \end{array}$

Note that these equations are independent of the radius. That parameter is here merely an overall scaling factor for the transformation to three dimensions, but will be included for consistency with later surfaces.

Numerical integration for an arbitrary starting point and a spray of geodesics around it:

The geodesics visibly trace out a spherical surface from any starting point. All of these geodesics are closed great circles on the sphere. Closed geodesics do not occur in general on arbitrary surfaces, which will become evident with more complicated cases.

One curiosity of the code for this graphic is how the initial velocities are assigned in the two-dimensional angular space. Since the differential of v is multiplied by sin u in the metric, it must be divided by this quantity to achieve parity with the variable u when assigning values equally spaced around a unit circle. Otherwise there is a noticeable and misleading artifact from the coordinate system in the resulting geodesics. Try altering the underlying code to see it!

Further, the regions around u = 0 and u = π are avoided due to the singularity in the second geodesic equation. The same applies to the initial velocities, where a small offset is used to avoid passing through these poles. This will no longer be necessary when adaptive Runge-Kutta integration is added to Math.

The embedding of this surface in three-dimensional space is easily shown by parametrizing an arbitrary point in three-dimensional space as

$\begin{array}{l} x = r sin u cos v \\ y = r sin u sin v \\ z = r cos u \end{array}$

and evaluating the corresponding metric:

$d s^{2} = d x^{2} + d y^{2} + d z^{2} = d r^{2} + r^{2} d u^{2} + r^{2} {sin}^{2} u d v^{2}$

This is the metric of a flat three-dimensional space expressed in spherical coordinates. The metric of the spherical surface by contrast has only the two dimensions parametrized by angles. The restriction of the radius to a constant value is what produces the curvature of the surface in its ambient three-dimensional space.

A generalization of the simple spherical surface is an ellipsoid of revolution parametrized with

$\begin{array}{l} x = a sin u cos v \\ y = a sin u sin v \\ z = c cos u \end{array} 0 \leq u \leq π, 0 \leq v \leq 2 π$

with an independent scaling factor in the z direction. The metric formed by summing squared differentials is

$d s^{2} = (a^{2} {cos}^{2} u + c^{2} {sin}^{2} u) d u^{2} + a^{2} {sin}^{2} u d v^{2}$

This metric is diagonal like the last but with functional dependence for the first metric function. The geodesic equations are a bit more complicated:

$\begin{array}{l} \overset{\cdot \cdot}{u} = - \frac{E_{1}}{2 E} {\overset{\cdot}{u}}^{2} + \frac{G_{1}}{2 E} {\overset{\cdot}{v}}^{2} = \frac{sin u cos u}{a^{2} {cos}^{2} u + c^{2} {sin}^{2} u} [(a^{2} - c^{2}) {\overset{\cdot}{u}}^{2} + a^{2} {\overset{\cdot}{v}}^{2}] \\ \overset{\cdot \cdot}{v} = - \frac{G_{1}}{G} \overset{\cdot}{u} \overset{\cdot}{v} = - 2 cot u \overset{\cdot}{u} \overset{\cdot}{v} \end{array}$

Simplification is possible for the first equation, but it will be left in this form to avoid introducing additional singularities. It reduces to the first equation for the spherical surface when $a = c$ as expected.

Numerical integration for an arbitrary starting point and a spray of geodesics around it:

The geodesics now visibly trace out an ellipsoid of rotation from any starting point. One major difference from the spherical surface is that these geodesics do not in general close upon themselves, but cover the surface densely as the maximum interval is increased. The initial velocities are again divided by their factors in the metric to reduce artifacts from the coordinate system.

A further generalization of the spherical surface is an ellipsoid parametrized with

$\begin{array}{l} x = a sin u cos v \\ y = b sin u sin v \\ z = c cos u \end{array} 0 \leq u \leq π, 0 \leq v \leq 2 π$

where there are now three independent constants scaling each axial direction. The metric formed by summing squared differentials is

$\begin{array}{l} d s^{2} = [{cos}^{2} u (a^{2} {cos}^{2} v + b^{2} {sin}^{2} v) + c^{2} {sin}^{2} u] d u^{2} \\ + 2 (b^{2} - a^{2}) sin u cos u sin v cos v d u d v \\ + {sin}^{2} u (a^{2} {sin}^{2} v + b^{2} {cos}^{2} v) d φ^{2} \end{array}$

All metric functions here are nonzero and dependent on both variables of parametrization. The corresponding geodesic equations are so complicated that there is not much point in writing them out in full, apart from in the code itself.

Numerical integration for an arbitrary starting point and a spray of geodesics around it:

The geodesics now trace out a general ellipsoid from any starting point. Simplifications with double-angle identities are used in the code to reduce function calls, but this beast still can take some time to evaluate. The equations are mathematically less stable under numerical integration than the previous two examples, so input parameter ranges have been restricted.

This last example indicates how one can visualize geodesics for any complicated two-dimensional surface. The underlying code already has the full geodesic equations for the general case. All that would need to change for other surfaces is the set of metric function definitions and their partial derivatives in the code.

The remainder of this presentation will visualize a couple less complicated compact surfaces for which the geodesic equations are easily expressed. While one can visualize noncompact surfaces similarly, the results are less striking than for compact surfaces because geodesics mostly run off to infinity.

The surface of a torus can be parametrized with

$\begin{array}{l} x = (a + b sin u) cos v \\ y = (a + b sin u) sin v \\ z = b cos u \end{array} 0 \leq u \leq 2 π, 0 \leq v \leq 2 π$

where $a > b$ and the first angular variable runs from the top of the torus around the circumference of the tube. The metric formed by summing squared differentials is

$d s^{2} = b^{2} d u^{2} + (a + b sin u)^{2} d v^{2}$

The metric is diagonal and depends on only one angular variable. The geodesic equations are

$\begin{array}{l} \overset{\cdot \cdot}{u} = \frac{G_{1}}{2 E} {\overset{\cdot}{v}}^{2} = \frac{1}{b} (a + b sin u) cos u {\overset{\cdot}{v}}^{2} \\ \overset{\cdot \cdot}{v} = - \frac{G_{1}}{G} \overset{\cdot}{u} \overset{\cdot}{v} = - \frac{2 b cos u}{a + b sin u} \overset{\cdot}{u} \overset{\cdot}{v} \end{array}$

and do not have singularities in the first angular variable, unlike the surfaces of the sphere and ellipsoid of revolution, as long as a is nonzero.

Numerical integration for an arbitrary starting point and a spray of geodesics around it:

Like the sphere, a toroidal surface can have closed geodesics, but they are special cases. One is visible with the default settings: experiment a bit to find others.

A surface similar to an ellipsoid can be generated by revolution of the ovals of Cassini. Like ellipses these ovals are defined with respect to two centers or foci. For an ellipse the sum of distances to the foci is constant, but for the ovals of Cassini it is the product of the distances. If the foci lie on the x-axis equal distances from the origin, the ovals have the implicit equation

$[(x - a)^{2} + y^{2}] [(x + a)^{2} + y^{2}] = b^{4}$

A parametrization of the curves is found by setting

$x = R (u) cos u y = R (u) sin u$

and solving for the introduced radial variable as a function of angle:

$\begin{array}{c} R^{4} - 2 a^{2} R^{2} cos 2 u + a^{4} = b^{4} \\ R (u) = \pm \sqrt{a^{2} cos 2 u \pm \sqrt{b^{4} - a^{4} {sin}^{2} 2 u}} \end{array}$

The domain of the angular variable u and signs of radicals are chosen to keep the expression real, and that depends on the relative values of the two parameters. When $b > a$ the entire angular domain $0 \leq u \leq 2 π$ is allowed, along with the positive signs of the radicals. When $b \leq a$ the angular variable is restricted to

$- \frac{1}{2} {sin}^{- 1} \frac{b^{2}}{a^{2}} \leq u \leq \frac{1}{2} {sin}^{- 1} \frac{b^{2}}{a^{2}}, b \leq a$

to keep the inner radical real, and all four combinations of signs are necessary to draw the entire curve, except for the the special case $b = a$ where only the inner positive sign can be used. Here is how the curve looks as the parameters are varied:

The ovals becomes circles in the limit $b \to \infty$ for fixed a, as well as $a \to 0$ for nonzero b.

A surface of revolution can be formed by rotation about either the z-axis as for the ellipsoid above or about the x-axis. In the interest of possible future application, the latter will be chosen, with the parametrization

$\begin{array}{l} x = R (u) cos u \\ y = R (u) sin u sin v \\ z = R (u) sin u cos v \end{array} 0 \leq u \leq 2 π, 0 \leq v \leq 2 π$

with the restriction $b > a$ for simplicity. The metric formed by summing squared differentials is

$d s^{2} = [R^{2} + {R^{'}}^{2}] d u^{2} + R^{2} {sin}^{2} u d v^{2}$

Like the ellipsoid of rotation and the torus, the metric is diagonal and depends on only one angular variable. As an aside, for revolution about the z-axis the sine here would be replaced by a cosine.

The geodesic equations for the chosen axis of revolution are

$\begin{array}{l} \overset{\cdot \cdot}{u} = - \frac{E_{1}}{2 E} {\overset{\cdot}{u}}^{2} + \frac{G_{1}}{2 E} {\overset{\cdot}{v}}^{2} = - \frac{1}{2} \frac{d ln [R^{2} + {R^{'}}^{2}]}{d u} {\overset{\cdot}{u}}^{2} + \frac{R R^{'} {sin}^{2} u + R^{2} sin u cos u}{R^{2} + {R^{'}}^{2}} {\overset{\cdot}{v}}^{2} \\ \overset{\cdot \cdot}{v} = - \frac{G_{1}}{G} \overset{\cdot}{u} \overset{\cdot}{v} = - 2 [\frac{R^{'}}{R} + cot u] \overset{\cdot}{u} \overset{\cdot}{v} \end{array}$

The second equation has singularities in the first angular variable like the surfaces of the sphere and ellipsoid of revolution. Differentiating the defining equation for $R (u)$ gives

$R^{'} = - \frac{a^{2} sin 2 u}{R^{2} - a^{2} cos 2 u} R = - \frac{a^{2} sin 2 u}{\sqrt{b^{4} - a^{4} {sin}^{2} 2 u}} R$

from which one can find

$R^{2} + {R^{'}}^{2} = \frac{b^{4}}{b^{4} - a^{4} {sin}^{2} 2 u} R^{2}$

The geodesic equations then can be written fairly compactly as

$\begin{array}{l} \overset{\cdot \cdot}{u} = - [\frac{R^{'}}{R} + 2 (\frac{R^{'}}{R})^{2} cot 2 u] {\overset{\cdot}{u}}^{2} \\ + (1 - \frac{a^{4}}{b^{4}} {sin}^{2} 2 u) [sin u cos u + \frac{R^{'}}{R} {sin}^{2} u] {\overset{\cdot}{v}}^{2} \\ \overset{\cdot \cdot}{v} = - 2 [\frac{R^{'}}{R} + cot u] \overset{\cdot}{u} \overset{\cdot}{v} \end{array}$

with the ratio $\frac{R^{'}}{R}$ defined above. As a check, since $\frac{R^{'}}{R} \to 0$ for either $b \to \infty$ or $a \to 0$ , these equations become those for the spherical surface in both limits as expected.

Numerical integration for an arbitrary starting point and a spray of geodesics around it:

The geodesics for $b < a$ will be left as a possible future addition.

Uploaded 2018.04.03 — Updated 2019.11.24 analyticphysics.com