Allgemeine Relativitätstheorie mit dem Computer

General Theory of Relativity on the Computer

Vorlesung gehalten an der J.W.Goethe-Universität in Frankfurt am Main (Sommersemester 2021)

von Dr.phil.nat. Dr.rer.pol. Matthias Hanauske

Frankfurt am Main 04.04.2021

Erster Vorlesungsteil: Allgemeine Relativitätstheorie mit Python

Bewegung eines Probekörpers um ein nicht-rotierendes schwarzes Loch in der Ebene

Teil III: Der ISCO und die Photonensphäre

Im Folgenden werden die beiden Konzepte

a) die innerste stabile kreisförmige Bahnbewegung (Innermost Stable Circular Orbit (ISCO)) eines Probekörpers um ein nicht-rotierendes schwarzes Loch

b) Die Photonensphäre eines nicht-rotierendes schwarzes Loches

vorgestellt.

a) Die innerste stabile kreisförmige Bahnbewegung (Innermost Stable Circular Orbit (ISCO)) eines Probekörpers

Wir beginnen mit dem ISCO eines schwarzen Loches und betrachten wieder die Geodätengleichung in vorgegebener Schwarzschild Raumzeit.

Zunächst wird das Python Modul "EinsteinPy" eingebunden, welches auf dem Modul SymPy basiert und symbolische Berechnungen in der Allgemeinen Relativitätstheorie relativ einfach möglich macht. Dann definieren wir die zugrundeliegenden Koordinaten und die kovariante Form der Raumzeit-Metrik der Schwarzschildmetrik.

from sympy import *
init_printing()
from einsteinpy.symbolic import *

t, r, theta, phi, M = symbols('t, r, theta, phi, M')
Metric = diag((1-2*M/r), -1/(1-2*M/r), -r**2, -r**2*sin(theta)**2).tolist()
g = MetricTensor(Metric, [t, r, theta, phi])

g.tensor()

$$\left[\begin{matrix}- \frac{2 M}{r} + 1 & 0 & 0 & 0\\0 & - \frac{1}{- \frac{2 M}{r} + 1} & 0 & 0\\0 & 0 & - r^{2} & 0\\0 & 0 & 0 & - r^{2} \sin^{2}{\left (\theta \right )}\end{matrix}\right]$$

Die Chistoffel Symbole (zweiter Art): $$ \Gamma^{\sigma}_{\mu \nu} = \frac{1}{2}g^{\sigma \alpha} \left( g_{\nu \alpha| \mu} + g_{\mu \alpha| \nu} - g_{\mu \nu| \alpha} \right)$$

Hier speziell $$ \Gamma^{1}_{1 1} = \Gamma^{r}_{r r}$$

chr = ChristoffelSymbols.from_metric(g)
print(chr.config)
chr.tensor()[1,1,1]

ull

$$\frac{2 M \left(\frac{M}{r} - \frac{1}{2}\right)}{r^{2} \left(- \frac{2 M}{r} + 1\right)^{2}}$$

Geodätische Bewegung eines Probekörpers

Die Geodätengleichung stellt ein System gekoppelter nichtlinearer Differentialgleichungen dar

$$ \begin{eqnarray} && \frac{d^2 t}{d\lambda^2} = - \Gamma^0_{\nu \rho} \frac{d x^\nu}{d\lambda} \frac{d x^\rho}{d\lambda} \\ && \frac{d^2 r}{d\lambda^2} = - \Gamma^1_{\nu \rho} \frac{d x^\nu}{d\lambda} \frac{d x^\rho}{d\lambda}\\ && \frac{d^2 \theta}{d\lambda^2} = - \Gamma^2_{\nu \rho} \frac{d x^\nu}{d\lambda} \frac{d x^\rho}{d\lambda}\\ && \frac{d^2 \phi}{d\lambda^2} = - \Gamma^3_{\nu \rho} \frac{d x^\nu}{d\lambda} \frac{d x^\rho}{d\lambda} \quad , \end{eqnarray} $$

wobei $\lambda$ ein affiner Parameter (z.B. die Eigenzeit $\tau$), t, r, $\theta$ und $\phi$ die Koordinaten und $\Gamma^\mu_{\nu \rho}$ die Christoffel Symbole zweiter Art darstellen.

Wir definieren die Geodätengleichung:

tau = symbols('tau')
x0 = Function('t')(tau)
x1 = Function('r')(tau)
x2 = Function('\\theta')(tau)
x3 = Function('\\phi')(tau)

xmu=GenericVector([x0, x1, x2, x3],[t, r, theta, phi], config='u',parent_metric=Metric)
Geod1=diff(xmu.tensor(),tau,tau)
Geod2T=tensorproduct(chr.tensor(),diff(xmu.tensor(),tau),diff(xmu.tensor(),tau))
Geod2=tensorcontraction(Geod2T, (1, 3),(2,4))
Geod=Geod1+Geod2
GeodEq0=Eq(Geod[0],0)
GeodEq1=Eq(Geod[1],0)
GeodEq2=Eq(Geod[2],0)
GeodEq3=Eq(Geod[3],0)

Wir lassen im Folgenden wieder nur ebene Bewegungen zu ($\theta=\frac{\pi}{2}, {\rm cos}(\theta)=0, {\rm sin}(\theta)=1, \frac{d \theta}{d\tau}=0, \frac{d^2 \theta}{d\tau^2}=0$). Unter diesen Annahmen sind nur drei der vier Geodätengleichungen relevant und diese lauten dann:

GeodEq0=GeodEq0.subs({(theta,pi/2),(diff(x2,tau),0)})
GeodEq0

$$\frac{2 M \frac{d}{d \tau} r{\left (\tau \right )} \frac{d}{d \tau} t{\left (\tau \right )}}{r^{2} \left(- \frac{2 M}{r} + 1\right)} + \frac{d^{2}}{d \tau^{2}} t{\left (\tau \right )} = 0$$

GeodEq1=GeodEq1.subs({(theta,pi/2),(diff(x2,tau),0)})
GeodEq1

$$- \frac{2 M \left(\frac{M}{r} - \frac{1}{2}\right) \left(\frac{d}{d \tau} t{\left (\tau \right )}\right)^{2}}{r^{2}} + \frac{2 M \left(\frac{M}{r} - \frac{1}{2}\right) \left(\frac{d}{d \tau} r{\left (\tau \right )}\right)^{2}}{r^{2} \left(- \frac{2 M}{r} + 1\right)^{2}} + 2 r \left(\frac{M}{r} - \frac{1}{2}\right) \left(\frac{d}{d \tau} \phi{\left (\tau \right )}\right)^{2} + \frac{d^{2}}{d \tau^{2}} r{\left (\tau \right )} = 0$$

GeodEq3=GeodEq3.subs({(theta,pi/2),(diff(x2,tau),0)})
GeodEq3

$$\frac{d^{2}}{d \tau^{2}} \phi{\left (\tau \right )} + \frac{2 \frac{d}{d \tau} \phi{\left (\tau \right )} \frac{d}{d \tau} r{\left (\tau \right )}}{r} = 0$$

