Geodesic Congruences in FRW, Schwarzschild and Kerr Spacetimes

Introduction

The theory of geodesic congruences is extensively covered in many textbooks (see References); what follows in the introduction is a brief summary. Consider a 1-parameter family of timelike geodesics $\gamma_s(\lambda)$, where $s$ labels each geodesic in the family whilst $\lambda$ is an affine parameter along each $\gamma_s$. Then the vector field $\xi \equiv \partial / \partial s$ is tangent to curves of constant $\lambda$ and is interpreted as a deviation vector between neighbouring geodesics.

In some neighbourhood of the family, $(s,\lambda, x^2, x^3)$ is a coordinate chart satisfying $\xi = \partial/\partial s$ and $u = \partial/\partial \lambda$. By the equality of mixed partial derivatives, the commutator of $u$ and $\xi$ is zero (i.e. $\xi$ is Lie transported along $u$),\begin{align*}
0 = [u, \xi]^a = (L_{u} \xi)^a = \xi^b \nabla_b u^a – u^b \nabla_b \xi^a
\end{align*}which implies that $\dfrac{D\xi^a}{d\lambda} = u^b \nabla_b \xi^a = (\nabla_b u^a) \xi^b \equiv {B^a}_b \xi^b$. In other words, the tensor ${B^a}_b \equiv \nabla_b u^a$ measures the failure of $\xi$ to be parallel transported along the geodesic.

A family of timelike geodesics is a congruence in a region $U \subseteq M$ if through each point $p \in U$ passes exactly one geodesic. Define a tensor $h_{ab} \equiv g_{ab} + u_a u_b$; this can be shown to satisfy $h^2 = h$,\begin{align*}
h^a_b h^b_c &= (\delta^a_b + u^a u_b)(\delta^b_c + u^b u_c) \\
&= \delta^a_c + u^a u_c + u^a u_c + u^a (u^b u_b)u_c \\
&= \delta^a_c + u^a u_c \\
&= h^a_c
\end{align*}assuming $u$ to be normalised as $u\cdot u = -1$. The tensor $h$ is therefore a projection onto the spacelike hypersurfaces $\Sigma$ which are orthogonal to the vector field $u$ (that is, $h^a_b u^b = 0$). It is also observed that $h^a_a = \delta^a_a + u^a u_a = 4-1 = 3$.

The tensor $B$ can be written as\begin{align*}
{B^a}_b = \dfrac{1}{3} \theta h^a_b + {\sigma^a}_b + {\omega^a}_b
\end{align*}where the expansion $\theta \equiv {B^a}_a$ is the trace of $B$, the shear tensor $\sigma_{ab} \equiv B_{(ab)} – \dfrac{1}{3}\theta h_{ab}$ is the symmetric-tracefree part of $B$, and the rotation tensor  $\omega_{ab} = B_{[ab]}$ is the antisymmetric part of $B$.

The evolution of the expansion $\theta$ can be shown (see Appendix) to satisfy\begin{align*}
\dfrac{d\theta}{d\lambda} = \dfrac{1}{3}\theta^2 – \sigma^2 + \omega^2 – R_{ab}u^a u^b
\end{align*}which is called Raychaudhuri’s equation.

To conclude the introduction, it is worth touching briefly upon the physical interpretation of $\theta$, $\sigma$ and $\omega$. Consider two small, spacelike hypersurfaces $\delta \Sigma(\lambda)$ and $\delta \Sigma(\lambda + \delta\lambda)$ in a small neighbourhood of some geodesic $\gamma_0$ in the congruence; both of these hypersurfaces are constructed to be orthogonal to $\gamma_0$ at the points of intersection. Consider a small vector $\mathbf{r}(\lambda)$ lying in $\delta \Sigma(\lambda)$. In affine parameter distance $\delta \lambda$, this vector changes by $\delta \mathbf{r} = \mathbf{B}\mathbf{r} \delta \lambda$.

The expansion matrix takes $\mathbf{r} \mapsto \mathbf{r} + \dfrac{\theta}{3} \mathbf{r} \delta \lambda$. Under such a transformation, the volume of a small sphere changes as $\delta V = 4\pi r^2 \delta r = 4\pi r^3 \dfrac{\delta r}{r} = \dfrac{4}{3}\pi r^3 \theta \delta \lambda = V \theta \delta \lambda$. Put differently, the expansion $\theta = \dfrac{1}{V} \dfrac{\delta V}{\delta \lambda}$ is the fractional change of volume per unit worldline parameter.

A shear is a volume-preserving transformation, that is, one of unit determinant; under the action of a shear, a sphere is mapped to an ellipsoid. Finally, since any antisymmetric tensor is dual to a (pseudo-)vector, the action of the rotation matrix on $\mathbf{r}$ is equivalent to $\mathbf{r} \mapsto \mathbf{r} + (\boldsymbol{\omega} \times \mathbf{r})\delta \lambda$, where $\boldsymbol{\omega} = (\omega_1, \omega_2, \omega_3)$. This is a rotation by angle $\omega \delta \lambda$ about an axis parallel to $\boldsymbol{\omega}$.

