Kerr_spher_null_geod.ipynb

Kernel: SageMath 9.2

Spherical null geodesics in Kerr spacetime

Computation with kerrgeodesic_gw

This Jupyter/SageMath notebook is relative to the lectures Geometry and physics of black holes.

It requires SageMath (version $\geq$ 8.2), with the package kerrgeodesic_gw (version $\geq$ 0.3.2). To install the latter, simply run

sage -pip install kerrgeodesic_gw

in a terminal.

version()

'SageMath version 9.2, Release Date: 2020-10-24'

First, we set up the notebook to use LaTeX-formatted display:

%display latex

and we ask for CPU demanding computations to be performed in parallel on 8 processes:

Parallelism().set(nproc=8)

A Kerr black bole is entirely defined by two parameters $(m, a)$, where $m$ is the black hole mass and $a$ is the black hole angular momentum divided by $m$. In this notebook, we shall set $m=1$ and we denote the angular momentum parameter $a$ by the symbolic variable a, using a0 for a specific numerical value:

a = var('a')
a0 = 0.95

The spacetime object is created as an instance of the class KerrBH:

from kerrgeodesic_gw import KerrBH
M = KerrBH(a)
print(M)

Kerr spacetime M

The Boyer-Lindquist coordinate $r$ of the event horizon:

rH = M.event_horizon_radius()
rH

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\sqrt{-a^{2} + 1} + 1$$

rH0 = rH.subs({a: a0})
show(LatexExpr(r'r_+ = '), rH0)