Um dieses System von drei Differentialgleichungen zweiter Ordnung mit Python nummerisch lösen zu können, müssen wir es zunächst in ein System von sechs Differentialgleichungen erster Ordnung umschreiben. Wir definieren dazu formal: $y_1=\frac{dt}{d\tau}$, $y_2=\frac{dr}{d\tau}$, $y_3=\frac{d\phi}{d\tau}$, $y_4=\frac{dy_1}{d\tau}=\frac{d^2t}{d\tau^2}$, $y_5=\frac{dy_1}{d\tau}=\frac{dy_2}{d\tau}=\frac{d^2r}{d\tau^2}$ und $y_6=\frac{dy_3}{d\tau}=\frac{d^2\phi}{d\tau^2}$. Wir legen des Weiteren die Masse des schwarzen Loches auf $M=1$ fest.

y1, y2, y3 = symbols('y_1, y_2, y_3')

F0=solve(GeodEq0,diff(x0,tau,tau))[0]
F1=solve(GeodEq1,diff(x1,tau,tau))[0]
F3=solve(GeodEq3,diff(x3,tau,tau))[0]

Eq(diff(x0,tau,tau),F0.subs({(diff(x0,tau),y1),(diff(x1,tau),y2),(diff(x3,tau),y3),(M,1)}))

$$\frac{d^{2}}{d \tau^{2}} t{\left (\tau \right )} = \frac{2 y_{1} y_{2}}{r \left(- r + 2\right)}$$

Eq(diff(x1,tau,tau),F1.subs({(diff(x0,tau),y1),(diff(x1,tau),y2),(diff(x3,tau),y3),(M,1)}))

$$\frac{d^{2}}{d \tau^{2}} r{\left (\tau \right )} = \frac{y_{2}^{2}}{\left(- r + 2\right)^{2}} + \frac{y_{3}^{2} \left(r - 2\right)^{3}}{\left(- r + 2\right)^{2}} - \frac{2 y_{2}^{2}}{r \left(- r + 2\right)^{2}} - \frac{y_{1}^{2}}{r^{2}} + \frac{2 y_{1}^{2}}{r^{3}}$$

Eq(diff(x3,tau,tau),F3.subs({(diff(x0,tau),y1),(diff(x1,tau),y2),(diff(x3,tau),y3),(M,1)}))

$$\frac{d^{2}}{d \tau^{2}} \phi{\left (\tau \right )} = - \frac{2 y_{2} y_{3}}{r}$$

y4=lambdify((r,y1,y2,y3), F0.subs({(diff(x0,tau),y1),(diff(x1,tau),y2),(diff(x3,tau),y3),(M,1)}))

y5=lambdify((r,y1,y2,y3), F1.subs({(diff(x0,tau),y1),(diff(x1,tau),y2),(diff(x3,tau),y3),(M,1)}))

y6=lambdify((r,y1,y2,y3), F3.subs({(diff(x0,tau),y1),(diff(x1,tau),y2),(diff(x3,tau),y3),(M,1)}))

Wir definieren das System der gekoppelten sechs Differentialgleichungen:

def DGLsys(vy, time):
    t, r, phi, y1, y2, y3 = vy
    dt = y1
    dr = y2
    dphi = y3
    dy1 = y4(r,y1,y2,y3)
    dy2 = y5(r,y1,y2,y3)
    dy3 = y6(r,y1,y2,y3)
    return np.array([dt,dr,dphi,dy1,dy2,dy3])

import numpy as np
import matplotlib.pyplot as plt 
import matplotlib
from scipy import integrate

Literatur:

Literatur:

1. General relativity : An introduction for physicists von M. P. Hobson, G. P. Efstathiou und A. N. Lasenby
2. Gravity : An introduction to Einstein's general relativity von James B. Hartle
3. Allgemeine Relativitätstheorie von Torsten Fließbach
4. Relativistic hydrodynamics von Luciano Rezzolla und Olindo Zanotti

Man kann zeigen (siehe z.B. Seite 206 in General relativity : An introduction for physicists by M. P. Hobson, G. P. Efstathiou and A. N. Lasenby), dass sich die erste und vierte Gleichung des Systems der Differentialgleichungen der Geodätengleichung in die folgenden Gleichungen umschreiben lässt: $$ \begin{eqnarray} \hbox{1. Gleichung:}&& \frac{d}{d\tau} \left[ \left( 1 - \frac{2 M}{r} \right) \frac{dt}{d\tau} \right] = 0 \\ \rightarrow && \left( 1 - \frac{2 M}{r} \right) \frac{dt}{d\tau} = e = \hbox{const} \\ \hbox{4. Gleichung:}&& \frac{d}{d\tau} \left( r^2 \hbox{sin}^2(\theta) \frac{d\phi}{d\tau} \right) = 0 \\ \rightarrow && r^2 \hbox{sin}^2(\theta) \frac{d\phi}{d\tau} = l = \hbox{const} \quad , \end{eqnarray} $$ wobei die während der Bewegungen erhaltenen Größen e (Teilchenenergie pro Masse, $e=\frac{E}{m}$) und l (Drehimpuls pro Masse m, $l=\frac{L}{m}$) sich aus der Definition des Viererimpulses $p_\mu = m \, u_\mu = m \, g_{\mu\nu} \frac{dx^\nu}{d\tau} = m \, \left( g_{tt} \, \frac{dt}{d\tau}, g_{rr} \, \frac{dr}{d\tau}, g_{\theta\theta} \, \frac{d\theta}{d\tau}, g_{\phi\phi} \, \frac{d\phi}{d\tau} \right) $ ergeben (Voraussetzung: diagonale Form der Metrik). Für die Schwarzschildmetrik erhält man:

$$ \begin{eqnarray} && p_t = p_0 = m \frac{dx_0}{d\tau} = m g_{00} \frac{dx^0}{d\tau} = m g_{00} \frac{dt}{d\tau} = m \left( 1 - \frac{2 M}{r} \right) \frac{dt}{d\tau} = m \, e \\ && p_\phi = p_3 = m \frac{dx_3}{d\tau} = m g_{33} \frac{dx^3}{d\tau} = m g_{33} \frac{d\phi}{d\tau} = - m \left( r^2 \hbox{sin}^2(\theta) \right) \frac{d\phi}{d\tau} = - m \, l \end{eqnarray} $$

Die Klassifikation möglicher Bahnen von Probekörpern in vorgegebener Schwarzschild-Raumzeit kann mittels eines definierten effektiven Potentials illustriert werden. Dieses Potential hängt von dem, bei der Bewegung erhaltenen Drehimpuls pro Masse m ab. Die im Zentralfeld möglichen Bewegungen werden mittels zweier erhaltener Größen (l: Drehimpuls pro Masse m und e: Energie pro Masse) charakterisiert. Die Definition des effektiven Potential erfolgt mittels der radialen, 2. Geodätengleichung. So definieren z.B. die Literaturangaben 1-3 (siehe Angaben oben) das effektive Potentail wie folgt: $$ \begin{eqnarray} \hbox{2. Gleichung:}&\rightarrow&\frac{1}{2} \left( \frac{dr}{d\tau} \right)^2 + V(r,M,l) = \frac{1}{2} \left( e^2 -1 \right) \\ \hbox{wobei:} && V(r,M,l) = \frac{l^2}{2 r^2} \left( 1 - \frac{2 M}{r} \right) - \frac{M}{r} \quad . \end{eqnarray} $$ In der Literaturangabe 4 wird das Potential V(r,M,l) hingegen wie folgt definiert:

