Solution (source code)

= Solution

Reflection across $\theta=\pi/2$ is an <isometry>; a <geodesic> initially tangent to its fixed equatorial plane remains in that plane by uniqueness of the <geodesic equation>. Thus take $\theta=\pi/2$, $\dot\theta=0$, and $\Sigma=r^2$. Define the algebraic abbreviations
$$
A=1-\frac{2M}{r},\qquad
B=\frac{2Ma}{r},\qquad
C=r^2+a^2+\frac{2Ma^2}{r}.
$$
The <conserved quantities> of part (b) become
$$
E=A\dot t+B\dot\phi,\qquad
L=-B\dot t+C\dot\phi,
\qquad
AC+B^2=\Delta.
$$
For $\Delta\ne0$, inversion gives
$$
\dot t=\frac{CE-BL}{\Delta},\qquad
\dot\phi=\frac{BE+AL}{\Delta}.
$$
The <null condition> can be written without expanding a quadratic square:
$$
0=-E\dot t+L\dot\phi+\frac{r^2}{\Delta}\dot r^2.
$$
Substituting the preceding expressions yields the <equatorial Kerr null geodesic radial equation>
$$
\boxed{
\dot r^2
=E^2+\frac{a^2E^2-L^2}{r^2}
+\frac{2M(L-aE)^2}{r^3}.
}
$$
An equivalent useful polynomial form is
$$
\boxed{
r^4\dot r^2=\mathcal R(r),\qquad
\mathcal R(r)=\left[E(r^2+a^2)-aL\right]^2
-\Delta(L-aE)^2.
}
$$
These expressions use an <affine parameter> $s$; a common rescaling of $E$ and $L$ accompanies its rescaling. The inversion was performed in the exterior $r>r_+$, where <Boyer-Lindquist coordinates> are regular. The polynomial radial identity has a smooth continuation to a horizon, but the formulas for $\dot t$ and $\dot\phi$ there must be replaced by regular-coordinate expressions. In the limit $a=0$, the result reduces to the <Schwarzschild spacetime> equation $\dot r^2=E^2-(1-2M/r)L^2/r^2$.