# Plotting Geodesics in Kerr

I'm interested in plotting the trajectories of null geodesics near a rotating black hole (given by the Kerr solution) which involves a system of first order differential equations.

Some Context (not entirely necessary):

Kerr spacetimes are stationary and axially symmetric, that is, the metric (denote it by $g$) is independent of $t$ and $\phi$ (with the remaining two coordinates being $r$ and $\theta$). Because of the existence of these two Killing fields I'm really only interested in "projection" of the trajectories in the $(r, \theta)$-plane.

Now by using conservation laws the equations of motion for a null geodesic are given by

$$\rho^4 \dot{r}^2 = E^2r^4 + (a^2E^2 - L_z^2 - \mathscr{L})r^2 + 2Mr(\mathscr{L} + (L_z - aE)^2) - a^2\mathscr{L} \tag{1}$$

$$\rho^4\dot{\theta}^2 = \mathscr{L} + (a^2 E^2 - L_z^2\csc^2\theta)\cos^2\theta \tag{2}$$

$$\rho^2\dot{\phi} = \frac{2aMrE + (\rho^2 - 2Mr)L_z\csc^2\theta}{\Delta}$$

$$\rho^2\dot{t} = \frac{\Sigma^2E - 2aMrL_z}{\Delta}$$

Where $\dot{x^\mu}$ denotes differentiation w.r.t. some affine parameter $\lambda$. We have that $E, L_z$ and $\mathscr{L}$ are constants, $a$ and $M$ are the Kerr parameters such that $0 < a^2 < M^2, M> 0$ ($M$ is usually taken to be 1) and

$$\rho^2 = r^2 + a^2\cos^2\theta, \quad \Delta = r^2 + a^2 - 2Mr$$

(at the moment I can't quite recall what $\Sigma$ is, but it's not particularly important since I'm really only interested in equations (1) and (2)).

Now as disgusting as all of these equations, together, may look, as mentioned above I'm really only interested in the $(r, \theta)$ motion, i.e. in equations (1) and (2). These two equations can be rewritten slightly by dividing both of them by $E^2$, yielding

$$\frac{\rho^4}{E^2} \dot{r}^2 = r^4 + (a^2 - \xi^2 - \eta)r^2 + 2M(\eta + (\xi - a)^2))r - a^2\eta \tag{3}$$

$$\frac{\rho^4}{E^2}\dot{\theta}^2 = \eta + a^2\cos^2\theta - \xi^2\cot^2\theta \tag{4}$$

The only restriction placed on the constants $E, \eta$ and $\xi$ is that $\eta > 0$. This is the system of coupled ODE's whose solutions I'd like to numerically approximate and plot.

My Problem:

I wish to plot the solutions of the above system of coupled ODE's (given by (3) and (4)). My thought regarding how to go about plotting these curves was to use NDSolve in the following way

M := 1
a := 0.8
e := 1
rminreal := M + Sqrt[M^2 - a^2]
\[Eta] = 3M
\[Xi] = 3.5

R1[t_] := (r[t]^2 + a^2 Cos[\[Theta][t]]^2)^2/e^2 r'[t]^2 - r[t]^4 - (a^2 - \[Xi]^2 - \[Eta]) r[t]^2  + 2 M ((a - \[Xi])^2 - \[Eta]) r[t] - a^2 \[Eta]

T1[t_] := (r[t]^2 + a^2 Cos[\[Theta][t]]^2)^2/e^2 \[Theta]'[t]^2 - (\[Xi]^2/Sin[\[Theta][t]]^2 - a^2) Cos[\[Theta][t]]^2

sol = NDSolve[{R1[t] == 0, T1[t] == 0, r[0] == 2, \[Theta][0] == Pi/4}, {r, \[Theta]}, {t, 0, 5}, MaxSteps -> 500000]

Piece = 4
EndPt = sol[[Piece]][[1]][[2]][[1]][[1]][[2]]
ParametricPlot[{{rs Cos[t \[Pi]/EndPt], rs Sin[t \[Pi]/EndPt]} , Evaluate[{r[t] Cos[\[Theta][t]], r[t] Sin[\[Theta][t]]} /. sol[[Piece]]]}, {t, 0, EndPt}, AspectRatio -> 1]


Here rminreal is the event horizon of the black hole.

There seems to be a few things that tend to go wrong when I try and run the code.

(1) The first thing is that the output tends to come in pieces, so instead of simply getting back two interpolatingfuncitons (one for r and one for \[Theta]) Mathematica returns a list of four pairs of interpolating functions, for example:

This seems like it can be remedied by increasing MaxSteps, but MaxSteps is already rather large and the program takes a while to evaluate anyway.

(2) Putting aside the chopping off of the domain and only getting back pieces of a solution, more often than not when I try to plot the solutions I get nothing back, that is, I get the plot of the event horizon but nothing from the interpolating function, specifically, the parameters and initial conditions used above yield nothing when plotted.

(3) Also from time to time (that is, for certain choices of constant $\eta$ and $\xi$) Mathematica throws up and error saying that the step size is essentially zero and its stopping the calculation.

I'm wondering if anybody knows why these things are happening and if there's is a way aroung them or if there's something i've done incorrectly (or inefficiently)?

I'm pretty new when it comes to using Mathematica so I definitely welcome any constructive criticism regarding my writing of the code. Moreover, I'm wondering have I gone about setting up the equations correctly? It seems as though there may be intricacies regarding NDSolve that I'm unaware of, and if there are better ways of going about plotting these solutions; I would appreciate feedback on how to do it.

It seems as though the the system of ODE's, itself, is pretty benign which is why I fear that I've probably botched something here.

Any feed back is much appreciated

-