$$ \begin{eqnarray} \hbox{2. Gleichung:}&\rightarrow& \left( \frac{dr}{d\tau} \right)^2 + \left( V(r,M,l) \right)^2 = e^2 \\ \hbox{wobei:} && V(r,M,l) = \sqrt{ \left( 1 - \frac{2 M}{r} \right) \left( 1 + \frac{l^2}{r^2} \right) } \quad . \end{eqnarray} $$

Im Folgenden betrachten wir uns die Form des effektiven Potentials (Literaturangabe 4) bei Variation des Parameters $l$ ($l \in [3,4.5]$ und $M=1$).

l = symbols('l')
VeffRez=sqrt(1-2*M/r)*sqrt(1+l**2/r**2)
VeffNewton=1-M/r + l**2/(2*r**2)
VeffRezl=lambdify((r,l),VeffRez.subs({(M,1)}))
VeffNewtonl=lambdify((r,l),VeffNewton.subs({(M,1)}))
rval = np.linspace(2.2,30, 1001)

params = {
    'figure.figsize'    : [10,8],
#    'text.usetex'       : True,
    'axes.titlesize' : 14,
    'axes.labelsize' : 16,  
    'xtick.labelsize' : 14 ,
    'ytick.labelsize' : 14 
}
matplotlib.rcParams.update(params) 

plt.cla()
plt.xlim(2.3, 20)
plt.ylim(0.8, 1.1)
plt.ylabel(r"$\rm V_{eff}(r)$")
plt.xlabel(r"$\rm Radius \, r \,[km]$")
plt.plot(rval,VeffRezl(rval,4.5), linewidth=1, color="yellow",label=r"$\rm l=4.5$");
plt.scatter(16.587632977355,VeffRezl(16.587632977355,4.5), c="yellow", s=60, marker='o');
plt.plot(rval,VeffRezl(rval,4), linewidth=1, color="green",label=r"$\rm l=4$");
plt.scatter(12,VeffRezl(12,4), c="green", s=60, marker='o');
plt.plot(rval,VeffRezl(rval,3.7), linewidth=1, color="red",label=r"$\rm l=3.7$");
plt.scatter(9.25,VeffRezl(9.25,3.7), c="red", s=60, marker='o');
plt.plot(rval,VeffRezl(rval,float(2*sqrt(3))), linewidth=1, color="black",label=r"$\rm l=2 \, \sqrt{3}$");
plt.scatter(6,VeffRezl(6,float(2*sqrt(3))), c="black", s=60, marker='o');
plt.plot(rval,VeffRezl(rval,3), linewidth=1, color="blue",label=r"$\rm l=3$");
plt.legend(loc='lower right', fontsize=18);

In der oberen Abbildung erkennt man, dass bei kleinen Drehimpuls Werten mit $l < 2 \, \sqrt{3} \approx 3.4641$ keine gebundenen stabilen Orbits des Probekörpers mehr möglich sind (fehlendes Minimum in der blauen Kurve). Für $l > 2 \, \sqrt{3}$ sind kreisförmige und elliptische gebundene Orbits möglich, wobei der jeweilige Punkt auf den Kurven den Radiuswert der kreisförmigen Bahn des Probekörpers kennzeichnet. Der Punkt auf der schwarzen Kurve stellt die Grenze der innersten stabilen kreisförmigen Bahnbewegungen dar, der sogenannte "Innermost Stable Circular Orbit (ISCO)" eines Probekörpers um ein nicht-rotierendes schwarzes Loch. Er befindet sich bei $r_{\rm ISCO}=6M$ und $l_{\rm ISCO}=2 \, \sqrt{3}\,M$, was wir weiter unten zeigen werden.

Die Existenz eines solchen ISCOs ist ein allgemein-relativistischer Effekt der in der Newtonschen Theorie nicht vorhanden ist. Um dies zu verdeutlichen, vergleichen wir das effektive Potential der allgemeinen Relativitätstheorie $V_{eff}(r,l)$ mit dem effektiven Newtonschen Potential $V_{Newton}(r,l)$. Für große Abstände geht das effektive Potential in das Newtonsche über und es gilt (für das effektive Potential der Literaturangabe 4):

$$ V_{eff}(r,l) \approx 1 - \frac{M}{r} + \frac{l^2}{2 \, r^2} = V_{Newton}(r,l) $$

Wir betrachten uns nun die Form dieses effektiven Newtonschen Potentials bei Variation des Parameters $l$ ($l \in [1,4]$ und $M=1$).

rval = np.linspace(0.1,30, 1001)
VeffNewton=1-M/r + l**2/(2*r**2)
VeffNewtonl=lambdify((r,l),VeffNewton.subs({(M,1)}))

plt.cla()
plt.xlim(0.1, 20)
plt.ylim(0.45, 1.1)
plt.ylabel(r"$\rm V_{Newton}(r)$")
plt.xlabel(r"$\rm Radius \, r \,[km]$")
plt.plot(rval,VeffNewtonl(rval,4), linewidth=1, color="green",label=r"$\rm l=4$");
rmin=solve(Eq(diff(VeffNewton.subs({(M,1),(l,4)}),r),0),r)[0]
plt.scatter(rmin,VeffNewtonl(rmin,4), c="green", s=60, marker='o');
plt.plot(rval,VeffNewtonl(rval,3.7), linewidth=1, color="red",label=r"$\rm l=3.7$");
rmin=solve(Eq(diff(VeffNewton.subs({(M,1),(l,3.7)}),r),0),r)[0]
plt.scatter(rmin,VeffNewtonl(rmin,3.7), c="red", s=60, marker='o');
plt.plot(rval,VeffNewtonl(rval,float(2*sqrt(3))), linewidth=1, color="black",label=r"$\rm l=2 \,\sqrt{3}$");
rmin=solve(Eq(diff(VeffNewton.subs({(M,1),(l,float(2*sqrt(3)))}),r),0),r)[0]
plt.scatter(rmin,VeffNewtonl(rmin,float(2*sqrt(3))), c="black", s=60, marker='o');
plt.plot(rval,VeffNewtonl(rval,3), linewidth=1, color="blue",label=r"$\rm l=3$");
rmin=solve(Eq(diff(VeffNewton.subs({(M,1),(l,3)}),r),0),r)[0]
plt.scatter(rmin,VeffNewtonl(rmin,3), c="blue", s=60, marker='o');

plt.plot(rval,VeffNewtonl(rval,2.5), linewidth=1, color="orange",label=r"$\rm l=2.5$");
rmin=solve(Eq(diff(VeffNewton.subs({(M,1),(l,2.5)}),r),0),r)[0]
plt.scatter(rmin,VeffNewtonl(rmin,2.5), c="orange", s=60, marker='o');

plt.plot(rval,VeffNewtonl(rval,2), linewidth=1, color="brown",label=r"$\rm l=2$");
rmin=solve(Eq(diff(VeffNewton.subs({(M,1),(l,2)}),r),0),r)[0]
plt.scatter(rmin,VeffNewtonl(rmin,2), c="brown", s=60, marker='o');

plt.plot(rval,VeffNewtonl(rval,1), linewidth=1, color="grey",label=r"$\rm l=1$");
rmin=solve(Eq(diff(VeffNewton.subs({(M,1),(l,1)}),r),0),r)[0]
plt.scatter(rmin,VeffNewtonl(rmin,1), c="grey", s=60, marker='o');

