When thinking about renders of things like warp drives and black holes we usually just expect to see a simple approximation or an artist rendition, assuming the math required to pull off something accurate would require someone with at least a PhD in Mathematical Physics. Which I won't tell that its completely untrue, but in this blog post I'll try to explain a way to do quite accurate visualizations within a 100 or so lines of code, for basically any kind of space-time for which you can write its metric as code. Although, the detailed mathematical derivation of this approach might be somewhat math heavy.

The main ingredient of any GR render is figuring out how the rays of light move around. Knowing how light moves we can trace rays from the camera into the scene and see where the light came from. So, to render the simplest scene without objects we simply trace a ray for each pixel and assign the color of the pixel to the color of the skybox in the direction in which the ray ends up pointing to.

• What are geodesics?
• Mathematical description of shortest path
• Lagrangian description of a geodesic
• Hamiltonian description of a geodesic
• Writing a Hamiltonian geodesic tracer in GLSL
• Conclusions
• References

What are geodesics?

How exactly do we trace rays in curves space? Any object inside a curved space follows something called a geodesic.

A geodesic is just a fancy word for, in some sense, a path of shortest length between 2 points inside a space.

I should note that there could be multiple of such paths, which are locally minimal, in the sense that you can't nudge the path to make it shorter, while globally there might be a shorter path. Also, in Minkowski space-time the definition is a bit more complicated, because of the imaginary distances (when \( ds^2 < 0 \) ).

In our case, instead of paths between 2 points, we are only interested in finding how a ray moves given an initial point and direction, but the definition above will still prove useful when deriving the equations describing a geodesic, which we will use here.

Mathematical description of shortest path

I'll try to quickly go through the derivation. If you wish to skip the math part, jump to Writing a Hamiltonian geodesic tracer in GLSL.

Mathematically speaking we have some coordinate system, a path, and a way to compute distances between 2 points.

A coordinate system being a set of several numbers labeling each point in the space. A path is a function that takes in the path parameter and outputs a coordinate. In General Relativity the path parameter is usually proper time (like a clock moving with the object, labeling each point), but it can be anything really. The way to compute distances is called a metric (and it's the main source of scary math here).

In physics, or more generally differential geometry, a metric is defined as an integral ("sum") of something called the metric tensor. A metric tensor is a bilinear form \( g(a, b) \), it essentially maps pairs of vectors to real numbers, and is a generalization of dot product for curved spaces. So using a metric tensor we can find the length of a vector in space, and also distances \( ds \) between infinitely close points in space.

\begin{equation}
ds^2 = g(dx, dx)
\end{equation}

In our case, where we describe vectors as a set of numbers, a metric is simply a matrix product of some matrix \( g_{\mu \nu} \) times the vectors. For our infinitesimal distance \( ds \) we get this expression:

\begin{equation}
ds^2 = \sum_{\mu \nu}^N g_{\mu \nu} dx^\mu dx^\nu
\end{equation}

Usually the sum is just implicitly assumed by Einstein notation [1].

\begin{equation}
ds^2 = g_{\mu \nu} dx^\mu dx^\nu
\end{equation}

Here we can actually see that for some simple choice of \( g_{\mu \nu} \) we can get the distances by Pythagoras' theorem. Specifically for the case when the metric tensor matrix is a unit matrix.

\begin{equation}
ds^2 = dx_1^2 + dx_2^2 + dx_3^2
\end{equation}

For a flat space-time like space, we get something similar but with the exception that the time coordinate component is with a negative sign.

\begin{equation}
ds^2 = - dx_0^2 + dx_1^2 + dx_2^2 + dx_3^2
\label{flat}
\end{equation}

Here I used the \( (- + + +) \) signature, but signs can actually be flipped without changing the geodesics, and in some cases, like for particle physics, it makes more sense to use the opposite \( (+ - - -) \) signature.

Going back to the main question of computing distances, to compute the length between 2 points, A and B, along some path \( x^i(t) \) we simply need to sum the infinitesimal distances together using an integral:

\begin{equation}
l = \int_A^B \sqrt{g_{\mu \nu} dx^\mu dx^\nu} = \int_A^B \sqrt{g_{\mu \nu} dx^\mu dx^\nu} \frac{dt}{dt} = \int_A^B \sqrt{g_{\mu \nu} \frac{dx^\mu}{dt} \frac{dx^\nu}{dt}} dt
\end{equation}

Where \( \frac{dx^i}{dt} \) is simply how fast the coordinate x changes with respect to the path parameter ("clock"), in some sense can be interpreted as the velocity.

Now our main question is how do we minimize the path length? Here is where we introduce a thing called calculus of variations, which is roughly speaking a way to find how a functional (distance) changes by infinitesimally small variations of its input function (path). Such derivatives have similar properties to normal function derivatives. And in fact, similarly to calculus, to find the extremum of a function (min, max or stationary point), we simply need to equate the variation to 0.

Lagrangian description of a geodesic

There is an entire branch of physics related to variational principles, which states that any kind of physical system has a value it likes to minimize (or more generally make unchanging under small variations of path, i.e. "stationary"). That value is called action, and the function under the integral is called the Lagrangian function of the system. The branch of physics studying Lagrangian functions of systems is called Lagrangian mechanics.

In our case the Lagrangian can be written like this:

\begin{equation}
L = \sqrt{g_{\mu \nu} \frac{dx^\mu}{dt} \frac{dx^\nu}{dt}}
\end{equation}

Turns out we don't need the square root for the minimum of the functional to be a geodesic, and we can simply remove it from our geodesic Lagrangian:

\begin{equation}
L = g_{\mu \nu} \frac{dx^\mu}{dt} \frac{dx^\nu}{dt}
\end{equation}

The proof of this you can find here [2]. The only difference such simplification makes is that the parametrization of the path might be different, but the path itself will be the same. Also notably, the equations for this specific case turn out to be the same, meaning the parametrization is also the same.

Additionally, we want this with a 1/2 factor, to simplify the equations down the line.

\begin{equation}
L = \frac{1}{2} g_{\mu \nu} \frac{dx^\mu}{dt} \frac{dx^\nu}{dt}
\end{equation}

So, our goal right now is to minimize this functional:

\begin{equation}
S = \int_A^B \frac{1}{2} g_{\mu \nu} \frac{dx^\mu}{dt} \frac{dx^\nu}{dt} dt
\end{equation}

In general the minimum of a functional like this can be found by applying the Euler-Lagrange equations [3]:

\begin{equation}
\frac{\partial L}{\partial x^i} - \frac{d}{dt} \frac{\partial L}{\partial \frac{dx^i}{dt}} = 0
\label{el}
\end{equation}

You can find a more detailed derivation here [4]

Let's derive the Euler-Lagrange equations for our geodesic Lagrangian (keep in mind that there is an equation for each coordinate \( x^i \) ):

\begin{equation}
\frac{\partial L}{\partial \frac{dx^i}{dt}} = \frac{1}{2} \frac{\partial }{\partial \frac{dx^i}{dt}} g_{\mu \nu} \frac{dx^\mu}{dt} \frac{dx^\nu}{dt} = \frac{1}{2} g_{i \nu} \frac{dx^\nu}{dt} + \frac{1}{2} g_{\mu i} \frac{dx^\mu}{dt} = g_{i \nu} \frac{dx^\nu}{dt} \label{el0}
\end{equation}

Then we take the derivative with respect to the path parameter:

\begin{equation}
\frac{d}{dt} \left( g_{i \nu} \frac{dx^\nu}{dt} \right) = \frac{d g_{i \nu} }{dt} \frac{dx^\nu}{dt} + g_{i \nu} \frac{d^2x^\nu}{dt^2} = \frac{d g_{i \nu} }{dx^\mu} \frac{dx^\mu}{dt} \frac{dx^\nu}{dt} + g_{i \nu} \frac{d^2x^\nu}{dt^2} \label{el1}
\end{equation}

And lastly:

\begin{equation}
\frac{\partial L}{\partial x^i} = \frac{1}{2} \frac{d g_{\mu \nu} }{dx^i} \frac{dx^\mu}{dt} \frac{dx^\nu}{dt} \label{el2}
\end{equation}

Substituting \eqref{el1} and \eqref{el2} into Euler-Lagrange equations \eqref{el} leads us to the equation of a geodesic:

\begin{equation}
\frac{1}{2} \frac{d g_{\mu \nu} }{dx^i} \frac{dx^\mu}{dt} \frac{dx^\nu}{dt} - \frac{d g_{i \nu} }{dx^\mu} \frac{dx^\mu}{dt} \frac{dx^\nu}{dt} - g_{i \nu} \frac{d^2x^\nu}{dt^2} = 0
\end{equation}

Multiplying by the metric tensor inverse \( - g^{i \nu} \) we get:

\begin{equation}
\frac{d^2x^i}{dt^2} + g^{i \nu} \left( \frac{d g_{i \nu} }{dx^\mu} - \frac{1}{2} \frac{d g_{\mu \nu} }{dx^i} \right) \frac{dx^\mu}{dt} \frac{dx^\nu}{dt} = 0
\end{equation}

And that's our final system of equations for a geodesic, we could of course also substitute the Christoffel symbols here, but for our application there is no difference. Of course, we could just use these equations for tracing geodesic rays and call it a day, but unfortunately this would require computing a whole lot of derivatives (in 4d space-time it's 64 of them to be specific), either manually, or by using numerical differentiation. Thankfully there is a way to avoid this, and in fact simplify the entire algorithm! (At a slight performance cost)

