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 II: Klassifizierung unterschiedlicher Bahnbewegungen

Im Folgenden wird die Geodätengleichung in vorgegebener Schwarzschild Raumzeit betrachtet. Die Geodätengleichung beschreibt wie sich ein Probekörper (Masse Probekörper $m << M$ Masse schwarzes Loch) im Raum bewegt und sagt voraus, dass diese Bewegung sich stets entlang der kürzesten Kurve, in der durch die Metrik beschriebenen gekrümmten Raumzeit, vollzieht.

Sie lässt sich demnach durch folgendes Variationsprinzip herleiten: $$\int_A^B ds= \int_A^B \sqrt{g_{\mu\nu} dx^\mu dx^\nu}=\int_A^B \sqrt{g_{\mu\nu} \frac{dx^\mu}{d\lambda} \frac{dx^\nu}{d\lambda}} d\lambda \rightarrow \hbox{Extremal}$$ , wobei sich dann die Geodätengleichung mittels der Euler-Lagrange Gleichungen $L = \sqrt{g_{\mu\nu} \frac{dx^\mu}{d\lambda} \frac{dx^\nu}{d\lambda}}$, bzw. alternativ $L = g_{\mu\nu} \frac{dx^\mu}{d\lambda} \frac{dx^\nu}{d\lambda}$ ergibt: $$ \frac{d}{d\lambda} \left( \frac{\partial L}{\partial \frac{\partial x^\mu}{\partial \lambda}} \right) - \frac{\partial L}{\partial x^\mu} = 0 \quad \rightarrow \quad \frac{d^2 x^\mu}{d\lambda^2} + \Gamma^\mu_{\nu \rho} \frac{d x^\nu}{d\lambda} \frac{d x^\rho}{d\lambda} ~=~ 0 $$ $\Gamma^\mu_{\nu \rho}$ sind die Christoffel Symbole zweiter Art und $\lambda$ ein affiner Parameter (z.B. die Eigenzeit $\tau$ des Probekörpers).

Wir betrachten im Folgenden 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 $$ in vorgegebener Schwarzschild-Raumzeit $$ g_{\mu\nu}=\left( \begin{array}{ccc} 1-\frac{2\,M}{r} & 0 & 0 & 0\\ 0& -\frac{1}{1-\frac{2\,M}{r}}& 0&0 \\ 0& 0& -r^2& 0\\ 0& 0& 0& -r^2 \hbox{sin}^2(\theta)\\ \end{array} \right) $$ , wobei wir wieder ein sphärisches Koordinatensystem benutzen ($x^\mu=\left(t,r,\theta,\phi \right)$). 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.

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

Wie schon oben erwähnt beschreibt die Geodätengleichung die Bewegung eines Probekörpers der sich entlang der kürzesten Kurve, in der durch die Metrik beschriebenen gekrümmten Raumzeit, vollzieht und ist ein System gekoppelter Differentialgleichungen, das wir uns im folgenden näher betrachten werden.

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

Wir definieren zunächst den ersten Term der Geodätengleichung $\frac{d^2 x^\mu}{d\tau^2}$:

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)
xmu.tensor()

$$\left[\begin{matrix}t{\left (\tau \right )} & r{\left (\tau \right )} & \theta{\left (\tau \right )} & \phi{\left (\tau \right )}\end{matrix}\right]$$

Geod1=diff(xmu.tensor(),tau,tau)
Geod1

$$\left[\begin{matrix}\frac{d^{2}}{d \tau^{2}} t{\left (\tau \right )} & \frac{d^{2}}{d \tau^{2}} r{\left (\tau \right )} & \frac{d^{2}}{d \tau^{2}} \theta{\left (\tau \right )} & \frac{d^{2}}{d \tau^{2}} \phi{\left (\tau \right )}\end{matrix}\right]$$

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

$$\left[\begin{matrix}\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{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) \sin^{2}{\left (\theta \right )} \left(\frac{d}{d \tau} \phi{\left (\tau \right )}\right)^{2} + 2 r \left(\frac{M}{r} - \frac{1}{2}\right) \left(\frac{d}{d \tau} \theta{\left (\tau \right )}\right)^{2} & - \sin{\left (\theta \right )} \cos{\left (\theta \right )} \left(\frac{d}{d \tau} \phi{\left (\tau \right )}\right)^{2} + \frac{2 \frac{d}{d \tau} \theta{\left (\tau \right )} \frac{d}{d \tau} r{\left (\tau \right )}}{r} & \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}\end{matrix}\right]$$

Die Geodätengleichung ist somit ein System von nichtlinearen gekoppelten Differentialgleichungen zweiter Ordnung. Die entsprechenden vier Gleichungen lauten:

Geod=Geod1+Geod2
GeodEq0=Eq(Geod[0],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=Eq(Geod[1],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) \sin^{2}{\left (\theta \right )} \left(\frac{d}{d \tau} \phi{\left (\tau \right )}\right)^{2} + 2 r \left(\frac{M}{r} - \frac{1}{2}\right) \left(\frac{d}{d \tau} \theta{\left (\tau \right )}\right)^{2} + \frac{d^{2}}{d \tau^{2}} r{\left (\tau \right )} = 0$$

GeodEq2=Eq(Geod[2],0)
GeodEq2

$$- \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} \theta{\left (\tau \right )} \frac{d}{d \tau} r{\left (\tau \right )}}{r} = 0$$

GeodEq3=Eq(Geod[3],0)
GeodEq3

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

Beispiel einer gebundenen elliptischen Bahn

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

In Abhängigkeit von den Anfangswerten können unterschiedliche Bahnen der Bewegung entstehen. Wir wählen zunächst als Beispiel einer gebundenen elliptischen Bahn. Zur Zeit $t=\tau=0$ befindet sich der Probekörper bei einem Radius von $r_0=10=5 R_S$ und Winkel $\phi_0=0$. Seine Anfangsgeschwindigkeit in radialer Richtung sei $\left. \frac{dr}{d\tau} \right|_{\tau=0}=0$ und seine Anfangsgeschwindigkeit in $\phi$Richtung sei $\left. \frac{dr}{d\tau} \right|_{\tau=0}=0.041$. 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 Anfangsbedingung $\left. \frac{dt}{d\tau} \right|_{\tau=0}=1.20835632162041$.

r0=10.0
t0=0.0
phi0=0.0
dr0=0.0
dphi0=0.041
dt0=LoesEq_ds[0].subs({(r,r0),(dr,dr0),(dphi,dphi0),(M,1)})
dt0

$$1.20835632162041$$

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

r0=10.0
t0=0.0
phi0=0.0
dr0=0.0
dphi0=0.041
dt0=1/np.sqrt(1-2/r0)*np.sqrt(1+r0**2*dphi0**2)

initialval = np.array([t0,r0,phi0,dt0,dr0,dphi0])

Loes = integrate.odeint(DGLsys, initialval, 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) 

plt.cla()
plt.ylabel(r"$\rm Radius \,\,r \,[km]$")
plt.xlabel(r"$\rm Eigenzeit \,\,\tau$")
plt.plot(tauval,Loes[:, 1],c="blue", linewidth=1.5, linestyle='-');

Das obere Diagramm zeigt den Abstand des Probekörpers (Radius r) als Funktion der Eigenzeit $\tau$. Die Trajektorie des Körpers wird in dem unteren Bild dargestellt, wobei der rote Punkt den Anfangsort des Probekörpers zeigt:

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

Wir werden in Kürze sehen, dass es Erhaltungsgrößen während der Bewegung des Probekörpers gibt, wobei die erhaltenen Größen l (Drehimpuls pro Masse m), e (Teilchenenergie pro Masse) und das Weglängenelement $ds^2$ in der numerischen Simulation in Näherung auch erhalten sind wie man in den folgenden Abbildungen sieht.

import matplotlib.gridspec as gridspec

fig = plt.figure(figsize=(16,5))
gs = gridspec.GridSpec(1, 3, width_ratios=[1,1,1], wspace=0.25)
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1])
ax3 = plt.subplot(gs[2])

ax1.set_ylabel(r"$\rm Drehimpuls \,\,l$")
ax1.set_xlabel(r"$\rm Eigenzeit \,\,\tau$")
ax2.set_ylabel(r"$\rm Energie \,\,e $")
ax2.set_xlabel(r"$\rm Eigenzeit \,\,\tau$")
ax3.set_ylabel(r"$\rm Weglängenelement \,\,ds^2$")
ax3.set_xlabel(r"$\rm Eigenzeit \,\,\tau$")

datal=Loes[:, 1]**2*Loes[:, 5]
dataE=(1-2/Loes[:, 1])*Loes[:, 3]
datads2=(1-2/Loes[:, 1])*Loes[:, 3]**2-Loes[:, 4]**2/(1-2/Loes[:, 1])-Loes[:, 5]**2*Loes[:, 1]**2

ax1.plot(tauval,datal,c="blue", linewidth=1.5, linestyle='-');
ax2.plot(tauval,dataE,c="blue", linewidth=1.5, linestyle='-');
ax3.plot(tauval,datads2,c="blue", linewidth=1.5, linestyle='-');

Klassifizierung unterschiedlicher Bahnbewegungen

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äßt: $$ \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} $$

l = symbols('l')
VeffFb=-M/r+l**2/(2*r**2)-M*l**2/r**3
VeffRez=sqrt(1-2*M/r)*sqrt(1+l**2/r**2)

setl=r0**2*dphi0
setE=dt0*((1-2/r0))
print("Drehimpuls l")
print(setl)
print("Energie e")
print(setE)
VeffFbl=lambdify((r),VeffFb.subs({(l,setl),(M,1)}))
VeffRezl=lambdify((r),VeffRez.subs({(l,setl),(M,1)}))
rval = np.linspace(2.85,30, 1001)

Drehimpuls l
4.1000000000000005
Energie e
0.9666850572963255

fig = plt.figure(figsize=(16,8))
gs = gridspec.GridSpec(1, 2, width_ratios=[1,1], wspace=0.21)
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1])
ax1.set_ylabel(r"$\rm V_1(r)$")
ax1.set_xlabel(r"$\rm Radius \, r \,[km]$")
ax2.set_ylabel(r"$\rm V_2(r)$")
ax2.set_xlabel(r"$\rm Radius \, r \,[km]$")
ax1.plot(rval,VeffRezl(rval),c="blue", linewidth=1.5, linestyle='-');
ax2.plot(rval,VeffFbl(rval),c="red", linewidth=1.5, linestyle='-');

ax1.scatter(r0,VeffRezl(r0),c="black", s=80, marker='o');
ax2.scatter(r0,VeffFbl(r0),c="black", s=80, marker='o');

Die obere Abbildung zeigt (in der Nomenklatur des 4.Buches (Diagramm links) bzw. des 1. bis 3. Buches (Diagramm rechts) ) das effektive Potential als Funktion des Radius bei den zuvor gewählten Anfangswerten (schwarze Punkte auf den Kurven).

Wir wollen uns nun die Bewegung des Probekörpers in einer Animation darstellen:

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

step=100
anzframes=100
xylim=18
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(Loes[0, 1]*np.cos(Loes[0, 2]),Loes[0, 1]*np.sin(Loes[0, 2]),c="red", s=60, marker='o');
    BH=plt.Circle((0, 0), 2, color='black')
    ax1.add_patch(BH);
    ax2.plot(rval,VeffRezl(rval),c="blue", linewidth=1.5, linestyle='-');
    ax2.scatter(r0,VeffRezl(r0),c="red", 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)
    return fig,

def animate(i):
    ax1.cla() 
    ax2.cla() 
    ax1.plot(Loes[:step*i, 1]*np.cos(Loes[:step*i, 2]),Loes[:step*i, 1]*np.sin(Loes[:step*i, 2]),c="blue", linewidth=0.5, linestyle='-');
    ax1.scatter(Loes[step*i, 1]*np.cos(Loes[step*i, 2]),Loes[step*i, 1]*np.sin(Loes[step*i, 2]),c="red", s=60, marker='o');
    ax2.plot(rval,VeffRezl(rval),c="blue", linewidth=1.5, linestyle='-');
    ax2.scatter(Loes[step*i, 1],VeffRezl(r0),c="red", s=80, marker='o');
    BH=plt.Circle((0, 0), 2, color='black')
    ax1.add_patch(BH);
    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)
    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 gleiche Animation in einem eingebetteten Diagramm der räumlichen Hypersphäre der Mannigfaltigkeit:

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

step=100
anzframes=100
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(Loes[0, 1]*np.cos(Loes[0, 2]),Loes[0, 1]*np.sin(Loes[0, 2]),np.real(zlam(Loes[0, 1]*np.cos(Loes[0, 2]),Loes[0, 1]*np.sin(Loes[0, 2]))),c="black", s=60, marker='o');
    ax1.contourf(X, Y, ZF,100, zdir='z', offset=0, cmap=cm.nipy_spectral,alpha=0.5)
    ax1.scatter(Loes[0, 1]*np.cos(Loes[0, 2]),Loes[0, 1]*np.sin(Loes[0, 2]),0.1,c="black", s=40, marker='o');
    ax2.plot(rval,VeffRezl(rval),c="black", linewidth=1.5, linestyle='-');
    ax2.scatter(r0,VeffRezl(r0),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)
    ax1.grid(False)
    return fig,

def animate(i):
    ax1.cla() 
    ax2.cla()
    ax1.plot_surface(X, Y, ZF, cmap=cm.nipy_spectral, linewidth=0, alpha=0.5)
    ax1.plot(Loes[:step*i, 1]*np.cos(Loes[:step*i, 2]),Loes[:step*i, 1]*np.sin(Loes[:step*i, 2]),np.real(zlam(Loes[:step*i, 1]*np.cos(Loes[:step*i, 2]),Loes[:step*i, 1]*np.sin(Loes[:step*i, 2]))),c="black", linewidth=1, linestyle='-');
    ax1.scatter(Loes[step*i, 1]*np.cos(Loes[step*i, 2]),Loes[step*i, 1]*np.sin(Loes[step*i, 2]),np.real(zlam(Loes[step*i, 1]*np.cos(Loes[step*i, 2]),Loes[step*i, 1]*np.sin(Loes[step*i, 2]))),c="black", s=60, marker='o');
    ax1.scatter(Loes[step*i, 1]*np.cos(Loes[step*i, 2]),Loes[step*i, 1]*np.sin(Loes[step*i, 2]),0.1,c="black", s=40, marker='o');
    ax1.contourf(X, Y, ZF,100, zdir='z', offset=0, cmap=cm.nipy_spectral,alpha=0.5)
    ax1.plot(Loes[:step*i, 1]*np.cos(Loes[:step*i, 2]),Loes[:step*i, 1]*np.sin(Loes[:step*i, 2]),0.1,c="black", linewidth=2.5, linestyle='-');
    ax2.plot(rval,VeffRezl(rval),c="black", linewidth=1.5, linestyle='-');
    ax2.scatter(Loes[step*i, 1],VeffRezl(r0),c="black", s=80, marker='o');
    ax1.set_ylabel(r"$\rm y \,[km]$")
    ax1.set_xlabel(r"$\rm x \,[km]$")
    ax1.set_zlabel("Einbettung 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)  
    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())

Im Folgenden werden die unterschiedlichen Bahnbewegungen um ein schwarzes Loch klassifiziert.

t0=0
phi0=0
dr0=0
setl=4.1
setVeffRez=VeffRez.subs({(l,setl),(M,1)})
EqCircOrb=Eq(diff(setVeffRez,r),0)
r0A=np.float(solve(EqCircOrb)[1])
dphi0A=setl/r0A**2
dt0A=np.float(LoesEq_ds[0].subs({(r,r0A),(dr,dr0),(dphi,dphi0A),(M,1)}))
setEA=dt0A*((1-2/r0A))
initialvalA = np.array([t0,r0A,phi0,dt0A,dr0,dphi0A])

r0B=8
dphi0B=setl/r0B**2
dt0B=np.float(LoesEq_ds[0].subs({(r,r0B),(dr,dr0),(dphi,dphi0B),(M,1)}))
setEB=dt0B*(1-2/r0B)

setEC=np.float(limit_seq(setVeffRez, r))
EqE1=Eq(setVeffRez,setEC)
r0C=np.float(solve(EqE1,r)[1])
dphi0C=setl/r0C**2
dt0C=np.float(LoesEq_ds[0].subs({(r,r0C),(dr,dr0),(dphi,dphi0C),(M,1)}))

setED=1.005
EqE1005=Eq(setVeffRez,setED)
r0D=solve(EqE1005,r)[2]
r0D=float(re(r0D))
dphi0D=setl/r0D**2
dt0D=np.float(LoesEq_ds[0].subs({(r,r0D),(dr,dr0),(dphi,dphi0D),(M,1)}))


r0E=225
dphi0E=setl/r0E**2
dr0E=-0.185388257
dr0E=-0.185388257
dt0E=np.float(LoesEq_ds[0].subs({(r,r0E),(dr,dr0E),(dphi,dphi0E),(M,1)}))
#setEE=1.01273
setEE=dt0E*((1-2/r0E))

tauval = np.linspace(0, 1280, 50001)

initialvalA = np.array([t0,r0A,phi0,dt0A,dr0,dphi0A])
initialvalB = np.array([t0,r0B,phi0,dt0B,dr0,dphi0B])
initialvalC = np.array([t0,r0C,phi0,dt0C,dr0,dphi0C])
initialvalD = np.array([t0,r0D,phi0,dt0D,dr0,dphi0D])
initialvalE = np.array([t0,r0E,phi0,dt0E,dr0E,dphi0E])

LoesA = integrate.odeint(DGLsys, initialvalA, tauval)
LoesB = integrate.odeint(DGLsys, initialvalB, tauval)
LoesC = integrate.odeint(DGLsys, initialvalC, tauval)
LoesD = integrate.odeint(DGLsys, initialvalD, tauval)
LoesE = integrate.odeint(DGLsys, initialvalE, tauval)

/home/hanauske/anaconda3/lib/python3.7/site-packages/scipy/integrate/odepack.py:236: ODEintWarning: Excess work done on this call (perhaps wrong Dfun type). Run with full_output = 1 to get quantitative information.
  warnings.warn(warning_msg, ODEintWarning)

zBett=sqrt(8*M*(r-2*M))
R = np.linspace(2.00001, 40, 1000)
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.85,40, 1001)

step=500
anzframes=100
xylim=35
zlim=20
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.contourf(X, Y, ZF,100, zdir='z', offset=0, cmap=cm.nipy_spectral,alpha=0.5)
    ax1.scatter(LoesA[0, 1]*np.cos(LoesA[0, 2]),LoesA[0, 1]*np.sin(LoesA[0, 2]),np.real(zlam(LoesA[0, 1]*np.cos(LoesA[0, 2]),LoesA[0, 1]*np.sin(LoesA[0, 2]))),c="blue", s=60, marker='o');
    ax1.scatter(LoesA[0, 1]*np.cos(LoesA[0, 2]),LoesA[0, 1]*np.sin(LoesA[0, 2]),0.1,c="blue", s=40, marker='o');
    ax1.scatter(LoesB[0, 1]*np.cos(LoesB[0, 2]),LoesB[0, 1]*np.sin(LoesB[0, 2]),np.real(zlam(LoesB[0, 1]*np.cos(LoesB[0, 2]),LoesB[0, 1]*np.sin(LoesB[0, 2]))),c="red", s=60, marker='o');
    ax1.scatter(LoesB[0, 1]*np.cos(LoesB[0, 2]),LoesB[0, 1]*np.sin(LoesB[0, 2]),0.1,c="red", s=40, marker='o');
    ax1.scatter(LoesC[0, 1]*np.cos(LoesC[0, 2]),LoesC[0, 1]*np.sin(LoesC[0, 2]),np.real(zlam(LoesC[0, 1]*np.cos(LoesC[0, 2]),LoesC[0, 1]*np.sin(LoesC[0, 2]))),c="green", s=60, marker='o');
    ax1.scatter(LoesC[0, 1]*np.cos(LoesC[0, 2]),LoesC[0, 1]*np.sin(LoesC[0, 2]),0.1,c="green", s=40, marker='o');
    ax1.scatter(LoesD[0, 1]*np.cos(LoesD[0, 2]),LoesD[0, 1]*np.sin(LoesD[0, 2]),np.real(zlam(LoesD[0, 1]*np.cos(LoesD[0, 2]),LoesD[0, 1]*np.sin(LoesD[0, 2]))),c="grey", s=60, marker='o');
    ax1.scatter(LoesD[0, 1]*np.cos(LoesD[0, 2]),LoesD[0, 1]*np.sin(LoesD[0, 2]),0.1,c="grey", s=40, marker='o');
    ax1.scatter(LoesE[0, 1]*np.cos(LoesE[0, 2]),LoesE[0, 1]*np.sin(LoesE[0, 2]),np.real(zlam(LoesE[0, 1]*np.cos(LoesE[0, 2]),LoesE[0, 1]*np.sin(LoesE[0, 2]))),c="black", s=60, marker='o');
    ax1.scatter(LoesE[0, 1]*np.cos(LoesE[0, 2]),LoesE[0, 1]*np.sin(LoesE[0, 2]),0.1,c="black", s=40, marker='o');
    ax2.plot(rval,VeffRezl(rval),c="black", linewidth=1.5, linestyle='-');
    ax2.scatter(r0,VeffRezl(r0),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=-130, elev=35)
    ax1.set_xlim(-xylim, xylim)
    ax1.set_ylim(-xylim, xylim)
    ax1.set_zlim(0, zlim)
    ax2.set_xlim(2, xylim)
    ax1.grid(False)
    return fig,