plt.legend(loc='lower right', fontsize=18);

Man erkennt, dass stabile kreisförmige Bahnen in der Newtonschen Theorie auch für kleine Abstände (kleine Drehimpulse) möglich sind und eine Grenze der innersten stabilen kreisförmigen Bahnbewegung (ein ISCO) nicht existiert. Zusätzlich sehen wir, dass Probekörper nie in das Zentrum bei $r=0$ fallen können, da es eine unendlich starke Potentialbarriere für kleine Radien gibt.

In der allgemeinen Relativitätstheorie sind jedoch nur stabile Bahnen möglich, falls der Wert des Drehimpulses l oberhalb einer gewissen Grenze ist, da nur dann ein Minimum im Potential vorhanden ist. Kreisförmige Bahnbewegungen sind dadurch charakterisiert, dass der Wert des Radius sich im Laufe der Zeit nicht verändert und somit sich der radiale Abstand des Probekörpers vom schwarzen Loch gerade im Minimum des effektiven Potentials befindet. Für kreisförmige Bahnbewegungen muss somit gelten: $\frac{dV(r,M,l)}{dr}=0$. Dies führt zu folgender Gleichung:

dr_VeffRez=diff(VeffRez,r)
EqCircular=Eq(dr_VeffRez,0)
EqCircular

$$\frac{M \sqrt{\frac{l^{2}}{r^{2}} + 1}}{r^{2} \sqrt{- \frac{2 M}{r} + 1}} - \frac{l^{2} \sqrt{- \frac{2 M}{r} + 1}}{r^{3} \sqrt{\frac{l^{2}}{r^{2}} + 1}} = 0$$

Löst man diese Gleichung nach r auf, so erhält man zwei Lösungen, wobei die zweite (positives Vorzeichen) dem stabilen Minimum und die erste (negatives Vorzeichen) dem instabilen Maximum entspricht:

$$ \begin{eqnarray} \frac{dV(r,M,l)}{dr}=0 \quad \Rightarrow \, r_\pm = \frac{l}{2 M} \left( l \pm \sqrt{l^2 -12 M^2} \right) \end{eqnarray} $$

LoesCircular=solve(EqCircular,r)
LoesCircular[0].simplify(),LoesCircular[1].simplify()

$$\left ( \frac{l \left(l - \sqrt{- 12 M^{2} + l^{2}}\right)}{2 M}, \quad \frac{l \left(l + \sqrt{- 12 M^{2} + l^{2}}\right)}{2 M}\right )$$

Die innerste kreisförmige Bahnbewegung (ISCO) hat zusätzlich die Eigenschaft, dass das Maximum und Minimum zu einem Sattelpunkt zusammen fällt - es muss somit zusätzlich gelten: $\frac{d^2V(r,M,l)}{d^2r}=0$. Dies führt zu folgender Gleichung:

drdr_VeffRez=diff(VeffRez,r,r)
EqSattel=Eq(drdr_VeffRez,0)
EqSattel

$$\frac{- \frac{M^{2} \sqrt{\frac{l^{2}}{r^{2}} + 1}}{r \left(- \frac{2 M}{r} + 1\right)^{\frac{3}{2}}} - \frac{2 M l^{2}}{r^{2} \sqrt{- \frac{2 M}{r} + 1} \sqrt{\frac{l^{2}}{r^{2}} + 1}} - \frac{2 M \sqrt{\frac{l^{2}}{r^{2}} + 1}}{\sqrt{- \frac{2 M}{r} + 1}} - \frac{l^{4} \sqrt{- \frac{2 M}{r} + 1}}{r^{3} \left(\frac{l^{2}}{r^{2}} + 1\right)^{\frac{3}{2}}} + \frac{3 l^{2} \sqrt{- \frac{2 M}{r} + 1}}{r \sqrt{\frac{l^{2}}{r^{2}} + 1}}}{r^{3}} = 0$$

Wir setzen in diese Gleichung für r die Lösung von "LoesCircular" ein:

EqSattel1=EqSattel.subs({(r,LoesCircular[1])}).simplify()
EqSattel1

$$\frac{16 M^{4} \left(4 M^{2} + \left(l + \sqrt{- 12 M^{2} + l^{2}}\right)^{2}\right) \left(- 4 M^{2} + l^{2} + l \sqrt{- 12 M^{2} + l^{2}}\right) \left(- 4 M^{2} \left(4 M^{2} + \left(l + \sqrt{- 12 M^{2} + l^{2}}\right)^{2}\right) \left(- 4 M^{2} + l^{2} + l \sqrt{- 12 M^{2} + l^{2}}\right) - 4 M^{2} \left(- 4 M^{2} + l^{2} + l \sqrt{- 12 M^{2} + l^{2}}\right)^{2} + \left(4 M^{2} + \left(l + \sqrt{- 12 M^{2} + l^{2}}\right)^{2}\right)^{2} \left(4 M^{2} - l^{2} - l \sqrt{- 12 M^{2} + l^{2}}\right) - \left(4 M^{2} + \left(l + \sqrt{- 12 M^{2} + l^{2}}\right)^{2}\right) \left(M^{2} \left(4 M^{2} + \left(l + \sqrt{- 12 M^{2} + l^{2}}\right)^{2}\right) - 3 \left(- 4 M^{2} + l^{2} + l \sqrt{- 12 M^{2} + l^{2}}\right)^{2}\right)\right)}{l^{5} \left(\frac{4 M^{2} + \left(l + \sqrt{- 12 M^{2} + l^{2}}\right)^{2}}{\left(l + \sqrt{- 12 M^{2} + l^{2}}\right)^{2}}\right)^{\frac{5}{2}} \left(\frac{- 4 M^{2} + l^{2} + l \sqrt{- 12 M^{2} + l^{2}}}{l \left(l + \sqrt{- 12 M^{2} + l^{2}}\right)}\right)^{\frac{5}{2}} \left(l + \sqrt{- 12 M^{2} + l^{2}}\right)^{11}} = 0$$

Lösen wir diese Gleichung nach l auf erhalten wir zwei Lösungen ($l_{\rm ISCO}=+/- \, 2\sqrt{3}\,M$), wobei diese einem rechts bzw. links rotierenden Probekörper am ISCO entsprechen.

LoesEqSattel1=solve(EqSattel1,l)
LoesEqSattel1

$$\left [ - 2 \sqrt{3} M, \quad 2 \sqrt{3} M\right ]$$

Setzen wir eine dieser Lösungen in "LoesCircular" ein, erhalten wir den Radiuswert der innersten kreisförmigen Bahnbewegung ($r_{\rm ISCO}=6M$):

LoesCircular[1].subs({(l,LoesEqSattel1[1])})

$$6 M$$

Zusammenfassend erhalten wir somit:

$$ \begin{eqnarray} \frac{d V(r,M,l)}{dr}=0 \,\,,\,\,\, \frac{d^2 V(r,M,l)}{dr^2}=0 \quad \underbrace{\Rightarrow}_{\hbox{ISCO}} r=6 M \,\,,\,\,\, l = 2 \sqrt{3} M \end{eqnarray} $$