Friedmann-Robertson-Walker spacetime

The Friedmann-Robertson-Walker (FRW) spacetime has metric\begin{align*}
g= -dt^2 + a^2(t)\left( \dfrac{dr^2}{1-kr^2} + r^2 d\Omega^2 \right)
\end{align*}where $d\Omega^2 = d\theta^2 + \sin^2{\theta}d\phi^2$.

Consider a congruence of comoving worldlines with tangent vector $u = \partial / \partial t$. In the coordinates $(t,r,\theta, \phi)$, the tangent vector to the congruence is simply $u^{\mu} = \delta^{\mu}_t = (1,0,0,0)$. This congruence is geodesic, since $u^{\nu} \nabla_{\nu} u^{\mu}$ vanishes:\begin{align*}
u^{\nu} \nabla_{\nu} u^{\mu} &= u^{\nu} \partial_{\nu} u^{\mu} + u^{\nu} \Gamma^{\mu}_{\rho \nu} u^{\rho} \\
&= \partial_t \delta^{\mu}_t + \delta^{\nu}_t \Gamma^{\mu}_{\rho \nu} \delta^{\rho}_t \\
&= \Gamma^{\mu}_{tt} \\
&= \dfrac{1}{2} g^{\mu \sigma} (2\partial_t g_{\sigma t} – \partial_{\sigma} g_{tt}) = 0
\end{align*}The expansion $\theta$ of the congruence is given by\begin{align*}
\theta = \nabla_{\mu} u^{\mu} &= \partial_{\mu} u^{\mu} + \Gamma^{\mu}_{\nu \mu} u^{\nu} = \Gamma^{\mu}_{t{\mu}} = \dfrac{1}{\sqrt{-g}} \dfrac{\partial}{\partial t}\left( \sqrt{-g} \right)
\end{align*}and since the metric determinant $g \equiv \mathrm{det}(g)$ satisfies $\sqrt{-g} = \dfrac{a^3 r^2 \sin{\theta}}{\sqrt{1-kr^2}}$ it follows that\begin{align*}
\theta = \dfrac{\sqrt{1-kr^2}}{a^3 r^2 \sin{\theta}} \left(\dfrac{3a^2 \dot{a} r^2 \sin{\theta}}{\sqrt{1-kr^2}}\right) = 3\dfrac{\dot{a}}{a}
\end{align*}Note that $3\dfrac{\dot{a}}{a} = \dfrac{1}{a^3} \dfrac{d}{dt} \left(a^3\right)$, which is consistent with the notion that the expansion $\theta$ measures the fractional rate of increase of the volume of the geodesic congruence. The shear tensor is given by\begin{align*}
\sigma_{ab} = \nabla_{(b} u_{a)} – \dfrac{1}{3} \theta h_{ab} &= -\Gamma^c_{ab} u_c – \dfrac{\dot{a}}{a} h_{ab}
\end{align*}Since $h_{\mu \nu}(t) = a^2(t) \tilde{h}_{\mu \nu}$ where $\tilde{h}_{\mu \nu}$ is time-independent, its time derivative satisfies $\dfrac{\partial h_{\mu \nu}}{\partial t} = 2a \dot{a} \tilde{h}_{\mu \nu} = 2\dfrac{\dot{a}}{a} h_{\mu \nu}(t)$. Therefore, evaluating $\sigma$ in coordinate basis,\begin{align*}
\sigma_{\mu \nu} &= \Gamma^t_{\mu \nu} – \dfrac{\dot{a}}{a}h_{\mu \nu} = – \dfrac{1}{2} g^{tt} \dfrac{\partial h_{\mu \nu}}{\partial t} – \dfrac{\dot{a}}{a}h_{\mu \nu} = \dfrac{\dot{a}}{a} h_{\mu \nu} – \dfrac{\dot{a}}{a}h_{\mu \nu} = 0
\end{align*}Finally the rotation tensor $\omega_{ab}$ satisfies\begin{align*}
\omega_{ab} = \nabla_{[b} u_{a]} = \dfrac{1}{2}(- \Gamma^c_{ab} u_c + \Gamma^c_{ba} u_c) = 0
\end{align*}and thus also vanishes due to the symmetry of the lower indices of the Christoffel symbol. It can in fact be shown from Frobenius’ theorem that the rotation tensor vanishes if and only if the geodesic congruence is hypersurface orthgonal; this is indeed true in the FRW model, since it follows from the assumptions of isotropy and homogeneity that the comoving observers are orthogonal to the hypersurfaces of spatial homogeneity. Collecting the previous results, Raychaudhuri’s equation reduces to\begin{align*}
\dfrac{d\theta}{d\tau} &= – \dfrac{1}{3}\theta^2 – R_{\mu \nu}u^{\mu} u^{\nu} \\
3\left( \dfrac{a\ddot{a} – \dot{a}^2}{a^2}\right)&= – 3\dfrac{\dot{a}^2}{a^2} – R_{tt} \\
\implies \dfrac{\ddot{a}}{a} &= – \dfrac{1}{3}R_{tt}
\end{align*}From Einstein’s equation it follows that $R_{tt} = 8\pi \left(T_{tt} – \dfrac{1}{2}Tg_{tt} \right)$, where for a perfect fluid $T_{tt} = \rho$ and $T = {T^{\mu}}_{\mu} = – \rho + 3p$. We have therefore arrived at the equation of motion\begin{align*}
\dfrac{\ddot{a}}{a} = – \dfrac{8\pi}{3}\left( \rho + \dfrac{1}{2} (-\rho + 3p)\right) = – \dfrac{4\pi}{3}(\rho + 3p)
\end{align*}This equation shows that if the strong energy condition is satisfied (that is if $\rho + 3p \geq 0$, such that $\ddot{a} \leq 0$) then there was a singularity at a previous time $t_*$ when $a(t_*) = 0$. If the time now is denoted by $t_0$, then an upper bound on $t_0 – t_*$ is obtained by setting $\ddot{a} = 0$ (along with the initial conditions $\dot{a}(t_0) \equiv H_0$ and $a(t_0) = 1$):\begin{align*}
\dot{a} = H_0 \implies a(t) = H_0 (t-t_0) + 1
\end{align*}Putting $a(t_*) = 0$ gives an upper bound of $t_0 – t_* = H_0^{-1}$ for the age of the universe.

