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 (Wintersemester 2026/27)

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

Frankfurt am Main 28.09.2026

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

Die Bewegung eines Probekörpers (eines Teilchens mit vernachlässigbarer Masse, das die Raumzeit selbst nicht deformiert) um ein nicht-rotierendes, ungeladenes schwarzes Loch wird durch die Schwarzschild-Metrik der Allgemeinen Relativitätstheorie bestimmt. In der Vorlesung Allgemeine Relativitätstheorie mit dem Computer behandelten wir die Eigenschaften der Schwarzschild-Metrik. Da das System sphärisch symmetrisch ist, bleibt der Drehimpuls erhalten. Die Bewegung eines Teilchens findet daher immer in einer konstanten Bahnebene statt, die man äquatorial bei $\theta = \pi/2$ ansetzen kann. In dem Jupyter Notebook Teil I: Radialer Fall eines Probekörpers in ein schwarzes Loch betrachteten wir zunächst den einfachsten Fall eines radialen Einfallens eines Testteilchens, wohingegen wir im Jupyter Notebook Teil II: Klassifizierung unterschiedlicher Bahnbewegungen betrachteten wir allgemeine geodätische Bewegungen um ein nicht-rotierendes, ungeladenes schwarzes Loch und klassifizierten diese mittels eines effektiven Potenzials.

Teil III: Der ISCO und die Photonensphäre

In diesem Jupyter Notebook werden wir im Folgenden das Konzept der Innersten stabilen kreisförmige Bahnbewegung (Innermost Stable Circular Orbit (ISCO)) vorstellen und den Begriff der Photonensphäre eines nicht-rotierenden schwarzen Loches diskutieren.

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

Der ISCO markiert die innerste stabile Kreisbahn, auf der sich ein Körper mit Masse dauerhaft bewegen kann. Der ISCO befindet sich bei drei Schwarzschild-Radien $r_{\rm ISCO} = 3 \, R_s = 6 \, M $. Für Radien $r > r_{\rm ISCO}$ existieren stabile kreisförmige Bahnen. Fällt ein Probekörper unter diesen Radius, verliert das effektive Potenzial sein lokales Minimum. Es gibt dort keine stabilen Kreisbahnen mehr. Erreicht Materie (z. B. in einer Akkretionsscheibe um ein schwarzes Loch) den ISCO, beginnt sie ohne Zufuhr weiterer Energie spiralförmig und rapide in das schwarze Loch zu stürzen.

Geodätische Bewegung eines Probekörpers

Um die Geodätengleichung

$$ \frac{d^2 x^\mu}{d\tau^2} + \Gamma^\mu_{\nu \rho} \frac{d x^\nu}{d\tau} \frac{d x^\rho}{d\tau} ~=~ 0 $$

zu berechnen, definieren zunächst den ersten Term der Gleichung $\frac{d^2 x^\mu}{d\tau^2}$:

Wir beginnen mit der Berechnung des ISCOs 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])
chr = ChristoffelSymbols.from_metric(g)

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)

und nun den zweiten Term $\Gamma^\mu_{\nu \rho} \frac{d x^\nu}{d\tau} \frac{d x^\rho}{d\tau}$:

Geod2T = tensorproduct(chr.tensor(),diff(xmu.tensor(),tau),diff(xmu.tensor(),tau))
Geod2 = tensorcontraction(Geod2T, (1, 3),(2,4))

Die Geodätengleichung ist ein System von nichtlinearen gekoppelten Differenzialgleichungen zweiter Ordnung. Wir definieren die entsprechenden vier Gleichungen:

Geod = Geod1 + Geod2
GeodEq0 = Eq(Geod[0],0)
GeodEq1 = Eq(Geod[1],0)
GeodEq2 = Eq(Geod[2],0)
GeodEq3 = Eq(Geod[3],0)

GeodEq0

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

GeodEq1

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

GeodEq2

$\displaystyle - \sin{\left(\theta \right)} \cos{\left(\theta \right)} \left(\frac{d}{d \tau} \phi{\left(\tau \right)}\right)^{2} + \frac{d^{2}}{d \tau^{2}} \theta{\left(\tau \right)} + \frac{2 \frac{d}{d \tau} r{\left(\tau \right)} \frac{d}{d \tau} \theta{\left(\tau \right)}}{r} = 0$

GeodEq3

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

Wir lassen im Folgenden 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)})
GeodEq1 = GeodEq1.subs({(theta,pi/2),(diff(x2,tau),0)})
GeodEq3 = GeodEq3.subs({(theta,pi/2),(diff(x2,tau),0)})

GeodEq0_a = Eq(diff(x0,tau,tau),solve(GeodEq0,diff(x0,tau,tau))[0])
GeodEq1_a = Eq(diff(x1,tau,tau),solve(GeodEq1,diff(x1,tau,tau))[0])
GeodEq3_a = Eq(diff(x3,tau,tau),solve(GeodEq3,diff(x3,tau,tau))[0])

GeodEq0_a

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

GeodEq1_a

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

GeodEq3_a

$\displaystyle \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}$

Um dieses System von drei gekoppelten Differentialgleichungen zweiter Ordnung mit Python nummerisch lösen zu können, müssen wir es zunächst in ein System von sechs Differenzialgleichungen erster Ordnung umschreiben. Details zum numerischen Lösung von Differenzialgleichungen finden Sie in der Vorlesung Einführung in die Programmierung für Studierende der Physik in den Unterpunkten Differenzialgleichungen: Numerische Lösung von Anfangswertproblemen und Systeme von gekoppelten Differenzialgleichungen und Differenzialgleichungen zweiter Ordnung.

Wir definieren dazu formal

$$ t(\tau) = y_1(\tau) \, , \quad r(\tau) = y_2(\tau), \quad \phi(\tau) = y_3(\tau) \, , \quad \frac{d t(\tau)}{d\tau} = y_4(\tau) \, , \quad \frac{d r(\tau)}{d\tau} = y_5(\tau), \quad \frac{d \phi(\tau)}{d\tau} = y_6 $$

und schreiben unser DGL-System wie folgt um:

y1, y2, y3, y4, y5, y6 = symbols('y_1, y_2, y_3, y_4, y_5, y_6')

SysDGL_0 = Eq(diff(x0,tau,tau),GeodEq0_a.rhs.subs({(r,y2),(diff(x0,tau),y4),(diff(x1,tau),y5)}))
SysDGL_0

$\displaystyle \frac{d^{2}}{d \tau^{2}} t{\left(\tau \right)} = \frac{2 M y_{4} y_{5}}{y_{2} \left(2 M - y_{2}\right)}$

SysDGL_1 = Eq(diff(x1,tau,tau),GeodEq1_a.rhs.subs({(r,y2),(diff(x0,tau),y4),(diff(x1,tau),y5),(diff(x3,tau),y6)}))
SysDGL_1

$\displaystyle \frac{d^{2}}{d \tau^{2}} r{\left(\tau \right)} = - \frac{2 M^{2} y_{5}^{2}}{y_{2} \left(2 M - y_{2}\right)^{2}} + \frac{2 M^{2} y_{4}^{2}}{y_{2}^{3}} + \frac{M y_{5}^{2}}{\left(2 M - y_{2}\right)^{2}} - \frac{M y_{4}^{2}}{y_{2}^{2}} + \frac{y_{6}^{2} \left(- 2 M + y_{2}\right)^{3}}{\left(2 M - y_{2}\right)^{2}}$

SysDGL_3 = Eq(diff(x1,tau,tau),GeodEq3_a.rhs.subs({(r,y2),(diff(x0,tau),y4),(diff(x1,tau),y5),(diff(x3,tau),y6)}))
SysDGL_3

$\displaystyle \frac{d^{2}}{d \tau^{2}} r{\left(\tau \right)} = - \frac{2 y_{5} y_{6}}{y_{2}}$

Wir werden später die DGLs mittels Python numerisch lösen und setzen im Folgenden $M=1$. Wir benutzen hierfür das Python-Modul SciPy, welches eine breite Kollektion von mathematischen Algorithmen und Funktionen bereitstellt. Im Speziellen werden wir in diesem Semester die Funktion 'solve_ivp(...)' verwenden, die sich im Untermodul 'scipy.integrate' befindet, welches Funktionen zum Lösen von gewöhnlichen Differenzialgleichungen bereitstellt. Zusätzlich berechnen wir die numerische Lösung auch mittels der älteren Methode "integrate.odeint()", die ebenfalls in dem Python-Modul "scipy" definiert ist.

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

