Circular equatorial 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.2, Release Date: 2020-10-24'
In [2]:
%display latex

Effective potential $\mathcal{V}(r)$

In [3]:
a = var('a')
r = var('r')
eps = var('eps', latex_name=r'\varepsilon')
l = var('l', latex_name=r'\ell') 
In [4]:
V(a, eps, l, r) = (1 - eps^2 - 2/r + (l^2 + a^2*(1 - eps^2))/r^2
                    - 2*(l - a*eps)^2/r^3)
V
Out[4]:

Constraint $r^2 - 3 m r \pm 2 a \sqrt{m r} > 0$

In [5]:
f(r, a) = r*(r-3) + 2*a*sqrt(r)
f
Out[5]:
In [6]:
g1 = plot(f(r, 0), (r, 0, 4), legend_label='a = 0', axes_labels=[r'$r/m$', ''])
g2 = plot(f(r, 0.5), (r, 0, 4), color='pink', legend_label='a = 0.5')
g3 = plot(f(r, 0.8), (r, 0, 4), color='red', legend_label='a = 0.8')
g4 = plot(f(r, 0.98), (r, 0, 4), color='brown', legend_label='a = 0.98')
g1 + g2 + g3 + g4
Out[6]:
In [7]:
f(r, a) = r*(r-3) - 2*a*sqrt(r)
f
Out[7]:
In [8]:
g1 = plot(f(r, 0), (r, 0, 6), legend_label='a = 0', axes_labels=[r'$r/m$', ''])
g2 = plot(f(r, 0.5), (r, 0, 6), color='pink', legend_label='a = 0.5')
g3 = plot(f(r, 0.8), (r, 0, 6), color='red', legend_label='a = 0.8')
g4 = plot(f(r, 0.98), (r, 0, 6), color='brown', legend_label='a = 0.98')
g1 + g2 + g3 + g4
Out[8]:

Existence of circular orbits: boundaries $r_{\rm ph}^\pm$ and $r_{\rm ph}^*$

In [9]:
r_min_p(a) = 4*cos(acos(-a)/3)^2
r_min_p
Out[9]:
In [10]:
r_min_m(a) = 4*cos(acos(a)/3)^2
r_min_m
Out[10]:
In [11]:
rp(a) = 1 + sqrt(1 - a^2)
rm(a) = 1 - sqrt(1 - a^2)
In [12]:
plot(r_min_p(a) - rp(a), (a, 0, 1), 
     axes_labels=[r'$a/m$', r'$(r_{\rm ph}^+ - r_+)/m$']) 
Out[12]:
In [13]:
rs(a) = 4*cos(acos(-a)/3 + 4*pi/3)^2
rs
Out[13]:
In [14]:
plot(rs, (0, 1), axes_labels=[r'$a/m$', r'$r_{\rm ph}^*/m$']) 
Out[14]:

Discrepancy between $r_{\rm ph}^*$ and $r_-$:

In [15]:
plot(rm(a) - rs(a), (a, 0, 1), 
     axes_labels=[r'$a/m$', r'$(r_- - r_{\rm ph}^*)/m$'], 
     frame=True, gridlines=True, axes=False) 
Out[15]:
In [16]:
graph = plot(r_min_p, (0, 1), axes_labels=[r'$a/m$', r'$r_0/m$'], 
             color='green', thickness=1.5, legend_label=r'$r_{\rm ph}^+$', 
             fill=4.1, fillcolor='lightgrey',
             frame=True, gridlines=True, axes=False) \
+ plot(r_min_m, (0, 1), linestyle='--', color='green', thickness=1.5, 
       legend_label=r'$r_{\rm ph}^-$', fill=4.1, fillcolor='grey') \
+ plot(rs, (0, 1), color='green', linestyle=':', thickness=2, 
       legend_label=r'$r_{\rm ph}^*$', fill=0, fillcolor='lightgrey') \
+ plot(lambda x: 2, (0, 1), color='orange', thickness=1.5, legend_label=r'$r_{\rm ergo}$') \
+ plot(rp, (0, 1), color='black', thickness=1.5, legend_label=r'$r_+$') \
+ plot(rm, (0, 1), color='peru', thickness=2, legend_label=r'$r_-$') 
graph.set_legend_options(handlelength=2, loc=(1.02, 0.3))
graph
Out[16]:
In [17]:
graph.save('gek_circ_orb_lim.pdf')

ISCO

In [18]:
Z1 = 1 + (1 - a^2)^(1/3)*((1 + a)^(1/3) + (1 - a)^(1/3))
Z1
Out[18]:
In [19]:
g = plot(Z1, (a, 0, 1), axes_labels=[r"$a/m$", r"$Z_1$"],
         thickness=1.5, frame=True, gridlines=True, axes=True)
g
Out[19]:
In [20]:
g.save("gek_Z1.pdf")
In [21]:
Z2 = sqrt(Z1^2 + 3*a^2)
Z2
Out[21]:
In [22]:
r_isco_p(a) = 3 + Z2 - sqrt((3 - Z1)*(3 + Z1 + 2*Z2))
r_isco_m(a) = 3 + Z2 + sqrt((3 - Z1)*(3 + Z1 + 2*Z2))

Energy and angular momentum as functions of $r_0$

In [23]:
sAp = sqrt(r^2 - 3*r + 2*a*sqrt(r))
sAm = sqrt(r^2 - 3*r - 2*a*sqrt(r))
eps_p(a, r) = (r^2 - 2*r + a*sqrt(r))/(r*sAp)
eps_m(a, r) = (r^2 - 2*r - a*sqrt(r))/(r*sAm)
In [24]:
eps_p(a, r)
Out[24]:
In [25]:
eps_m(a, r)
Out[25]:
In [26]:
l_p(a, r) = (r^2 + a^2 - 2*a*sqrt(r))/(sqrt(r)*sAp)
l_m(a, r) = -(r^2 + a^2 + 2*a*sqrt(r))/(sqrt(r)*sAm)
In [27]:
l_p(a, r)
Out[27]:
In [28]:
l_m(a, r)
Out[28]:
In [29]:
num_l_p(a, r) = r^2 + a^2 - 2*a*sqrt(r)
In [30]:
al = [0., 0.5, 0.9, 0.98, 1]
ah = -2
g = Graphics()
for a1 in al[1:]:
    g += plot(eps_p(a1, r), (r, 0.01, 0.999*rs(a1)), color=hue(a1/ah), 
              linestyle='-', thickness=1.5,
              axes_labels=[r"$r_0/m$", r"$\varepsilon_\pm$"], 
              frame=True, gridlines=True, axes=True)
g
Out[30]:
In [31]:
r_max = 10
for a1 in al:
    g += plot(eps_p(a1, r), (r, 1.001*r_min_p(a1), r_max), color=hue(a1/ah), 
              thickness=1.5, legend_label=r"$a = {}\, m$".format(float(a1)))
    ri = 1.00001*r_isco_p(a1)  # 1.00001 to avoid 0/0 for a1=1
    g += point((ri, eps_p(a1, ri)), color=hue(a1/ah), size=60)
for a1 in al:
    g += plot(eps_m(a1, r), (r, 1.001*r_min_m(a1), r_max), thickness=1.5,
              linestyle='--', color=hue(a1/ah))
    ri = r_isco_m(a1)
    if ri < r_max:
        g += point((ri, eps_m(a1, ri)), color=hue(a1/ah), marker=r'$\circ$', size=60)
g += line([(0,1), (r_max, 1)], color='grey')
show(g, ymin=-2, ymax=2.5, xmax=8)
In [32]:
g.save('gek_eps_circ_orb.pdf', ymin=-2, ymax=2.5, xmax=8)
In [33]:
show(g, xmin=4, xmax=r_max, ymin=0.91, ymax=1.03)
In [34]:
g.save('gek_eps_circ_orb_zoom.pdf', xmin=4, xmax=r_max, ymin=0.91, ymax=1.03)
In [35]:
g = Graphics()
for a1 in al[1:]:
    g += plot(l_p(a1, r), (r, 0.001, rs(a1)), color=hue(a1/ah), 
              linestyle='-', thickness=1.5, 
              axes_labels=[r"$r_0/m$", r"$\ell_\pm/m$"],
              frame=True, gridlines=True, axes=True)
show(g, ymin=-10, ymax=10)
In [36]:
for a1 in al:
    g += plot(l_p(a1, r), (r, 1.001*r_min_p(a1), r_max), color=hue(a1/ah), 
              thickness=1.5, legend_label=r"$a = {}\, m$".format(float(a1)))
    ri = 1.00001*r_isco_p(a1)  # 1.00001 to avoid 0/0 for a1=1
    g += point((ri, l_p(a1, ri)), color=hue(a1/ah), size=60)
for a1 in al:
    g += plot(l_m(a1, r), (r, 1.001*r_min_m(a1), r_max), thickness=1.5,
              linestyle='--', color=hue(a1/ah)) 
    ri = r_isco_m(a1)
    if ri < r_max:
        g += point((ri, l_m(a1, ri)), color=hue(a1/ah), marker=r'$\circ$', size=60)
show(g, xmax=8, ymin=-10, ymax=10)
In [37]:
g.save('gek_ell_circ_orb.pdf', xmax=8, ymin=-10, ymax=10)
In [38]:
show(g, xmin=4, xmax=r_max, ymin=-4.5, ymax=4)
In [39]:
g.save('gek_ell_circ_orb_zoom.pdf', xmin=4, xmax=r_max, ymin=-4.5, ymax=4)
In [40]:
r_max = 30
g = Graphics()
for a1 in al[1:]:
    g += parametric_plot((l_p(a1, r), eps_p(a1, r)), (r, 0.01, 0.99*rs(a1)), 
                         color=hue(a1/ah), linestyle=':', thickness=1.5, 
                         axes_labels=[r"$\ell/m$", r"$\varepsilon$"],
                         frame=True, gridlines=True, axes=True)
for a1 in reversed(al):
    ri = 1.00001*r_isco_p(a1)  # 1.00001 to avoid 0/0 for a1=1
    g += parametric_plot((l_p(a1, r), eps_p(a1, r)), (r, 1.01*r_min_p(a1), ri), 
                         color=hue(a1/ah), thickness=1.5, linestyle=':')
    g += parametric_plot((l_p(a1, r), eps_p(a1, r)), (r, ri, r_max), 
                         color=hue(a1/ah), thickness=1.5)
for a1 in al:
    ri = 1.00001*r_isco_m(a1)  # 1.00001 to avoid 0/0 for a1=1
    g += parametric_plot((l_m(a1, r), eps_m(a1, r)), (r, 1.01*r_min_m(a1), ri), 
                         thickness=1.5, linestyle=':', color=hue(a1/ah)) 
    g += parametric_plot((l_m(a1, r), eps_m(a1, r)), (r, ri, r_max), 
                         thickness=1.5, color=hue(a1/ah), 
                         legend_label=r"$a = {}\, m$".format(float(a1))) 
g += line([(-6,1), (6, 1)], color='grey')
show(g, xmin=-5, xmax=5, ymin=-0.4, ymax=1.2, aspect_ratio='automatic')
In [41]:
g.save("gek_ell_eps_circ_orb.pdf", xmin=-5, xmax=5, ymin=-0.4, 
       ymax=1.2, aspect_ratio='automatic')
In [42]:
show(g, xmin=1, xmax=5, ymin=0.5, ymax=1.2, aspect_ratio='automatic')

Orbital angular velocity

In [43]:
omega_p(a, r) = 1/(r^(3/2) + a)
omega_m(a, r) = -1/(r^(3/2) - a)
omega_p(a, r), omega_m(a, r)
Out[43]:
In [44]:
r_max = 10
g = Graphics()
for a1 in al[1:]:
    g += plot(omega_p(a1, r), (r, 0.001, rs(a1)), color=hue(a1/ah), 
              linestyle='-', thickness=1.5, 
              axes_labels=[r"$r_0/m$", r"$m\, \Omega_\pm$"],
              frame=True, gridlines=True, axes=True)
for a1 in al:
    g += plot(omega_p(a1, r), (r, 1.001*r_min_p(a1), r_max), color=hue(a1/ah), 
              thickness=1.5, legend_label=r"$a = {}\, m$".format(float(a1)))
    ri = 1.00001*r_isco_p(a1)  # 1.00001 to avoid 0/0 for a1=1
    g += point((ri, omega_p(a1, ri)), color=hue(a1/ah), size=60)
for a1 in al:
    g += plot(omega_m(a1, r), (r, 1.001*r_min_m(a1), r_max), thickness=1.5,
              linestyle='--', color=hue(a1/ah)) 
    ri = r_isco_m(a1)  
    g += point((ri, omega_m(a1, ri)), marker=r'$\circ$', color=hue(a1/ah), size=60)
show(g, xmax=8)
In [45]:
g.save('gek_omega_circ_orb.pdf', xmax=8)
In [46]:
show(g, xmin=4, xmax=r_max, ymin=-0.15, ymax=0.15)
In [47]:
g.save('gek_omega_circ_orb_zoom.pdf', xmin=4, xmax=r_max, ymin=-0.15, ymax=0.15)

Linear velocity with respect to the ZAMO

In [48]:
V_ZAMO_p(a, r) = (r^2 - 2*a*sqrt(r) + a^2)/((r^(3/2) + a)*sqrt(r^2 - 2*r + a^2))
V_ZAMO_p(a, r)
Out[48]:
In [49]:
V_ZAMO_m(a, r) = -(r^2 + 2*a*sqrt(r) + a^2)/((r^(3/2) - a)*sqrt(r^2 - 2*r + a^2))
V_ZAMO_m(a, r)
Out[49]:
In [50]:
r_max = 10
g = Graphics()
for a1 in al[1:]:
    g += plot(V_ZAMO_p(a1, r), (r, 0.001, rs(a1)), color=hue(a1/ah), 
              linestyle='-', thickness=1.5, 
              axes_labels=[r"$r_0/m$", r"$V_\pm^{(\varphi)}$"],
              frame=True, gridlines=True, axes=True)
for a1 in al:
    g += plot(V_ZAMO_p(a1, r), (r, 1.001*r_min_p(a1), r_max), color=hue(a1/ah), 
              thickness=1.5, legend_label=r"$a = {}\, m$".format(float(a1)))
    ri = 1.00001*r_isco_p(a1)  # 1.00001 to avoid 0/0 for a1=1
    g += point((ri, V_ZAMO_p(a1, ri)), color=hue(a1/ah), size=60)
for a1 in al:
    g += plot(V_ZAMO_m(a1, r), (r, 1.001*r_min_m(a1), r_max), thickness=1.5,
              linestyle='--', color=hue(a1/ah)) 
    ri = r_isco_m(a1)  
    g += point((ri, V_ZAMO_m(a1, ri)), marker=r'$\circ$', color=hue(a1/ah), size=60)
show(g, xmax=8)
In [51]:
g.save('gek_v_zamo_circ_orb.pdf', xmax=8)
In [52]:
a1 = 0.9
plot(r^2 - 2*a1*sqrt(r) + a1^2, (r, 0.98, 1.02))
Out[52]:
In [53]:
plot(r^2 - 6*r - 3*a1^2 + 8*a1*sqrt(r), (r, 0, 3))
Out[53]:
In [54]:
plot(r^2 - 6*r - 3*a1^2 - 8*a1*sqrt(r), (r, 0, 10))
Out[54]:

Behaviour of the function $\mathcal{V}(r)$ near a circular orbit

In [55]:
a1 = 0.9
r_isco_p(a1)
Out[55]:
In [56]:
r1 = 2.
eps1 = eps_p(a1, r1)
l1 = l_p(a1, r1)
V1(r) = V(a1, eps1, l1, r)
g = plot(V1, (1, 5), color="red", thickness=1.5, legend_label=r"$r_0 = 2m$",
         axes_labels=[r"$r/m$", r"$\mathcal{V}(r)$"],
         frame=True, gridlines=True, axes=True)
g += point((r1, 0), color="red", size=60)
r1 = 3.
eps1 = eps_p(a1, r1)
l1 = l_p(a1, r1)
V1(r) = V(a1, eps1, l1, r)
g += plot(V1, (1, 5), color="blue", thickness=1.5, legend_label=r"$r_0 = 3m$")
g += point((r1, 0), color="blue", size=60)
show(g, ymin=-0.04)