Schwarzschild spacetime

The Schwarzschild metric is\begin{align*}
g = -fdt^2 + f^{-1} dr^2 + r^2 d\Omega^2
\end{align*}where $f = 1- \dfrac{2M}{r}$. Consider the vector field\begin{align*}
u = \dfrac{1}{\sqrt{1-\dfrac{3M}{r}}} \left( \dfrac{\partial}{\partial t} + \sqrt{\dfrac{M}{r^3}} \dfrac{\partial}{\partial \theta}\right)
\end{align*}which represents observers in circular orbits (notice that the radial coordinate $r$ remains fixed along the geodesics). It is timelike and normalised as $u \cdot u = -1$, since\begin{align*}
g_{\mu \nu}u^{\mu} u^{\nu} &= g_{tt} u^t u^t + g_{\theta \theta} u^{\theta} u^{\theta} \\
&= \dfrac{1}{1- \dfrac{3M}{r}} \left( -f + \dfrac{M}{r} \right) = -1
\end{align*}It can also be verified that $u$ is geodesic;
\begin{align*}
u^{\nu} \nabla_{\nu} u^{\mu} &= u^{\nu}(\partial_{\nu} u^{\mu} + \Gamma^{\mu}_{\rho \nu} u^{\rho}) \\
&= u^t \Gamma^{\mu}_{\rho t} u^{\rho} + u^{\theta} \Gamma^{\mu}_{\rho \theta} u^{\rho}\end{align*}To evaluate this, the following Christoffel symbols are required:\begin{align*}
\Gamma^{\mu}_{\rho t} &= \dfrac{1}{2} g^{\mu \sigma} \left(\partial_{\rho} g_{\sigma t} – \partial_{\sigma} g_{\rho t}\right) \\
\Gamma^{\mu}_{\rho \theta} &= \dfrac{1}{2} g^{\mu \sigma} \left(\partial_{\rho} g_{\sigma \theta} – \partial_{\sigma} g_{\rho \theta}\right)
\end{align*}It is observed that $\Gamma^{\mu}_{\theta t} = \Gamma^{\mu}_{t\theta} = 0$ vanish, whilst\begin{align*}
\Gamma^{\mu}_{tt} &= – \dfrac{1}{2} g^{\mu r} \partial_r g_{tt} = \dfrac{M}{r^2} g^{\mu r} \\
\Gamma^{\mu}_{\theta \theta} &= – \dfrac{1}{2} g^{\mu r} \partial_r g_{\theta \theta} = -rg^{\mu r}
\end{align*}therefore\begin{align*}
u^{\nu} \nabla_{\nu} u^{\mu} &= u^t \Gamma^{\mu}_{tt} u^t + u^{\theta} \Gamma^{\mu}_{\theta \theta} u^{\theta} \\
&= \dfrac{1}{1- \dfrac{3M}{r}} \left( \dfrac{M}{r^2} g^{\mu r} + \dfrac{M}{r^3}(-rg^{\mu r})\right) = 0
\end{align*}The expansion $\Theta$ of this congruence is\begin{align*}
\nabla_{\mu} u^{\mu} &= \Gamma^{\mu}_{t\mu} u^t + \Gamma^{\mu}_{\theta \mu} u^{\theta}
\end{align*}where, since $\sqrt{-g} = \sqrt{r^4 \sin^2 \theta} = r^2 \sin{\theta}$, the relevant Christoffel symbols are given by\begin{align*}
\Gamma^{\mu}_{t\mu} &= \dfrac{1}{r^2 \sin{\theta}} \dfrac{\partial}{\partial t} (r^2 \sin{\theta}) = 0 \\
\Gamma^{\mu}_{\theta \mu} &= \dfrac{1}{r^2 \sin{\theta}} \dfrac{\partial}{\partial \theta} (r^2 \sin{\theta}) = \cot{\theta}
\end{align*}which implies that\begin{align*}
\Theta \equiv \nabla_{\mu} u^{\mu} = \cot{\theta} \sqrt{\dfrac{M}{r^3(1 – 3M/r)}}
\end{align*}In the Northern hemisphere $\theta \in \left(0, \dfrac{\pi}{2} \right) \implies \cot{\theta} > 0$ which implies that the geodesics are diverging, whilst in the Southern hemisphere $\theta \in \left(\dfrac{\pi}{2}, \pi \right) \implies \cot{\theta} < 0$ which implies that the geodesics are converging. At the North and South poles ($\theta = 0, \theta = \pi$ respectively), $\cot{\theta}$ is undefined.