Wir definieren das System der gekoppelten sechs Differenzialgleichungen:

set_M = 1
f_4 = lambdify((y1,y2,y3,y4,y5,y6), SysDGL_0.rhs.subs(M,set_M))
f_5 = lambdify((y1,y2,y3,y4,y5,y6), SysDGL_1.rhs.subs(M,set_M))
f_6 = lambdify((y1,y2,y3,y4,y5,y6), SysDGL_3.rhs.subs(M,set_M))

def DGLsys(vy, time):
    y1,y2,y3,y4,y5,y6 = vy
    dy1 = y4
    dy2 = y5
    dy3 = y6
    dy4 = f_4(y1,y2,y3,y4,y5,y6)
    dy5 = f_5(y1,y2,y3,y4,y5,y6)
    dy6 = f_6(y1,y2,y3,y4,y5,y6)
    return np.array([dy1,dy2,dy3,dy4,dy5,dy6])

def DGLsys_ivp(time, vy):
    y1,y2,y3,y4,y5,y6 = vy
    dy1 = y4
    dy2 = y5
    dy3 = y6
    dy4 = f_4(y1,y2,y3,y4,y5,y6)
    dy5 = f_5(y1,y2,y3,y4,y5,y6)
    dy6 = f_6(y1,y2,y3,y4,y5,y6)
    return np.array([dy1,dy2,dy3,dy4,dy5,dy6])

Bevor wir die Geodätengleichung numerisch lösen und die innerste stabile kreisförmige Bahnbewegung eines Probekörpers visualisieren können, benötigen wir die Anfangsbedingungen des ISCOSs.

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 Differenzialgleichungen 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 Potenzials illustriert werden. Dieses Potenzial 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 Potenzial erfolgt mittels der radialen, 2. Geodätengleichung. So definieren z.B. die Literaturangaben 1-3 (siehe Angaben oben) das effektive Potenzial 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 Pozential 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} $$

Wir betrachten nun kreisförmige gebundene Orbits, bei denen der Abstand zum schwarzen Loch (der Radiuswert $r$) während der gesamten Trajektorie konstant ist. Solche Trajektorien sind nur bei den Extrema des effektiven Potenzial $V(r,M,l)$ möglich, wobei stabile Orbits nur im lokalen Minimum möglich sind.

$$ \frac{d V}{dr} ~=~ 0 \,\, , \quad \frac{d^2 V}{dr^2} ~<~ 0 $$

Die Ableitungen berechnen sich zu:

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,set_M)}))
VeffNewtonl = lambdify((r,l),VeffNewton.subs({(M,set_M)}))
rval = np.linspace(2.2,30, 1001)

d_VeffRez = diff(VeffRez,r)
d_VeffRez

$\displaystyle \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}}$

Die Extrema von $V(r,M,l)$ liegen bei den folgenden Radiuswerten.

Extrem_V = solve(d_VeffRez,r)
Extrem_V

$\displaystyle \left[ \frac{l^{2}}{2 M} - \frac{l \sqrt{- 12 M^{2} + l^{2}}}{2 M}, \ \frac{l^{2}}{2 M} + \frac{l \sqrt{- 12 M^{2} + l^{2}}}{2 M}\right]$

Setzt man die freien Parameterwerte z.B. auf $l=4.5$ und $M=1$, so erhält man für die Extrema von $V(r,M=1,l=4.5)$ die folgenden Radiuswerte.

set_l = 4.5
Extrem_V_set = solve(d_VeffRez.subs({(M,set_M),(l,set_l)}),r)
Extrem_V_set

$\displaystyle \left[ 3.66236702264472, \ 16.5876329773553\right]$

Anhand der Werte für die zweite Ableitung sieht man welcher dieser Werte ein Maximum bzw. Minimum ist.

dd_VeffRez = diff(VeffRez,r,r)
dd_VeffRez

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

[dd_VeffRez.subs({(M,set_M),(l,set_l),(r,Extrem_V_set[0])}),dd_VeffRez.subs({(M,set_M),(l,set_l),(r,Extrem_V_set[1])})]

