Numerical integration of light paths in a Schwarzschild metric
Differential equations of orbit
The Schwarzschild metric is one of the most famous solutions to the Einstein field equations, and the line element in this metric (in natural units ) is given by:
We are interested in the trajectory of a light ray in such a metric. Since the metric is spherically symmetric, any light ray that starts with a certain must stay in the same plane, hence we can arbitrarily set and do away with all the terms.
Light follows a null (lightlike) trajectory given by . In the absence of external forces, it should also travel along a geodesic. These are governed by the geodesic equations, which can be derived using Euler-Lagrange equations. Due to the symmetry of the metric, applying the Euler-Lagrange equations to the metric gives us two conserved quantities:
where refers to derivative with respect to an affine parameter.
Using the null condition, we have
This can be expressed in terms of , and differentiating again gives the second-order differential equation for :
This can be easily converted into a first-order differential equation to be solved numerically by setting a variable . So we have these 3 differential equations to compute numerically:
(This can of course be solved analytically in the weak gravity limit, which gives the light bending equation .)
In principle, we need the initial values of , , and to start the numerical simulation. However, if we fix the incoming velocity to be horizontal, then we would only need to specify the initial and coordinates.
The initial conditions then can be given as follows:
Then, the only free parameters to specify and , in addition to mass.
For a mass of (corresponding to Schwarzschild black hole radius of 2), this is a plot of the trajectories with different impact parameters :
And they do fit quite well with the theoretical deflection angle, for large impact parameters: