ges_effective_potential.ipynb

Kernel: SageMath 8.8

Effective potential for timelike geodesic motion in Schwarzschild spactime

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

Click here to download the worksheet file (ipynb format). To run it, you must start SageMath with the Jupyter notebook, with the command sage -n jupyter

%display latex

version()

Effective potential

var('r L')
V(r,L) = - 1/r + L^2/(2*r^2) - L^2/r^3
V(r,L)

l_values = [0, 2, 3.0, 2*sqrt(3), 3.8, 4.2, 4.6, 5.0]

graph = Graphics()
for l in l_values:
    graph += plot(V(r,l), (r, 0.1, 10), plot_points=300, 
                  legend_label="$\\ell/m = {:.4f}$".format(float(l)),
                  color=hue((1+l)/7), thickness=2,
                  xmin=0, xmax=10, ymin=-1.5, ymax=0.2,
                  axes_labels=[r'$r/m$', r'$V_{\ell}(r)$'])
graph += line([(2,-1.6), (2,0.25)], color='black', thickness=1.5, linestyle='--')
graph

graph.save("ges_eff_pot.pdf", legend_loc="lower right")

graph = Graphics()
for l in l_values:
    graph += plot(V(r,l), (r, 0.1, 25), plot_points=300, 
                  legend_label="$\\ell/m = {:.4f}$".format(float(l)),
                  color=hue((1+l)/7), thickness=2,
                  xmin=0, xmax=25, ymin=-0.16, ymax=0.16,
                  axes_labels=[r'$r/m$', r'$V_{\ell}(r)$'])
graph += line([(2,-0.17), (2,0.17)], color='black', thickness=1.5, linestyle='--')
graph

V(2,L)

Circular orbits

Circular orbits are obtained by imposing $\mathrm{d}r/\mathrm{d}\tau = 0$ and $\mathrm{d}^2r/\mathrm{d}\tau^2 = 0$, which imply $\mathrm{d}V/\mathrm{d}r=0$:

diff(V(r,L),r)

s = solve(diff(V(r,L),r)==0, r)
s

for l in l_values:
    r1 = s[0].subs(L=l).rhs()
    r2 = s[1].subs(L=l).rhs()
    print("l={}: r1={}, r2={}".format(l, r1, r2))

l=0: r1=0, r2=0
l=2: r1=-2*sqrt(-2) + 2, r2=2*sqrt(-2) + 2
l=3.00000000000000: r1=4.50000000000000 - 2.59807621135332*I, r2=4.50000000000000 + 2.59807621135332*I
l=2*sqrt(3): r1=6, r2=6
l=3.80000000000000: r1=4.25210512315547, r2=10.1878948768445
l=4.20000000000000: r1=3.83277632344408, r2=13.8072236765559
l=4.60000000000000: r1=3.61893686280608, r2=17.5410631371939
l=5.00000000000000: r1=3.48612181134003, r2=21.5138781886600

for l in l_values:
    if l>=2*sqrt(3):
        r2 = s[1].subs(L=l).rhs()
        V2 = V(r2,l)
        if l== 2*sqrt(3):
            graph += point((r2,V2), color=hue((1+l)/7), size=80, 
                           marker=r'$\circ$')
        else:
            graph += point((r2,V2), color=hue((1+l)/7), size=60)
graph

graph.save("ges_eff_pot_zoom.pdf")

Coordinate $r$ on a circular orbit of angular momentum $L$:

r_circ1 = s[0].rhs()
r_circ2 = s[1].rhs()
r_circ1, r_circ2

graph = plot(r_circ1, (L, 2*sqrt(3), 5), linestyle='--', thickness=2)
graph += plot(r_circ2, (L, 2*sqrt(3), 5), thickness=2)
show(graph, xmin=3.4, xmax=5, ymin=0, ymax=22, gridlines='minor',
     axes_labels=[r'$\ell/m$', r'$r/m$'])

L_circ = r / sqrt(r-3)
L_circ

graph = plot(L_circ, (r, 3.2, 6), linestyle='--', thickness=2)
graph += plot(L_circ, (r, 6, 12), thickness=2, plot_points=100)
graph += line([(3,0), (3, 8)], color='black', linestyle='--')
graph += line([(6,0), (6, 2*sqrt(3))], color='black', linestyle='--')
show(graph, xmin=2, xmax=12, ymin=3.2, ymax=6, gridlines=True,
     axes_labels=[r'$r/m$', r'$\ell/m$'])