$\displaystyle \left[ -0.067312449891573, \ 0.000175703779112795\right]$

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

rcParams.update({
    'figure.figsize'    : [10,8],
    'axes.titlesize' : 14,
    'axes.labelsize' : 16,  
    'xtick.labelsize' : 14 ,
    'ytick.labelsize' : 14 
})

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

set_data = [
    (4.5, "yellow", r"$\rm l=4.5$", 1),
    (4, "green", r"$\rm l=4$", 1),
    (3.7, "red", r"$\rm l=3.7$", 1),
    (float(2 * sqrt(3)), "black", r"$\rm l=2 \, \sqrt{3}$", 0),
    (3, "blue", r"$\rm l=3$", None)  
]

for l_val, color, label_text, i in set_data:
    plt.plot(rval, VeffRezl(rval, l_val), linewidth=1, color=color, label=label_text)
    if i is not None:
        Extrem_V_set = solve(d_VeffRez.subs({(M, set_M), (l, l_val)}), r)
        x = Extrem_V_set[i]
        y = float(VeffRez.subs({(M, set_M), (l, l_val), (r, x)}))
        plt.scatter(x, y, c=color, s=60, marker="o")

plt.legend(loc="lower right", fontsize=18)
plt.show()

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 Potenzial der allgemeinen Relativitätstheorie $V(r,M=1,l) = V(r,l)$ mit dem effektiven Newtonschen Potenzial $V_{Newton}(r,l)$. Für große Abstände geht das effektive Potenzial in das Newtonsche über und es gilt (für das effektive Potenzial der Literaturangabe 4):

$$ \lim\limits_{r \to \infty}V(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 Potenzials 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,set_M)}))

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

set_data = [
    (4, "green", r"$\rm l=4$"),
    (3.7, "red", r"$\rm l=3.7$"),
    (float(2 * sqrt(3)), "black", r"$\rm l=2 \,\sqrt{3}$"),
    (3, "blue", r"$\rm l=3$"),
    (2.5, "orange", r"$\rm l=2.5$"),
    (2, "brown", r"$\rm l=2$"),
    (1, "grey", r"$\rm l=1$")
]

for l_val, color, label_text in set_data:
    plt.plot(rval, VeffNewtonl(rval, l_val), linewidth=1, color=color, label=label_text)
    rmin_set = solve(Eq(diff(VeffNewton.subs({(M, set_M), (l, l_val)}), r), 0), r)
    if rmin_set:
        rmin = rmin_set[0]
        plt.scatter(rmin, VeffNewtonl(rmin, l_val), c=color, s=60, marker='o')

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

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 Potenzialbarriere für kleine Radien gibt.

Berechnung des ISCOs

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 Potenzial 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 Potenzials befindet. Für kreisförmige Bahnbewegungen muss somit gelten: $\frac{dV(r,M,l)}{dr}=0$. Dies führt zu folgender Gleichung:

EqCircular = Eq(d_VeffRez,0)
EqCircular

$\displaystyle \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()

$\displaystyle \left( \frac{l \left(l - \sqrt{- 12 M^{2} + l^{2}}\right)}{2 M}, \ \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:

EqSattel = Eq(dd_VeffRez,0)
EqSattel

$\displaystyle \frac{- \frac{2 M l^{2}}{r^{2} \sqrt{- \frac{2 M}{r} + 1} \sqrt{\frac{l^{2}}{r^{2}} + 1}} + \frac{M \sqrt{\frac{l^{2}}{r^{2}} + 1} \left(\frac{M}{r \left(\frac{2 M}{r} - 1\right)} - 2\right)}{\sqrt{- \frac{2 M}{r} + 1}} - \frac{l^{2} \sqrt{- \frac{2 M}{r} + 1} \left(\frac{l^{2}}{r^{2} \left(\frac{l^{2}}{r^{2}} + 1\right)} - 3\right)}{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

$\displaystyle \frac{16 M^{4} \sqrt{\frac{- 4 M^{2} + l^{2} + l \sqrt{- 12 M^{2} + l^{2}}}{l \left(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) - \left(4 M^{2} + \left(l + \sqrt{- 12 M^{2} + l^{2}}\right)^{2}\right)^{2} \left(- 3 M^{2} + l^{2} + l \sqrt{- 12 M^{2} + l^{2}}\right) + \left(8 M^{2} + 3 \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)^{2}\right)}{l^{2} \sqrt{\frac{4 M^{2} + \left(l + \sqrt{- 12 M^{2} + l^{2}}\right)^{2}}{\left(l + \sqrt{- 12 M^{2} + l^{2}}\right)^{2}}} \left(4 M^{2} + \left(l + \sqrt{- 12 M^{2} + l^{2}}\right)^{2}\right) \left(l + \sqrt{- 12 M^{2} + l^{2}}\right)^{4} \left(4 M^{2} - l^{2} - l \sqrt{- 12 M^{2} + l^{2}}\right)^{2}} = 0$

Lösen wir diese Gleichung nach l auf erhalten wir zwei reellwertige 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

$\displaystyle \left[ - 2 \sqrt{3} M, \ 2 \sqrt{3} M, \ 2 i 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])})

$\displaystyle 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} $$