Hamiltonian description of a geodesic

So here comes the star of the show - Hamiltonian mechanics. Hamiltonian equations of motion have a really nice form which allows to easily write a computer program that integrates them by using Euler integration.

\begin{equation}
\frac{dp^i}{dt} = - \frac{\partial H}{\partial x^i}
\end{equation}
\begin{equation}
\frac{dx^i}{dt} = \frac{\partial H}{\partial p^i}
\end{equation}

Where \( p \) is the so called generalized momentum, it's the derivative of the Lagrangian with respect to the coordinate path parameter derivative:

\begin{equation}
p_i = \frac{\partial L}{\partial \frac{dx^i}{dt} }
\label{momentumdef}
\end{equation}

And to get the Hamiltonian itself you need to apply the Legendre Transform [6] on the Lagrangian:

\begin{equation}
H = \sum_{i}^N p^i \frac{dx^i}{dt} - L
\label{legandre}
\end{equation}

A different derivation of Hamilton's equations of motion can be found here [5].

And for our case the momentum would be the following, which we already computed when writing down the Euler-Lagrange equations \eqref{el0}, and given the definition of generalized momentum \eqref{momentumdef}:

\begin{equation}
p_i = g_{i j} \frac{dx^j}{dt}
\label{momentum}
\end{equation}

To get the "time" derivatives you simply need to multiply both sides by the metric tensor inverse:

\begin{equation}
\frac{dx^i}{dt} = g^{i j} p_j
\label{dxdt}
\end{equation}

And the Hamiltonian itself:

\begin{equation}
H = \sum_{i}^N \frac{dx^i}{dt} p i - L = g{i j} \frac{dx^i}{dt} \frac{dx^j}{dt} - \frac{1}{2} g{i j} \frac{dx^i}{dt} \frac{dx^j}{dt} = \frac{1}{2} g{i j} \frac{dx^i}{dt} \frac{dx^j}{dt} = L
\end{equation}

Turns out that for this simple choice of a geodesic Lagrangian, the Hamiltonian is equal to the Lagrangian!

Also, we want to know the Hamiltonian as a function of the generalized momentum by substituting \eqref{dxdt} into the Hamiltonian equation:

\begin{equation}
H = \frac{1}{2} g_{i j} \frac{dx^i}{dt} \frac{dx^j}{dt} = \frac{1}{2} g^{i j} p_i p_j
\label{hamiltonian}
\end{equation}

While the equations of motion will simply be:

\begin{equation}
\frac{dp_i}{dt} = - \frac{\partial H}{\partial x^i}
\label{eqmotion1}
\end{equation}

\begin{equation}
\frac{dx^i}{dt} = g^{i j} p_j
\label{eqmotion2}
\end{equation}

This is all we need to write a numerical geodesic integrator!

Writing a Hamiltonian geodesic tracer in GLSL

You might have noticed that in the final Hamilton's equations of motion I didn't write out \( \frac{\partial H}{\partial x^i} \), this is actually important! We want to keep the derivative of the Hamiltonian as is, because then instead of computing the 64 derivatives of the metric tensor, we only need 4 to find the Hamiltonian gradient. This is the main simplification of the geodesic tracing algorithm.

Here we will use the GLSL shading language, since it has variables and functions which map quite well to the mathematical operations we will perform here. On top of that we can easily then make a real time GR visualization shader.

First of all, we need a function that evaluates the metric tensor at a 4d point in space and time. Let's use the Alcubierre warp drive [7] metric as an example, since it is quite simple.

mat4 Metric(vec4 x)
{
  //Alcubierre metric
  const float R = 1.0;
  const float sigma = 35.0; 
  const float v = 1.1;

  float x = v*x.x;
  float r = sqrt(sqr(x.y - x) + x.z*x.z + x.w*x.w);
  float f = 0.5*(tanh(sigma*(r + R)) - tanh(sigma*(r - R)))/tanh(sigma*R);
  float gtt = v*v*f*f - 1.0;
  float gxt = -v*f;

  return mat4(gtt, gxt,  0,  0,
              gxt,   1,  0,  0,
                0,   0,  1,  0,
                0,   0,  0,  1);
}

In our case x is a 4D vector representing position. The first component x.x or x[0] being time. As an output we get a 4 by 4 matrix represented by mat4 in GLSL.

Then we need to write down the Hamiltonian \eqref{hamiltonian}. The Hamiltonian is a function that takes 2 things, the position in space-time, and the 4d momentum vector, and outputs a scalar.

float Hamiltonian(vec4 x, vec4 p)
{
  mat4 g_inv = inverse(Metric(x));
  return 0.5*dot(g_inv*p,p);
}

As a bonus here is the Lagrangian

float Lagrangian(vec4 x, vec4 dxdt)
{
  return 0.5*dot(Metric(x)*dxdt,dxdt);
}

Surprisingly enough that's it, GLSL already has a matrix inverse function inverse(), on top of it the Hamiltonian is just the dot product (in GLSL sense) of g_inv*p and p, which are the contravariant and covariant momentum vectors respectively. The contravariant momentum actually just being the time derivative of the coordinate dxdt, i.e. dot(dxdt,p).

After this we need to compute the 4D gradient of the Hamiltonian. We can do this by using a forward numerical difference in all 4 spacial directions, using some small value eps:

vec4 HamiltonianGradient(vec4 x, vec4 p)
{
  const float eps = 0.001;
  return (vec4(Hamiltonian(x + vec4(eps,0,0,0), p),
               Hamiltonian(x + vec4(0,eps,0,0), p),
               Hamiltonian(x + vec4(0,0,eps,0), p),
               Hamiltonian(x + vec4(0,0,0,eps), p)) - Hamiltonian(x,p))/eps;
}

Now that we have the Hamiltonian gradient, we can finally write down the equation of motion \eqref{eqmotion1} \eqref{eqmotion2} integration code

vec4 IntegrationStep(inout vec4 x, inout vec4 p)
{
  const float TimeStep = 0.1;
  p = p - TimeStep * HamiltonianGradient(x, p);
  x = x + TimeStep * inverse(Metric(x)) * p;
}

You might ask, "wait, that's it?", and indeed that is all you need to integrate the geodesic. Of course, it is quite slow since we do a whopping 6 matrix inverse evaluations, which can be optimized down to 1, by replacing most Hamiltonians with Lagrangians which don’t have inverses, since they are equal. Even better is to have the metric inverse already computed analytically, but it’s not possible for every metric, especially for an implicitly defined one.

There is of course the last problem, while initializing the space-time position is easy, how do we initialize the value of the momentum vector p when starting to trace?

Before tracing the geodesic, you can use the equation \eqref{momentum}

p = Metric(x) * dxdt;

But what is dxdt? It's nothing more than the 4D direction the ray moves inside space-time. There are 3 categories the directions can fall into:

• Time-like, when \( A < 0 \)
• Null, when \( A = 0 \)
• Space-like, when \( A > 0 \)

Where \(A\) is
\begin{equation}
A = g_{\mu \nu} \frac{dx^\mu}{dt} \frac{dx^\nu}{dt}
\end{equation}

or in GLSL, A = dot(Metric(x) * dxdt, dxdt)

When simulating how light travels we just want null directions, which lead to null geodesic solutions. On the other hand, if you want to model an object moving slower than light you need a time-like geodesic. And space-like geodesics for tachyonic stuff, which doesn't happen in real life, so we ignore it.

So, assuming a flat space metric \eqref{flat} and some 3D direction in space our p for a light ray would be

vec4 GetNullMomentum(vec3 dir)
{
  return Metric(x) * vec4(1.0, normalize(dir));
}

And the inverse of this operation

vec3 GetDirection(vec4 p)
{
  vec4 dxdt = inverse(Metric(x)) * p;
  return normalize(dxdt.yzw);
}

So, in the end the final simple algorithm will look like this:

void TraceGeodesic(inout vec3 pos, inout vec3 dir, inout float time)
{
  vec4 x = vec4(time, pos);
  vec4 p = GetNullMomentum(dir);

  const int steps = 256;
  for(int i = 0; i < steps; i++)
  {
    IntegrationStep(x, p);
    //you can add a stop condition here when x is below the event horizon for example
  }

  pos = x.yzw;
  time = x.x;
  dir = GetDirection(p);
}

Essentially this is just a 4D ray marching algorithm where the direction of the ray changes every step. In this specific case the size of the step also changes, which can be avoided by normalizing the momentum p = normalize(p). This only changes the step length, and doesn't change the geodesic path, i.e., it works just like a dynamic reparameterization of the path. The time step of the integration can also be varied depending on the metric used. For example, in the case of black holes I change the time step proportionally to the distance to the event horizon, so that the accuracy of the geodesic is roughly proportional to the curvature of space. This is an important optimization to get accurate results, while keeping the computational cost relatively small.

Example shadertoy code:

mat4 diag(vec4 a)
{
    return mat4(a.x,0,0,0,
                0,a.y,0,0,
                0,0,a.z,0,
                0,0,0,a.w);
}

mat4 Metric(vec4 x)
{
    //Kerr-Newman metric in Kerr-Schild coordinates 
    const float a = 0.8;
    const float m = 1.0;
    const float Q = 0.0;
    vec3 p = x.yzw;
    float rho = dot(p,p) - a*a;
    float r2 = 0.5*(rho + sqrt(rho*rho + 4.0*a*a*p.z*p.z));
    float r = sqrt(r2);
    vec4 k = vec4(1, (r*p.x + a*p.y)/(r2 + a*a), (r*p.y - a*p.x)/(r2 + a*a), p.z/r);
    float f = r2*(2.0*m*r - Q*Q)/(r2*r2 + a*a*p.z*p.z);
    return f*mat4(k.x*k, k.y*k, k.z*k, k.w*k)+diag(vec4(-1,1,1,1));
}

float Hamiltonian(vec4 x, vec4 p)
{
    mat4 g_inv = inverse(Metric(x));
    return 0.5*dot(g_inv*p,p);
}

/*
float Lagrangian(vec4 x, vec4 dxdt)
{
    return 0.5*dot(Metric(x)*dxdt,dxdt);
}
*/

vec4 HamiltonianGradient(vec4 x, vec4 p)
{
    const float eps = 0.001;
    return (vec4(Hamiltonian(x + vec4(eps,0,0,0), p),
                 Hamiltonian(x + vec4(0,eps,0,0), p),
                 Hamiltonian(x + vec4(0,0,eps,0), p),
                 Hamiltonian(x + vec4(0,0,0,eps), p)) - Hamiltonian(x,p))/eps;
}

void IntegrationStep(inout vec4 x, inout vec4 p)
{
    const float TimeStep = 0.15;
    p = p - TimeStep * HamiltonianGradient(x, p);
    x = x + TimeStep * inverse(Metric(x)) * p;
}

vec4 GetNullMomentum(vec4 x, vec3 dir)
{
    return Metric(x) * vec4(1.0, normalize(dir));
}

vec3 GetDirection(vec4 x, vec4 p)
{
    vec4 dxdt = inverse(Metric(x)) * p;
    return normalize(dxdt.yzw);
}

void TraceGeodesic(inout vec3 pos, inout vec3 dir, inout float time)
{
    vec4 x = vec4(time, pos);
    vec4 p = GetNullMomentum(x, dir);

    const int steps = 256;
    for(int i = 0; i < steps; i++)
    {
        IntegrationStep(x, p);
        //you can add a stop condition here when x is below the event horizon for example
    }

    pos = x.yzw;
    time = x.x;
    dir = GetDirection(x, p);
}

void mainImage(out vec4 fragColor, in vec2 fragCoord) {
    vec2 uv = 2.0 * (fragCoord - 0.5 * iResolution.xy) / max(iResolution.x, iResolution.y);

    vec3 RayPos = vec3( 0.0,  0.0,  32.0);
    vec3 RayDir = vec3(uv.x, uv.y, -1.0);
    float  Time = 0.0;

    RayDir = normalize(RayDir);

    TraceGeodesic(RayPos, RayDir, Time);

    fragColor = vec4(texture(iChannel0, RayDir).rgb, 1.0);
}

You can check out this Shadertoy implementation to see some of the optimizations, like variable timestep, replacing Hamiltonians with Lagrangians, using a symmetric matrix inversion function (a bit faster), reusing some of the computed values (restart if the Shadertoy is black):

mat4 diag(vec4 a)
{
    return mat4(a.x,0,0,0,
                0,a.y,0,0,
                0,0,a.z,0,
                0,0,0,a.w);
}

mat4 Metric(vec4 x)
{
  //Kerr-Newman metric in Kerr-Schild coordinates 
  const float a = 0.8;
  const float m = 1.0;
  const float Q = 0.0;
  vec3 p = x.yzw;
  float rho = dot(p,p) - a*a;
  float r2 = 0.5*(rho + sqrt(rho*rho + 4.0*a*a*p.z*p.z));
  float r = sqrt(r2);
  vec4 k = vec4(1, (r*p.x + a*p.y)/(r2 + a*a), (r*p.y - a*p.x)/(r2 + a*a), p.z/r);
  float f = r2*(2.0*m*r - Q*Q)/(r2*r2 + a*a*p.z*p.z);
  return f*mat4(k.x*k, k.y*k, k.z*k, k.w*k)+diag(vec4(-1,1,1,1));    
}