The rotation tensor is given by\begin{align*}
\omega_{\mu \nu} = \nabla_{[\nu} u_{\mu]} &= \dfrac{1}{2}(\partial_{\nu}u_{\mu} – \Gamma^{\rho}_{\mu \nu} u_{\rho} – \partial_{\mu} u_{\nu} + \Gamma^{\rho}_{\nu \mu} u_{\rho}) \\
&= \dfrac{1}{2}(\partial_{\nu}u_{\mu} – \partial_{\mu}u_{\nu})
\end{align*}The only non-zero components are $\omega_{tr} = – \omega_{rt}$ and $\omega_{\theta r} = -\omega_{r\theta}$, which are explicitly\begin{align*}
\omega_{tr} &= \dfrac{1}{2}\dfrac{\partial}{\partial r}(g_{tt} u^t) = \dfrac{1}{2} \dfrac{\partial}{\partial r} \dfrac{-f}{\sqrt{1-\dfrac{3M}{r}}} = -\dfrac{M}{4r^2} \dfrac{(1- 6M/r)}{\left(1-3M/r\right)^{3/2}} \\ \\
\omega_{\theta r} &= \dfrac{1}{2} \dfrac{\partial}{\partial r}(g_{\theta \theta} u^{\theta}) = \dfrac{1}{2} \dfrac{\partial}{\partial r} \sqrt{\dfrac{Mr}{1-\dfrac{3M}{r}}} = \dfrac{1}{4} \sqrt{\dfrac{M}{r}} \dfrac{(1- 6M/r)}{\left(1-3M/r\right)^{3/2}}
\end{align*}Its square $\omega^2 = g^{ac} g^{bd} \omega_{cd} \omega_{ab}$ satisfies\begin{align*}
\omega^2 &= 2g^{tt} g^{rr} \omega_{tr} \omega_{tr} + 2g^{\theta \theta} g^{rr} \omega_{\theta r} \omega_{\theta r} \\
&= \dfrac{M}{8r^3} \dfrac{\left(1-6M/r \right)^2}{\left(1-3M/r\right)^3}\left[f – M/r \right] \\
&=\dfrac{M}{8r^3} \dfrac{\left(1-6M/r\right)^2}{\left(1-3M/r \right)^2}
\end{align*}The shear tensor is given by\begin{align*}
\sigma_{\mu \nu} = \nabla_{(\nu} u_{\mu)} &= \dfrac{1}{2}(\nabla_{\nu}u_{\mu}-\nabla_{\mu}u_{\nu})\\
&= \dfrac{1}{2}(\partial_{\nu}u_{\mu}+\partial_{\mu}u_{\nu}) – \Gamma^{\rho}_{\mu \nu} u_{\rho}
\end{align*}The non-zero components of the shear tensor are $\sigma_{tt}, \sigma_{rr}, \sigma_{\theta \theta}, \sigma_{\phi \phi}, \sigma_{t\theta}, \sigma_{tr}$ and $\sigma_{r\theta}$. For example, $\sigma_{tt}$ is given by\begin{align*}
\sigma_{tt} &= -\Gamma^{\rho}_{tt} u_{\rho} – \dfrac{1}{3}\theta h_{tt} \\
&= – \dfrac{M}{r^2} g^{\rho r}u_{\rho} – \dfrac{1}{3}\theta h_{tt} \\
&=-\dfrac{1}{3}\theta \left(g_{tt} + (g_{tt}u^t)^2 \right) \\
&= \dfrac{\theta f}{3}\left( 1- \dfrac{f}{\left( 1- 3M/r \right)} \right) \\
&= -\dfrac{\cot{\theta}}{3} \sqrt{\dfrac{M^3}{r^5}} \dfrac{(1- 2M/r)}{(1-3M/r)^{3/2}}
\end{align*}To calculate $\sigma^2$ directly would take too long, so Raychaudhuri’s equation will be used instead. The derivative of the expansion $\Theta$ is\begin{align*}
\dfrac{d\Theta}{d\tau} &= u^{\mu} \partial_{\mu} \Theta \\
&= u^{\theta} \dfrac{\partial \Theta}{\partial \theta} \\
&= – \sqrt{\dfrac{M}{r^3 – 3Mr^2}} \cdot \csc^2{\theta} \sqrt{\dfrac{M}{r^3 – 3Mr^2}} \\
&= -\dfrac{M \csc^2{\theta}}{r^3(1-3M/r)}
\end{align*}Since the Ricci tensor $R_{ab}$ is zero for the (vacuum) Schwarzschild solution, Raychaudhuri’s equation reads\begin{align*}
\sigma^2 &= -\dfrac{d\Theta}{d\lambda} – \dfrac{1}{3}\Theta^2 + \omega^2 \\
&= \dfrac{M \csc^2{\theta}}{r^3(1-3M/r)} – \dfrac{M\cot^2{\theta}}{3r^3(1-3M/r)} + \dfrac{M}{8r^3} \dfrac{\left(1-6M/r\right)^2}{\left(1-3M/r \right)^2} \\
&= \dfrac{M}{r^3(1-3M/r)} \left( 1 + \dfrac{2}{3} \cot^2{\theta} + \dfrac{1}{8}\dfrac{\left(1-6M/r\right)^2}{\left(1-3M/r \right)} \right)
\end{align*}This timelike vector field is an example of a geodesic congruence with all three quantities $\theta$, $\sigma$ and $\omega$ are non-zero. It is also observed that the quantity $\sigma^2 – \omega^2$ is positive for $r > 3M$.