Numerische Simulation

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

$\displaystyle - 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

$\displaystyle \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}, \ \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,set_M)}))

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

#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,set_M)}))

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)
CLoes_ivp = integrate.solve_ivp(DGLsys_ivp, [tauval[0], tauval[-1]], Cinitialval, t_eval=tauval, rtol=10**(-10), atol=10**(-10))

rcParams.update({
    'figure.figsize'    : [8,8],
    'axes.titlesize' : 14,
    'axes.labelsize' : 16,  
    'xtick.labelsize' : 14 ,
    'ytick.labelsize' : 14 
})

ax = plt.gca()
ax.cla() 
ax.set_ylabel(r"$\rm y \,[km]$")
ax.set_xlabel(r"$\rm x \,[km]$")
ax.plot(CLoes[:800, 1]*np.cos(CLoes[:800, 2]),CLoes[:800, 1]*np.sin(CLoes[:800, 2]),c="blue", linewidth=1.5, linestyle='-');
ax.plot(CLoes_ivp.y[1][:800]*np.cos(CLoes_ivp.y[2][:800]),CLoes_ivp.y[1][:800]*np.sin(CLoes_ivp.y[2][:800]),c="red", linewidth=2, linestyle=':');
BH = plt.Circle((0, 0), 2, color='black')
ax.add_patch(BH);

Das obere Diagramm zeigt die Trajektorie des Probekörpers auf dem ISCO, wobei die blaue Kurve die Ergebnisse der odeint-Methode und die gepunktete rote Kurve die solve_ivp-Simulation zeigt. Beide Simulationen stimmen gut miteinander überein.

Wir stellen nun alle berechneten Trajektorien in einer Animation dar.

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=(11, 5))
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(fig)
HTML(ani.to_html5_video())

bzw. ohne Verwendung von ax1.cla() in jedem Animationsschritt und mit einer Trennung von Animationsobjekten und statischen Elementen.

rval = np.linspace(2.2, 30, 1001)
step = 100
anzframes = 100
xylim = 15

DATA_CONFIG = [
    {"label": "A", "data": ALoes, "r0": Ar0, "param": 4.0, "color": "green"},
    {"label": "B", "data": BLoes, "r0": Br0, "param": 3.7, "color": "red"},
    {"label": "C", "data": CLoes, "r0": Cr0, "param": float(2 * np.sqrt(3)), "color": "black"}
]

fig = plt.figure(figsize=(11, 5))
gs = gridspec.GridSpec(1, 2, width_ratios=[1, 1], wspace=0.21)
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1])

for cfg in DATA_CONFIG:
    ax2.plot(rval, VeffRezl(rval, cfg["param"]), c=cfg["color"], linewidth=1.5, linestyle='-')
    ax2.scatter(cfg["r0"], VeffRezl(cfg["r0"], cfg["param"]), c=cfg["color"], 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)

bh_patch = plt.Circle((0, 0), 2, color='black')
ax1.add_patch(bh_patch)

lines = {}
scatters = {}
for cfg in DATA_CONFIG:
    lbl = cfg["label"]
    lines[lbl], = ax1.plot([], [], c=cfg["color"], linewidth=0.5, linestyle='-')
    scatters[lbl] = ax1.scatter([], [], c=cfg["color"], s=60, marker='o')