This is, of course, not the limit for optimization. The main other optimization is computing analytical inverses of the metric. For a large class of metrics you can use the Sherman–Morrison formula [9]

Some useful things to keep in mind

Since we used numerical finite differences, the results can actually depend quite a lot on the relative values of the float numbers. For example, the accuracy of the numerical derivatives is a lot lower far from the coordinate system center, so you'd need to vary the size of eps to avoid excessive numerical noise. Also, metrics quite often have numerical singularities which you should avoid, unless you want to get NaN results.

I usually avoid metrics in spherical coordinates due to their polar axis singularity, which has strong visual effects which are extremely hard to avoid even with a tiny varying timestep, although such metrics are usually mathematically simpler and allow for larger timesteps without breaking the look of the Black hole. For spherically symmetric metrics, like non-spinning black holes and wormholes there is a trick to avoid the polar singularity altogether! The thing about spherical symmetry is that the geodesic is always moving inside a 2d plane, which can be mapped to the equatorial plane of the coordinate system, basically reducing the 3d + time problem to 2d + time. (Scott Manley has a video explaining how he rendered wormholes, in there he used this trick + precomputing a lookup table to simplify the computation by a lot)

I've also used the dimensionality reduction trick in my shadertoy wormhole:

And finally you could always derive the equations analytically, while this is the most annoying method it is usually the fastest performance-wise. A compromise solution would be derive the equations automatically, this approach is used by geodesic_raytracing made by James Berrow (you should follow him on Twitter, he has a lot of cool stuff on this topic).

Figuring out if the ray has fallen inside the event horizon is actually not trivial, and there is no universal method, and while you could just set the color to 0 if the ray is below the event horizon surface, this is incorrect when viewing things from inside the black hole. Tracing the rays should also be done backwards in time, since we trace the rays from the camera, not to the camera, this has a noticeble effect on the resulting render, if not done this also results in completely dark renders inside of black holes, even though light does exist under the event horizon, and can reach from the outside.

Redshift in General Relativity can be computed simply from the ratio of the dot products of the velocity of the object at emission/absorption times the momentum of the photon, where the momentum of the photon is parallel to the photons direction of movement. The momentum needs to be parallel transported along the geodesic from emission to absorption, and the good news is that we can just use the geodesic 4-velocity instead, since it also is technically parallel transported along the geodesic and is pointed in the direction of movement. (Just in case, this also means we need to avoid renormalizing the generalized momentum when integrating the equations)

Also notably, the dot product in General Relativity is simply defined from the metric tensor \( g(u,v) \)

\begin{equation}
u \cdot v = g_{\mu \nu} u^\mu v^\nu
\end{equation}

Finally, if you only want to compute geodesics for Schwarzschild black holes, you can simply use the equation for a particle with mass 1 in a certain classical force field, details are in
this blog post.

Conclusions

Kerr black holes

Alcubierre warp fields and wormholes

Volumetric accretion disks around a Kerr black hole

Fast volumetric ray tracing with geodesics is quite difficult, and we needed to separate the ray marching loop into 2 loops, main loop being the geodesic steps, and the second loop being the volumetric substeps. Since we also use blue noise, it was necessary to keep the steps uniform along the geodesic, otherwise there would be clear artifacts in the volume, which required a few tricks with having a variable number of substeps per geodesic step.

Combining this with SDF's is somewhat easier, you need to vary the geodesic step to be the min() between the current step size and the SDF. Using this I've also tried to make a really simple path tracer in Unity with a Kerr black hole, naturally it was quite slow.

Note that rendering of moving objects is waaay harder, and requires either to have a space-time SDF, or some insane acceleration structure for triangles. And on top of that the entire history of the scene's past needs to be kept in memory, the only simple cases are when the moving objects can be represented as analytical functions you can sample in space and time, like the volumetric accretion disk in SpaceEngine.

References

• [1] Einstein notation
• [2] Equivalence of squared Lagrangian to Lagrangian
• [3] Euler-Lagrange equations
• [4] Euler-Lagrange equations derivation
• [5] Hamiltonian equations derivation
• [6] Legendre Transform
• [7] Alcubierre metric
• [8] Kerr–Newman metric in Kerr-Schild coordinates
• [9] Sherman–Morrison formula
• [10] Space-time path tracer

Author: FastFourierTransform
Original post on the forum: link

I want to try to address this. As Watsisname said; establishing a non-inertial co-rotating frame of reference has many advantages in terms of understanding the actual situation but it also has some disadvantages (non-inertial forces as Coriolis and centrifugal have to be invoked and even then are difficult or impossible to represent in a potential dwell map). Inertial frames of reference are great because of one thing: they are what we are used to picture in our minds frequently and they are the ones where physical laws are coherent without any extra tweaking needed (no fictitious forces needed). The explanation from the point of view of an inertial observer should be equivalent to Watsisname's but it should be given in other terms. I also feel that the actual inertial perspective on the subject should be extremely complex to wrap our heads around without doing many alienating calculations.

I haven't been able to come with a simple discourse for the stability of L4 and L5 nor for Trojan orbits, but I think the fact that Lagrange points exist as "dynamical" equilibrium points is feasible within an inertial framework, so here it goes:
Since we are not co-rotating is useful to remember that Lagrange points are not fixed in space and move themselves. Also very important to remember is Kepler's third law, which here will be expressed by the fact that objects farther from the Sun have lower orbital speeds (the angular speed will be dependent as the -3/2 power of the distance to the Sun). We place ourselves at rest relative to the Sun-Earth barycenter.

So. Lets say we want our probe to orbit the Sun but always been at the same distance from Earth without any additional effort, profiting from the gravitational landscape of the Solar System to accomplish that. If the Earth was a massless object then any point in Earth's orbit would be a suitable place to put our probe since it would travel just as Earth does but leading the march or lagging always the same constant amount. But the Earth has mass and some influence so your probe would be pulled back or ahead by the planet and it would be unable to maintain its heliocentric Earth-like orbit. This is interesting because Watsisname did the reverse, he started with a huge mass and then slowly went to the small Earth mass to show how nearly the entire orbit of Earth was covered with a low-inclination plain potential field. If we went a step further, to a massless Earth, we indeed would have equilibrium points all across Earth's orbit.

L1 point

So the L1 point is the easiest to understand. But first you have to abandon the misconception that the gravitational tug of the Sun equals that of Earth in that point. This is not at all about that (in fact as we saw in previous posts, even the Moon is farther from Earth than the equal-force point). So no, L1 is not the place where the Earth screens the Sun's gravity and you can move freely as if there was no Sun nor Earth until they move and lose that particular vantage point, no! L1 is located closer to the Sun than Earth so an object at that distance in a circular heliocentric orbit should get ahead of the Earth since it moves faster in its orbit, but the Earth here is creating his magic; having the Earth right in the opposite direction of the Sun what your probe feels is that the Sun-Earth system is equivalent to a Sun without Earth but this Sun having a little less gravitational attraction. So the Earth is orbiting the Sun but our probe orbits as if it orbited a lighter Sun, this means that even if it is closer to it our probe can still move as slow as Earth due to the lesser influence of this lighter Sun. So the the probe moves at the same angular velocity around the Sun as the Earth because it feels the Sun less and it is closer (both effects balance out) so the L1 - Sun - Earth configuration is maintained constant. Bare in mind that the vector sum of the gravitational force of Earth and that of the Sun on an object placed here at L1 will have the same direction as to point to the Sun-Earth barycenter (so in fact L1 is orbiting the Sun-Earth barycenter also).

