Kerr_circular_orbits.ipynb

Kernel: SageMath 9.2

Circular equatorial orbits in Kerr spacetime

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

version()

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

%display latex

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

a = var('a')
r = var('r')
eps = var('eps', latex_name=r'\varepsilon')
l = var('l', latex_name=r'\ell') 

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

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

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

f(r, a) = r*(r-3) + 2*a*sqrt(r)
f

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

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

f(r, a) = r*(r-3) - 2*a*sqrt(r)
f

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

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

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

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

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

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

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

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

plot(r_min_p(a) - rp(a), (a, 0, 1), 
     axes_labels=[r'$a/m$', r'$(r_{\rm ph}^+ - r_+)/m$']) 

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

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

plot(rs, (0, 1), axes_labels=[r'$a/m$', r'$r_{\rm ph}^*/m$']) 

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

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) 

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

graph.save('gek_circ_orb_lim.pdf')

ISCO

Z1 = 1 + (1 - a^2)^(1/3)*((1 + a)^(1/3) + (1 - a)^(1/3))
Z1

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

g = plot(Z1, (a, 0, 1), axes_labels=[r"$a/m$", r"$Z_1$"],
         thickness=1.5, frame=True, gridlines=True, axes=True)
g

g.save("gek_Z1.pdf")

Z2 = sqrt(Z1^2 + 3*a^2)
Z2

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

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$

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)

eps_p(a, r)

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

eps_m(a, r)

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

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)

l_p(a, r)

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

l_m(a, r)

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

num_l_p(a, r) = r^2 + a^2 - 2*a*sqrt(r)

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

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)

g.save('gek_eps_circ_orb.pdf', ymin=-2, ymax=2.5, xmax=8)

show(g, xmin=4, xmax=r_max, ymin=0.91, ymax=1.03)

g.save('gek_eps_circ_orb_zoom.pdf', xmin=4, xmax=r_max, ymin=0.91, ymax=1.03)

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)

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)

g.save('gek_ell_circ_orb.pdf', xmax=8, ymin=-10, ymax=10)

show(g, xmin=4, xmax=r_max, ymin=-4.5, ymax=4)

g.save('gek_ell_circ_orb_zoom.pdf', xmin=4, xmax=r_max, ymin=-4.5, ymax=4)

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

g.save("gek_ell_eps_circ_orb.pdf", xmin=-5, xmax=5, ymin=-0.4, 
       ymax=1.2, aspect_ratio='automatic')

show(g, xmin=1, xmax=5, ymin=0.5, ymax=1.2, aspect_ratio='automatic')

Orbital angular velocity

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)

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

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)

g.save('gek_omega_circ_orb.pdf', xmax=8)

show(g, xmin=4, xmax=r_max, ymin=-0.15, ymax=0.15)

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

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)

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

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)

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

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)

g.save('gek_v_zamo_circ_orb.pdf', xmax=8)

a1 = 0.9
plot(r^2 - 2*a1*sqrt(r) + a1^2, (r, 0.98, 1.02))

plot(r^2 - 6*r - 3*a1^2 + 8*a1*sqrt(r), (r, 0, 3))

plot(r^2 - 6*r - 3*a1^2 - 8*a1*sqrt(r), (r, 0, 10))

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

a1 = 0.9
r_isco_p(a1)

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

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)

g.save('gek_V_stability.pdf', ymin=-0.04)

Stability of orbits with $r_0 < r_*$:

g = Graphics()
for r1, color in zip([0.4, 0.5], ['blue', 'turquoise']):
    eps1 = eps_p(a1, r1)
    l1 = l_p(a1, r1)
    V1(r) = V(a1, eps1, l1, r)
    g += plot(V1, (0.3, 0.6), color=color, axes_labels=[r"$r/m$", r"$\mathcal{V}(r)$"],
         frame=True, gridlines=True, axes=True)
    g += point((r1, 0), color=color, size=60)
show(g, ymin=-0.01)

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}^+$',
             frame=True, gridlines=True, axes=False) \
        + plot(r_min_m, (0, 1), linestyle='--', color='green', thickness=1.5, 
               legend_label=r'$r_{\rm ph}^-$') \
        + plot(rs, (0, 1), color='green', linestyle=':', thickness=2, 
               legend_label=r'$r_{\rm ph}^*$') \
        + 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 += plot(r_isco_p, (0, 1), color='red', thickness=1.5, legend_label=r'$r_{\rm ISCO}^+$') \
         + plot(r_isco_m, (0, 1), color='red', linestyle='--', thickness=1.5, 
                legend_label=r'$r_{\rm ISCO}^-$')
graph.set_legend_options(handlelength=2, loc=(1.02, 0.36))
graph

Check that the ISCO correspond to a minimum of $\varepsilon_\pm(r_0)$ and $\ell_\pm(r_0)$

a1 = 0.9
g = plot(eps_p(a1, r), (r, 2.3, 2.35), color=hue(a1/ah), 
              linestyle='-', thickness=1.5,
              axes_labels=[r"$r_0/m$", r"$\varepsilon_+$"], 
              frame=True, gridlines=True, axes=True)
ri = r_isco_p(a1)
g += point((ri , eps_p(a1, ri)), color="blue", size=60)
g

a1 = 0.9
g = plot(l_p(a1, r), (r, 2.3, 2.35), color=hue(a1/ah), 
              linestyle='-', thickness=1.5,
              axes_labels=[r"$r_0/m$", r"$\ell_+/m$"], 
              frame=True, gridlines=True, axes=True)
ri = r_isco_p(a1)
g += point((ri , l_p(a1, ri)), color="blue", size=60)
g

Marginally bound circular orbit

r_mb_p(a) = 2 - a + 2*sqrt(1 - a)
r_mb_p

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

r_mb_m(a) = 2 + a + 2*sqrt(1 + a)
r_mb_m

$$\renewcommand{\Bold}[1]{\mathbf{#1}}a \ {\mapsto}\ a + 2 \, \sqrt{a + 1} + 2$$

graph += plot(r_mb_p, (0, 1), color='burlywood', thickness=1.5, 
              legend_label=r'$r_{\rm mb}^+$') \
         + plot(r_mb_m, (0, 1), color='burlywood', linestyle='--', thickness=1.5, 
                legend_label=r'$r_{\rm mb}^-$')
graph

graph.save("gek_circ_orb_isco.pdf")