Using matplotlib, I am path tracing some 2D photons bending due to a non-rotating standard black hole defined by Schwarzschild metric. I set my initial velocities in terms of r (radius), phi (angle), and t (time) with respect to the affine parameter lambda and then iteratively update the space-time vector based on the Christoffel symbols at the respective point.

The main concerns are is why the null geodesic statements don't print something closer to zero.

Everything scaled down by 1e6 to plot

import numpy as np
import matplotlib.pyplot as plt
from math import sin, cos, sqrt, atan2, pi

# Constants
M = 9e32  # Mass of the black hole in kg
c = 299792458  # Speed of light in m/s
G = 6.6743e-11  # Gravitational constant in N m^2/kg^2

# Simulation parameters
curve_num = 2000  # Resolution of the curve
d_lambda = 1e-4  # Step size for integration

# Black hole parameters
horizon = 2 * G * M / c**2  # Event horizon radius
photon_sphere = 3 * G * M / c**2  # Photon sphere radius

# Test point parameters
num = 30  # Number of photons
x_i = 7.5
y_i = 2.5
y_f = 4

# Initial velocity in Cartesian coordinates
v_x = -c  # Speed of light in negative x-direction
v_y = 0   # No velocity in the y-direction

# Prepare for matplotlib plotting
trajectories = []

for i in range(num):
    trajectory = []

    # Define initial test point
    if num > 1:
        test_point = np.array([x_i, y_i + i / (num - 1) * (y_f - y_i), 0]) # linear interpolation
    else:
        test_point = np.array([x_i, y_i, 0])

    # Convert to polar coordinates
    r = np.linalg.norm(test_point) * 1e6  # Radius in meters
    p = atan2(test_point[1], test_point[0])  # Initial planar angle

    # Metric coefficients
    g_tt = -(1 - 2 * G * M / (r * c**2))
    g_rr = 1 / (1 - 2 * G * M / (r * c**2))
    g_pp = r**2

    # Initial velocities
    v_r = (v_x * cos(p) + v_y * sin(p))  # Radial velocity
    v_p = (-v_x * sin(p) + v_y * cos(p)) / r  # Angular velocity
    v_t = sqrt(-(g_rr * v_r**2 + g_pp * v_p**2) / g_tt)

    # Check the null geodesic condition for the initial point
    #print(g_tt * v_t**2 + g_rr * v_r**2 + g_pp * v_p**2)

    # Integrate geodesics
    for j in range(curve_num):
        if r > horizon + 10000:
            # Precompute common terms
            term1 = G * M / (r**2 * c**2)  # Common term: GM / r^2c^2
            term2 = 1 - 2 * G * M / (r * c**2)  # Common term: 1 - 2GM / rc^2

            # Christoffel symbols using common terms
            Γ_r_tt = term1 * term2
            Γ_r_rr = -term1 / term2
            Γ_r_pp = -r * term2
            Γ_p_rp = 1 / r
            #Γ_t_rt = term1 / term2 # ignoring time marching

            # Update change in velocities
            dv_r = -Γ_r_tt * v_t**2 - Γ_r_rr * v_r**2 - Γ_r_pp * v_p**2
            dv_p = -2 * Γ_p_rp * v_r * v_p
            #dv_t = -2 * Γ_t_rt * v_r * v_t # ignoring time marching

            # Update velocities
            v_r += dv_r * d_lambda
            v_p += dv_p * d_lambda
            #v_t += dv_t * d_lambda # ignoring time marching

            # Update positions
            r += v_r * d_lambda
            p += v_p * d_lambda

            # Metric tensor components (update for new r)
            g_tt = -term2
            g_rr = 1 / term2
            g_pp = r**2

            # Recalculate v_t from the metric components
            v_t = sqrt(-(g_rr * v_r**2 + g_pp * v_p**2) / g_tt)

            # Check the null geodesic condition at each step
            #print(g_tt * v_t**2 + g_rr * v_r**2 + g_pp * v_p**2)

            # Store Cartesian coordinates
            x = (r / 1e6) * cos(p)
            y = (r / 1e6) * sin(p)

            # Only store points within the -10 to 10 range
            if -10 <= x <= 10 and -10 <= y <= 10:
                trajectory.append((x, y))
        else:
            break

    trajectories.append(trajectory)

# Plot using matplotlib
plt.figure(figsize=(8, 8))

# Plot each trajectory
for trajectory in trajectories:
    trajectory = np.array(trajectory)
    if len(trajectory) > 0:  # Plot only if there are points
        plt.plot(trajectory[:, 0], trajectory[:, 1])

# Plot the event horizon
circle = plt.Circle((0, 0), horizon / 1e6, color='black', fill=True, label="Event Horizon")
plt.gca().add_artist(circle)

# Plot the photon sphere
photon_sphere_circle = plt.Circle((0, 0), photon_sphere / 1e6, color='red', fill=False, linestyle='--', linewidth=1.5, label="Photon Sphere")
plt.gca().add_artist(photon_sphere_circle)