Null geodesics in Kerr spacetime

The Kerr metric, expressed in oblate spheroidal coordinates, is\begin{align*}
g &= -dt^2 + \dfrac{r^2 + a^2 \cos^2{\theta}}{r^2 + a^2}dr^2 + (r^2 + a^2 \cos^2{\theta})d\theta^2 + (r^2+a^2)\sin^2{\theta}d\phi^2 \\
&= -dt^2 + \dfrac{\Sigma}{r^2 + a^2} dr^2 + \Sigma d\theta^2 + (r^2 + a^2) \sin^2{\theta} d\phi^2
\end{align*}where $\Sigma \equiv r^2 + a^2 \cos^2{\theta}$. Consider the vector field\begin{align*}
k = \dfrac{\partial}{\partial t} + \dfrac{\partial}{\partial r} + \dfrac{a}{r^2 + a^2} \dfrac{\partial}{\partial \phi}
\end{align*}which is null, that is, $k \cdot k = 0$:\begin{align*}
g_{\mu \nu}k^{\mu}k^{\nu} &= g_{tt} k^t k^t + g_{rr}k^r k^r + g_{\phi \phi}k^{\phi}k^{\phi} \\
&= -1 + \dfrac{\Sigma}{r^2 + a^2} + \dfrac{a^2 \sin^2{\theta}}{r^2 + a^2} \\
&= -1 + \dfrac{r^2 + a^2}{r^2 + a^2} \\
&= 0
\end{align*}$k$ is also geodesic; write\begin{align*}
k^{\nu} \nabla_{\nu} k^{\mu} &= k^{\nu} \partial_{\nu} k^{\mu} + k^{\nu} \Gamma^{\mu}_{\rho \nu} k^{\rho} \\
&= k^r \partial_r k^{\mu} + k^{\nu} \Gamma^{\mu}_{\rho \nu} k^{\rho}
\end{align*}Since only the components $k^t$, $k^r$ and $k^{\phi}$ are non-zero, the dummy indices $\rho$ and $\nu$ run through $t$, $r$ and $\phi$. Furthermore, because the metric is independent of $t$ and because $g_{tt}$ is a constant, it can be seen that $\Gamma^{\mu}_{tt} = \Gamma^{\mu}_{t\phi} = \Gamma^{\mu}_{tr} = 0$. The only required Christoffel symbols are therefore\begin{align*}
\Gamma^{\mu}_{\phi \phi} &= – \dfrac{1}{2} g^{\mu \sigma} \dfrac{\partial g_{\phi \phi}}{\partial x^{\sigma}} \\
&= – \dfrac{1}{2} g^{\mu \theta} \dfrac{\partial g_{\phi \phi}}{\partial \theta} – \dfrac{1}{2} g^{\mu r} \dfrac{\partial g_{\theta \theta}}{\partial r} \\
&= -g^{\mu \theta}(r^2 + a^2)\sin{\theta} \cos{\theta} – g^{\mu r}r\sin^2{\theta} \\ \\
\Gamma^{\mu}_{\phi r} &= \dfrac{1}{2} g^{\mu \sigma} \dfrac{\partial g_{\sigma \phi}}{\partial r} \\
&= \dfrac{1}{2} g^{\mu \phi} \dfrac{\partial g_{\phi \phi}}{\partial r} \\
&= g^{\mu \phi}r \sin^2{\theta} \\ \\
\Gamma^{\mu}_{rr} &= \dfrac{1}{2}g^{\mu \sigma} \left( \dfrac{2\partial g_{\sigma r}}{\partial r} – \dfrac{\partial g_{rr}}{\partial x^{\sigma}} \right) \\
&= g^{\mu r} \dfrac{\partial g_{rr}}{\partial r} – \dfrac{1}{2} g^{\mu \sigma} \dfrac{\partial g_{rr}}{\partial x^{\sigma}} \\
&= \dfrac{1}{2} g^{\mu r} \dfrac{\partial g_{rr}}{\partial r} – \dfrac{1}{2} g^{\mu \theta} \dfrac{\partial g_{rr}}{\partial \theta} \\
&= g^{\mu r} \dfrac{ra^2 \sin^2{\theta}}{(r^2 + a^2)} + g^{\mu \theta} \dfrac{a^2 \sin{\theta}\cos{\theta}}{r^2 + a^2}
\end{align*}After writing out the sum $k^{\nu} \Gamma^{\mu}_{\rho \nu} k^{\rho} = k^{\phi} \Gamma^{\mu}_{\phi \phi} k^{\phi} + 2k^{\phi} \Gamma^{\mu}_{r \phi} k^{r} + k^{r} \Gamma^{\mu}_{rr} k^{r}$, one arrives at\begin{align*}
k^{\nu} \nabla_{\nu} k^{\mu} &= k^r \partial_r k^{\mu} + \dfrac{a^2}{(r^2 + a^2)^2}(-g^{\mu \theta}(r^2 + a^2)\sin{\theta} \cos{\theta} – g^{\mu r}r\sin^2{\theta}) \\
&\quad\quad\quad + \dfrac{2a}{r^2 + a^2} g^{\mu \phi} r\sin^2{\theta} \\
&\quad\quad\quad + g^{\mu r} \dfrac{ra^2 \sin^2{\theta}}{(r^2 + a^2)^2} + g^{\mu \theta} \dfrac{a^2 \sin{\theta}\cos{\theta}}{r^2 + a^2} \\ \\
&= k^r \partial_r k^{\mu} + g^{\mu \phi} \dfrac{2ar \sin^2{\theta}}{r^2 + a^2}
\end{align*}Because $k^{\phi}$ is the only component of $k$ which depends on $r$ and because the metric is diagonal, if $\mu \neq \phi$ then this expression is zero. If $\mu = \phi$, then noting that $g^{\phi \phi} = \dfrac{1}{(r^2 + a^2)\sin^2{\theta}}$ it follows that\begin{align*}
k^{\nu}\nabla_{\nu} k^{\phi} = \dfrac{-2ar}{(r^2 + a^2)^2} + \dfrac{1}{(r^2 + a^2)\sin^2{\theta}} \cdot \dfrac{2ar\sin^2{\theta}}{r^2 + a^2} = 0
\end{align*}which is likewise zero. Therefore, $k$ is geodesic.