L2 point

L2 obeys a similar interaction. Instead of the tug of Earth screening a small part of the Sun's influence we get the Earth adding to the pull of the Sun (since L2 the Sun and the Earth are both in the same direction). A probe in L2 would, in principle, lag behind Earth because of Kepler's third law since it is farther from the Sun, but at the same time the Sun "feels stronger" (thanks to the extra pull of Earth) so you end up going faster than expected. Both effects balance out in the L2 point. And since you move at the same rate as the Earth here (same period or same angular velocity) you will always maintain the same configuration with respect to the Sun and to Earth and that in turn means that you will maintain this behaviour and the L2 point wouldn't go away from you (you will move with it). Again, as in L1, bare in mind that the vector sum of the gravitational force of Earth and that of the Sun on an object placed here (at L2) will have the same direction as to point to the Sun-Earth barycenter.

L3 point

L3 is similar to L2 but symmetrical to the Sun. At the opposite side of Earth's orbit you notice the pull of the Sun as if Earth nearly didn't existed (you are ~2 AU from it) but Earth still contributes to make you feel the combined pulls of the Earth and the Sun in the same direction as if the Sun had gained some mass and was a little stronger than the Sun we feel from Earth. This means that you can orbit this "more massive Sun" quicker if you are in Earth's orbit but if you get a little farther from the Sun than Earth's orbit you would reach a state of "dynamical equilibrium" again were you are moving at the same speed as Earth does. This is interesting because we often see pictures of L3 been located in Earth's orbital path but in fact it is outside of it, like L2. The difference with L2 is that for you to beat the angular speed increase of the strong gravitational pull of the Sun and Earth you (now in L3) just need to get a little farther but since in L2 the Earth is quite close the influence that has to be overcome is greater and thus you move considerably farther from Earth's orbit than in the case of L3 to archive the same balance. Again, as in L1 and L2, bare in mind that the vector sum of the gravitational force of Earth and that of the Sun on an object placed here (at L3) will have the same direction as to point to the Sun-Earth barycenter.

L4 and L5 points

Since the gravitational field of the Sun-Earth system is mirror-symmetric to the Sun-Earth line then both L4 and L5 points have the exact same underlying explanation. Remember that L1, L2 and L3 all were orbiting the Sun-Earth barycenter because the net force (the result of adding Earth's tug and the Sun's tug) always pointed there? Well this fact was simple to grasp because in the 3 cases we had only vectors in the Sun-Earth line (so the net result would always lie in the same line as the Sun's and Earth's pull). But what if we could have a net force still pointed to the Sun-Earth barycenter but with the Earth's pull and the Sun's not been parallel from our position? Well, that's exactly the L4 and L5 points. Here you will experience the Sun pulling in one direction strongly and the Earth pulling in an other (not parallel to that of the Sun) more weakly but enough to have your net force dragging you to the Sun-Earth barycenter instead of just the Sun. Since this geometrical situation can be archived in many places we have to also include the restriction that you maintain the same orbital speed as Earth (if not you would feel the net force pointing to the Sun-Earth barycenter just for a brief moment and then you would be pulled slowly to a different direction as the distance to the Sun or to the Earth would vary). With both restrictions (having your net force pointing to the barycenter and moving at the same orbital period as Earth as to maintain that net force pointing to the barycenter) you end up with only two places in the entire system, and those are the L4 and L5 positions.

L1 stability