# Set plot limits explicitly
plt.xlim(-10, 10)
plt.ylim(-10, 10)

# Configure plot
plt.title("Photon Trajectories Around a Black Hole")
plt.xlabel("x (scaled by 1e6 meters)")
plt.ylabel("y (scaled by 1e6 meters)")
plt.axhline(0, color="gray", linewidth=0.5)
plt.axvline(0, color="gray", linewidth=0.5)
plt.axis('equal')
plt.grid()
plt.legend()
plt.show()

I expect my null geodesic print statement inside the for loop to be either 0 or a relatively small integer. For example, the first null geodesic print statement may be 0, 16, -8, etc. due to a small amount of imprecision with the large float addition (e17 magnitudes).

I have tried to debug by replacing my "else" statement with "elif j == 1" and looking at the very next iteration, it can be seen the null geodesic prints a much larger float.

I believe figuring out the null geodesic error will reveal if my trajectories are incorrect.

Using matplotlib, I am path tracing some 2D photons bending due to a non-rotating standard black hole defined by Schwarzschild metric. I set my initial velocities in terms of r (radius), phi (angle), and t (time) with respect to the affine parameter lambda and then iteratively update the space-time vector based on the Christoffel symbols at the respective point.

The main concerns are is why the null geodesic statements don't print something closer to zero.

Everything scaled down by 1e6 to plot

import numpy as np
import matplotlib.pyplot as plt
from math import sin, cos, sqrt, atan2, pi

# Constants
M = 9e32  # Mass of the black hole in kg
c = 299792458  # Speed of light in m/s
G = 6.6743e-11  # Gravitational constant in N m^2/kg^2

# Simulation parameters
curve_num = 2000  # Resolution of the curve
d_lambda = 1e-4  # Step size for integration

# Black hole parameters
horizon = 2 * G * M / c**2  # Event horizon radius
photon_sphere = 3 * G * M / c**2  # Photon sphere radius

# Test point parameters
num = 30  # Number of photons
x_i = 7.5
y_i = 2.5
y_f = 4

# Initial velocity in Cartesian coordinates
v_x = -c  # Speed of light in negative x-direction
v_y = 0   # No velocity in the y-direction

# Prepare for matplotlib plotting
trajectories = []

for i in range(num):
    trajectory = []

    # Define initial test point
    if num > 1:
        test_point = np.array([x_i, y_i + i / (num - 1) * (y_f - y_i), 0]) # linear interpolation
    else:
        test_point = np.array([x_i, y_i, 0])

    # Convert to polar coordinates
    r = np.linalg.norm(test_point) * 1e6  # Radius in meters
    p = atan2(test_point[1], test_point[0])  # Initial planar angle

    # Metric coefficients
    g_tt = -(1 - 2 * G * M / (r * c**2))
    g_rr = 1 / (1 - 2 * G * M / (r * c**2))
    g_pp = r**2

    # Initial velocities
    v_r = (v_x * cos(p) + v_y * sin(p))  # Radial velocity
    v_p = (-v_x * sin(p) + v_y * cos(p)) / r  # Angular velocity
    v_t = sqrt(-(g_rr * v_r**2 + g_pp * v_p**2) / g_tt)

    # Check the null geodesic condition for the initial point
    #print(g_tt * v_t**2 + g_rr * v_r**2 + g_pp * v_p**2)

    # Integrate geodesics
    for j in range(curve_num):
        if r > horizon + 10000:
            # Precompute common terms
            term1 = G * M / (r**2 * c**2)  # Common term: GM / r^2c^2
            term2 = 1 - 2 * G * M / (r * c**2)  # Common term: 1 - 2GM / rc^2

            # Christoffel symbols using common terms
            Γ_r_tt = term1 * term2
            Γ_r_rr = -term1 / term2
            Γ_r_pp = -r * term2
            Γ_p_rp = 1 / r
            #Γ_t_rt = term1 / term2 # ignoring time marching

            # Update change in velocities
            dv_r = -Γ_r_tt * v_t**2 - Γ_r_rr * v_r**2 - Γ_r_pp * v_p**2
            dv_p = -2 * Γ_p_rp * v_r * v_p
            #dv_t = -2 * Γ_t_rt * v_r * v_t # ignoring time marching

            # Update velocities
            v_r += dv_r * d_lambda
            v_p += dv_p * d_lambda
            #v_t += dv_t * d_lambda # ignoring time marching

            # Update positions
            r += v_r * d_lambda
            p += v_p * d_lambda

            # Metric tensor components (update for new r)
            g_tt = -term2
            g_rr = 1 / term2
            g_pp = r**2

            # Recalculate v_t from the metric components
            v_t = sqrt(-(g_rr * v_r**2 + g_pp * v_p**2) / g_tt)

            # Check the null geodesic condition at each step
            #print(g_tt * v_t**2 + g_rr * v_r**2 + g_pp * v_p**2)

            # Store Cartesian coordinates
            x = (r / 1e6) * cos(p)
            y = (r / 1e6) * sin(p)

            # Only store points within the -10 to 10 range
            if -10 <= x <= 10 and -10 <= y <= 10:
                trajectory.append((x, y))
        else:
            break

    trajectories.append(trajectory)