To evaluate the expansion, shear and rotation of this null congruence, it is first required to obtain an auxiliary null vector $N$ satisfying $N \cdot N = 0$ and $N \cdot k = -1$ on some initial hypersurface $\Sigma$. The choice of $N$ is not unique; consider a vector of the form\begin{align*}
N \equiv N^t \dfrac{\partial}{\partial t} + N^r\dfrac{\partial}{\partial r}\end{align*}with $N^{\theta} = N^{\phi} = 0$. The first condition implies that\begin{align*}
g_{tt}N^t N^t + g_{rr} N^r N^r &= -(N^t)^2 + \dfrac{\Sigma}{r^2 + a^2} (N^r)^2 = 0
\end{align*}or otherwise that $N^t = N^r \sqrt{\dfrac{\Sigma}{r^2 + a^2}}$. The second condition implies that\begin{align*}
N \cdot k &= g_{tt} N^t k^t + g_{rr} N^r k^r \\
&= -N^t + \dfrac{\Sigma}{r^2 + a^2} N^r \overset{!}{=} -1
\end{align*}The solution to this system is\begin{align*}
N = \dfrac{\sqrt{r^2 + a^2}}{\sqrt{r^2 + a^2} – \sqrt{\Sigma}} \left( \dfrac{\partial}{\partial t} + \sqrt{\dfrac{r^2 + a^2}{\Sigma}} \dfrac{\partial}{\partial r}\right)
\end{align*}It is then demanded that $N$ be parallel transported along the null geodesics, that is, $k^b \nabla_b N^a = 0$. It follows that $N \cdot N = 0$ and $N \cdot k =-1$ everywhere, since\begin{align*}
\dfrac{d(N\cdot N)}{d\lambda} &= k^b \nabla_b (N^a N_a) = 2N_a k^b \nabla_b N^a = 0\\ \\
\dfrac{d(N \cdot k)}{d\lambda} &= k^b \nabla_b (N^a k_a) = N_a k^b \nabla_b k^a + k_a k^b \nabla_b N^a = 0
\end{align*}where $k^b \nabla_b k^a = 0$ from the geodesic equation. It is also demanded that the deviation vector $\eta$ be orthogonal to the tangent vector $k$, in analogy to the timelike case. The vector $\eta$ can be written uniquely as\begin{align*}
\eta = aN + bk + \hat{\eta}
\end{align*}where $\hat{\eta}$ occupies the two dimensional vector subspace orthogonal to both $k$ and $N$ (and with $a$ and $b$ given by $a = – \eta \cdot k$ and $b = – \eta \cdot N$). If a tensor $p_{ab}$ is defined by\begin{align*}
p_{ab} = g_{ab} + k_a N_b + N_a k_b
\end{align*}then one can verify that $\hat{\eta}^a = p^a_b \eta^b$, since\begin{align*}
p^a_b \eta^b &= (\delta^a_b + k^a N_b + N^a k_b) \eta^b \\
&= \eta^a + k^a(N_b \eta^b) + N^a (k_b \eta^b) \\
&= \eta^a -bk^a – aN^a \\
&= \hat{\eta}^a
\end{align*}Furthermore, the tensor $p$ satisfies\begin{align*}
p^a_b p^b_c &= (\delta^a_b + k^a N_b + N^a k_b)(\delta^b_c + k^b N_c + N^b k_c) \\
&= \delta^a_c + k^a N_c + N^a k_c + k^a N_c + N^a k_c\\
&\quad\quad\quad+ k^a (N_b k^b)N_c + k^a (N_b N^b) k_c \\
&\quad\quad\quad+ N^a (k_b k^b) N_c + N^a (k_b N^b) k_c \\ \\
&= \delta^a_c + k^a N_c + N^a k_c \\
&= p^a_c
\end{align*}and is therefore a projection onto the vector subspace orthogonal to $k$ and $N$. Defining ${\hat{B}^a}_b = p^a_c p^d_b {B^c_d}$, it can be shown that an expression for $\dfrac{D\hat{\eta}^a}{d\lambda}$ analogous to the timelike case is obtained:\begin{align*}
\dfrac{D\hat{\eta}^a}{d\lambda} = {\hat{B}^a}_b \hat{\eta}^b
\end{align*}The expansion, shear and rotation are given by $\theta = {\hat{B}}^a_a$, $\hat{\sigma}_{ab} = \hat{B}_{(ab)} – \dfrac{1}{2}\theta p_{ab}$ and $\hat{\omega}_{ab} = \hat{B}_{[ab]}$. The tensor $\hat{B}$ can be written in a more explicit form using the definition of $p_{ab}$,\begin{align*}
\hat{B}_{ab} = B_{ab} + k_a N^c B_{cb} + k_b B_{ac} N^c + k_a k_b B_{cd} N^c N^d
\end{align*}It’s worth noting that $B$ and $\hat{B}$ have equal trace,\begin{align*}
{\hat{B}^a}_a &= {B^a}_a + k^a N^c B_{ca} + k_a {B^a}_c N^c + (k^a k_a) B_{cd} N^c N^d \\
&= {B^a}_a + k^a (\nabla_a k_c) N^c + k^a (\nabla_c k_a) N^c \\
&= {B^a}_a
\end{align*}which follows since $k^a \nabla_c k_a = \dfrac{1}{2} \nabla_c(k^a k_a) = 0$.