This is difficult to address in an inertial frame but let's try. To understand if L1 is a stable or unstable equilibrium point we are going to imagine our probe placed in 4 different places (all very close to L1). Let's take a slightly higher orbit from L1 (this means that we are in a Sun - system barycenter - L1 - our probe - Earth - L2 configuration). In this situation we have a greater pull from Earth so our net force will be pointing the same way (to the barycenter) but with less intensity (since here the Earth is now screening more of the Sun's influence). Been in a higher orbit alone would mean we are moving slower due to Kepler's 3th law but also we move slower because of the "lighter Sun" we have now. Because of that our probe would lag behind Earth. As this situation continues to evolve the net force would start to have a component pointing outside of the barycenter of the system (as the angle between the pull from the Sun and the pull from Earth gets lower than 180º). This in principle should be nice because it should impart some tangential motion on our probe so that it reaches L1 again, but sadly we are in a to high orbit for this to work. Our probe will gain a little acceleration but would end of still losing the race to Earth. But! this could happen if instead of moving to a higher orbit from L1 we maintain the same orbit as L1 but place our probe behind L1 as it moves. Then the orbit would have a smaller circumference than in the previous case and this new tangential pull would be enough to make our probe reach L1. Once the probe reached L1 it would probably have to much speed to maintain that position so it would come ahead of L1 and the reverse would hold. The probe would be oscillating like in a valley in this case. Lastly let's think about the Sun - system barycenter - our probe - L1 - Earth - L2 configuration (our probe is in a closer orbit than L1). In this situation we are faster than Earth since the pull of Earth is not enough (since we are farther from it) to overcome the effect of a smaller quicker orbit around the Sun. The probe would surpass Earth and it would end pulling us backwards, but not enough strongly. As the tangential acceleration performs its magic the probe would slow down and lower the orbit even farther from L1 making the probe lose its place of near equilibrium.
What does all of this means? It means that any displacement from L1 performed in the Sun-Earth line would end up been amplified and any deviation in the perpendicular direction of that line will be corrected and the equilibrium restored. So here we have for L1 a saddle point as was expected (stable equilibrium in one direction but unstable in other direction).

L2 and L3 stability

Same thing as before but here the situation is more simple. Any small deviation from L2 in the Sun-Earth line can be explained in the same way for higher or lower orbits from L2; we are now considering a "heavier Sun" (due to Earth's contributing to the net pull in the same direction as the Sun). If we are in a lower orbit (with respect to L2) the probe will move faster and will overcome the Earth, but then a tangential force would try to slow down the probe because Earth would be a little behind it. If it slows down the orbit shrinks again so a positive feedback loop is generated where the probe ends been separated more and more from L2. For a higher orbit the the probe would lag behind Earth so a tangential pull would appear that would speed up the probe throwing it to an even higher orbit and escaping L2 more and more. So L2 is an unstable equilibrium point for any direction you choose to displace your probe from it. L3 behaves exactly in the same manner due to the same reasoning but since it is farther away from Earth the instability is less abrupt and we still have a lot of margin to move the spacecraft around L3 without loosing control so quickly (this is a peak with a less inclined slope).

L4 and L5 stability:

Here is where I get lost. The geometrical configuration here is way more complex than in the other cases since it ceases to be a 1-dimensional problem and becomes 2D. I don't know exactly hot to picture this in my mind but any small deviation from L4 has to be diminished since as Watsisname explained this is a stable equilibrium. My guess is that this problem is equivalent to the problem where we only have the Sun and the probe, but this Sun changes mass slightly and harmonically and wobbles from one place to the other in such a way that it pulls the probe as if it was Earth. An animation should make this more visually clear. Let's say, for example, that we are ahead of L4 (still behind Earth). In this situation Earth's pull will increase in magnitude and the direction of the net force will change slightly (from the barycenter towards the Sun). The speed then would increase so the probe would archive a higher orbit and slow down, but this somehow would happen in such a way that L4 would surpass our probe because of that slowdown and we would end up behind L4 quickly, where the reverse would hold. We would be somehow oscillating around L4. How is this the case for displacements from L4 in any direction I don't really know but sure an animated diagram with the forces involved would clarify this a lot (maybe one day I try to do it).

Author: Watsisname
Original post on the forum: link

Introduction

When two masses move in a circular orbit about each other, five points of gravitational equilibrium are generated, which move along with them.  A particle placed at one of these points will have a very small acceleration, as if gravity was in some way cancelled out.  More strangely still, it can be possible to place a particle in a stable orbit around some of these points, despite there being no mass there to attract it!  

These five special locations are known as the Lagrange Points, discovered by Leonhard Euler and Joseph-Louis Lagrange in the 18th century.  Not only are they of special interest in celestial mechanics, they also have important applications for space missions and the placing of satellites.  We even find asteroids caught around some of the Lagrange points: so called "Trojan asteroids".

Why do these points exist?  To investigate, we will again use the concept of a potential, as we saw in the discussion of Mercury's orbital precession.  As a reminder, in physics a "potential" is an imagined landscape of hills and valleys that describes how an object reacts to forces.  We can think of an orbit as being the motion of a particle trapped in a valley of a potential.  This time however, we will generate the potential in a somewhat unusual way: we will adopt a frame of reference which is rotating.  There are pros and cons to this approach.  One one hand, motion in a rotating frame of reference is not very intuitive, and we will need to invoke "non-inertial" forces to understand them.  On the other hand, doing it this way will make life a lot simpler, as then the masses and the Lagrange points themselves will appear stationary.  Motions around the Lagrange points will also be much easier to visualize and interpret.

But, seeing as how much I myself struggled when I was learning about non-inertial forces in rotating frames of reference, I think a quick digression into that subject would be a good idea.

Frames of Reference and Non-Inertial Forces

Most of us are comfortable with physics in "inertial" frames.  An inertial frame is one in which objects subjected to no forces will move in straight lines with constant speed due to inertia.  In other words, an inertial frame is one in which the "law of inertia", or Newton's First law, is valid.  But what if we put ourselves on a turn table, and try rolling a ball back and forth to each another?  We will discover it is quite hard -- the ball will seem to be deflected by some mysterious force.  

Someone watching from the side will still say the ball moves in a straight line relative to the ground.  It is us on the spinning turn-table that accounts for the weird trajectory we see.  However, we can be just as valid in saying that the ball is deflected by something.  We may invoke new forces, unique to our rotating frame.  There is a centrifugal force which tends to fling things away from the axis of rotation, and a Coriolis force which tends to make moving objects turn and follow curling paths.  We all have felt centrifugal force whenever riding in a car going around a turn.  And many are probably familiar with how Coriolis forces on the spinning Earth causes air currents to deflect, and explains the rotation of hurricanes.  These forces are sometimes said to be "fictitious", in the sense that there is no real physical interaction responsible for them.  They are simply consequences of being in a strange frame of reference.  But when we are in that reference frame, they seem totally real.

I have found no more wonderful a demonstration and explanation of these fictitious or non-inertial forces than in this film put together by University of Toronto physicists in 1960.  If you have not seen this, it is an absolute must watch.  Or if you're in a hurry:

-The first 13 minutes build the ideas of different frames of reference.
-13:25 to 17:00 looks at non-inertial (accelerated) frames.  
-17:04 is a beautiful example of motion in a rotating frame.  This segment is extremely relevant to the physics we will be looking at here.

The Generalized Potential for Two Masses

Here we go!  Let's imagine two masses in orbit around each other, and examine the forces and the potential they produce.  For starters, and for symmetry, let's first let the two masses be equal.  We will also "lock them in place", so that (for now) there are no non-inertial forces.  What does the potential look like for this system?

For any point in space there will be a gravitational force.  Actually two of them: one from each of the masses.  The gravitational force from each is given by

where r with a hat over it means "in the direction from the mass to the point in question", since force is a vector, like an arrow, which has magnitude and direction.  The minus sign in front flips it around so that it points to the mass.

In terms of the potential (the landscape whose steepness describes the force), this becomes

This expression lacks any vectors, since each point is assigned only a magnitude (a height in the landscape).  What does this landscape look like?  Below I plot it for the two equal masses, 1 unit apart along the x axis, fixed in place and centered around the origin.  Darker/redder colors indicate more negative values of the potential.  Around each mass is a bottomless hole.  Drop a ball into this landscape, and it will fall into the hole associated with the nearest mass.

As FFT has said earlier, the steepness of the potential (how closely packed the contour lines are) relates to the magnitude of the force at that point.  On a topographic map, close lines means steep terrain.  A ball on steeper terrain will accelerate down the hill faster, and likewise a particle on a steep region of a potential will experience greater force.

To help see how the potential and the force/acceleration are connected, below I'll plot the exact same situation in terms of the acceleration. I'll plot on a log scale, so that details at very small accelerations are visible (whereas on a linear plot they would be overwhelmed by the very large accelerations close to the masses).  Here, darker blues represent greater acceleration.

Unsurprisingly, near the origin the acceleration drops to zero.  There the gravitational force due to the mass at the left exactly cancels the force due to the mass on the right.  We can also see a hint of this from the potential.  This location is a "saddle point", where the landscape rises up if we move along the y axis, but drops down if we move along the x axis.  Trying to balance an object here would be very difficult, with a tendency for it to "roll off".

The Co-Rotating Frame: how do non-inertial forces modify the potential?

Now let's set these two equal masses into a circular orbit around each other, and place ourselves into a rotating frame of reference, such that the masses appear to stay still.  What new forces appear, and how do they affect the potential and the motions?

The first new force this will introduce is centrifugal.  This force acts to fling things away from the center of rotation of the system, about which our frame is rotating with angular velocity Ω (meaning the speed at which the angle is changing).  We can find the angular velocity of the frame by using Kepler's Third Law:

where R is the separation distance between the two masses M₁ and M_2.  With the frame rotating at this angular speed, the centrifugal force is

where r is the distance from the point of interest to the center of mass of the system about which the frame is rotating.  The r-hat in this equation tells us the centrifugal force is directed away from the axis of rotation, as it should be.  The strength of the force also grows in proportion to the distance from the center of the frame, and also with the square of the angular speed that the frame is rotating.  

In terms of the potential, the centrifugal force adds a new term to it, with magnitude

Let's see what this looks like now, again for the two equal masses.  Now the masses are orbiting in a circle, and we are in a frame which rotates with them.

There they are!  Five points where the net force is zero!  One is in the exact center as before, where the gravitational forces cancel and the centrifugal force is zero since it is at the center of rotation.  This is the L1 Lagrange point, in the case of the two masses being equal.  

Two new points are added on either side of the masses.  These are L2 and L3.  Here the centrifugal force is exactly equal and opposite to the sum of the gravitational forces.  In the potential they also look like saddles, and are unstable.

Finally, two new points are added above and below.  These are L4 and L5, where again the centrifugal and combined gravitational forces cancel.

You might notice that the inclusion of centrifugal force makes the potential drop down again at large distances, whereas with gravitation alone, it rose upward.  This is because if you start out far away from the center, and match your speed to rotate along with the system, then you are moving quite fast, and have far too much angular momentum to be in a circular orbit.  Your side-ways speed will carry you off to larger distances.  On the other hand, if you start out too close, then you have too little speed to be in orbit and so you drop closer in.

Let's do something fun!  Earlier I said that L4 and L5 are stable points about which we can place something in orbit.  So let me just drop in a particle very close to the L4 point, and....

Wait, what?  Why did it spiral away?

Well, I had lied a little bit.  L4 and L5 are not always stable.  And you might have noticed that the potential at these points looks like the top of a hill.  Not the bottom of a valley! (Lighter colors instead of darker).  How could there possibly be stable orbits around them?  (Actually I did not understand this until now.  I had suffered the misconception that L4 and L5 are shaped like valleys, until I started working through this in detail!)

What is going on?

L4 and L5 are the tops of hills, not the bottoms of valleys.  So of course an object placed near them will slide away.  But once it does, a new non-inertial or "fictitious" force takes hold: 

The Coriolis Force:  

This is the force which causes the deflection of the puck seen in the demonstration in the above video with the turn table.  It is a velocity dependent force, and it always acts at right angles to the velocity.  This means it will cause a moving object to turn, but not change speed.  The Coriolis force is given by

The cross (X) product between Ω and v means that the direction of the force is given by the right hand rule, and then we flip the direction 180° because of the minus sign in front.  So if the frame of reference is rotating counter-clockwise, then the force on a particle moving upward will be to the right.  That is, it will deflect trajectories clockwise.  This is also why high pressure systems (air flow away from center) spin clockwise in the northern hemisphere, and why low pressure systems (air flow toward the center) spin counter-clockwise.

Because the Coriolis force is velocity dependent, it is impossible to define a potential for it.  If we move faster, the force will be stronger, and we'll take a different path, with a different change in potential energy.  So there is not a unique landscape that we can use to model motions subject to this force.  Instead, we can model the general potential with just gravitation and centrifugal force, and to simulate motions we can include the effect of Coriolis force by direct computation on the computer.  This is how I generated the above plot.

In the case of the masses of M₁ and M₂ being equal, the particle dropped near L4 slides down the potential hill, and then as it speeds up, the Coriolis force whisks it around in a growing spiral.  There are no stable orbits around any of the Lagrange points in this scenario.  

In order for L4 and L5 to be stable, one of the masses must be at least 24.96 times the mass of the other.  Understanding why this is the case (or why this number in particular) is hard, and requires a careful analysis of the equations of motion.  If any readers are interested in the mathematical details, I can refer you to some suggested reading at the end.  But for our purposes, let's just take it as an established fact.  What we can do is look at more examples of the potential where we vary the mass ratio, check what happens visually, and see if we can qualitatively make sense of it.

The generalized potential for a 10:1 mass ratio:

Let's set the mass of one object to be 1/10th the mass of the other.

The acceleration plot is particularly illuminating.  The Lagrange points are starting to look a lot like their commonly-shown configuration, with L4 and L5 leading/lagging the secondary mass by 60°.  This 10:1 ratio is often used in figures that portray the Lagrange points for the Sun and Earth.  Which is an exaggeration of truth, because Earth is of course much less than 1/10th the mass of the Sun!  But it strikes a nice balance, giving an idea of the shape of the potential for the Sun and planet.  Using the true mass ratio of Sun and Earth would be quite unhelpful for visualizing it, as the Sun's influence dominates the landscape.

While we're at the 10:1 ratio, let's drop in another particle near L4:

Once again it spirals away.  10:1 is still too small of a ratio of masses, when we need ~25:1 at least.  But we can see that we're going in the right direction!  The spiral is a bit tighter, and the particle wraps around the L4 point completely before dropping down and hitting the smaller mass.  (Oops.)  

Finally, let's increase the ratio to match the Sun and Jupiter.  That is 1047:1 -- definitely enough for L4 and L5 to be stable.  And since this is a real situation with specific masses and distances, I'll add in units.  The potential is in terms of energy, here shown in gigajoules per kilogram.  Distances are in AU.

It's pretty hard to see the details here (again why most figures online show 10:1 instead of the real thing).  Even the log-scale acceleration plot isn't terribly useful.  But we can at least see that the potential in this case has a crescent ridge along Jupiter's orbit, with Jupiter making a small depression around itself.  The L1, L2, and L3 points still exist and are saddle-shaped, but they are not visible in this plot scale.  Similarly, L4 and L5 are still separate hills along the crescent, but we can't distinguish them from one another.

Let's zoom in to the area around L4, and drop in a particle again.

Here I said "WOW!"  I had no idea that the orbits around Lagrange points look like this.  But they do -- spiraling about in these crazy elongated paths.  I was extremely surprised and excited when the computer output this trajectory.  These orbits are known as "Trojan orbits".  As mentioned before, they have enormous astrophysical significance, and in fact many asteroids are caught in them, especially around Jupiter's L4 and L5 points.  One has also been found in a Trojan orbit of Earth!

Here is an even larger one (starting farther away from the L4 point), with an overview including the Sun and Jupiter but without the mess of the potential.  

As a sanity check, I went to the literature. Check out Figure 3 from this paper "Horseshoe and Trojan orbits associated with Jupiter and Saturn", by Everhart, 1973.  A very close match!

How can we explain this behaviour, qualitatively?

As before, when the particle slides down the potential hill from the Lagrange point, Coriolis force whisks it aside, turning it clockwise.  But soon it is turned so much that it rises back up the potential hill, at a location farther along the orbit -- either closer to or farther from Jupiter.  As it rises up the hill it slows down again, the Coriolis force weakens, and the arc shrinks.  Then the process repeats -- falling down the well, getting deflected to the side by Coriolis, and so on.

But the shape of this potential hill is a like a very broad crescent.  Farther away from the Lagrange point, it gets thinner, and a little less high.  This means the particle has an easier and easier time rising back up to the top of the "ridge".  Finally it reaches a point where it is able to cross over the ridge, and the process repeats on the other side!  It is the Coriolis force that keeps the particle circulating around near the top of this potential hill -- with the hill itself being generated by gravity and centrifugal force.

For larger deviations from the L4/L5 points, the particle will spiral along more and more elongated paths.  Once the departure gets too large, it can even cross over to the other side of the Sun-Jupiter line.  Now we have Horseshoe orbits!   Some asteroids are known to be caught in these orbits as well.

One other thing to add real quick, and then this will be the last for a while.  I wanted to show that it is also possible to orbit around the unstable Lagrange points (L1 for example), but the orbit is extremely delicate and will fall apart if not carefully maintained.

Here it is For Jupiter's L1 point.  The particle starts out a little bit closer to Jupiter than the L1 point, and with an initial velocity downward.  Coriolis force deflects its path clockwise, and with the initial position and velocity tuned "just right", it ends up circulating in an oval.  For the values I chose, it completed about an orbit and a half before falling apart.

To give an idea of the timing involved, this path is plotted over 12.5 years, and the first orbit took about 5.5 years.  Plenty of time for a spacecraft trying to follow this course to make necessary adjustments to stay on it.  This has actually been done with missions to Earth's L1 and L2 points.  Here is the path the LISA Pathfinder used to get to and orbit around the Earth-Sun L1 point.  Again some maintenance is required to follow this orbit.

Concluding Remarks & Suggested Reading

I think I will end it here, as I think this has covered all (and more!) that I had hoped to go through.  I myself have learned an enormous amount while working through this, and came across many surprises and previous misconceptions.  I did not know for example that L4 and L5 are actually repulsive, and that the stability of orbits around them depends on the Coriolis force.  Nor did I really have a sense of what Trojan orbits really look like or why they behave as they do.  This was extremely exciting to work through and discover, and I hope you all find this an interesting and enjoyable exercise as I did.

Cheers!

References & Further reading:

An excellent, more technical treatment of Lagrange points, including deriving them and calculating their stability, can be found here. (written by Neil Cornish of Montana State University).

Another good one can be found here from Tipler's Modern Physics text.

For those interested in a broader scale of classical or celestial mechanics, including deriving equations of motion in non-inertial frames of reference, I can highly recommend the text "Classical Mechanics" by John Taylor.

Background

Often in discussions of Mercury's perihelion advance, the effect is shown greatly exaggerated.  Sometimes even the orbit is shown much more eccentric than it really is.  Before we dive in, let's take a moment to get some perspective on the true shape of the orbit and appreciate the scale of the solar system.

 
All orbits are ellipses, but most planetary orbits are close enough to circular that it can be hard to tell at a glance.  Mercury is the biggest exception with an eccentricity of 0.21.  (Or Pluto at 0.25).  What's interesting, and what this post is all about, is that Mercury's orbit is not stationary.  The ellipse slowly shifts around, it's perihelion point advancing by about 574 arcseconds (or 0.159 degrees) per century.  Much of this (531 arcseconds) can be explained by the perturbations from the other planets.  However, the remaining 43 arcseconds per century are "anomalous", unaccounted for by Newtonian mechanics.  

The puzzle of this additional observed precession was first noted by LeVerrier in 1859.  Following the successful prediction of the existence and location of Neptune from similar perturbations on Uranus' orbit, LeVerrier proposed that Mercury's orbital precession could be caused by another undiscovered planet inside of Mercury's orbit, which he named Vulcan.  Astronomers searched for this missing inner planet during solar eclipses.  Of course, Vulcan was never found... 

The mystery was not solved until Einstein began developing his general theory of relativity.  When he applied his general relativistic equations to the problem in 1915, he found that they exactly predicted the additional 43 arcseconds per century.  He considered this one of his greatest achievements.

Imagine my joy at the recognition of the feasibility of general covariance and at the result that the equations correctly yield the perihelion motion of Mercury. I was beside myself for several days in joyous excitement.

Physics

So, why does it happen? Is there an extra energy involved?  (Spoiler: Yes, in a way!)

For an elliptical orbit, we can think of the orbital motion as being made of two parts: an oscillation in azimuth (angle around the Sun), and an oscillation radially (in and out).  In Newtonian mechanics, the periods of these two oscillations are exactly equal.  That is, each time the planet completes one cycle around the Sun, it also exactly completes one cycle radially in and out.  Therefore in Newtonian mechanics there is no precession of the orbit (besides that which is caused by the influences of the other planets).

In general relativity this is no longer true.  The periods of the two oscillations are different!  To see why, it is easiest to examine the shape of the "effective potential", which you can think of as being just like the potential energy well around a spherical mass, except also taking into account the angular momentum of the orbiting body.  

If you are unfamiliar with the concept of a potential energy well, imagine a rolling landscape of hills and valleys.  If you place a ball at the bottom of a valley, it will just sit there.  It's in a stable equilibrium.  On the other hand, the top of a hill is an unstable equilibrium.  If you displace the ball slightly from the top of the hill and let it go, it will continue to roll down, exchanging gravitational potential energy for kinetic energy as it drops in height.  Finally, if you have some bowl-shaped valley and drop the ball from rest at some initial height, then it will roll down to the bottom, and (ignoring friction) keep going, rolling back up the other side until it reaches its former height, and then roll back again, oscillating back and forth indefinitely.  

Because the gravitational potential energy near Earth's surface is simply proportional to your altitude, if you make a plot of the potential energy of a ball in this landscape as a function of its position, it will be exactly the shape of the landscape.  This form of modelling of the potential is extremely useful.  In general, we can use "potential energy wells" as a way to understand motions of objects subjected to different kinds of attractive and repulsive forces, from balls rolling down hills to the vibrations of atoms bound together in molecules.

An elliptical orbit is an oscillation in an effective potential well, as well.  There is a hill to climb as you move outward due to the gravitational attraction of the star, and there is also a hill to climb as you move inward!  This is because the sideways motion is associated with angular momentum, and angular momentum acts like a repulsive centrifugal force.  By conservation of angular momentum, the sideways velocity increases as your orbit swings inward, increasing the centrifugal repulsion and making it harder to approach the center.  If the angular momentum is large enough, then the inward fall can be halted, and you swing back out again, making an orbit!  

So if we can compute the effective potential well for a planet orbiting the Sun, we can figure out some things about its motion.  The anomalous precession of an orbit can also be understood this way, by comparing the motion from the Newtonian potential with what happens in the general relativistic version of the potential.  Let's try it!

Math

The goal here is to try to understand how things work conceptually, and I will try to emphasize those concepts (especially physical concepts) as we go along.  The hard truth is that much of it will still require going through the math in order to access it, but I'll do my best to turn that math into something visualizable and comprehensible.

In Newtonian mechanics, the gravitational potential energy of a mass m a distance r around a spherical mass M is

This describes a simple curve which plummets downward indefinitely as you move to smaller radii.  As it should.  If you drop something near a massive object and give it no sideways motion (no angular momentum), then it will fall straight into it.

If we instead give the object some angular momentum L, then the potential is modified:

where μ is the reduced mass, which for a planet orbiting the Sun (m << M), reduces to approximately m.

Let's see what this looks like for Mercury:

Again the way to think about this is that the potential describes a landscape that a ball will frictionlessly roll across.  To visualize that, I added the large red dot to represent Mercury's position as a function of time.  The acceleration is determined by the slope of the potential, and I plotted the motion with 2 days per frame.  With 44 frames in all, this traces out one complete period of Mercury's orbit in the radial motion.  I also added a horizontal line to show the total energy of Mercury (its gravitational potential energy + kinetic), which is constant along the orbit.  Where this line is above the effective potential defines the range of distances from the Sun that Mercury's orbit will cover.  The vertical lines represent the extremes (Mercury's perihelion and aphelion distances).

Why is the Newtonian effective potential shaped this way?  
Because of Mercury's amount of angular momentum, the shape of the effective potential that it sees near the Sun is a valley with hills on either side.  The hill at large radii is due to the Sun's gravitational attraction, while the hill at small radii is due to centrifugal repulsion.  Even though the gravitational force grows stronger at closer distances, for an orbit the centrifugal force gets stronger more quickly, due to the conservation of angular momentum which increases your sideways speed.

A useful concept to keep in mind here is that the period of the radial oscillation depends on the curvature of the potential well.  If the well opens up more sharply, then the average acceleration is greater, and the period is shorter. 

And here's one other useful trick.  If the amplitude of the oscillation is not too large (does not reach too far away from the minimum of the well), then we can approximate the well as a parabola around the minimum.  Then for a parabolic well the motion is described very simply as a simple harmonic oscillator.  For Mercury this approximation isn't particularly good (the well is noticeably asymmetric over the region Mercury covers), but it's not terrible either.  

For a simple harmonic oscillator, the frequency is given by

which for the radial motion leads to

The frequency for the angular oscillation on the other hand is given by

Which leads to the exact same expression:

So we see in Newtonian mechanics the periods are exactly equal and there is no anomalous precession.  Now let's see how this gets modified when we move to General Relativity.

The Effective Potential in General Relativity

Mercury's Orbit in General Relativity

Going Further: The Bizarre Orbits in Strongly Curved Space-Time

Don't miss discussion of the original post on the forum: it has a lot of interesting questions and answers.