Wir wollen uns nun drei kreisförmige Bahnbewegung um ein nicht-rotierendes schwarzes Loch (Masse $M=1$) darstellen, wobei eine davon der ISCO ist. Wir wählen zur Zeit $t=\tau=0$ die Anfangswerte der Radien zu $r_0=12, 9.5$ und $6$ km und starten jeweils bei einem Winkel $\phi_0=0$. Die Anfangsgeschwindigkeit in radialer Richtung sei $\left. \frac{dr}{d\tau} \right|_{\tau=0}=0$ und die Anfangsgeschwindigkeit in $\phi$-Richtung ergibt sich durch die festgelegten, konstanten Drehimpuls Werte der kreisförmigen Bahnbewegungen mittels der Gleichung $\left. \frac{d\phi}{d\tau} \right|_{\tau=0} = \frac{l}{r_0^2}$ ($l=4, 3.7$ und $2\sqrt{3}$). Die weitere nötige Anfangsbedingung ( $\left. \frac{dt}{d\tau} \right|_{\tau=0}$ ) erhalten wir mittels der Bedingung $\left( \frac{ds}{d\tau} \right)^2= g_{\mu\nu} \frac{dx^\mu}{d\tau} \frac{dx^\nu}{d\tau} = g_{\mu\nu} u^\mu u^\nu = 1$, die für massive Probekörper immer erfüllt sein muss ($u^\mu$ ist die Vierergeschwindigkeit des Probekörpers). Die Gleichung für das infinitesimale Weglängenelement $ds^2=1$ in der Ebene lautet:

dt, dr, dtheta, dphi = symbols('dt, dr, d\\theta, d\\phi')
dx=GenericVector([dt, dr, dtheta, dphi],[t, r, theta, phi], config='u',parent_metric=Metric)
dxdx=tensorproduct(dx.tensor(),dx.tensor())
dxdxT=Tensor(dxdx,config="uu")
ds2a=tensorproduct(g.tensor(),dxdxT.tensor())
ds2aT=Tensor(ds2a,config="lluu")
ds2=tensorcontraction(ds2aT.tensor(), (0, 2),(1, 3))
Eq_ds=Eq(ds2.subs({(theta,pi/2),(dtheta,0)}),1)
Eq_ds

$$- d\phi^{2} r^{2} - \frac{dr^{2}}{- \frac{2 M}{r} + 1} + dt^{2} \left(- \frac{2 M}{r} + 1\right) = 1$$

Wir lösen diese Gleichung nach $dt$ auf:

LoesEq_ds=solve(Eq_ds,dt)
LoesEq_ds

$$\left [ \frac{\sqrt{r \left(- 2 M d\phi^{2} r^{2} - 2 M + d\phi^{2} r^{3} + dr^{2} r + r\right)}}{- 2 M + r}, \quad \frac{\sqrt{r \left(- 2 M d\phi^{2} r^{2} - 2 M + d\phi^{2} r^{3} + dr^{2} r + r\right)}}{2 M - r}\right ]$$

Wir benutzen die erste Lösung (positive Lösung) für unsere Anfangsbedingungen $\left. \frac{dt}{d\tau} \right|_{\tau=0}$. Die Anfangsbedingungen lauten somit:

t0=0.0
phi0=0.0
dr0=0.0

#Grüne Kurve
Ar0=12.0
Adphi0=4.0/Ar0**2
Adt0=float(LoesEq_ds[0].subs({(r,Ar0),(dr,dr0),(dphi,Adphi0),(M,1)}))

#Rote Kurve
Br0=9.25
Bdphi0=3.7/Br0**2
Bdt0=float(LoesEq_ds[0].subs({(r,Br0),(dr,dr0),(dphi,Bdphi0),(M,1)}))

#Schwarze Kurve (ISCO)
Cr0=6.0
Cdphi0=float(2*sqrt(3)/Cr0**2)
Cdt0=float(LoesEq_ds[0].subs({(r,Cr0),(dr,dr0),(dphi,Cdphi0),(M,1)}))

tauval = np.linspace(0, 800, 10001)

Ainitialval = np.array([t0,Ar0,phi0,Adt0,dr0,Adphi0])
Binitialval = np.array([t0,Br0,phi0,Bdt0,dr0,Bdphi0])
Cinitialval = np.array([t0,Cr0,phi0,Cdt0,dr0,Cdphi0])

ALoes = integrate.odeint(DGLsys, Ainitialval, tauval)
BLoes = integrate.odeint(DGLsys, Binitialval, tauval)
CLoes = integrate.odeint(DGLsys, Cinitialval, tauval)

params = {
    'figure.figsize'    : [8,8],
#    'text.usetex'       : True,
    'axes.titlesize' : 14,
    'axes.labelsize' : 16,  
    'xtick.labelsize' : 14 ,
    'ytick.labelsize' : 14 
}
matplotlib.rcParams.update(params) 

import matplotlib.gridspec as gridspec
import matplotlib.animation as animation
from IPython.display import HTML

rval = np.linspace(2.2,30, 1001)
step=100
anzframes=100
xylim=15
fig = plt.figure(figsize=(15,7))
gs = gridspec.GridSpec(1, 2, width_ratios=[1,1], wspace=0.21)
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1])

def init():
    ax1.scatter(ALoes[0, 1]*np.cos(ALoes[0, 2]),ALoes[0, 1]*np.sin(ALoes[0, 2]),c="green", s=60, marker='o');
    ax1.scatter(BLoes[0, 1]*np.cos(BLoes[0, 2]),BLoes[0, 1]*np.sin(BLoes[0, 2]),c="red", s=60, marker='o');
    ax1.scatter(CLoes[0, 1]*np.cos(CLoes[0, 2]),CLoes[0, 1]*np.sin(CLoes[0, 2]),c="black", s=60, marker='o');
    BH=plt.Circle((0, 0), 2, color='black')
    ax1.add_patch(BH);
    ax2.plot(rval,VeffRezl(rval,4),c="green", linewidth=1.5, linestyle='-');
    ax2.scatter(Ar0,VeffRezl(Ar0,4),c="green", s=80, marker='o');
    ax2.plot(rval,VeffRezl(rval,3.7),c="red", linewidth=1.5, linestyle='-');
    ax2.scatter(Br0,VeffRezl(Br0,3.7),c="red", s=80, marker='o');
    ax2.plot(rval,VeffRezl(rval,float(2*sqrt(3))),c="black", linewidth=1.5, linestyle='-');
    ax2.scatter(Cr0,VeffRezl(Cr0,float(2*sqrt(3))),c="black", s=80, marker='o');
    ax1.set_ylabel(r"$\rm y \,[km]$")
    ax1.set_xlabel(r"$\rm x \,[km]$")
    ax2.set_ylabel(r"$\rm V(r)$")
    ax2.set_xlabel(r"$\rm Radius \, r \,[km]$")
    ax1.set_xlim(-xylim, xylim)
    ax1.set_ylim(-xylim, xylim)
    ax2.set_xlim(2.3, xylim)
    ax2.set_ylim(0.9, 1.01)
    return fig,