def init():
    for cfg in DATA_CONFIG:
        lbl = cfg["label"]
        d = cfg["data"]
        x0 = d[0, 1] * np.cos(d[0, 2])
        y0 = d[0, 1] * np.sin(d[0, 2])
        lines[lbl].set_data([], [])
        scatters[lbl].set_offsets(np.c_[x0, y0])  
    return [lines[lbl] for lbl in lines] + [scatters[lbl] for lbl in scatters]


def animate(i):
    idx = step * i
    
    for cfg in DATA_CONFIG:
        lbl = cfg["label"]
        d = cfg["data"]
        x_path = d[:idx, 1] * np.cos(d[:idx, 2])
        y_path = d[:idx, 1] * np.sin(d[:idx, 2])
        lines[lbl].set_data(x_path, y_path)
        x_head = d[idx, 1] * np.cos(d[idx, 2])
        y_head = d[idx, 1] * np.sin(d[idx, 2])
        scatters[lbl].set_offsets(np.c_[x_head, y_head])    
    return [lines[lbl] for lbl in lines] + [scatters[lbl] for lbl in scatters]

ani = animation.FuncAnimation(
    fig, animate, init_func=init, frames=anzframes, interval=100, blit=True
)
plt.close(fig)

HTML(ani.to_html5_video())

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

Die Photonensphäre ist der Bereich, in dem masselose Teilchen (Photonen) auf einer Kreisbahn wandern können. Sie liegt exakt bei $r_{PhSph}=3 \, M$. Die Bahn der Photonen ist jedoch instabil. Winzigste Störungen führen dazu, dass das Photon entweder ins schwarze Loch stürzt oder ins Unendliche entweicht. Da Licht hier kreisen kann, markiert die Photonensphäre die theoretische Grenze für optische Phänomene (wie den sogenannten „Photonenring“ bei astronomischen Aufnahmen).

In der Literatur wird die Bewegung von Licht um ein schwarzes Loch mittels eines definierten, effektiven Potenzials 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 Potenzial 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 Potenzial als Funktion des Radius bei einer Masse von M=1:

VeffPhoton = (1-2*M/r)/r**2
VeffPhotonl = lambdify((r),VeffPhoton.subs({(M,set_M)}))
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 es keine stabilen gebundenen Photonen-Trajektorien gibt (fehlendes Minimum beim oberen Potenzial). Es existiert jedoch eine quasi-stabile gebundenen Photonen-Trajektorie bei r=3 km, die sogenannte "Photonensphäre" (oranger Punkt in der oberen Abbildung, allgemein bei $r_{PhSph}=3 \, M$). 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 für masselose 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

$\displaystyle \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)}}, \ \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 die Kreisbahn des Photons in einer Animation 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,set_M)}))

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=(11, 5))
gs = gridspec.GridSpec(1, 2, width_ratios=[1,1], wspace=0.25)
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(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=(11, 5))
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(fig)
HTML(ani.to_html5_video())

bzw. ohne Verwendung von ax1.cla() in jedem Animationsschritt und mit einer Trennung von Animationsobjekten und statischen Elementen.

step, anzframes, xylim, zlim = 20, 150, 18, 15

DATA_CONFIG = [
    {"data": ALoes, "r0": Ar0, "param": 4.0, "color": "green"},
    {"data": BLoes, "r0": Br0, "param": 3.7, "color": "red"},
    {"data": CLoes, "r0": Cr0, "param": float(2 * np.sqrt(3)), "color": "black"},
    {"data": DLoes, "r0": None, "param": None, "color": "orange"}
]

fig = plt.figure(figsize=(11, 5))
gs = gridspec.GridSpec(1, 2, width_ratios=[2, 1], wspace=0.21)
ax1 = plt.subplot(gs[0], projection='3d')
ax2 = plt.subplot(gs[1])

for axis in [ax1.xaxis, ax1.yaxis, ax1.zaxis]:
    axis.pane.fill = False
    axis.pane.set_edgecolor('white')