# Plot using matplotlib
plt.figure(figsize=(8, 8))

# Plot each trajectory
for trajectory in trajectories:
    trajectory = np.array(trajectory)
    if len(trajectory) > 0:  # Plot only if there are points
        plt.plot(trajectory[:, 0], trajectory[:, 1])

# Plot the event horizon
circle = plt.Circle((0, 0), horizon / 1e6, color='black', fill=True, label="Event Horizon")
plt.gca().add_artist(circle)

# Plot the photon sphere
photon_sphere_circle = plt.Circle((0, 0), photon_sphere / 1e6, color='red', fill=False, linestyle='--', linewidth=1.5, label="Photon Sphere")
plt.gca().add_artist(photon_sphere_circle)

# Set plot limits explicitly
plt.xlim(-10, 10)
plt.ylim(-10, 10)

# Configure plot
plt.title("Photon Trajectories Around a Black Hole")
plt.xlabel("x (scaled by 1e6 meters)")
plt.ylabel("y (scaled by 1e6 meters)")
plt.axhline(0, color="gray", linewidth=0.5)
plt.axvline(0, color="gray", linewidth=0.5)
plt.axis('equal')
plt.grid()
plt.legend()
plt.show()

I expect my null geodesic print statement inside the for loop to be either 0 or a relatively small integer. For example, the first null geodesic print statement may be 0, 16, -8, etc. due to a small amount of imprecision with the large float addition (e17 magnitudes).

I have tried to debug by replacing my "else" statement with "elif j == 1" and looking at the very next iteration, it can be seen the null geodesic prints a much larger float.

I believe figuring out the null geodesic error will reveal if my trajectories are incorrect.

• python
• math
• simulation
• physics

2 Answers 2

You are missing a factor of 2 in the expression for the acceleration component dv_p. Change that line to

dv_p = - 2 * Γ_p_rp * v_r * v_p

This is because Γ_p_pr = Γ_p_rp and you need to include that term as well in the summation.

Neil Butcher's post also picks out another non-zero Christoffel symbol Γ_t_rt that I missed and would be necessary to get v_t right and then affect the trajectories. BUT ...

I believe the metric check is a red herring - you initialise v_t at the start with that metric and I think you can (and should) do the same again each time the metric components are recalculated from a new r. That is then bound to keep you on a geodesic.

For completeness, I believe that you are solving systems of the form

or, in terms of your quasi-velocities,

If you do have to print the null-geodesic error then I think it should be normalised. You do have to see it in the context of c², which is very large. Anyway, making the factor-of-2 change does reduce this error by 10 orders of magnitude.

You will have to tell me the outcome of the plotting. I don't have, or know how to use, Blender.

As you define a variable horizon it would be useful to use that (or the more common r_s), rather than repeating 2 * G * M / c**2

In summary:

• put the extra factor of 2 in the expression for dv_p
• update v_t from the metric coefficients every time you recalculate the metric from a new r.

I know nothing about the plotting but there are a few omissions with the terms of the accelerations

You need to double any a asymmetric terms

You calculate Γ_p_rp and don't bother to calculate Γ_p_pr (due to the symmetry). However you do need to compensate for that in your accelerations. eg 2.0*Γ_p_rp * v_r * v_p

You need to include the 'acceleration' of time

(I don't know what to call dv_t, but you need to include it)

There is a symbol totally omitted

Γ_t_rt = (G * M / (r**2 * c**2)) / (1 - 2 * G * M / (r * c**2))

this needs to be included and it causes to updates in v_t

dv_t = 2.0*Γ_t_rt * v_t * v_r
v_t += dv_t * d_lambda

Putting this together gives the updated code

        # Relevant Christoffel symbols

        Γ_t_rt = (G * M / (r**2 * c**2)) / (1 - 2 * G * M / (r * c**2))
        Γ_r_tt = (G * M / (r**2 * c**2)) * (1 - 2 * G * M / (r * c**2))
        Γ_r_rr = -(G * M / (r**2 * c**2)) / (1 - 2 * G * M / (r * c**2))
        Γ_r_pp = -r * (1 - 2 * G * M / (r * c**2))
        Γ_p_rp = 1 / r

        # Get accelerations
        dv_t = -2.0*Γ_t_rt * v_t * v_r
        dv_r = -Γ_r_tt * v_t**2 - Γ_r_rr * v_r**2 - Γ_r_pp * v_p**2
        dv_p = -2.0*Γ_p_rp * v_r * v_p

        # Update velocities
        v_t += dv_t * d_lambda
        v_r += dv_r * d_lambda
        v_p += dv_p * d_lambda