def animate(i):
    ax1.cla() 
    BH=plt.Circle((0, 0), 2, color='black')
    ax1.add_patch(BH);
    ax1.plot(ALoes[:step*i, 1]*np.cos(ALoes[:step*i, 2]),ALoes[:step*i, 1]*np.sin(ALoes[:step*i, 2]),c="green", linewidth=0.5, linestyle='-');
    ax1.scatter(ALoes[step*i, 1]*np.cos(ALoes[step*i, 2]),ALoes[step*i, 1]*np.sin(ALoes[step*i, 2]),c="green", s=60, marker='o');
    ax1.plot(BLoes[:step*i, 1]*np.cos(BLoes[:step*i, 2]),BLoes[:step*i, 1]*np.sin(BLoes[:step*i, 2]),c="red", linewidth=0.5, linestyle='-');
    ax1.scatter(BLoes[step*i, 1]*np.cos(BLoes[step*i, 2]),BLoes[step*i, 1]*np.sin(BLoes[step*i, 2]),c="red", s=60, marker='o');
    ax1.plot(CLoes[:step*i, 1]*np.cos(CLoes[:step*i, 2]),CLoes[:step*i, 1]*np.sin(CLoes[:step*i, 2]),c="black", linewidth=0.5, linestyle='-');
    ax1.scatter(CLoes[step*i, 1]*np.cos(CLoes[step*i, 2]),CLoes[step*i, 1]*np.sin(CLoes[step*i, 2]),c="black", s=60, marker='o');   
    ax1.set_ylabel(r"$\rm y \,[km]$")
    ax1.set_xlabel(r"$\rm x \,[km]$")
    ax2.set_ylabel(r"$\rm V(r)$")
    ax2.set_xlabel(r"$\rm Radius \, r \,[km]$")
    ax1.set_xlim(-xylim, xylim)
    ax1.set_ylim(-xylim, xylim)
    ax2.set_xlim(2.3, xylim)
    ax2.set_ylim(0.9, 1.01)
    return fig,

ani = animation.FuncAnimation(fig,animate,init_func=init,frames=anzframes,interval=100)
plt.close(ani._fig)
HTML(ani.to_html5_video())

b) Bewegung von Licht (Photonen) um ein schwarzes Loch in der Ebene (Die Photonensphäre)

In der Literatur wird die Bewegung von Licht um ein schwarzes Loch mittels eines definierten, effektiven Potentials illustriert:

1. General relativity : An introduction for physicists by M. P. Hobson, G. P. Efstathiou and A. N. Lasenby

2. Gravity : An introduction to Einstein's general relativity by James B. Hartle

3. Allgemeine Relativitätstheorie by Torsten Fließbach

4. Relativistic hydrodynamics by Luciano Rezzolla and Olindo Zanotti

Dieses Potential hängt allein von der Masse des schwarzen Lochs ab:

$$ \begin{equation} V_{\rm Photon} = \frac{1-2M/r}{r^2} \end{equation} $$

Die folgende Abbildung zeigt das effektive Potential als Funktion des Radius bei einer Masse von M=1:

VeffPhoton=(1-2*M/r)/r**2
VeffPhotonl=lambdify((r),VeffPhoton.subs({(M,1)}))
rval = np.linspace(1.8, 10, 1001)

plt.cla()
plt.xlim(1.8, 10)
plt.ylim(0, 0.04)
plt.ylabel(r"$\rm V_{\rm Photon}(r)$")
plt.xlabel(r"$\rm Radius \, r \,[km]$")
plt.plot(rval,VeffPhotonl(rval), linewidth=1, color="black");
plt.scatter(3,VeffPhotonl(3),c="orange", s=80, marker='o');

Wir erkennen, dass keine stabilen gebundenen Photonen-Trajektorien gibt (fehlendes Minimum beim oberen Potential). Es existiert jedoch eine quasi-stabile gebundenen Photonen-Trajektorie bei r=3 km, die sogenannte "Photonensphäre" (roter Punkt in der oberen Abbildung, allgemein bei $r_{PhSph}=3M$). Wir wollen uns im Folgenden diesen Photonenorbit näher betrachten und in einer Animation darstellen. Wir benutzen dazu wieder die Geodätengleichung, berücksichtigen jedoch bei den Anfangsbedingungen, dass das infinitesimale Weglängenelement bei Photonen Null ist ($ds^2=0$).Für die Photonensphäre gelten somit die folgenden Anfangswerte:

EqPh_ds=Eq(ds2.subs({(theta,pi/2),(dtheta,0)}),0)
LoesEqPh_ds=solve(EqPh_ds,dphi)
LoesEqPh_ds

$$\left [ - \sqrt{\frac{\left(- 2 M dt + dr r + dt r\right) \left(2 M dt + dr r - dt r\right)}{r^{3} \left(2 M - r\right)}}, \quad \sqrt{\frac{\left(- 2 M dt + dr r + dt r\right) \left(2 M dt + dr r - dt r\right)}{r^{3} \left(2 M - r\right)}}\right ]$$

Wir stellen wieder in einer Animation die Kreisbahn des Photons dar:

t0=0.0
phi0=0.0
dr0=0.0
#Photonensphaere
Dr0=3.0
Ddt0=1.0
Ddphi0=float(LoesEqPh_ds[1].subs({(r,3),(dr,0),(dt,1),(M,1)}))

Dinitialval = np.array([t0,Dr0,phi0,Ddt0,dr0,Ddphi0])
DLoes = integrate.odeint(DGLsys, Dinitialval, tauval)

step=15
anzframes=100
xylim=5
fig = plt.figure(figsize=(15,7))
gs = gridspec.GridSpec(1, 2, width_ratios=[1,1], wspace=0.21)
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1])

def init():
    ax1.scatter(DLoes[0, 1]*np.cos(DLoes[0, 2]),DLoes[0, 1]*np.sin(DLoes[0, 2]),c="red", s=60, marker='o');
    BH=plt.Circle((0, 0), 2, color='black')
    ax1.add_patch(BH);
    ax2.plot(rval,VeffPhotonl(rval), linewidth=1, color="black");
    ax2.scatter(3,VeffPhotonl(3),c="orange", s=80, marker='o');
    ax1.set_ylabel(r"$\rm y \,[km]$")
    ax1.set_xlabel(r"$\rm x \,[km]$")
    ax2.set_ylabel(r"$\rm V(r)$")
    ax2.set_xlabel(r"$\rm Radius \, r \,[km]$")
    ax1.set_xlim(-xylim, xylim)
    ax1.set_ylim(-xylim, xylim)
    ax2.set_xlim(1.8, 10)
    ax2.set_ylim(0, 0.04)
    return fig,

def animate(i):
    ax1.cla() 
    BH=plt.Circle((0, 0), 2, color='black')
    ax1.add_patch(BH);
    ax1.plot(DLoes[:step*i, 1]*np.cos(DLoes[:step*i, 2]),DLoes[:step*i, 1]*np.sin(DLoes[:step*i, 2]),c="orange", linewidth=0.5, linestyle='-');
    ax1.scatter(DLoes[step*i, 1]*np.cos(DLoes[step*i, 2]),DLoes[step*i, 1]*np.sin(DLoes[step*i, 2]),c="orange", s=60, marker='o');
    ax1.set_ylabel(r"$\rm y \,[km]$")
    ax1.set_xlabel(r"$\rm x \,[km]$")
    ax2.set_ylabel(r"$\rm V(r)$")
    ax2.set_xlabel(r"$\rm Radius \, r \,[km]$")
    ax1.set_xlim(-xylim, xylim)
    ax1.set_ylim(-xylim, xylim)
    ax2.set_xlim(1.8, 10)
    ax2.set_ylim(0, 0.04)
    return fig,

ani = animation.FuncAnimation(fig,animate,init_func=init,frames=anzframes,interval=100)
plt.close(ani._fig)
HTML(ani.to_html5_video())

