Kerr_spher_photon_existence.ipynb

Kernel: SageMath 9.3.beta5

Existence of spherical photon orbits in Kerr spacetime

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

version()

'SageMath version 9.3.beta5, Release Date: 2020-12-27'

%display latex

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

We use $m=1$ and denote $r_0$ simply by $r$.

a, r = var('a 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}}$$

The radii of the two horizons:

rp(a) = 1 + sqrt(1 - a^2)
rm(a) = 1 - sqrt(1 - a^2)

Critical radii $r_{\rm ph}^{**}$, $r_{\rm ph}^*$, $r_{\rm ph}^+$, $r_{\rm ph}^-$, $r_{\rm ph}^{\rm ms}$ and $r_{\rm ph}^{\rm pol}$

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

$$\renewcommand{\Bold}[1]{\mathbf{#1}}a \ {\mapsto}\ \cos\left(\frac{2}{3} \, \pi + \frac{2}{3} \, \arcsin\left(a\right)\right) + \frac{1}{2}$$

rph_s(a) = 4*cos(acos(-a)/3 + 4*pi/3)^2
rph_s

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

rph_p(a) = 4*cos(acos(-a)/3)^2
rph_p

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

rph_m(a) = 4*cos(acos(a)/3)^2
rph_m

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

We add the radius of the marginally stable orbit:

rph_ms(a) = 1 - (1 - a^2)^(1/3)
rph_ms

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

as well as the radius of outer and inner polar orbits:

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$$

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

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

Plots for $a=0.95 \, m$

a0 = 0.95
show(LatexExpr(r'r_{\rm ph}^{**} = '), n(rph_ss(a0)))
show(LatexExpr(r'r_{\rm ph}^{\rm ms} = '), n(rph_ms(a0)))
show(LatexExpr(r'r_{\rm ph}^{*} = '), n(rph_s(a0)))
show(LatexExpr(r'r_- = '), n(rm(a0)))
show(LatexExpr(r'r_+ = '), n(rp(a0)))
show(LatexExpr(r'r_{\rm ph}^+ = '), n(rph_p(a0)))
show(LatexExpr(r'r_{\rm ph}^- = '), n(rph_m(a0)))

$$\renewcommand{\Bold}[1]{\mathbf{#1}}r_{\rm ph}^{**} = -0.477673658836338$$

$$\renewcommand{\Bold}[1]{\mathbf{#1}}r_{\rm ph}^{\rm ms} = 0.539741795874205$$

$$\renewcommand{\Bold}[1]{\mathbf{#1}}r_{\rm ph}^{*} = 0.658372153864346$$

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

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

$$\renewcommand{\Bold}[1]{\mathbf{#1}}r_{\rm ph}^+ = 1.38628052846298$$

$$\renewcommand{\Bold}[1]{\mathbf{#1}}r_{\rm ph}^- = 3.95534731767268$$

gq = plot(qsph(a0, r), (r, -2, 0.9), thickness=2, legend_label=r'$q$', 
          axes_labels=[r'$r_0/m$', r'$\ell/m,\ q/m^2$'],
          frame=True, gridlines=True) \
     + plot(qsph(a0, r), (r, 1.1, 5), thickness=2)
show(gq, ymin=-4)

gl = plot(lsph(a0, r), (r, -2, 0.99), color='red', thickness=2, legend_label=r'$\ell$') \
     + plot(lsph(a0, r), (r, 1.002, 5), color='red', thickness=2)

show(gq+gl, ymin=-5, ymax=10)

def qmin(a1, r):
    l0 = abs(lsph(a1, r))
    if l0 < a1:
        return -(a1 - l0)^2
    return 0

gqmin = plot(lambda r: qmin(a0, r), (r, -2, 5), color='grey', thickness=2, fill=-8)

hor = line([(rp(a0), -10), (rp(a0), 30)], color='black', thickness=2) \
      + line([(rm(a0), -10), (rm(a0), 30)], color='peru', thickness=2)
llim = line([(-2, a0), (5, a0)], color='red', thickness=1.5, 
            linestyle=':', legend_label=r'$|\ell| = a$') \
       + line([(-2, -a0), (5, -a0)], color='red', thickness=1.5, 
              linestyle=':')

graph = gq + gl + gqmin + hor + llim
for r1 in [rph_ss(a0), rph_s(a0), rph_p(a0), rph_m(a0)]:
    graph += point((r1, qsph(a0, r1)), size=60, color='blue', zorder=100)
graph += text(r'$r_{\rm ph}^{**}$', (rph_ss(a0) + 0.1, 2), fontsize=14, zorder=101)
graph += text(r'$r_{\rm ph}^*$', (rph_s(a0) + 0.23, 2), fontsize=14, zorder=101)
graph += text(r'$r_{\rm ph}^+$', (rph_p(a0) + 0.25, 2), fontsize=14, zorder=101)
graph += text(r'$r_{\rm ph}^-$', (rph_m(a0) + 0.2, 2), fontsize=14, zorder=101)
graph.set_legend_options(handlelength=2)
show(graph, xmin=-1.5, xmax=4.5, ymin=-7, ymax=28)
graph.save('gik_spher_orb_exist.pdf', xmin=-1.5, xmax=4.5, ymin=-7, ymax=28)

Zoom:

gq_stable = plot(qsph(a0, r), (r, 0, rph_ms(a0)), thickness=5) 
gl_stable = plot(lsph(a0, r), (r, 0, rph_ms(a0)), thickness=5, color='red') 
graph = gq + gq_stable + gl + gl_stable + gqmin + hor + llim
for r1 in [rph_ss(a0), rph_s(a0), rph_p(a0), rph_m(a0)]:
    graph += point((r1, qsph(a0, r1)), size=60, color='blue', zorder=100)
graph += text(r'$r_{\rm ph}^*$', (rph_s(a0) + 0.12, 0.15), fontsize=14, zorder=101)
graph += text(r'$r_{\rm ph}^+$', (rph_p(a0) + 0.09, 0.15), fontsize=14, zorder=101)
r1 = rph_ss(a0)
graph += text(r'$r_{\rm ph}^{**}$', (r1, 0.15), fontsize=14, zorder=101)
graph += line([(r1, 0), (r1, qsph(a0, r1))], linestyle='--', color='blue')
r1 = rph_ms(a0)
graph += text(r'$r_{\rm ph}^{\rm ms}$', (r1, -0.2), color='black',
              fontsize=14, zorder=101)
graph += line([(r1, 0), (r1, lsph(a0, r1))], linestyle='--', color='black')
graph.set_legend_options(handlelength=2)
show(graph, xmin=-1, xmax=1.5, ymin=-1.5, ymax=2)
graph.save('gik_spher_orb_exist_zoom.pdf', xmin=-1, xmax=1.5, ymin=-1.5, ymax=2)

Plots for $a=0.5 \, m$

a0 = 0.5

gq = plot(qsph(a0, r), (r, -2, 0.9), thickness=2, legend_label=r'$q$', 
          axes_labels=[r'$r_0/m$', r'$\ell/m,\ q/m^2$'],
          frame=True, gridlines=True) \
     + plot(qsph(a0, r), (r, 1.1, 5), thickness=2)
gl = plot(lsph(a0, r), (r, -2, 0.99), color='red', thickness=2, legend_label=r'$\ell$') \
     + plot(lsph(a0, r), (r, 1.002, 5), color='red', thickness=2)
gqmin = plot(lambda r: qmin(a0, r), (r, -2, 5), color='grey', thickness=2, fill=-8)
hor = line([(rp(a0), -10), (rp(a0), 30)], color='black', thickness=2) \
      + line([(rm(a0), -10), (rm(a0), 30)], color='peru', thickness=2)
llim = line([(-2, a0), (5, a0)], color='red', thickness=1.5, 
            linestyle=':', legend_label=r'$|\ell| = a$') \
       + line([(-2, -a0), (5, -a0)], color='red', thickness=1.5, 
              linestyle=':')
graph = gq + gl + gqmin + hor + llim
for r1 in [rph_ss(a0), rph_s(a0), rph_p(a0), rph_m(a0)]:
    graph += point((r1, qsph(a0, r1)), size=60, color='blue', zorder=100)
graph.set_legend_options(handlelength=2)
show(graph, xmin=-1.5, xmax=4.5, ymin=-7, ymax=28)

gq_stable = plot(qsph(a0, r), (r, 0, rph_ms(a0)), thickness=5) 
gl_stable = plot(lsph(a0, r), (r, 0, rph_ms(a0)), thickness=5, color='red') 
graph = gq + gq_stable + gl + gl_stable + gqmin + hor + llim
show(graph, xmin=-1, xmax=1.5, ymin=-0.2, ymax=0.6)

show(graph, xmin=-1, xmax=1.5, ymin=-0.2, ymax=0.02)

show(graph, xmin=-1, xmax=1.5, ymin=-0.01, ymax=0.01)

Plots for $a=0.998\, m$

a0 = 0.998

gq = plot(qsph(a0, r), (r, -2, 0.999), thickness=2, legend_label=r'$q$', 
          axes_labels=[r'$r_0/m$', r'$\ell/m,\ q/m^2$'],
          frame=True, gridlines=True) \
     + plot(qsph(a0, r), (r, 1.001, 5), thickness=2)
gl = plot(lsph(a0, r), (r, -2, 0.999), color='red', thickness=2, legend_label=r'$\ell$') \
     + plot(lsph(a0, r), (r, 1.001, 5), color='red', thickness=2)
gqmin = plot(lambda r: qmin(a0, r), (r, -2, 5), color='grey', thickness=2, fill=-8)
hor = line([(rp(a0), -10), (rp(a0), 30)], color='black', thickness=2) \
      + line([(rm(a0), -10), (rm(a0), 30)], color='peru', thickness=2)
llim = line([(-2, a0), (5, a0)], color='red', thickness=1.5, 
            linestyle=':', legend_label=r'$|\ell| = a$') \
       + line([(-2, -a0), (5, -a0)], color='red', thickness=1.5, 
              linestyle=':')
graph = gq + gl + gqmin + hor + llim
for r1 in [rph_ss(a0), rph_s(a0), rph_p(a0), rph_m(a0)]:
    graph += point((r1, qsph(a0, r1)), size=60, color='blue', zorder=100)
graph.set_legend_options(handlelength=2)
show(graph, xmin=-1.5, xmax=4.5, ymin=-7, ymax=28)

gq_stable = plot(qsph(a0, r), (r, 0, rph_ms(a0)), thickness=5) 
gl_stable = plot(lsph(a0, r), (r, 0, rph_ms(a0)), thickness=5, color='red') 
graph = gq + gq_stable + gl + gl_stable + gqmin + hor + llim
show(graph, xmin=-1, xmax=1.5, ymin=-1.5, ymax=2.5)

Solutions of the cubic equation $2 r^3 - 3 m r^2 + a^2 m = 0$

P = 2*r^3 - 3*r^2 + a^2 
P

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

a0 = 0.95
plot(P.subs({a: a0}), (r, -0.7, 1.5))

Pdep = (P/2).subs({r: x + 1/2}).simplify_full()
Pdep

$$\renewcommand{\Bold}[1]{\mathbf{#1}}x^{3} + \frac{1}{2} \, a^{2} - \frac{3}{4} \, x - \frac{1}{4}$$

pp = -3/4
qq = 1/2*(a^2 - 1/2)

sqrt(-pp/3)

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\frac{1}{2}$$

3*qq/(2*pp)

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

3*qq/(2*pp)*sqrt(-3/pp)

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

g = Graphics()
for k in range(3):
    rk = 1/2 + cos(2/3*asin(a) + 2*k*pi/3)
    g += plot(rk, (a, 0, 1), color=hue((3-k)/3), thickness=2)
g

Checking that $x_1 = \cos\left(\frac{2}{3}\arcsin(a) + \frac{2\pi}{3}\right)$ is a solution:

s = Pdep.subs({x: cos(2/3*asin(a) + 2*pi/3)}) 
s

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

s.expand_trig().simplify_trig().reduce_trig().expand_trig()

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

Plot of the critical radii as functions of $a$

graph = plot(rph_p, (0, 1), axes_labels=[r'$a/m$', r'$r_0/m$'], 
             color='green', thickness=1.5, legend_label=r'$r_{\rm ph}^+$', 
             fill=rph_m, fillcolor='palegreen',
             frame=True, gridlines=True, axes=True) \
+ plot(rph_m, (0, 1), linestyle='--', color='green', thickness=1.5, 
       legend_label=r'$r_{\rm ph}^-$') \
+ plot(rph_s, (0, 1), color='green', linestyle=':', thickness=2, 
       legend_label=r'$r_{\rm ph}^*$', fill=rph_ss, fillcolor='palegreen') \
+ plot(rph_ss, (0, 1), color='green', linestyle='-.', thickness=1.5, 
       legend_label=r'$r_{\rm ph}^{**}$') \
+ plot(rph_ms, (a, 0, 1), color='blue', linestyle='-', thickness=1.5, 
       legend_label=r'$r_{\rm ph}^{\rm ms}$', fill=0, fillcolor='darkseagreen') \
+ plot(rp, (0, 1), color='black', thickness=1.5, legend_label=r'$r_+$') \
+ plot(rm, (0, 1), color='darkgoldenrod', thickness=2, legend_label=r'$r_-$') 
graph.set_legend_options(handlelength=2, loc=(1.02, 0.36))
show(graph)
graph.save('gik_spher_orb_range.pdf')

Properties of the spherical orbit at $r=r_{\rm ph}^{**}$

graph = plot(lambda a: lsph(a, rph_ss(a)), (0, 1), color='red', thickness=2, 
             legend_label=r'$\ell_{\rm ph}^{**}$',
             axes_labels=[r'$a/m$', r'$\ell_{\rm ph}^{**}/m,\ \  q_{\rm ph}^{**}/m^2$'],
             frame=True, gridlines=True, axes=False) \
        + plot(lambda a: qsph(a, rph_ss(a)), (0.001, 0.999), color='blue', thickness=2,
               legend_label=r'$q_{\rm ph}^{**}$')
graph.set_legend_options(handlelength=2)
graph.save("gik_ell_q_rss.pdf")
graph

lsph(1, rph_ss(1))

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

qsph(1, rph_ss(1))

$$\renewcommand{\Bold}[1]{\mathbf{#1}}-\frac{9}{16}$$

n(_)

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

th_s = lambda a: arcsin(sqrt(abs(lsph(a, rph_ss(a)))/a))

graph = plot(th_s, (0.0001, 0.9999), color='blue', thickness=2, 
             axes_labels=[r'$a/m$', r'$\theta^{**}_{\rm ph}$'],
             frame=True, gridlines=True, axes=False)
graph.save("gik_theta_ss.pdf")
graph

n(pi/6)

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

Angular momentum of the circular photon orbits in the equatorial plane

graph = plot(lambda a: lsph(a, rph_s(a)), (0, 1), color='red', thickness=2, 
             linestyle=':', legend_label=r'$\ell_{\rm ph}^{*}$',
             axes_labels=[r'$a/m$', r'$\ell/m$'],
             frame=True, gridlines=True, axes=False) \
        + plot(lambda a: lsph(a, rph_p(a)), (0, 1), color='red', thickness=2, 
               legend_label=r'$\ell_{\rm ph}^+$') \
        + plot(lambda a: lsph(a, rph_m(a)), (0, 1), color='red', thickness=2, 
               linestyle='--', legend_label=r'$\ell_{\rm ph}^-$')
graph.set_legend_options(handlelength=2)
graph.save("gik_ell_circ_equat.pdf", ymin=-8, ymax=6)
show(graph, ymin=-8, ymax=6)

Stability of the photon spherical orbits

Generic function $\mathcal{R}(r)$:

l = var('l', latex_name=r'\ell')
q = var('q')
R(a, l, q, r) = r^4 + (a^2 - l^2 - q)*r^2 + 2*(q + (l - a)^2)*r - a^2*q
R

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

Function $\mathcal{R}(r)$ for a spherical orbit:

r0 = var('r_0')

R0(a, r0, r) = R(a, lsph(a, r0), qsph(a, r0), r).simplify_full()
R0(a, r0, r)

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\frac{r_{0}^{6} - {\left(2 \, r^{2} - 8 \, r - 9\right)} r_{0}^{4} - 6 \, r_{0}^{5} + r^{4} - 4 \, {\left(a^{2} + 4 \, r\right)} r_{0}^{3} + {\left(r^{4} + 8 \, a^{2} r + 6 \, r^{2}\right)} r_{0}^{2} - 2 \, {\left(2 \, a^{2} r^{2} + r^{4}\right)} r_{0}}{r_{0}^{2} - 2 \, r_{0} + 1}$$

We check that $r_0$ is a double root of $\mathcal{R}(r)$:

R0(a, r0, r).factor()

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

s = (R0(a, r0, r)*(r0 - 1)^2/(r - r0)^2).simplify_full()
s

$$\renewcommand{\Bold}[1]{\mathbf{#1}}2 \, {\left(r - 3\right)} r_{0}^{3} + r_{0}^{4} + {\left(r^{2} - 4 \, r + 9\right)} r_{0}^{2} + r^{2} - 2 \, {\left(2 \, a^{2} + r^{2} - r\right)} r_{0}$$

bool(s*(r - r0)^2 == R0(a, r0, r).factor().numerator())

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\mathrm{True}$$

s.collect(r)

$$\renewcommand{\Bold}[1]{\mathbf{#1}}r_{0}^{4} + {\left(r_{0}^{2} - 2 \, r_{0} + 1\right)} r^{2} - 4 \, a^{2} r_{0} - 6 \, r_{0}^{3} + 2 \, {\left(r_{0}^{3} - 2 \, r_{0}^{2} + r_{0}\right)} r + 9 \, r_{0}^{2}$$

b = r0^3 - 6*r0^2 + 9*r0 - 4*a^2
b.factor()

$$\renewcommand{\Bold}[1]{\mathbf{#1}}r_{0}^{3} - 4 \, a^{2} - 6 \, r_{0}^{2} + 9 \, r_{0}$$

bool(R0(a, r0, r) == (r - r0)^2*(r^2 + 2*r0*r + r0*b/(r0 - 1)^2))

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\mathrm{True}$$

Check of the formula $$\mathcal{R}(r) = (r - r_0)^2 \left( r^2 + 2 r_0 r - a^2 \frac{q}{r_0^2} \right)$$

bool(R0(a, r0, r) == (r - r0)^2*(r^2 + 2*r0*r - a^2/r0^2*qsph(a, r0)))

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\mathrm{True}$$

Plot of $\mathcal{R}(r)$ for various values of $r_0$

Rss(a, r) = R0(a, rph_ss(a), r)
Rs(a, r) = R0(a, rph_s(a), r)
Rp(a, r) = R0(a, rph_p(a), r)
Rm(a, r) = R0(a, rph_m(a), r)
Rms(a, r) = R0(a, rph_ms(a), r)

a0 = 0.95
# show(LatexExpr(r'r_{\rm ph}^{**} = ') + latex(n(rph_ss(a0))) + L)
show(LatexExpr(r'r_{\rm ph}^{**} = ' + '{}  m'.format(n(rph_ss(a0)))))
plot(lambda r: Rss(a0, r), (-1, 1), axes_labels=[r'$r/m$', r'$\mathcal{R}(r)$'])

$$\renewcommand{\Bold}[1]{\mathbf{#1}}r_{\rm ph}^{**} = -0.477673658836338 m$$

A stablle inner orbit:

r0 = 0.4
show(LatexExpr(r'r_0 = ' + '{}  m'.format(float(r0))))
plot(lambda r: R0(a0, r0, r), (-0.2, 1), axes_labels=[r'$r/m$', r'$\mathcal{R}(r)$'])

$$\renewcommand{\Bold}[1]{\mathbf{#1}}r_0 = 0.4 m$$

The marginally stable orbit:

show(LatexExpr(r'r_{\rm ph}^{\rm ms} = ' + '{}  m'.format(n(rph_ms(a0)))))
plot(lambda r: Rms(a0, r), (-1, 1), axes_labels=[r'$r/m$', r'$\mathcal{R}(r)$'])

$$\renewcommand{\Bold}[1]{\mathbf{#1}}r_{\rm ph}^{\rm ms} = 0.539741795874205 m$$

show(LatexExpr(r'r_{\rm ph}^{*} = ' + '{}  m'.format(n(rph_s(a0)))))
plot(lambda r: Rs(a0, r), (-1, 1), axes_labels=[r'$r/m$', r'$\mathcal{R}(r)$'])

$$\renewcommand{\Bold}[1]{\mathbf{#1}}r_{\rm ph}^{*} = 0.658372153864346 m$$

show(LatexExpr(r'r_{\rm ph}^{+} = ' + '{}  m'.format(n(rph_p(a0)))))
plot(lambda r: Rp(a0, r), (-1/2, 2), axes_labels=[r'$r/m$', r'$\mathcal{R}(r)$'])

$$\renewcommand{\Bold}[1]{\mathbf{#1}}r_{\rm ph}^{+} = 1.38628052846298 m$$

show(LatexExpr(r'r_{\rm ph}^{-} = ' + '{}  m'.format(n(rph_m(a0)))))
plot(lambda r: Rm(a0, r), (-1/2, 5), axes_labels=[r'$r/m$', r'$\mathcal{R}(r)$'])

$$\renewcommand{\Bold}[1]{\mathbf{#1}}r_{\rm ph}^{-} = 3.95534731767268 m$$

R2(a, l, q, r) = diff(R(a, l, q, r), r, 2).simplify_full()
R2

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

R2o(a, r) = R2(a, lsph(a, r), qsph(a, r), r).simplify_full().factor()
R2o

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

Computation of $\frac{\mathrm{d}q}{\mathrm{d}r_0}$

dqdr(a, r) = diff(qsph(a, r), r).simplify_full().factor()
dqdr(a, r)

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

Let us check that marginally stable orbits correspond to an extremum (actually a maximum) of $q(r_0)$:

dqdr(a, rph_ms(a)).simplify_full()

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

The maximum at $r_0 = 3 m$ does not depend on $a$:

qsph(a, 3)

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

while for $r_0 = r_{\rm ph}^{ms}$, we have

qms = qsph(a, rph_ms(a)).simplify_full()
qms

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

q1 = SR(qms._sympy_().simplify())
q1

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

(qms - q1).simplify_full()

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

qms = q1

taylor(qms, a, 0, 10)

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\frac{38}{729} \, a^{10} + \frac{4}{81} \, a^{8} + \frac{1}{27} \, a^{6}$$

qms.subs({a: 1})

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

plot(qms, (a, 0, 1), axes_labels=[r'$a/m$', r'$q_{\rm ms}/m^2$'],
     thickness=2, frame=True, gridlines=True)

dldr(a, r) = diff(lsph(a, r), r).simplify_full().factor()
dldr(a, r)

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

dldr(a, rph_ms(a)).simplify_full()

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

lms = lsph(a, rph_ms(a)).simplify_full()
lms

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

taylor(lms, a, 0, 6)

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\frac{4}{27} \, a^{5} + \frac{1}{3} \, a^{3} + a$$

lms.limit(a=1)

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

plot(lms, (a, 0, 1), color='red', thickness=2,
     axes_labels=[r'$a/m$', r'$\ell_{\rm ms}/m$'],
     frame=True, gridlines=True) \
+ plot(a, (a, 0, 1), color='black', linestyle='--')

Plot of spherical orbits in the meriodinal plane

Choice of a radial function $\hat r$

R2(r) = (r + sqrt(r^2 + 4))/2
R2

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

R4(r) = (r + (r^4 + 16)^(1/4))/2
R4

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

n(R2(2)), n(R4(2))

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

n(R2(6)), n(R4(6))

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

g = plot(R4, (-8, 8), color='green',
         axes_labels=[r'$r/m$', r'$\hat{r}/m$'])
show(g, aspect_ratio=1, frame=True, gridlines=True)

g = plot(R2, (-8, 8), axes_labels=[r'$r/m$', r'$\hat{r}/m$']) \
    + plot(r, (r, 0, 8), linestyle='--') \
    + plot(exp(r), (r, -8, 2.1), color='red') \
    + plot(R4, (-8, 8), color='green')
show(g, aspect_ratio=1, frame=True, gridlines=True)

$\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, {\ell}, q \right) \ {\mapsto} \ \arccos\left(\sqrt{\frac{1}{2} \, \sqrt{{\left(\frac{{\ell}^{2} + q}{a^{2}} - 1\right)}^{2} + \frac{4 \, q}{a^{2}}} - \frac{{\ell}^{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, {\ell}, q \right) \ {\mapsto} \ \arccos\left(\sqrt{-\frac{1}{2} \, \sqrt{{\left(\frac{{\ell}^{2} + q}{a^{2}} - 1\right)}^{2} + \frac{4 \, q}{a^{2}}} - \frac{{\ell}^{2} + q}{2 \, a^{2}} + \frac{1}{2}}\right)$$

Plots

a0 = 0.95

R = R2
th = var('th')
# We use parametric_plot() instead of arc() because of legend_label
g = parametric_plot((R(0)*sin(th), R(0)*cos(th)), 
                    (th, 0, pi), color='darkorange', thickness=3,
                    linestyle=':', legend_label=r'$r=0$',
                    axes_labels=[r'$\frac{\hat{r}}{m}\sin\theta$', 
                                 r'$\frac{\hat{r}}{m}\cos\theta$'])
g += parametric_plot((R(rm(a0))*sin(th), R(rm(a0))*cos(th)), 
                     (th, 0, pi), color='peru', thickness=3,
                     legend_label=r'$r=r_-$')
g += parametric_plot((R(rp(a0))*sin(th), R(rp(a0))*cos(th)), 
                     (th, 0, pi), color='black', thickness=3,
                     legend_label=r'$r=r_+$')

sing = point((R(0),0), size=60, color='orangered', zorder=100) # ring singularity
rminf = point((0,0), size=60, marker='o', color='white', 
              markeredgecolor='black', zorder=100)

g += sing + rminf
g.set_legend_options(handlelength=2.5, loc=(0.9, 0.9), 
                     font_size='large')

ns = 9
rmin, rmax = RR(1.00001*rph_p(a0)), RR(0.99999*rph_m(a0))

g += parametric_plot([R(r)*sin(theta0(a0, lsph(a0, r), qsph(a0, r))),
                      R(r)*cos(theta0(a0, lsph(a0, r), qsph(a0, r)))],
                     (r, rmin, rmax), color='darkgreen', thickness=2,
                     plot_points=200, fill=True, fillcolor='palegreen', 
                     fillalpha=0.2)
g += parametric_plot([R(r)*sin(theta0(a0, lsph(a0, r), qsph(a0, r))),
                      -R(r)*cos(theta0(a0, lsph(a0, r), qsph(a0, r)))],
                     (r, rmin, rmax), color='darkgreen', thickness=2,
                     plot_points=200, fill=True, fillcolor='palegreen',
                     fillalpha=0.2)
dr = (rmax - rmin)/(ns - 1)
for i in range(ns):
    r0 = rmin + i*dr
    l0 = lsph(a0, r0)
    q0 = qsph(a0, r0)
    th0 = theta0(a0, l0, q0)
    g += arc((0, 0), R(r0), sector=(th0 - pi/2, pi/2 - th0), 
             color='green')

rmin, rmax = RR(0), RR(rph_s(a0))
g += parametric_plot([R(r)*sin(theta0(a0, lsph(a0, r), qsph(a0, r))),
                      R(r)*cos(theta0(a0, lsph(a0, r), qsph(a0, r)))],
                     (r, rmin, rmax), color='darkgreen', thickness=2,
                     fill=True, fillcolor='palegreen', fillalpha=0.2)
g += parametric_plot([R(r)*sin(theta0(a0, lsph(a0, r), qsph(a0, r))),
                      -R(r)*cos(theta0(a0, lsph(a0, r), qsph(a0, r)))],
                     (r, rmin, rmax), color='darkgreen', thickness=2,
                     fill=True, fillcolor='palegreen', fillalpha=0.2)

# outer polar orbits:
r0 = rph_pol(a0)
print(lsph(a0, r0))
q0 = qsph(a0, r0)
th0 = 0
print(theta0(a0, 0, q0))
g += arc((0, 0), R(r0), sector=(th0 - pi/2, pi/2 - th0), 
         linestyle='--', color='green')
    
ns = 3
rmin, rmax = RR(rph_ms(a0)), 0.9999*RR(rph_s(a0))
dr = (rmax - rmin)/(ns - 1)
for i in range(ns):
    r0 = rmin + i*dr
    l0 = lsph(a0, r0)
    q0 = qsph(a0, r0)
    th0 = theta0(a0, l0, q0)
    g += arc((0, 0), R(r0), sector=(th0 - pi/2, pi/2 - th0), 
             color='green')

ns = 5
rmin, rmax = RR(0), RR(rph_ms(a0))
dr = (rmax - rmin)/(ns - 1)
for i in range(ns):
    r0 = rmin + i*dr
    l0 = lsph(a0, r0)
    q0 = qsph(a0, r0)
    th0 = theta0(a0, l0, q0)
    g += arc((0, 0), R(r0), sector=(th0 - pi/2, pi/2 - th0), 
             color='green')

r0 = rph_ms(a0)
l0 = lsph(a0, r0)
q0 = qsph(a0, r0)
th0 = theta0(a0, l0, q0)
g += arc((0, 0), R(r0), sector=(th0 - pi/2, pi/2 - th0), 
         thickness=2, color='blue', zorder=100)

    
rmin, rmax = RR(rph_ss(a0)), -1e-8

g += parametric_plot([R(r)*sin(theta0(a0, lsph(a0, r), qsph(a0, r))),
                      R(r)*cos(theta0(a0, lsph(a0, r), qsph(a0, r)))],
                     (r, rmin, rmax), color='darkgreen', thickness=2)
g += parametric_plot([R(r)*sin(theta1(a0, lsph(a0, r), qsph(a0, r))),
                      R(r)*cos(theta1(a0, lsph(a0, r), qsph(a0, r)))],
                     (r, rmin, rmax), color='darkgreen', thickness=2)
g += parametric_plot([R(r)*sin(theta0(a0, lsph(a0, r), qsph(a0, r))),
                      -R(r)*cos(theta0(a0, lsph(a0, r), qsph(a0, r)))],
                     (r, rmin, rmax), color='darkgreen', thickness=2)
g += parametric_plot([R(r)*sin(theta1(a0, lsph(a0, r), qsph(a0, r))),
                      -R(r)*cos(theta1(a0, lsph(a0, r), qsph(a0, r)))],
                     (r, rmin, rmax), color='darkgreen', thickness=2)
    
ns = 5
dr = (rmax - rmin)/(ns - 1)
for i in range(ns):
    r0 = rmin + i*dr
    l0 = lsph(a0, r0)
    q0 = qsph(a0, r0)
    th0 = theta0(a0, l0, q0)
    th1 = theta1(a0, l0, q0)
    g += arc((0, 0), R(r0), sector=(pi/2 - th1, pi/2 - th0), color='green')
    g += arc((0, 0), R(r0), sector=(th0 - pi/2, th1 - pi/2), color='green')
g += arc((0,0), R(rph_ss(a0)), sector=(- pi/2, pi/2), color='grey', 
         linestyle=':', thickness=1)

# inner polar orbits:
r0 = rph_pol_in(a0)
print(n(lsph(a0, r0)))
q0 = qsph(a0, r0)
th0 = 0
th1 = theta1(a0, 0, q0)
print(n(theta0(a0, 0, q0)))
g += arc((0, 0), R(r0), sector=(pi/2 - th1, pi/2 - th0), 
         linestyle='--', color='green')
g += arc((0, 0), R(r0), sector=(th0 - pi/2, th1 - pi/2), 
         linestyle='--', color='green')

g

-6.26333640524599e-16
0.000000000000000
6.68118973136174e-16
0.000000000000000

Adding the ergoregion:

x, y = var('x y')
Rxy = sqrt(x^2 + y^2)
r1 = Rxy - 1/Rxy
costh2 = y^2/(x^2 + y^2)
g += region_plot(r1^2 - 2*r1 + a0^2*costh2 < 0, (0, 3), (-3, 3), 
                 incol='whitesmoke', bordercol='grey')

Coloring the region of vortical spherical orbits:

def vortical_orb(x, y):
    R0 = sqrt(x*x + y*y)
    r0 = R0 - 1/R0
    if r0 < rph_ss(a0) or r0 > 0:
        return False
    phi = atan2(y, x)
    th = pi/2 - phi
    l0 = lsph(a0, r0)
    q0 = qsph(a0, r0)
    th0 = theta0(a0, l0, q0)
    th1 = theta1(a0, l0, q0)
    if (th >= th0 and th <= th1) or (th >= pi - th1 and th <= pi - th0):
        return True
    return False

g += region_plot(vortical_orb, (0, 1), (-1, 1), incol='palegreen', alpha=0.2)

Coloring the region of stable spherical orbits:

def stable_orb(x, y):
    R0 = sqrt(x*x + y*y)
    r0 = R0 - 1/R0
    if r0 <= 0 or r0 > rph_ms(a0):
        return False
    phi = atan2(y, x)
    th = pi/2 - phi
    l0 = lsph(a0, r0)
    q0 = qsph(a0, r0)
    th0 = theta0(a0, l0, q0)
    if th < th0 or th > pi - th0:
        return False
    return True

g += region_plot(stable_orb, (1, 1.5), (-0.7, 0.7), incol='green', alpha=0.5)

Adding circular orbits:

for r1 in [rph_s(a0), rph_p(a0), rph_m(a0)]:
    g += point((R(r1), 0), size=60, color='green', zorder=100)
r1 = rph_ss(a0)
sin_th_ss = sqrt( r1/(r1 - 1)*(1 - r1^2/a0^2) )
cos_th_ss = sqrt( 1 - sin_th_ss^2 ) 
g += point((R(r1)*sin_th_ss, R(r1)*cos_th_ss), size=60, color='green', zorder=100)
g += point((R(r1)*sin_th_ss, -R(r1)*cos_th_ss), size=60, color='green', zorder=100)
show(g, xmin=0, xmax=4.5, ticks=[1, 1], figsize=8)
g.save("gik_spo_meridional.pdf", xmin=0, xmax=4.5, ticks=[1, 1], figsize=8)

show(g, xmin=0, xmax=2.2, ymin=-1, ymax=1, show_legend=False)
g.save("gik_spo_meridional_zoom.pdf", xmin=0, xmax=2.2, ymin=-1, ymax=1, 
       show_legend=False)