def animate(i):
    ax1.cla() 
    ax2.cla()
    ax1.plot_surface(X, Y, ZF, cmap=cm.nipy_spectral, linewidth=0, alpha=0.5)
    ax1.contourf(X, Y, ZF,100, zdir='z', offset=0, cmap=cm.nipy_spectral,alpha=0.5)
    Loes=LoesA
    ax1.plot(Loes[:step*i, 1]*np.cos(Loes[:step*i, 2]),Loes[:step*i, 1]*np.sin(Loes[:step*i, 2]),np.real(zlam(Loes[:step*i, 1]*np.cos(Loes[:step*i, 2]),Loes[:step*i, 1]*np.sin(Loes[:step*i, 2]))),c="blue", linewidth=1, linestyle='-');
    ax1.plot(Loes[:step*i, 1]*np.cos(Loes[:step*i, 2]),Loes[:step*i, 1]*np.sin(Loes[:step*i, 2]),0.1,c="blue", linewidth=2.5, linestyle='-');
    ax1.scatter(Loes[step*i, 1]*np.cos(Loes[step*i, 2]),Loes[step*i, 1]*np.sin(Loes[step*i, 2]),np.real(zlam(Loes[step*i, 1]*np.cos(Loes[step*i, 2]),Loes[step*i, 1]*np.sin(Loes[step*i, 2]))),c="blue", s=60, marker='o');
    ax1.scatter(Loes[step*i, 1]*np.cos(Loes[step*i, 2]),Loes[step*i, 1]*np.sin(Loes[step*i, 2]),0.1,c="blue", s=40, marker='o');
    Loes=LoesB
    ax1.plot(Loes[:step*i, 1]*np.cos(Loes[:step*i, 2]),Loes[:step*i, 1]*np.sin(Loes[:step*i, 2]),np.real(zlam(Loes[:step*i, 1]*np.cos(Loes[:step*i, 2]),Loes[:step*i, 1]*np.sin(Loes[:step*i, 2]))),c="red", linewidth=1, linestyle='-');
    ax1.plot(Loes[:step*i, 1]*np.cos(Loes[:step*i, 2]),Loes[:step*i, 1]*np.sin(Loes[:step*i, 2]),0.1,c="red", linewidth=2.5, linestyle='-');
    ax1.scatter(Loes[step*i, 1]*np.cos(Loes[step*i, 2]),Loes[step*i, 1]*np.sin(Loes[step*i, 2]),np.real(zlam(Loes[step*i, 1]*np.cos(Loes[step*i, 2]),Loes[step*i, 1]*np.sin(Loes[step*i, 2]))),c="red", s=60, marker='o');
    ax1.scatter(Loes[step*i, 1]*np.cos(Loes[step*i, 2]),Loes[step*i, 1]*np.sin(Loes[step*i, 2]),0.1,c="red", s=40, marker='o');
    Loes=LoesC
    ax1.plot(Loes[:step*i, 1]*np.cos(Loes[:step*i, 2]),Loes[:step*i, 1]*np.sin(Loes[:step*i, 2]),np.real(zlam(Loes[:step*i, 1]*np.cos(Loes[:step*i, 2]),Loes[:step*i, 1]*np.sin(Loes[:step*i, 2]))),c="green", linewidth=1, linestyle='-');
    ax1.plot(Loes[:step*i, 1]*np.cos(Loes[:step*i, 2]),Loes[:step*i, 1]*np.sin(Loes[:step*i, 2]),0.1,c="green", linewidth=2.5, linestyle='-');
    ax1.scatter(Loes[step*i, 1]*np.cos(Loes[step*i, 2]),Loes[step*i, 1]*np.sin(Loes[step*i, 2]),np.real(zlam(Loes[step*i, 1]*np.cos(Loes[step*i, 2]),Loes[step*i, 1]*np.sin(Loes[step*i, 2]))),c="green", s=60, marker='o');
    ax1.scatter(Loes[step*i, 1]*np.cos(Loes[step*i, 2]),Loes[step*i, 1]*np.sin(Loes[step*i, 2]),0.1,c="green", s=40, marker='o');
    Loes=LoesD
    ax1.plot(Loes[:step*i, 1]*np.cos(Loes[:step*i, 2]),Loes[:step*i, 1]*np.sin(Loes[:step*i, 2]),np.real(zlam(Loes[:step*i, 1]*np.cos(Loes[:step*i, 2]),Loes[:step*i, 1]*np.sin(Loes[:step*i, 2]))),c="grey", linewidth=1, linestyle='-');
    ax1.plot(Loes[:step*i, 1]*np.cos(Loes[:step*i, 2]),Loes[:step*i, 1]*np.sin(Loes[:step*i, 2]),0.1,c="grey", linewidth=2.5, linestyle='-');
    ax1.scatter(Loes[step*i, 1]*np.cos(Loes[step*i, 2]),Loes[step*i, 1]*np.sin(Loes[step*i, 2]),np.real(zlam(Loes[step*i, 1]*np.cos(Loes[step*i, 2]),Loes[step*i, 1]*np.sin(Loes[step*i, 2]))),c="grey", s=60, marker='o');
    ax1.scatter(Loes[step*i, 1]*np.cos(Loes[step*i, 2]),Loes[step*i, 1]*np.sin(Loes[step*i, 2]),0.1,c="grey", s=40, marker='o');
    Loes=LoesE
    ax1.plot(Loes[:step*i, 1]*np.cos(Loes[:step*i, 2]),Loes[:step*i, 1]*np.sin(Loes[:step*i, 2]),np.real(zlam(Loes[:step*i, 1]*np.cos(Loes[:step*i, 2]),Loes[:step*i, 1]*np.sin(Loes[:step*i, 2]))),c="black", linewidth=1, linestyle='-');
    ax1.plot(Loes[:step*i, 1]*np.cos(Loes[:step*i, 2]),Loes[:step*i, 1]*np.sin(Loes[:step*i, 2]),0.1,c="black", linewidth=2.5, linestyle='-');
    ax1.scatter(Loes[step*i, 1]*np.cos(Loes[step*i, 2]),Loes[step*i, 1]*np.sin(Loes[step*i, 2]),np.real(zlam(Loes[step*i, 1]*np.cos(Loes[step*i, 2]),Loes[step*i, 1]*np.sin(Loes[step*i, 2]))),c="black", s=60, marker='o');
    ax1.scatter(Loes[step*i, 1]*np.cos(Loes[step*i, 2]),Loes[step*i, 1]*np.sin(Loes[step*i, 2]),0.1,c="black", s=40, marker='o');

    ax2.plot(rval,VeffRezl(rval),c="black", linewidth=1.5, linestyle='-');
    ax2.scatter(LoesA[step*i, 1],VeffRezl(LoesA[0, 1]),c="blue", s=80, marker='o');
    ax2.scatter(LoesB[step*i, 1],VeffRezl(LoesB[0, 1]),c="red", s=80, marker='o');
    ax2.plot(LoesB[:step*i, 1],VeffRezl(LoesB[0, 1]*LoesB[:step*i, 1]/LoesB[:step*i, 1]),c="red", linewidth=1, linestyle=':');
    ax2.scatter(LoesC[step*i, 1],VeffRezl(LoesC[0, 1]),c="green", s=80, marker='o');
    ax2.plot(LoesC[:step*i, 1],VeffRezl(LoesC[0, 1]*LoesC[:step*i, 1]/LoesC[:step*i, 1]),c="green", linewidth=1, linestyle=':');
    ax2.scatter(LoesD[step*i, 1],VeffRezl(LoesD[0, 1]),c="grey", s=80, marker='o');
    ax2.plot(LoesD[:step*i, 1],VeffRezl(LoesD[0, 1]*LoesD[:step*i, 1]/LoesD[:step*i, 1]),c="grey", linewidth=1, linestyle=':');
    ax2.scatter(LoesE[step*i, 1],setEE,c="black", s=80, marker='o');
    ax2.plot(LoesE[:step*i, 1],setEE*LoesE[:step*i, 1]/LoesE[:step*i, 1],c="black", linewidth=1, linestyle=':');
    
    
    ax1.set_ylabel(r"$\rm y \,[km]$")
    ax1.set_xlabel(r"$\rm x \,[km]$")
    ax1.set_zlabel("Einbettung z");
    ax2.set_ylabel(r"$\rm V(r)$")
    ax2.set_xlabel(r"$\rm Radius \, r \,[km]$")
    ax1.view_init(azim=-130, elev=35)
    ax1.set_xlim(-xylim, xylim)
    ax1.set_ylim(-xylim, xylim)
    ax1.set_zlim(0, zlim)  
    ax2.set_xlim(2, xylim)
    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())