Returning to the example at hand, the tensor ${B^a}_b = \nabla_b k^a$ is obtained by expanding in coordinate basis\begin{align*}
{B^{\mu}}_{\nu} = \nabla_{\nu} k^{\mu} &= \partial_{\nu} k^{\mu} + \Gamma^{\mu}_{\rho \nu} k^{\rho} \\
&= \partial_{\nu} k^{\mu} + \Gamma^{\mu}_{t \nu} + \Gamma^{\mu}_{r\nu} + \Gamma^{\mu}_{\phi\nu} \dfrac{a}{r^2 + a^2} \\
&= \Gamma^{\mu}_{r\nu} + \Gamma^{\mu}_{\phi \nu} \dfrac{a}{r^2 + a^2}
\end{align*}Since the traces of $B$ and $\hat{B}$ are equal, as explained, the expansion $\theta$ can be calculated directly from $B$:\begin{align*}
\theta &= \Gamma^{\mu}_{r\mu} + \Gamma^{\mu}_{\phi \mu} \dfrac{a}{r^2 + a^2}
\end{align*}Given that the metric determinant satisfies $\sqrt{-g} = \Sigma \sin{\theta}$, it follows that\begin{align*}
\Gamma^{\mu}_{r\mu} &= \dfrac{1}{\Sigma \sin{\theta}} \dfrac{\partial}{\partial r} (\Sigma \sin{\theta}) \\
&= \dfrac{2r}{\Sigma}
\end{align*}whilst $\dfrac{\partial}{\partial \phi}(\Sigma \sin{\theta}) = 0$ implies that $\Gamma^{\mu}_{\phi \mu} = 0$. Therefore the expansion of this congruence is $\theta = \dfrac{2r}{\Sigma}$. It can further be shown that the shear tensor of this congruence vanishes, whilst its rotation tensor satisfies $\omega^2 = \dfrac{2a^2 \cos^2{\theta}}{\Sigma^2}$; as mentioned earlier, the non-vanishing of $\omega$ indicates that the congruence is not hypersurface orthogonal.

