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.
Share Improve this question edited Jan 21 at 22:37 Adam Swearingen asked Jan 15 at 7:46 Adam SwearingenAdam Swearingen 2032 silver badges7 bronze badges 0
2 Answers
Reset to default 1You 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 c2, 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 rs), rather than repeating 2 * G * M / c**2
In summary:
- put the extra factor of 2 in the expression for
dv_p - update
v_tfrom 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