Existence of spherical photon orbits in Kerr spacetime

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

In [1]:
version()
Out[1]:
'SageMath version 9.3.beta5, Release Date: 2020-12-27'
In [2]:
%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$.

In [3]:
a, r = var('a r') 
In [4]:
lsph(a, r) = (r^2*(3 - r) - a^2*(r + 1))/(a*(r -1))
lsph
Out[4]:
In [5]:
qsph(a, r) = r^3 / (a^2*(r - 1)^2) * (4*a^2 - r*(r - 3)^2)
qsph
Out[5]:

The radii of the two horizons:

In [6]:
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}$

In [7]:
rph_ss(a) = 1/2 + cos(2/3*asin(a) + 2*pi/3)
rph_ss
Out[7]:
In [8]:
rph_s(a) = 4*cos(acos(-a)/3 + 4*pi/3)^2
rph_s
Out[8]:
In [9]:
rph_p(a) = 4*cos(acos(-a)/3)^2
rph_p
Out[9]:
In [10]:
rph_m(a) = 4*cos(acos(a)/3)^2
rph_m
Out[10]:

We add the radius of the marginally stable orbit:

In [11]:
rph_ms(a) = 1 - (1 - a^2)^(1/3)
rph_ms
Out[11]:

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

In [12]:
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
Out[12]:
In [13]:
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
Out[13]:

Plots for $a=0.95 \, m$

In [14]:
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)))
In [15]:
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)
In [16]:
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)
In [17]:
show(gq+gl, ymin=-5, ymax=10)
In [18]:
def qmin(a1, r):
    l0 = abs(lsph(a1, r))
    if l0 < a1:
        return -(a1 - l0)^2
    return 0
In [19]:
gqmin = plot(lambda r: qmin(a0, r), (r, -2, 5), color='grey', thickness=2, fill=-8)
In [20]:
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=':')
In [21]:
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:

In [22]:
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$

In [23]:
a0 = 0.5
In [24]:
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)
In [25]:
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)
In [26]:
show(graph, xmin=-1, xmax=1.5, ymin=-0.2, ymax=0.02)
In [27]:
show(graph, xmin=-1, xmax=1.5, ymin=-0.01, ymax=0.01)

Plots for $a=0.998\, m$

In [28]:
a0 = 0.998
In [29]:
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)
In [30]:
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$

In [31]:
P = 2*r^3 - 3*r^2 + a^2 
P
Out[31]:
In [32]:
a0 = 0.95
plot(P.subs({a: a0}), (r, -0.7, 1.5))
Out[32]:
In [33]:
Pdep = (P/2).subs({r: x + 1/2}).simplify_full()
Pdep
Out[33]:
In [34]:
pp = -3/4
qq = 1/2*(a^2 - 1/2)
In [35]:
sqrt(-pp/3)
Out[35]:
In [36]:
3*qq/(2*pp)
Out[36]:
In [37]:
3*qq/(2*pp)*sqrt(-3/pp)
Out[37]:
In [38]:
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
Out[38]:

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

In [39]:
s = Pdep.subs({x: cos(2/3*asin(a) + 2*pi/3)}) 
s
Out[39]:
In [40]:
s.expand_trig().simplify_trig().reduce_trig().expand_trig()
Out[40]:

Plot of the critical radii as functions of $a$

In [41]:
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}^{**}$

In [42]:
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
Out[42]:
In [43]:
lsph(1, rph_ss(1))
Out[43]:
In [44]:
qsph(1, rph_ss(1))
Out[44]:
In [45]:
n(_)
Out[45]:
In [46]:
th_s = lambda a: arcsin(sqrt(abs(lsph(a, rph_ss(a)))/a))
In [47]:
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
Out[47]:
In [48]:
n(pi/6)
Out[48]:

Angular momentum of the circular photon orbits in the equatorial plane

In [49]:
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)$:

In [50]:
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
Out[50]:

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

In [51]:
r0 = var('r_0')
In [52]:
R0(a, r0, r) = R(a, lsph(a, r0), qsph(a, r0), r).simplify_full()
R0(a, r0, r)
Out[52]:

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

In [53]:
R0(a, r0, r).factor()
Out[53]:
In [54]:
s = (R0(a, r0, r)*(r0 - 1)^2/(r - r0)^2).simplify_full()
s
Out[54]:
In [55]:
bool(s*(r - r0)^2 == R0(a, r0, r).factor().numerator())
Out[55]:
In [56]:
s.collect(r)
Out[56]:
In [57]:
b = r0^3 - 6*r0^2 + 9*r0 - 4*a^2
b.factor()
Out[57]:
In [58]:
bool(R0(a, r0, r) == (r - r0)^2*(r^2 + 2*r0*r + r0*b/(r0 - 1)^2))
Out[58]:

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

In [59]:
bool(R0(a, r0, r) == (r - r0)^2*(r^2 + 2*r0*r - a^2/r0^2*qsph(a, r0)))
Out[59]:

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

In [60]:
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)
In [61]:
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)$'])
Out[61]:

A stablle inner orbit:

In [62]:
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)$'])
Out[62]:

The marginally stable orbit:

In [63]:
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)$'])
Out[63]:
In [64]:
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)$'])
Out[64]:
In [65]:
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)$'])
Out[65]:
In [66]:
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)$'])
Out[66]:
In [67]:
R2(a, l, q, r) = diff(R(a, l, q, r), r, 2).simplify_full()
R2
Out[67]:
In [68]:
R2o(a, r) = R2(a, lsph(a, r), qsph(a, r), r).simplify_full().factor()
R2o
Out[68]:

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

In [69]:
dqdr(a, r) = diff(qsph(a, r), r).simplify_full().factor()
dqdr(a, r)
Out[69]:

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

In [70]:
dqdr(a, rph_ms(a)).simplify_full()
Out[70]:

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

In [71]:
qsph(a, 3)
Out[71]:

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

In [72]:
qms = qsph(a, rph_ms(a)).simplify_full()
qms
Out[72]:
In [73]:
q1 = SR(qms._sympy_().simplify())
q1
Out[73]:
In [74]:
(qms - q1).simplify_full()
Out[74]:
In [75]:
qms = q1
In [76]:
taylor(qms, a, 0, 10)
Out[76]:
In [77]:
qms.subs({a: 1})
Out[77]:
In [78]:
plot(qms, (a, 0, 1), axes_labels=[r'$a/m$', r'$q_{\rm ms}/m^2$'],
     thickness=2, frame=True, gridlines=True)
Out[78]:
In [79]:
dldr(a, r) = diff(lsph(a, r), r).simplify_full().factor()
dldr(a, r)
Out[79]:
In [80]:
dldr(a, rph_ms(a)).simplify_full()
Out[80]:
In [81]:
lms = lsph(