Wir betrachten uns die drei kreisförmigen Bahnbewegungen der Probekörper und die oben dargestellte Trajektorie des masselosen Lichtteilchens auf der Photonensphäre in einem eingebetteten Diagramm der räumlichen Hypersphäre der Mannigfaltigkeit der Schwarzschildmetrik:

from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm

x,y = symbols('x,y')
zBett=sqrt(8*M*(r-2*M))
R = np.linspace(2.00001, 20, 300)
PHI = np.linspace(0, 2*np.pi, 300)
R, PHI = np.meshgrid(R, PHI)
X = R*np.cos(PHI)
Y = R*np.sin(PHI)
zlam = lambdify((x,y), zBett.subs({(M,1),(r,sqrt(x**2+y**2))}))
ZF = np.array(list(np.real(zlam(X,Y))))
rval = np.linspace(2.2,30, 1001)

step=20
anzframes=150
xylim=18
zlim=15
fig = plt.figure(figsize=(14,6))
gs = gridspec.GridSpec(1, 2, width_ratios=[2,1], wspace=0.21)
ax1 = plt.subplot(gs[0],projection='3d')
ax2 = plt.subplot(gs[1])
ax1.xaxis.pane.fill = False
ax1.yaxis.pane.fill = False
ax1.zaxis.pane.fill = False
ax1.xaxis.pane.set_edgecolor('white')
ax1.yaxis.pane.set_edgecolor('white')
ax1.zaxis.pane.set_edgecolor('white')

def init():
    ax1.plot_surface(X, Y, ZF, cmap=cm.nipy_spectral, linewidth=0, alpha=0.5)
    ax1.scatter(ALoes[0, 1]*np.cos(ALoes[0, 2]),ALoes[0, 1]*np.sin(ALoes[0, 2]),np.real(zlam(ALoes[0, 1]*np.cos(ALoes[0, 2]),ALoes[0, 1]*np.sin(ALoes[0, 2]))),c="green", s=60, marker='o');
    ax1.contourf(X, Y, ZF,100, zdir='z', offset=0, cmap=cm.nipy_spectral,alpha=0.5)
    ax1.scatter(ALoes[0, 1]*np.cos(ALoes[0, 2]),ALoes[0, 1]*np.sin(ALoes[0, 2]),0.1,c="green", s=40, marker='o');
    ax1.scatter(BLoes[0, 1]*np.cos(BLoes[0, 2]),BLoes[0, 1]*np.sin(BLoes[0, 2]),np.real(zlam(BLoes[0, 1]*np.cos(BLoes[0, 2]),BLoes[0, 1]*np.sin(BLoes[0, 2]))),c="red", s=60, marker='o');
    ax1.scatter(BLoes[0, 1]*np.cos(BLoes[0, 2]),BLoes[0, 1]*np.sin(BLoes[0, 2]),0.1,c="red", s=40, marker='o');
    ax1.scatter(CLoes[0, 1]*np.cos(CLoes[0, 2]),CLoes[0, 1]*np.sin(CLoes[0, 2]),np.real(zlam(CLoes[0, 1]*np.cos(CLoes[0, 2]),CLoes[0, 1]*np.sin(CLoes[0, 2]))),c="black", s=60, marker='o');
    ax1.scatter(CLoes[0, 1]*np.cos(CLoes[0, 2]),CLoes[0, 1]*np.sin(CLoes[0, 2]),0.1,c="black", s=40, marker='o'); 
    ax1.scatter(DLoes[0, 1]*np.cos(DLoes[0, 2]),DLoes[0, 1]*np.sin(DLoes[0, 2]),np.real(zlam(DLoes[0, 1]*np.cos(DLoes[0, 2]),DLoes[0, 1]*np.sin(DLoes[0, 2]))),c="yellow", s=60, marker='o');
    ax1.scatter(DLoes[0, 1]*np.cos(DLoes[0, 2]),DLoes[0, 1]*np.sin(DLoes[0, 2]),0.1,c="yellow", s=40, marker='o'); 
    ax2.plot(rval,VeffRezl(rval,4),c="green", linewidth=1.5, linestyle='-');
    ax2.scatter(Ar0,VeffRezl(Ar0,4),c="green", s=80, marker='o');
    ax2.plot(rval,VeffRezl(rval,3.7),c="red", linewidth=1.5, linestyle='-');
    ax2.scatter(Br0,VeffRezl(Br0,3.7),c="red", s=80, marker='o');
    ax2.plot(rval,VeffRezl(rval,float(2*sqrt(3))),c="black", linewidth=1.5, linestyle='-');
    ax2.scatter(Cr0,VeffRezl(Cr0,float(2*sqrt(3))),c="black", s=80, marker='o');
    ax1.set_ylabel(r"$\rm y \,[km]$")
    ax1.set_xlabel(r"$\rm x \,[km]$")
    ax1.set_zlabel("Visualisierung z");
    ax2.set_ylabel(r"$\rm V(r)$")
    ax2.set_xlabel(r"$\rm Radius \, r \,[km]$")
    ax1.view_init(azim=35, elev=20)
    ax1.set_xlim(-xylim, xylim)
    ax1.set_ylim(-xylim, xylim)
    ax1.set_zlim(0, 15)
    ax2.set_xlim(2.3, xylim)
    ax2.set_ylim(0.9, 1.01)
    ax1.grid(False)
    return fig,