Appendix: Evolution of $\theta$ (Raychaudhuri’s equation)

The derivative of the tensor $B$ along the geodesics of a timelike congruence satisfies\begin{align*}
\dfrac{DB_{ab}}{d\lambda} = u^c \nabla_c B_{ab} &= u^c \nabla_c \nabla_b u_a \\
&= u^c (\nabla_b \nabla_c u_a – R_{adbc}u^d) \\
&= \nabla_b (u^c \nabla_c u_a) – (\nabla_b u^c)(\nabla_c u_a) – R_{adbc}u^c u^d \\
&= – {B^c}_b B_{ac} – R_{acbd} u^c u^d
\end{align*}The trace of this equation is\begin{align*}
\dfrac{d\theta}{d\lambda} &= -B^{ca}B_{ac} – R_{cd}u^c u^d \\
&= -\left(\dfrac{1}{3} \theta h^{ca} + \sigma^{ca} + \omega^{ca} \right) \left(\dfrac{1}{3} \theta h_{ac} + \sigma_{ac} + \omega_{ac} \right) – R_{cd} u^c u^d \\
&= – \dfrac{1}{9} \theta^2 h^a_a – \sigma^{ca}\sigma_{ac} – \dfrac{2}{3}\theta h^{ca}\sigma_{ac} – \omega^{ca}\omega_{ac} – R_{cd}u^c u^d
\end{align*}since the cross terms with the antisymmetric tensor $\omega_{ab}$ vanish. Furthermore,\begin{align*}
h^{ca}\sigma_{ac} &= h^{ca}(B_{(ac)} – \dfrac{1}{3}\theta h_{ac}) \\
&= h^{ca} B_{(ac)} – \theta \\
&= (g^{ca} + u^c u^a)B_{(ac)} – \theta \\
&= {B^a}_a + \dfrac{1}{2} u^c u^a(\nabla_c u_a + \nabla_a u_c) – \theta \\
&= 0
\end{align*}using the geodesic equation; the equation therefore reduces to\begin{align*}
\dfrac{d\theta}{d\lambda} &= – \dfrac{1}{3}\theta^2 – \sigma^{ab}\sigma_{ab} + \omega^{ab}\omega_{ab} – R_{ab} u^a u^b \\
&\equiv \dfrac{1}{3}\theta^2 – \sigma^2 + \omega^2 – R_{ab}u^a u^b
\end{align*}An important corollary is the focusing theorem. If the strong energy condition, $R_{ab} u^a u^b \geq 0$, holds, and if the geodesic congruence is hypersurface orthogonal, $\omega = 0$, then\begin{align*}
\dfrac{d\theta}{d\lambda} = – \dfrac{1}{3}\theta^2 -\sigma^2 -R_{ab}u^a u^b
\end{align*}which means that the condition $\dfrac{d\theta}{d\lambda} \leq -\dfrac{1}{3}\theta^2$ also holds. Integrating yields\begin{align*}
\int \theta^{-2} d\theta &\leq – \dfrac{1}{3} \int d\lambda \\
\theta(0)^{-1} – \theta^{-1} &\leq -\dfrac{1}{3}\lambda \implies \theta \leq \dfrac{1}{\theta(0)^{-1} + \dfrac{1}{3}\lambda}
\end{align*}If $\theta(0) < 0$, then the denominator on the RHS approaches $0$ from the left at affine parameter $3|\theta(0)|^{-1}$, therefore $\theta$ approaches $-\infty$, within finite affine parameter $\lambda \leq 3|\theta(0)|^{-1}$; this implies the geodesics eventually converge to a caustic.

References

[1] H. Reall, Part 3 Black Holes (2020)

[2] E. Poisson, A relativist’s toolkit: the mathematics of black hole mechanics, (Cambridge University Press, Cambridge, UK, 2004)

[3] R. M. Wald, General Relativity, (University of Chicago Press, Chicago, USA, 1984)

[4] S. W. Hawking and G. F. R. Ellis, The large scale structure of spacetime, (Cambridge University Press, Cambridge, UK, 1973)