graph.save("ges_ell_circ_orbit.pdf", xmin=2, xmax=12, ymin=3.2, ymax=6, 
           gridlines=True, axes_labels=[r'$r/m$', r'$\ell/m$'])

V_circ = V(r, L_circ).simplify_full()
V_circ

E_circ = sqrt(2*V_circ + 1).simplify_full()
E_circ

graph = plot(E_circ, (r, 3.2, 6), linestyle='--', thickness=2)
graph += plot(E_circ, (r, 6, 12), thickness=2)
graph += line([(3,0), (3, 8)], color='black', linestyle='--')
graph += line([(6,0), (6, 4/(3*sqrt(2)))], color='black', linestyle='--')
graph += line([(4,0), (4, 1)], color='black', linestyle='--')
graph += line([(2,1), (12, 1)], color='red')
show(graph, xmin=2, xmax=12, ymin=0.92, ymax=1.2, gridlines=True,
     axes_labels=[r'$r/m$', r'$\varepsilon$'])

graph.save("ges_ener_circ_orbit.pdf", xmin=2, xmax=12, ymin=0.92, 
           ymax=1.2, gridlines=True,
           axes_labels=[r'$r/m$', r'$\varepsilon$'])

diff(E_circ, r).simplify_full()

E_circ.subs(r=6)

n(_)

(E_circ / L_circ).simplify_full()

graph = parametric_plot((L_circ, E_circ), (r, 3.2, 6), linestyle='--', thickness=2)
graph += parametric_plot((L_circ, E_circ), (r, 6, 65), thickness=2)
show(graph, xmin=3, xmax=8, ymin=0.92, ymax=1.05, aspect_ratio=30, gridlines=True,
     axes_labels=[r'$\ell/m$', r'$\varepsilon$'])

graph.save("ges_circ_eps_ell.pdf", xmin=3, xmax=8, ymin=0.92, ymax=1.05, 
           aspect_ratio=30, gridlines=True,
           axes_labels=[r'$\ell/m$', r'$\varepsilon$'])

Non-circular orbits

l0 = 4.2
E0 = 0.973
V0 = (E0^2-1)/2
V0

sol = solve(V(r,l0) == V0, r)
for s in sol:
    show(CDF(s.rhs()))

r1 = CDF(sol[0].rhs()).real_part()
r2 = CDF(sol[2].rhs()).real_part()
r1, r2

r_circ = r_circ2.subs({L: l0})
r_circ

E_circ.subs(r=r_circ)

graph = plot(V(r,l0), (r, 2.7, 31), plot_points=300, color=hue((1+l0)/7), 
             thickness=2, axes_labels=[r'$r/m$', r'$V_{\ell}(r)$'])
graph += plot(V(r,l0), (r, r1, r2), color=hue((1+l0)/7), thickness=4)
graph += line([(2,-0.17), (2,0.17)], color='black', thickness=1.5, linestyle='--')
graph += line([(r1,V0), (r2,V0)], color='red', thickness=2)
graph += line([(r1,V0), (r1,2*V0)], color='black', thickness=1.5, linestyle=':') + \
         text(r'$r_{\rm per}$', (r1, -0.045), color='black', fontsize=16)
graph += line([(r2,V0), (r2,2*V0)], color='black', thickness=1.5, linestyle=':') + \
         text(r'$r_{\rm apo}$', (r2, -0.046), color='black', fontsize=16)
graph += line([(r1,V0), (-1,V0)], color='black', thickness=1.5, linestyle=':') + \
         text(r'$\frac{1}{2}(\varepsilon^2-1)$', (-3, V0), color='black', 
              fontsize=16)
# graph += point((r_circ, V(r_circ,l0)), color=hue((1+l0)/7), size=80)
show(graph, xmin=0, xmax=30, ymin=-0.04, ymax=0.03, 
     gridlines=True, frame=True, axes=False)

graph.save("ges_eff_pot_bound.pdf", xmin=0, xmax=30, ymin=-0.04, ymax=0.03, 
           gridlines=True, frame=True, axes=False)