def animate(i):
    ax1.cla() 
    ax1.plot_surface(X, Y, ZF, cmap=cm.nipy_spectral, linewidth=0, alpha=0.15)
    ax1.plot(ALoes[:step*i, 1]*np.cos(ALoes[:step*i, 2]),ALoes[:step*i, 1]*np.sin(ALoes[:step*i, 2]),np.real(zlam(ALoes[:step*i, 1]*np.cos(ALoes[:step*i, 2]),ALoes[:step*i, 1]*np.sin(ALoes[:step*i, 2]))),c="green", linewidth=1, linestyle='-');
    ax1.scatter(ALoes[step*i, 1]*np.cos(ALoes[step*i, 2]),ALoes[step*i, 1]*np.sin(ALoes[step*i, 2]),np.real(zlam(ALoes[step*i, 1]*np.cos(ALoes[step*i, 2]),ALoes[step*i, 1]*np.sin(ALoes[step*i, 2]))),c="green", s=60, marker='o');
    ax1.scatter(ALoes[step*i, 1]*np.cos(ALoes[step*i, 2]),ALoes[step*i, 1]*np.sin(ALoes[step*i, 2]),0.1,c="green", s=40, marker='o');
    ax1.plot(ALoes[:step*i, 1]*np.cos(ALoes[:step*i, 2]),ALoes[:step*i, 1]*np.sin(ALoes[:step*i, 2]),0.1,c="green", linewidth=2.5, linestyle='-');
    ax1.plot(BLoes[:step*i, 1]*np.cos(BLoes[:step*i, 2]),BLoes[:step*i, 1]*np.sin(BLoes[:step*i, 2]),np.real(zlam(BLoes[:step*i, 1]*np.cos(BLoes[:step*i, 2]),BLoes[:step*i, 1]*np.sin(BLoes[:step*i, 2]))),c="red", linewidth=1, linestyle='-');
    ax1.scatter(BLoes[step*i, 1]*np.cos(BLoes[step*i, 2]),BLoes[step*i, 1]*np.sin(BLoes[step*i, 2]),np.real(zlam(BLoes[step*i, 1]*np.cos(BLoes[step*i, 2]),BLoes[step*i, 1]*np.sin(BLoes[step*i, 2]))),c="red", s=60, marker='o');
    ax1.scatter(BLoes[step*i, 1]*np.cos(BLoes[step*i, 2]),BLoes[step*i, 1]*np.sin(BLoes[step*i, 2]),0.1,c="red", s=40, marker='o');
    ax1.plot(BLoes[:step*i, 1]*np.cos(BLoes[:step*i, 2]),BLoes[:step*i, 1]*np.sin(BLoes[:step*i, 2]),0.1,c="red", linewidth=2.5, linestyle='-');
    ax1.plot(CLoes[:step*i, 1]*np.cos(CLoes[:step*i, 2]),CLoes[:step*i, 1]*np.sin(CLoes[:step*i, 2]),np.real(zlam(CLoes[:step*i, 1]*np.cos(CLoes[:step*i, 2]),CLoes[:step*i, 1]*np.sin(CLoes[:step*i, 2]))),c="black", linewidth=1, linestyle='-');
    ax1.scatter(CLoes[step*i, 1]*np.cos(CLoes[step*i, 2]),CLoes[step*i, 1]*np.sin(CLoes[step*i, 2]),np.real(zlam(CLoes[step*i, 1]*np.cos(CLoes[step*i, 2]),CLoes[step*i, 1]*np.sin(CLoes[step*i, 2]))),c="black", s=60, marker='o');
    ax1.scatter(CLoes[step*i, 1]*np.cos(CLoes[step*i, 2]),CLoes[step*i, 1]*np.sin(CLoes[step*i, 2]),0.1,c="black", s=40, marker='o');
    ax1.plot(CLoes[:step*i, 1]*np.cos(CLoes[:step*i, 2]),CLoes[:step*i, 1]*np.sin(CLoes[:step*i, 2]),0.1,c="black", linewidth=2.5, linestyle='-');
    ax1.plot(DLoes[:step*i, 1]*np.cos(DLoes[:step*i, 2]),DLoes[:step*i, 1]*np.sin(DLoes[:step*i, 2]),np.real(zlam(DLoes[:step*i, 1]*np.cos(DLoes[:step*i, 2]),DLoes[:step*i, 1]*np.sin(DLoes[:step*i, 2]))),c="orange", linewidth=1, linestyle='-');
    ax1.scatter(DLoes[step*i, 1]*np.cos(DLoes[step*i, 2]),DLoes[step*i, 1]*np.sin(DLoes[step*i, 2]),np.real(zlam(DLoes[step*i, 1]*np.cos(DLoes[step*i, 2]),DLoes[step*i, 1]*np.sin(DLoes[step*i, 2]))),c="orange", s=60, marker='o');
    ax1.scatter(DLoes[step*i, 1]*np.cos(DLoes[step*i, 2]),DLoes[step*i, 1]*np.sin(DLoes[step*i, 2]),0.1,c="orange", s=40, marker='o');
    ax1.plot(DLoes[:step*i, 1]*np.cos(DLoes[:step*i, 2]),DLoes[:step*i, 1]*np.sin(DLoes[:step*i, 2]),0.1,c="orange", linewidth=2.5, linestyle='-');
    ax1.contourf(X, Y, ZF,100, zdir='z', offset=0, cmap=cm.nipy_spectral,alpha=0.15)
    ax1.set_ylabel(r"$\rm y \,[km]$")
    ax1.set_xlabel(r"$\rm x \,[km]$")
    ax1.set_zlabel("Einbettung z");
    ax1.view_init(azim=35, elev=20)
    ax1.set_xlim(-xylim, xylim)
    ax1.set_ylim(-xylim, xylim)
    ax1.set_zlim(0, 15)  
    ax1.grid(False)
    return fig,

ani = animation.FuncAnimation(fig,animate,init_func=init,frames=anzframes,interval=100)
plt.close(ani._fig)
HTML(ani.to_html5_video())

Die Visualisierung kann man auch wieder mittels der "Plotly Python Open Source Graphing Library" erstellen ( siehe https://plotly.com/python/ ).

import plotly.graph_objects as go

trichter = go.Surface(x=X, y=Y, z=ZF, colorscale="spectral_r", showscale=True, name="Embedded Surface", hoverinfo='none')

Loes=ALoes
circular1_orbit =go.Scatter3d(
    x=Loes[:, 1]*np.cos(Loes[:, 2]),
    y=Loes[:, 1]*np.sin(Loes[:, 2]),
    z=np.real(zlam(Loes[:, 1]*np.cos(Loes[:, 2]),Loes[:, 1]*np.sin(Loes[:, 2]))),
    name="Circular orbit I",
    mode='lines',
    line=dict(color='green', width=8),
    hoverinfo='none'
    )

Loes=BLoes
circular2_orbit =go.Scatter3d(
    x=Loes[:, 1]*np.cos(Loes[:, 2]),
    y=Loes[:, 1]*np.sin(Loes[:, 2]),
    z=np.real(zlam(Loes[:, 1]*np.cos(Loes[:, 2]),Loes[:, 1]*np.sin(Loes[:, 2]))),
    name="Circular orbit II",
    mode='lines',
    line=dict(color='red', width=8),
    hoverinfo='none'
    )

Loes=CLoes
ISCO_orbit =go.Scatter3d(
    x=Loes[:, 1]*np.cos(Loes[:, 2]),
    y=Loes[:, 1]*np.sin(Loes[:, 2]),
    z=np.real(zlam(Loes[:, 1]*np.cos(Loes[:, 2]),Loes[:, 1]*np.sin(Loes[:, 2]))),
    name="ISCO",
    mode='lines',
    line=dict(color='black', width=8),
    hoverinfo='none'
    )

Loes=DLoes
Photon_orbit =go.Scatter3d(
    x=Loes[:, 1]*np.cos(Loes[:, 2]),
    y=Loes[:, 1]*np.sin(Loes[:, 2]),
    z=np.real(zlam(Loes[:, 1]*np.cos(Loes[:, 2]),Loes[:, 1]*np.sin(Loes[:, 2]))),
    name="Photon orbit",
    mode='lines',
    line=dict(color='orange', width=8),
    hoverinfo='none'
    )

layout = go.Layout(autosize=True,
                  width=700, height=700,
                  margin=dict(l=15, r=50, b=65, t=90),
                  scene_aspectmode='cube',
                  scene = dict(
                    xaxis_showbackground=False,
                    yaxis_showbackground=False,
                    zaxis_showbackground=False,
                    xaxis_gridcolor="lightgrey",
                    yaxis_gridcolor="lightgrey",
                    zaxis_gridcolor="lightgrey",
                    xaxis_range=[-22,22],
                    yaxis_range=[-22,22],
                    zaxis_range=[0,15],
                    xaxis_title='x',
                    yaxis_title='y',
                    zaxis_title='f'),
                  legend=dict(
                    yanchor="top",
                    y=0.99,
                    xanchor="left",
                    x=0.01)
            )

data=[trichter,circular1_orbit,circular2_orbit,ISCO_orbit,Photon_orbit]
fig=go.Figure(data=data, layout=layout)

fig.show()