/home/hanauske/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:84: RuntimeWarning: invalid value encountered in true_divide

In der oberen Animation erkennt man neben den gebundenen kreisförmigen (A: blau) und elliptischen (B: rot) Bahnen, den parabolischen (C: grün) und hyperbolischen (D: grau) Bahnverläufen auch eine durch das schwarze Loch eingefangene Bahn (E: capture orbit, schwarz). Diese ist in der unteren Abbildung nochmals separat veranschaulicht.

Loes=LoesE
ax = plt.gca()
ax.cla() 
ax.set_ylabel(r"$\rm y \,[km]$")
ax.set_xlabel(r"$\rm x \,[km]$")
ax.plot(Loes[:, 1]*np.cos(Loes[:, 2]),Loes[:, 1]*np.sin(Loes[:, 2]),c="black", linewidth=1.5, linestyle='-');
ax.scatter(Loes[0, 1]*np.cos(Loes[0, 2]),Loes[0, 1]*np.sin(Loes[0, 2]),c="black", s=60, marker='o');
ax.axis('equal');

Die Visualisierung kann man auch mittels der "Plotly Python Open Source Graphing Library" erstellen ( siehe https://plotly.com/python/ ). In der folgenden interaktiven Abbildung wird der "Capture orbit" und der "Elliptische Orbit" visualisiert. Zusätzlich wird die Eigen- und Koordinatenzeit ausgegeben, wenn man mit der Maus genau auf eine der beiden Probekörper-Trajektorien fokussiert.

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=LoesE
labels=[]
for i in range(len(Loes[:, 0])):
    labels.append('Eigenzeit '+str(round(tauval[i], 2))+', Koordinatenzeit = '+str(round(Loes[i, 0], 2)))
capture_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="Capture orbit",
    mode='lines',
    line=dict(color='black', width=8),
    text=labels,
    hoverinfo='text'
    )

Loes=LoesB
labels=[]
for i in range(len(Loes[:, 0])):
    labels.append('Eigenzeit '+str(round(tauval[i], 2))+', Koordinatenzeit = '+str(round(Loes[i, 0], 2)))
elliptic_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="Elliptic orbit",
    mode='lines',
    line=dict(color='red', width=8),
    text=labels,
    hoverinfo='text'
    )

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=[-40,40],
                    yaxis_range=[-40,40],
                    zaxis_range=[0,20],
                    xaxis_title='x',
                    yaxis_title='y',
                    zaxis_title='f'),
                  legend=dict(
                    yanchor="top",
                    y=0.99,
                    xanchor="left",
                    x=0.01)
            )

data=[trichter,elliptic_orbit,capture_orbit]
fig=go.Figure(data=data, layout=layout)

fig.show()