I'm interested in plotting the trajectories of null geodesics near an uncharged rotating black hole (described by the Kerr solution) which involves a system of first order differential equations. Kerr spacetimes are stationary and axially symmetric, and so because of the existence of these two Killing fields I'm really only interested in the "projection" of these trajectories in the $(r, \theta)$-plane.
By using conservation laws the equations describing these null geodesic are given by the following system of ODE's (which can be found, specifically, at the end of $\S$62 in Chandrasekhar's book 'The Mathematical Theory of Black Holes', but also in just about any other book describing the geometry of Kerr black holes)
$$ \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 of motion/conserved quantities, $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 I don't think it's particularly important since I feel as though I'm really only interested in equations (1) and (2)).
Equations (1) and (2) 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}$$
Explicitly,
$$\eta = \frac{\mathscr{L}}{E^2} \quad and \quad \xi = \frac{L_z}{E} $$
In this form it seems as though (according to $\S$63 of Chandrasekhar) the only restriction placed on the constants $E, \eta$ and $\xi$ is that $\eta > 0$.
I have a couple questions regarding how (theoretically) one would go about numerically approximating (and subsequently plotting) the $r$ and $\theta$ solutions of this system of differential equations.
First, to me, it seems as though since I'm only interested in the $(r, \theta)$ motion of any given null geodesic, would I then only have worry about numerically solving (3) and (4), that is, even though there are 4 equations governing the total motion of the geodesic, since (3) and (4) are independent of $t$ and $\phi$ it seems like equations (3) and (4) are kind of standalone equations governing the motion in the $(r, \theta)$-plane. Is this assertion of mine correct? Need I only deal with equations (3) and (4)?
Another concern of mine is that; I've read that the geodesics in and around a black hole (with $\eta > 0$) oscillate symmetrically about the equatorial plane ($\theta = \pi/2$) and the geodesics which do this 'take infinitely long' to fall into the black hole, that is, as they oscillate toward the black hole it takes infinite coordinate time to fall into the black hole. This fact concerns me because when plugging (3) and (4) into Mathematica, the equations need to be parameterized, so really the system looks like
$$ \frac{[r(\lambda)^2 + a^2\cos^2\theta(\lambda)]^2}{E^2} \dot{r}(\lambda)^2 = r(\lambda)^4 + (a^2 - \xi^2 - \eta)r(\lambda)^2 + 2M(\eta + (\xi - a)^2))r(\lambda) - a^2\eta$$
$$\frac{[r(\lambda)^2 + a^2\cos^2\theta(\lambda)]^2}{E^2}\dot{\theta}(\lambda)^2 = \eta + a^2\cos^2\theta(\lambda) - \xi^2\cot^2\theta(\lambda)$$
In which case, does this 'infinite amount of coordinate time' it takes to fall into the black hole affect how big of a domain one would have to allow for $\lambda$? if it does; theoretically this isn't an issue but practically it's not really feasible.
These last two questions have a little more to do with implementation and less to do with the theory.
Is dealing with these ((3) and (4)) equations actually any more useful than simply working directly from the geodesic equations, from a numerical stand point, that is?
And lastly, I've recently found out that one can decouple (3) and (4) and I was wondering (again from a numerical standpoint) if there's any benefit in working with decoupled equations as opposed to these coupled ones?
Any input on any of my above questions/concerns would be greatly appreciated.