ax1.plot_surface(X, Y, ZF, cmap=cm.nipy_spectral, linewidth=0, alpha=0.15)
ax1.contourf(X, Y, ZF, 100, zdir='z', offset=0, cmap=cm.nipy_spectral, alpha=0.15)

ax1.set_xlabel(r"$\rm x \,[km]$")
ax1.set_ylabel(r"$\rm y \,[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, zlim)
ax1.grid(False)

for cfg in DATA_CONFIG[:3]:
    ax2.plot(rval, VeffRezl(rval, cfg["param"]), c=cfg["color"], linewidth=1.5)
    ax2.scatter(cfg["r0"], VeffRezl(cfg["r0"], cfg["param"]), c=cfg["color"], s=80)
ax2.set_ylabel(r"$\rm V(r)$")
ax2.set_xlabel(r"$\rm Radius \, r \,[km]$")
ax2.set_xlim(2.3, xylim)
ax2.set_ylim(0.9, 1.01)

plots_3d = []
for cfg in DATA_CONFIG:
    color = cfg["color"]
    ln3d, = ax1.plot([], [], [], c=color, linewidth=1)
    ln2d, = ax1.plot([], [], [], c=color, linewidth=2.5) # Projektion auf z=0.1
    scat3d = ax1.scatter([], [], [], c=color, s=60)
    scat2d = ax1.scatter([], [], [], c=color, s=40)     # Projektion auf z=0.1
    
    plots_3d.append({
        "data": cfg["data"],
        "ln3d": ln3d, "ln2d": ln2d,
        "scat3d": scat3d, "scat2d": scat2d
    })

def init():
    for p in plots_3d:
        d = p["data"]
        x0 = d[0, 1] * np.cos(d[0, 2])
        y0 = d[0, 1] * np.sin(d[0, 2])
        z0 = np.real(zlam(x0, y0))
        p["ln3d"].set_data([], [])
        p["ln3d"].set_3d_properties([])
        p["ln2d"].set_data([], [])
        p["ln2d"].set_3d_properties([])
        p["scat3d"].set_offsets(np.array([[x0, y0]]))
        p["scat3d"].set_3d_properties([z0], 'z')
        p["scat2d"].set_offsets(np.array([[x0, y0]]))
        p["scat2d"].set_3d_properties([0.1], 'z')
    return fig,

def animate(i):
    idx = step * i
    if idx == 0:
        return init()

    for p in plots_3d:
        d = p["data"]
        x = d[:idx, 1] * np.cos(d[:idx, 2])
        y = d[:idx, 1] * np.sin(d[:idx, 2])
        z = np.real(zlam(x, y))
        p["ln3d"].set_data(x, y)
        p["ln3d"].set_3d_properties(z)
        p["ln2d"].set_data(x, y)
        p["ln2d"].set_3d_properties(np.full_like(x, 0.1))
        p["scat3d"].set_offsets(np.array([[x[-1], y[-1]]]))
        p["scat3d"].set_3d_properties([z[-1]], 'z')
        p["scat2d"].set_offsets(np.array([[x[-1], y[-1]]]))
        p["scat2d"].set_3d_properties([0.1], 'z')
    return fig,

ani = animation.FuncAnimation(fig, animate, init_func=init, frames=anzframes, interval=100, blit=True)
plt.close(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(renderer='iframe')

Zusammenfassung der Zonen um ein nicht-rotierendes, ungeladenes schwarzes Loch

• $r > 6 \, M$: Stabile Kreisbahnen für massereiche Teilchen sind möglich
• $r = r_{\rm ISCO} = 6 \, M$: Der ISCO. Die letzte Grenze für stabile Bahnen.
• $3 \, M < r < 6 \, M$: Nur noch instabile Kreisbahnen für Teilchen mit Masse möglich. Teilchen stürzen unaufhaltsam in das schwarze Loch, wenn sie nicht stark beschleunigen.
• $r = r_{PhSph} = 3 \, M$: Die Photonensphäre. Instabile Kreisbahn für Licht.
• $2 \, M < r < 3 \, M$: Keine Kreisbahnen mehr möglich, weder für Teilchen mit Masse noch für Licht.
• $r = R_{\rm s} = 2 \, M$: Der Ereignishorizont. Ab hier ist kein Entkommen mehr möglich.