$$\renewcommand{\Bold}[1]{\mathbf{#1}}r_+ = 1.31224989991992$$

show(LatexExpr(r'r_- = '),
     M.inner_horizon_radius().subs({a: a0}))

$$\renewcommand{\Bold}[1]{\mathbf{#1}}r_- = 0.687750100080080$$

The method boyer_lindquist_coordinates() returns the chart of Boyer-Lindquist coordinates BL and allows the user to instanciate the Python variables (t, r, th, ph) to the coordinates $(t,r,\theta,\phi)$:

BL.<t, r, th, ph> = M.boyer_lindquist_coordinates()
BL

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left(M,(t, r, {\theta}, {\phi})\right)$$

The metric tensor is naturally returned by the method metric():

g = M.metric()
g.display()

$$\renewcommand{\Bold}[1]{\mathbf{#1}}g = \left( -\frac{a^{2} \cos\left({\theta}\right)^{2} + r^{2} - 2 \, r}{a^{2} \cos\left({\theta}\right)^{2} + r^{2}} \right) \mathrm{d} t\otimes \mathrm{d} t + \left( -\frac{2 \, a r \sin\left({\theta}\right)^{2}}{a^{2} \cos\left({\theta}\right)^{2} + r^{2}} \right) \mathrm{d} t\otimes \mathrm{d} {\phi} + \left( \frac{a^{2} \cos\left({\theta}\right)^{2} + r^{2}}{a^{2} + r^{2} - 2 \, r} \right) \mathrm{d} r\otimes \mathrm{d} r + \left( a^{2} \cos\left({\theta}\right)^{2} + r^{2} \right) \mathrm{d} {\theta}\otimes \mathrm{d} {\theta} + \left( -\frac{2 \, a r \sin\left({\theta}\right)^{2}}{a^{2} \cos\left({\theta}\right)^{2} + r^{2}} \right) \mathrm{d} {\phi}\otimes \mathrm{d} t + \left( \frac{2 \, a^{2} r \sin\left({\theta}\right)^{4} + {\left(a^{2} r^{2} + r^{4} + {\left(a^{4} + a^{2} r^{2}\right)} \cos\left({\theta}\right)^{2}\right)} \sin\left({\theta}\right)^{2}}{a^{2} \cos\left({\theta}\right)^{2} + r^{2}} \right) \mathrm{d} {\phi}\otimes \mathrm{d} {\phi}$$

Spherical photon orbits

Functions $\ell(r_0)$ and $q(r_0)$ for spherical photon orbits:

r = var('r')
lsph(a, r) = (r^2*(3 - r) - a^2*(r + 1))/(a*(r -1))
lsph

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left( a, r \right) \ {\mapsto} \ -\frac{a^{2} {\left(r + 1\right)} + {\left(r - 3\right)} r^{2}}{a {\left(r - 1\right)}}$$

qsph(a, r) = r^3 / (a^2*(r - 1)^2) * (4*a^2 - r*(r - 3)^2)
qsph

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left( a, r \right) \ {\mapsto} \ -\frac{{\left({\left(r - 3\right)}^{2} r - 4 \, a^{2}\right)} r^{3}}{a^{2} {\left(r - 1\right)}^{2}}$$

$\theta$-turning points:

theta0(a, l, q) = arccos(sqrt(1/2*(1 - (l^2+q)/a^2 + sqrt((1 - (l^2+q)/a^2)^2 + 4*q/a^2))))
theta0

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left( a, l, q \right) \ {\mapsto} \ \arccos\left(\sqrt{\frac{1}{2} \, \sqrt{{\left(\frac{l^{2} + q}{a^{2}} - 1\right)}^{2} + \frac{4 \, q}{a^{2}}} - \frac{l^{2} + q}{2 \, a^{2}} + \frac{1}{2}}\right)$$

theta1(a, l, q) = arccos(sqrt(1/2*(1 - (l^2+q)/a^2 - sqrt((1 - (l^2+q)/a^2)^2 + 4*q/a^2))))
theta1

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left( a, l, q \right) \ {\mapsto} \ \arccos\left(\sqrt{-\frac{1}{2} \, \sqrt{{\left(\frac{l^{2} + q}{a^{2}} - 1\right)}^{2} + \frac{4 \, q}{a^{2}}} - \frac{l^{2} + q}{2 \, a^{2}} + \frac{1}{2}}\right)$$

Spherical photon orbit at $r_0 = 3 m$ ($q = q_{\rm max} = 27 m^2$)

r0 = 3.
E = 1
L = lsph(a0, r0)
Q = qsph(a0, r0)
L, Q

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left(-1.90000000000000, 27.0000000000000\right)$$

theta0(a0, L, Q)

$$\renewcommand{\Bold}[1]{\mathbf{#1}}0.345877348029357$$

P = M.point((0, r0, pi/2, 0), name='P')
print(P)

Point P on the Kerr spacetime M

A geodesic is constructed by providing the range $[\lambda_{\rm min},\lambda_{\rm max}]$ of the affine parameter $\lambda$, the initial point and either

• (i) the Boyer-Lindquist components $(p^t_0, p^r_0, p^\theta_0, p^\phi_0)$ of the initial 4-momentum vector $p_0 = \left. \frac{\mathrm{d}x}{\mathrm{d}\lambda}\right| _{\lambda_{\rm min}}$,

• (ii) the four integral of motions $(\mu, E, L, Q)$

• or (iii) some of the components of $p_0$ along with with some integrals of motion.

lmax = 100 # lambda_max

Li = M.geodesic([0, lmax], P, mu=0, E=E, L=L, Q=Q, a_num=a0,
                name='Li', latex_name=r'\mathcal{L}', verbose=True)

Initial tangent vector: 

$$\renewcommand{\Bold}[1]{\mathbf{#1}}p = 3.00000000000000 \frac{\partial}{\partial t } + \left( 9.36596633575423 \times 10^{-9} \right) \frac{\partial}{\partial r } + 0.577350269189626 \frac{\partial}{\partial {\theta} }$$

The curve was correctly set.
Parameters appearing in the differential system defining the curve are [a].

v0 = Li.initial_tangent_vector()
Li = M.geodesic([0, lmax], P, pt0=v0[0], pr0=0, pth0=v0[2], pph0=v0[3],
                a_num=a0, name='Li', latex_name=r'\mathcal{L}', verbose=True)

Initial tangent vector: 

$$\renewcommand{\Bold}[1]{\mathbf{#1}}p = 3.00000000000000 \frac{\partial}{\partial t } + 0.577350269189626 \frac{\partial}{\partial {\theta} }$$

The curve was correctly set.
Parameters appearing in the differential system defining the curve are [a].

The numerical integration of the geodesic equation is performed via integrate(), by providing the integration step $\delta\lambda$:

Li.integrate(step=0.001, method='dopri5')
Li.check_integrals_of_motion(0.999*lmax)

print("Final point: ")
show(BL[:], "=", BL(Li(0.999*lmax)))

Final point: 

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left(t, r, {\theta}, {\phi}\right) \verb|=| \left(291.1717434701057, 3.0, 2.144493861874287, -39.55667803565552\right)$$

lplot = 0.32*lmax
print("max lambda (plot): ", lplot)
Li.plot(prange=(0, lplot), plot_points=2000, thickness=3, color='green',
        display_tangent=True, scale=0.15, width_tangent=2, color_tangent='green', 
        plot_points_tangent=20, horizon_color='lightgrey') \
     + line([(0,0,-3), (0,0,3)], color='black', thickness=2)

max lambda (plot):  32.0000000000000

graph = Li.plot(coordinates='xy', prange=(0, lplot), plot_points=2000, 
                thickness=3, color='green', display_tangent=True, scale=0.15, 
                width_tangent=2, color_tangent='green', plot_points_tangent=20, 
                horizon_color='lightgrey', axes_labels=[r'$x/m$', r'$y/m$'])
graph.save("gik_spher_3d_r_30_xy.pdf")
graph

Li.plot(coordinates='xz', prange=(0, lplot), plot_points=2000, 
        thickness=3, color='green', display_tangent=True, scale=0.15, 
        width_tangent=2, color_tangent='green', plot_points_tangent=20, 
        horizon_color='lightgrey', axes_labels=[r'$x/m$', r'$z/m$'])

Prograde spherical photon orbit at $r_0 = 1.6m$

r0 = 1.6
L = lsph(a0, r0)
Q = qsph(a0, r0)
L, Q

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left(2.17105263157895, 5.97569713758080\right)$$

theta0(a0, L, Q)

$$\renewcommand{\Bold}[1]{\mathbf{#1}}0.705442812649839$$

P = M.point((0, r0, pi/2, 0), name='P')
lmax = 70
Li = M.geodesic([0, lmax], P, mu=0, E=E, L=L, Q=Q, a_num=a0,
                name='Li', latex_name=r'\mathcal{L}', verbose=True)

Initial tangent vector: 

$$\renewcommand{\Bold}[1]{\mathbf{#1}}p = 7.66666666666666 \frac{\partial}{\partial t } + \left( 1.74608999691187 \times 10^{-24} + 2.85158136717879 \times 10^{-8}i \right) \frac{\partial}{\partial r } + 0.954892151626195 \frac{\partial}{\partial {\theta} } + 2.45614035087719 \frac{\partial}{\partial {\phi} }$$

The curve was correctly set.
Parameters appearing in the differential system defining the curve are [a].

v0 = Li.initial_tangent_vector()
Li = M.geodesic([0, lmax], P, pt0=v0[0], pr0=0, pth0=v0[2], pph0=v0[3],
                a_num=a0, name='Li', latex_name=r'\mathcal{L}', verbose=True)

Initial tangent vector: 

$$\renewcommand{\Bold}[1]{\mathbf{#1}}p = 7.66666666666666 \frac{\partial}{\partial t } + 0.954892151626195 \frac{\partial}{\partial {\theta} } + 2.45614035087719 \frac{\partial}{\partial {\phi} }$$

The curve was correctly set.
Parameters appearing in the differential system defining the curve are [a].

Li.integrate(step=0.0004, method='dopri5')
Li.check_integrals_of_motion(0.999*lmax)

print("Final point: ")
show(BL[:], "=", BL(Li(0.999*lmax)))

Final point: 

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left(t, r, {\theta}, {\phi}\right) \verb|=| \left(492.9599102672748, 1.599974600261711, 0.8276182451833887, 185.0198306183178\right)$$

lplot = 0.11*lmax
print("max lambda (plot): ", lplot)
Li.plot(prange=(0, lplot), plot_points=2000, thickness=3, color='green',
        display_tangent=True, scale=0.03, width_tangent=1, color_tangent='green', 
        plot_points_tangent=20, horizon_color='lightgrey') \
    + line([(0,0,-1.5), (0,0,1.5)], color='black', thickness=2)

max lambda (plot):  7.70000000000000

Li.plot(prange=(0, lmax), plot_points=2000, thickness=3, color='green',
        display_tangent=True, scale=0.03, width_tangent=1, color_tangent='green', 
        plot_points_tangent=40, horizon_color='lightgrey') \
    + line([(0,0,-1.5), (0,0,1.5)], color='black', thickness=2)

Spherical photon orbit at $r_0=2.8m$

r0 = 2.8
L = lsph(a0, r0)
Q = qsph(a0, r0)
L, Q

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left(-1.08859649122807, 26.2604206422489\right)$$

theta0(a0, L, Q)

$$\renewcommand{\Bold}[1]{\mathbf{#1}}0.206050206829550$$

P = M.point((0, r0, pi/2, 0), name='P')
lmax = 100
Li = M.geodesic([0, lmax], P, mu=0, E=E, L=L, Q=Q, a_num=a0,
                name='Li', latex_name=r'\mathcal{L}', verbose=True)

Initial tangent vector: 

$$\renewcommand{\Bold}[1]{\mathbf{#1}}p = 3.22222222222222 \frac{\partial}{\partial t } + \left( 4.65527024480212 \times 10^{-25} + 7.60263326216718 \times 10^{-9}i \right) \frac{\partial}{\partial r } + 0.653634213345212 \frac{\partial}{\partial {\theta} } + 0.116959064327486 \frac{\partial}{\partial {\phi} }$$

The curve was correctly set.
Parameters appearing in the differential system defining the curve are [a].

v0 = Li.initial_tangent_vector()
Li = M.geodesic([0, lmax], P, pt0=v0[0], pr0=0, pth0=v0[2], pph0=v0[3],
                a_num=a0, name='Li', latex_name=r'\mathcal{L}', verbose=True)

Initial tangent vector: 

$$\renewcommand{\Bold}[1]{\mathbf{#1}}p = 3.22222222222222 \frac{\partial}{\partial t } + 0.653634213345212 \frac{\partial}{\partial {\theta} } + 0.116959064327486 \frac{\partial}{\partial {\phi} }$$

The curve was correctly set.
Parameters appearing in the differential system defining the curve are [a].

Li.integrate(step=0.001, method='dopri5')
Li.check_integrals_of_motion(0.999*lmax)

print("Final point: ")
show(BL[:], "=", BL(Li(0.999*lmax)))

Final point: 

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left(t, r, {\theta}, {\phi}\right) \verb|=| \left(310.4228769082733, 2.8, 2.462502125969338, -39.07005458983995\right)$$

lplot = 0.38*lmax
print("max lambda (plot): ", lplot)
Li.plot(prange=(0, lplot), plot_points=2000, thickness=3, color='green',
        display_tangent=True, scale=0.12, width_tangent=2, color_tangent='green', 
        plot_points_tangent=20, horizon_color='lightgrey')

max lambda (plot):  38.0000000000000

graph = Li.plot(coordinates='xy', prange=(0, lplot), plot_points=2000, 
                thickness=3, color='green', display_tangent=True, scale=0.01, 
                width_tangent=1, color_tangent='green', 
                plot_points_tangent=20, horizon_color='lightgrey',
                axes_labels=[r'$x/m$', r'$y/m$'])
graph.save("gik_spher_3d_r_28_xy.pdf")
graph

Li.plot(coordinates='txy', prange=(0, lplot), plot_points=2000, 
        thickness=3, color='green', display_tangent=True, scale=0.2, 
        width_tangent=2, color_tangent='green', 
        plot_points_tangent=20, horizon_color='lightgrey')

Retrograde spherical photon orbit at $r_0 = 3.9 m$

r0 = 3.9
L = lsph(a0, r0)
Q = qsph(a0, r0)
L, Q

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left(-6.57395644283122, 3.52474056409564\right)$$

theta0(a0, L, Q)

$$\renewcommand{\Bold}[1]{\mathbf{#1}}1.29003447185548$$

P = M.point((0, r0, pi/2, 0), name='P')

lmax = 100
Li = M.geodesic([0, lmax], P, mu=0, E=E, L=L, Q=Q, a_num=a0,
                name='Li', latex_name=r'\mathcal{L}', verbose=True)

Initial tangent vector: 

$$\renewcommand{\Bold}[1]{\mathbf{#1}}p = 2.37931034482759 \frac{\partial}{\partial t } + \left( 2.22112222707613 \times 10^{-8} \right) \frac{\partial}{\partial r } + 0.123433875307859 \frac{\partial}{\partial {\theta} } -0.326678765880218 \frac{\partial}{\partial {\phi} }$$

The curve was correctly set.
Parameters appearing in the differential system defining the curve are [a].

v0 = Li.initial_tangent_vector()
Li = M.geodesic([0, lmax], P, pt0=v0[0], pr0=0, pth0=v0[2], pph0=v0[3],
                a_num=a0, name='Li', latex_name=r'\mathcal{L}', verbose=True)

Initial tangent vector: 

$$\renewcommand{\Bold}[1]{\mathbf{#1}}p = 2.37931034482759 \frac{\partial}{\partial t } + 0.123433875307859 \frac{\partial}{\partial {\theta} } -0.326678765880218 \frac{\partial}{\partial {\phi} }$$

The curve was correctly set.
Parameters appearing in the differential system defining the curve are [a].

Li.integrate(step=0.001, method='dopri5')
Li.check_integrals_of_motion(0.999*lmax)

print("Final point: ")
show(BL[:], "=", BL(Li(0.999*lmax)))

Final point: 

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left(t, r, {\theta}, {\phi}\right) \verb|=| \left(237.38270015383975, 3.9, 1.6881755637536546, -34.301693666892874\right)$$

lplot = 0.55*lmax
print("max lambda (plot): ", lplot)
Li.plot(prange=(0, lplot), plot_points=2000, thickness=3, color='green',
        display_tangent=True, scale=0.15, width_tangent=2, color_tangent='green', 
        plot_points_tangent=20, horizon_color='lightgrey')

max lambda (plot):  55.0000000000000

graph = Li.plot(coordinates='xy', prange=(0, lplot), plot_points=500, 
                thickness=3, color='green', display_tangent=True, scale=0.01, 
                width_tangent=1, color_tangent='green', 
                plot_points_tangent=20, horizon_color='lightgrey',
                axes_labels=[r'$x/m$', r'$y/m$'])
graph.save("gik_spher_3d_r_39_xy.pdf")
graph

A polar spherical photon orbit

rph_pol(a) = 1 + 2*sqrt(1 - a^2/3)*cos(1/3*arccos((1 - a^2)/(1 - a^2/3)^(3/2)))
rph_pol

$$\renewcommand{\Bold}[1]{\mathbf{#1}}a \ {\mapsto}\ 2 \, \sqrt{-\frac{1}{3} \, a^{2} + 1} \cos\left(\frac{1}{3} \, \arccos\left(-\frac{a^{2} - 1}{{\left(-\frac{1}{3} \, a^{2} + 1\right)}^{\frac{3}{2}}}\right)\right) + 1$$

r0 = rph_pol(a0)
L = lsph(a0, r0)
Q = qsph(a0, r0)
r0, L, Q

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left(2.49269429554008, -6.26333640524599 \times 10^{-16}, 22.8640201857508\right)$$

theta0(a0, L, Q)

$$\renewcommand{\Bold}[1]{\mathbf{#1}}0.000000000000000$$

lmax = 30
P = M.point((0, r0, pi/2, 0), name='P')
Li = M.geodesic([0, lmax], P, mu=0, E=E, L=0, Q=Q, a_num=a0,
                name='Li', latex_name=r'\mathcal{L}', verbose=True)

Initial tangent vector: 

$$\renewcommand{\Bold}[1]{\mathbf{#1}}p = 3.67971815257238 \frac{\partial}{\partial t } + \left( 1.31343299672075 \times 10^{-24} + 2.14499886438313 \times 10^{-8}i \right) \frac{\partial}{\partial r } + 0.769552507817167 \frac{\partial}{\partial {\theta} } + 0.357746396090726 \frac{\partial}{\partial {\phi} }$$

The curve was correctly set.
Parameters appearing in the differential system defining the curve are [a].

v0 = Li.initial_tangent_vector()
Li = M.geodesic([0, lmax], P, pt0=v0[0], pr0=0, pth0=v0[2], pph0=v0[3],
                a_num=a0, name='Li', latex_name=r'\mathcal{L}', verbose=True)

Initial tangent vector: 

$$\renewcommand{\Bold}[1]{\mathbf{#1}}p = 3.67971815257238 \frac{\partial}{\partial t } + 0.769552507817167 \frac{\partial}{\partial {\theta} } + 0.357746396090726 \frac{\partial}{\partial {\phi} }$$

The curve was correctly set.
Parameters appearing in the differential system defining the curve are [a].

Li.integrate(step=0.002, method='dopri5')
#Li.check_integrals_of_motion(0.99*lmax)

/home/eric/sage/9.2/local/lib/python3.8/site-packages/scipy/integrate/_ode.py:1179: UserWarning: dopri5: step size becomes too small
  warnings.warn('{:s}: {:s}'.format(self.__class__.__name__,

lplot = 0.5*lmax
print("max lambda (plot): ", lplot)
Li.plot(prange=(0, lplot), plot_points=2000, thickness=3, color='green',
        display_tangent=True, scale=0.2, width_tangent=2, color_tangent='green', 
        plot_points_tangent=20, horizon_color='lightgrey')

max lambda (plot):  15.0000000000000

Inner spherical orbits

To plot the inner orbits, and in particular orbits with $r<0$, we use a map from Kerr spacetime to the Euclidean space based on the radial coordinate $$\hat{r} := \frac{1}{2} \left( r + \sqrt{r^2 + 4m^2} \right)$$ instead of the Boyer-Lindquist $r$.

E4 = M.map_to_Euclidean().codomain()
X4 = E4.cartesian_coordinates()
t, x, y, z = X4[:]
X4

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left(\mathbb{E}^{4},(t, x, y, z)\right)$$

rhat = 1/2*(r + sqrt(r^2 + 4))
rhat

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\frac{1}{2} \, r + \frac{1}{2} \, \sqrt{r^{2} + 4}$$

F = M.diff_map(E4, {(BL, X4): [t, rhat*sin(th)*cos(ph),
                               rhat*sin(th)*sin(ph),
                               rhat*cos(th)]}, name='F')
F.display()

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\begin{array}{llcl} F:& M & \longrightarrow & \mathbb{E}^{4} \\ & \left(t, r, {\theta}, {\phi}\right) & \longmapsto & \left(t, x, y, z\right) = \left(t, \frac{1}{2} \, {\left(r + \sqrt{r^{2} + 4}\right)} \cos\left({\phi}\right) \sin\left({\theta}\right), \frac{1}{2} \, {\left(r + \sqrt{r^{2} + 4}\right)} \sin\left({\phi}\right) \sin\left({\theta}\right), \frac{1}{2} \, {\left(r + \sqrt{r^{2} + 4}\right)} \cos\left({\theta}\right)\right) \end{array}$$

Marginally stable spherical photon orbit

r0 = 1 - (1 - a0^2)^(1/3)
r0

$$\renewcommand{\Bold}[1]{\mathbf{#1}}0.539741795874205$$

L = lsph(a0, r0)
Q = qsph(a0, r0)
L, Q

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left(1.53893384905757, 0.282109845505914\right)$$

theta0(a0, L, Q)

$$\renewcommand{\Bold}[1]{\mathbf{#1}}1.17336440458345$$

Since all orbits with $0< r_0 \leq r_{\rm ph}^*$ (such as the marginally stable spherical orbit) have $E<0$, we set E = -1 and L = -L:

lmax = 50
P = M.point((0, r0, pi/2, 0), name='P')
Li = M.geodesic([0, lmax], P, mu=0, E=-1, L=-L, Q=Q, a_num=a0,
                name='Li', latex_name=r'\mathcal{L}', verbose=True)

Initial tangent vector: 

$$\renewcommand{\Bold}[1]{\mathbf{#1}}p = 7.69077392677338 \frac{\partial}{\partial t } + \left( 3.13205131833926 \times 10^{-24} + 5.11502797462896 \times 10^{-8}i \right) \frac{\partial}{\partial r } + 1.82321137638829 \frac{\partial}{\partial {\theta} } + 5.62672311935440 \frac{\partial}{\partial {\phi} }$$

The curve was correctly set.
Parameters appearing in the differential system defining the curve are [a].

v0 = Li.initial_tangent_vector()
Li = M.geodesic([0, lmax], P, pt0=v0[0], pr0=0, pth0=v0[2], pph0=v0[3],
                a_num=a0, name='Li', latex_name=r'\mathcal{L}', verbose=True)

Initial tangent vector: 

$$\renewcommand{\Bold}[1]{\mathbf{#1}}p = 7.69077392677338 \frac{\partial}{\partial t } + 1.82321137638829 \frac{\partial}{\partial {\theta} } + 5.62672311935440 \frac{\partial}{\partial {\phi} }$$

The curve was correctly set.
Parameters appearing in the differential system defining the curve are [a].

Li.integrate(step=0.001, method='dopri5')
Li.check_integrals_of_motion(0.999*lmax)

print("Final point: ")
show(BL[:], "=", BL(Li(0.999*lmax)))

Final point: 

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left(t, r, {\theta}, {\phi}\right) \verb|=| \left(302.91453819062644, 0.5397417961010423, 1.5090634888129522, 210.47540919172215\right)$$

Plot with respect to $\hat{r}$ coordinate:

sing_ring = circle((0,0,0), 1, thickness=3, color='orangered')
r_minf = sphere((0,0,0), 0.03, color='black') \
         + text3d("r=\u2212\u221E", (0,0,0.07), fontfamily='serif', 
                  fontstyle='italic', fontsize=20)

lplot = 0.06*lmax
print("max lambda (plot): ", lplot)
graph = Li.plot_integrated(chart=X4, mapping=F,
                   ambient_coords=(x,y,z), prange=(0, lplot),
                   thickness=3, plot_points=2000, color='green',
                   display_tangent=True, scale=0.012, width_tangent=1, 
                   color_tangent='green', plot_points_tangent=20,
                   label_axes=False)
graph += sing_ring + r_minf \
         + line([(0,0,-0.53), (0,0,0.53)], color='black', thickness=2)
graph._extra_kwds['axes_labels'] = ['\u0302x', '\u0302y', '\u0302z']
graph

max lambda (plot):  3.00000000000000

lplot = 0.4*lmax
print("max lambda (plot): ", lplot)
graph = Li.plot_integrated(chart=X4, mapping=F,
                   ambient_coords=(x,y,z), prange=(0, lplot),
                   thickness=3, plot_points=2000, color='green',
                   display_tangent=True, scale=0.012, width_tangent=1, 
                   color_tangent='green', plot_points_tangent=20,
                   label_axes=False)
graph += sing_ring + r_minf \
         + line([(0,0,-0.53), (0,0,0.53)], color='black', thickness=2)
graph._extra_kwds['axes_labels'] = ['\u0302x', '\u0302y', '\u0302z']
graph

max lambda (plot):  20.0000000000000

Stable inner spherical photon orbit at $r_0 = 0.5 m$

r0 = 0.5
L = lsph(a0, r0)
Q = qsph(a0, r0)
L, Q

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left(1.53421052631579, 0.268698060941828\right)$$

theta0(a0, L, Q)

$$\renewcommand{\Bold}[1]{\mathbf{#1}}1.17991588772518$$

lmax = 30
P = M.point((0, r0, pi/2, 0), name='P')
Li = M.geodesic([0, lmax], P, mu=0, E=-1, L=-L, Q=Q, a_num=a0,
                name='Li', latex_name=r'\mathcal{L}', verbose=True)

Initial tangent vector: 

$$\renewcommand{\Bold}[1]{\mathbf{#1}}p = 7.00000000000000 \frac{\partial}{\partial t } + 2.07344374774655 \frac{\partial}{\partial {\theta} } + 5.26315789473685 \frac{\partial}{\partial {\phi} }$$

The curve was correctly set.
Parameters appearing in the differential system defining the curve are [a].

v0 = Li.initial_tangent_vector()
Li = M.geodesic([0, lmax], P, pt0=v0[0], pr0=0, pth0=v0[2], pph0=v0[3],
                a_num=a0, name='Li', latex_name=r'\mathcal{L}', verbose=True)

Initial tangent vector: 

$$\renewcommand{\Bold}[1]{\mathbf{#1}}p = 7.00000000000000 \frac{\partial}{\partial t } + 2.07344374774655 \frac{\partial}{\partial {\theta} } + 5.26315789473685 \frac{\partial}{\partial {\phi} }$$

The curve was correctly set.
Parameters appearing in the differential system defining the curve are [a].

Li.integrate(step=0.001, method='dopri5')
Li.check_integrals_of_motion(0.999*lmax)

print("Final point: ")
show(BL[:], "=", BL(Li(0.999*lmax)))

Final point: 

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left(t, r, {\theta}, {\phi}\right) \verb|=| \left(160.26568413055006, 0.499999999999837, 1.4610686498265664, 113.33060855013213\right)$$

Plot in (Cartesian) Boyer-Lindquist coordinates:

Li.plot(prange=(0, lmax), plot_points=2000, thickness=3, color='green',
        display_tangent=True, scale=0.01, width_tangent=1, color_tangent='green', 
        plot_points_tangent=20, plot_horizon=False)

Plot with respect to $\hat{r}$ coordinate:

lplot = 0.25*lmax
print("max lambda (plot): ", lplot)
graph1 = Li.plot_integrated(chart=X4, mapping=F,
                   ambient_coords=(x,y,z), prange=(0, lplot),
                   thickness=3, plot_points=2000, color='green',
                   display_tangent=True, scale=0.015, width_tangent=1, 
                   color_tangent='green', plot_points_tangent=20,
                   label_axes=False)
graph1 += sing_ring + r_minf
graph1._extra_kwds['axes_labels'] = ['\u0302x', '\u0302y', '\u0302z']
graph1

max lambda (plot):  7.50000000000000

Vortical circular photon orbit at $r_0 = r_{\rm ph}^{**}$

rph_ss = n(1/2 + cos(2/3*asin(a0) + 2*pi/3))
rph_ss

$$\renewcommand{\Bold}[1]{\mathbf{#1}}-0.477673658836338$$

L = lsph(a0, rph_ss)
Q = qsph(a0, rph_ss)
L, Q 

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left(-0.229456449433387, -0.519183008263141\right)$$

th0 = theta0(a0, L, Q)
th1 = theta1(a0, L, Q)
th0, th1

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left(0.513765577758655, 0.513765602376411\right)$$

th_ss = arcsin(sqrt(abs(L)/a0))
th_ss

$$\renewcommand{\Bold}[1]{\mathbf{#1}}0.513765590067533$$

rh_ss = rhat.subs({r: rph_ss})
rh_ss

$$\renewcommand{\Bold}[1]{\mathbf{#1}}0.789289150745023$$

lmax = 10
P = M.point((0, rph_ss, th_ss, 0), name='P')
Li = M.geodesic([0, lmax], P, mu=0, E=E, L=L, Q=Q, a_num=a0,
                name='Li', latex_name=r'\mathcal{L}', verbose=True)

Initial tangent vector: 

$$\renewcommand{\Bold}[1]{\mathbf{#1}}p = 0.323260590036169 \frac{\partial}{\partial t } + \left( 1.63266670244580 \times 10^{-8} \right) \frac{\partial}{\partial r } + \left( 1.99960017051231 \times 10^{-8} \right) \frac{\partial}{\partial {\theta} } -1.40881021577044 \frac{\partial}{\partial {\phi} }$$

The curve was correctly set.
Parameters appearing in the differential system defining the curve are [a].

v0 = Li.initial_tangent_vector()
Li = M.geodesic([0, lmax], P, pt0=v0[0], pr0=0, pth0=0, pph0=v0[3],
                a_num=a0, name='Li', latex_name=r'\mathcal{L}', verbose=True)

Initial tangent vector: 

$$\renewcommand{\Bold}[1]{\mathbf{#1}}p = 0.323260590036169 \frac{\partial}{\partial t } -1.40881021577044 \frac{\partial}{\partial {\phi} }$$

The curve was correctly set.
Parameters appearing in the differential system defining the curve are [a].

Li.integrate(step=0.001, method='dopri5')
Li.check_integrals_of_motion(0.999*lmax)

show("Final point: ", BL[:], "=", BL(Li(0.999*lmax)))

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\verb|Final|\phantom{\verb!x!}\verb|point:| \left(t, r, {\theta}, {\phi}\right) \verb|=| \left(3.229373294461296, -0.47767365883633794, 0.513765590067533, -14.074014055546579\right)$$

lplot = 0.444*lmax
print("max lambda (plot): ", lplot)
vort_circ = Li.plot_integrated(chart=X4, mapping=F,
                   ambient_coords=(x,y,z), prange=(0, lplot),
                   thickness=3, plot_points=2000, color='lightgreen',
                   display_tangent=True, scale=0.1, width_tangent=1, 
                   color_tangent='lightgreen', plot_points_tangent=5,
                   label_axes=False)
graph = vort_circ + sing_ring + r_minf \
        + line([(0,0,0), (0,0,1.3*rh_ss)], color='black', thickness=2)
graph._extra_kwds['axes_labels'] = ['\u0302x', '\u0302y', '\u0302z']
graph

max lambda (plot):  4.44000000000000

Vortical inner spherical photon orbit at $r_0 = -0.46 m$

r0 = -0.46
L = lsph(a0, r0)
Q = qsph(a0, r0)
L, Q

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left(-0.176485940879596, -0.461285155596124\right)$$

th0 = theta0(a0, L, Q)
th1 = theta1(a0, L, Q)
th0, th1

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left(0.281889824619659, 0.731306339672501\right)$$

lmax = 20
P = M.point((0, r0, (th0+th1)/2, 0), name='P')
Li = M.geodesic([0, lmax], P, mu=0, E=E, L=L, Q=Q, a_num=a0,
                name='Li', latex_name=r'\mathcal{L}', verbose=True)

Initial tangent vector: 

$$\renewcommand{\Bold}[1]{\mathbf{#1}}p = 0.357024161116205 \frac{\partial}{\partial t } + \left( 1.84776601065215 \times 10^{-8} \right) \frac{\partial}{\partial r } + 0.396167049160545 \frac{\partial}{\partial {\theta} } -1.22114489356954 \frac{\partial}{\partial {\phi} }$$

The curve was correctly set.
Parameters appearing in the differential system defining the curve are [a].

v0 = Li.initial_tangent_vector()
Li = M.geodesic([0, lmax], P, pt0=v0[0], pr0=0, pth0=v0[2], pph0=v0[3],
                a_num=a0, name='Li', latex_name=r'\mathcal{L}', verbose=True)

Initial tangent vector: 

$$\renewcommand{\Bold}[1]{\mathbf{#1}}p = 0.357024161116205 \frac{\partial}{\partial t } + 0.396167049160545 \frac{\partial}{\partial {\theta} } -1.22114489356954 \frac{\partial}{\partial {\phi} }$$

The curve was correctly set.
Parameters appearing in the differential system defining the curve are [a].

Li.integrate(step=0.001, method='dopri5')
Li.check_integrals_of_motion(0.999*lmax)

print("Final point: ")
show(BL[:], "=", BL(Li(0.999*lmax)))

Final point: 

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left(t, r, {\theta}, {\phi}\right) \verb|=| \left(6.6100475583006695, -0.45999868463861343, 0.39359482789056505, -28.895385729902912\right)$$

lplot = lmax
print("max lambda (plot): ", lplot)
graph = Li.plot_integrated(chart=X4, mapping=F,
                   ambient_coords=(x,y,z), prange=(0, lplot),
                   thickness=3, plot_points=2000, color='green',
                   display_tangent=True, scale=0.1, width_tangent=1, 
                   color_tangent='green', plot_points_tangent=20,
                   label_axes=False)

max lambda (plot):  20

graph += sing_ring + r_minf + vort_circ \
         + line([(0,0,0), (0,0,0.9)], color='black', thickness=2)
graph._extra_kwds['axes_labels'] = ['\u0302x', '\u0302y', '\u0302z']
graph

graph + graph1