I am writing a numerical solver for the null geodesics of an arbitrary spacetime metric, but currently mostly targetting black hole spacetimes. I am attempting to numerically integrate the geodesic equations:
$$ \frac{d^2 x^\mu}{ds^2} = -\Gamma^\mu_{\alpha \beta} \frac{dx^\alpha}{ds} \frac{dx^\beta}{ds} $$
Given $c = G = 1$, $M = 15000$, and initial conditions $x^\mu = \langle 0, 6M, \frac{\pi}{4}, 0\rangle$ and $v^\mu = \left \langle \frac{r_0^2 + a^2}{\Delta}, -1, 0, \frac{a}{\Delta} \right \rangle$, this is the code of my solver:
import numpy as np
import matplotlib.pyplot as plt
import torch
def calculate_christoffel(jacob, g_inv, dims):
# based on https://github.com/AndreaAntoniali/Riemann-tensor-calculator/blob/main/Riemann_Calculations.ipynb
gamma = np.zeros((dims, dims, dims))
for beta in range(dims):
for mu in range(dims):
for nu in range(dims):
for alpha in range(dims):
gamma[beta,mu,nu] = 1/2 * g_inv[alpha][beta] * (jacob[alpha][mu][nu] + jacob[alpha][nu][mu] - jacob[mu][nu][alpha])
return gamma
def christoffel_at_point_4d(metric, inverse_metric, t, r, theta, phi, dims):
coord = torch.tensor([t, r, theta, phi], requires_grad=True)
g_inv = inverse_metric(coord)
jacobian = torch.autograd.functional.jacobian(metric, coord, create_graph=True)
return calculate_christoffel(jacobian, g_inv, dims)
def kerr_metric(coords, M=15000, a=0.97):
t = coords[0]
r = coords[1]
theta = coords[2]
phi = coords[3]
sigma = r ** 2 + a ** 2 * torch.cos(theta) ** 2
delta = r ** 2 - 2 * M * r + a ** 2
return torch.tensor([
[-(1 - (2 * M * r) /sigma), 0., 0., -((2 * M * r * a * torch.sin(theta) ** 2) / sigma)],
[0., sigma / delta, 0., 0.],
[0., 0., sigma, 0.],
[-((2 * M * r * a * torch.sin(theta) ** 2) / sigma), 0., 0., (r ** 2 + a ** 2 + (2 * M * r * a ** 2)/sigma * torch.sin(theta) ** 2) * torch.sin(theta) ** 2]
])
def kerr_inverse_metric(coords, M=15000, a=0.97):
# Based off https://www.roma1.infn.it/teongrav/onde19_20/kerr.pdf
t = coords[0]
r = coords[1]
theta = coords[2]
phi = coords[3]
sigma = r ** 2 + a ** 2 * torch.cos(theta) ** 2
delta = r ** 2 - 2 * M * r + a ** 2
return torch.tensor([
[-1 / delta * (r ** 2 + a ** 2 + (2 * M * r * a ** 2) / sigma * torch.sin(theta) ** 2), 0., 0., -(2 * M * r * a) / (sigma * delta)],
[0., delta / sigma, 0., 0.],
[0., 0., 1 / sigma, 0.],
[-(2 * M * r * a) / (sigma * delta), 0., 0., (delta - a ** 2 * torch.sin(theta) ** 2) / (sigma * delta * torch.sin(theta) ** 2)]
])
def kerr_d_ds(X, s, metric=kerr_metric, inverse_metric=kerr_inverse_metric):
"""
The value of the first and second
derivatives with respect to an affine
parameter s
"""
# Create a new vector to hold the positions and velocities
u = np.zeros(X.shape)
# X is a vector with 4 components of position
# and 4 components of velocity
x = X[:4]
velocities = X[4:]
# Find christoffel symbols given the position and the metric
# here t is coordinate time, not the affine parameter s, and
# also not the proper time
x0, x1, x2, x3 = x
Gamma = christoffel_at_point_4d(metric, inverse_metric, x0, x1, x2, x3, 4)
# Given the christoffel symbols, calculate the next position
for mu in range(4):
for i in range(4):
for j in range(4):
# Solve for x components
# we sum due to the Einstein summation convention
u[mu] += -Gamma[mu][i][j] * velocities[i] * velocities[j]
# Solve for v components
u[4:] = velocities
return u
def rk4(f, u0, t0, tf, n):
"""
Expects a function in the form y' = f(y, t) to solve
f - function, MUST return numpy array of derivatives
u0 - initial values
t0 - initial time
tf - final time
n - number of samples
"""
t = np.linspace(t0, tf, n+1)
u = np.zeros((n+1, len(u0)))
u[0] = u0
h = t[1]-t[0]
for i in range(n):
k1 = h * f(u[i], t[i])
k2 = h * f(u[i] + 0.5 * k1, t[i] + 0.5h)
k3 = h f(u[i] + 0.5 * k2, t[i] + 0.5h)
k4 = h f(u[i] + k3, t[i] + h)
u[i+1] = u[i] + (k1 + 2*(k2 + k3 ) + k4) / 6
return u, t
c = 1 # as we use geometrized units
r0 = 6 * 15000 # orbital radius of 6M
theta0 = np.pi / 4 # if we use pi / 2 we get NaNs
phi0 = 0
v_r0 = c
v_theta0 = 0
v_phi0 = 0
u0 = np.array([0, r0, theta0, phi0, 0, v_r0, v_theta0, v_phi0])
tf = 1500
x_sol, t_sol = rk4(kerr_d_ds, u0, 0, tf, 50)
r = x_sol[:,1]
theta = x_sol[:,2]
phi = x_sol[:,3]
fig = plt.figure()
ax = plt.axes(projection="3d")
x = r * np.sin(theta) * np.cos(phi)
y = r * np.sin(theta) * np.sin(phi)
z = r * np.cos(theta)
ax.plot3D(x, y, z)
ax.set_title("Null geodesics 3D plot")
plt.show()
However, when I run the solver, this is the value of $r$:
array([90000., 90000., 90000., 90000., 90000., 90000., 90000., 90000.,
90000., 90000., 90000., 90000., 90000., 90000., 90000., 90000.,
90000., 90000., 90000., 90000., 90000., 90000., 90000., 90000.,
90000., 90000., 90000., 90000., 90000., 90000., 90000., 90000.,
90000., 90000., 90000., 90000., 90000., 90000., 90000., 90000.,
90000., 90000., 90000., 90000., 90000., 90000., 90000., 90000.,
90000., 90000., 90000.])
I'm wondering why the radial, as well as azimuthal and polar velocities, are constant. I would not expect that a photon placed at an orbit of $6M$ from a Kerr black hole would suddenly stop moving. I've already factored out issues relating to the photon being too close to the event horizon/ergospheres or numerical precision errors, so why does the solver show this result?
