Einführung 

Ein wenig mehr als hundert Jahre nachdem Albert Einstein seine Feldgleichungen der Allgemeine Relativiätstheorie der Öffentlichkeit präsentierte, und er damit die Grundlage für Gravitationswellen und schwarzer Löcher formulierte, ist seit einigen Wochen ein Meilenstein in der Geschichte der Astronomie in aller Munde: 

Abb. 1: Erstes Bild eines schwarzen Lochs im Vergleich mit den Simulationsergebnissen des Maple- Worksheets. 

Das Bild zeigt das schwarze Loch im Zentrum unserer Nachbargalaxie M87; bzw. ein wenig genauer, die um ein schwarze Loch entstehende Radiostrahlung (das Bild wurde mittels eines weltweiten Verbunds von Radiowellenteleskopen (EHT: Event Horizon Teleskop) sichtbar gemacht). In Kürze (voraussichtlich im Sommer 2019) werden die aufgenommenen Bilder von Sagittarius A*, das schwarze Loch im Zentrum unserer Galaxie, ebenfalls ausgewertet sein, und weitere Bilder in verbesserter Auflösung werden wohl in einigen Jahren folgen. In beiden Systemen (SL in M87 und Sagittarius A*) kreist Materie in einer Akretionsscheibe um ein supermassives SL und emittiert beim Einfallen Radiowellenstrahlung.Die Bewegung eines Probekörpers (Masse verschwindend klein gegenüber der Masse des schwarzen Lochs) um ein schwarzes Loch ist somit ein astrophysikalisch sehr relevantes Problem. Schon Jahre bevor die ersten Bilder eines SLs entstanden. galt es schon als so gut wie bestätigt, dass im Zentrum unserer Galaxie ein superschweres schwarzes Loch existiert und man verfolgte schon seit Jahrzenten (siehe z.B. R.Genzel et,al.) die Bewegung einzelner sogenannter S-Sterne um dieses Zentrum. 

In diesem Maple-Worksheet wird die Geodätengleichung in vorgegebener Schwarzschild-Raumzeit behandelt und die möglichen Bahnen von Probekörpern mittels eines definierten effektiven Potentials klassifiziert. Neben den gebundenen kreisförmigen und elliptischen Bahnen, den parabolischen und hyperbolischen Bahnverläufen ist auch eine durch das schwarze Loch eingefangene Bahn (capture orbit) möglich. Diese letzteren Bahnen beschreiben die in das schwarze Loch einfallenden Partikel, welche, aufgrund der extremem gravitativen Beschleinigungskräfte, sich aufheizen und elektromagnetische Strahlung emittieren. Ein Verständnis des Erscheinungsbildes des schwarzen Lochs ist möglich indem man sich grundlegenden Begriffe, wie Ereignishorizont, Photonensphäre und der Orbit der letzten stabilen Kreisbahn um ein SL veranschaulicht.  

In diesem Maple-Worksheet werden die, mathematisch anspruchsvollen Gleichungen der Allgemeinen Relativitätstheorie (ART) am Beispiel des schwarzen Lochs in M87 behandelt. Im ersten Teil (Grundlegende Größen der Allgemeinen Relativitätstheorie) erlernen die Studierenden die Verwendung von Computeralgebra-Systemen (Maple und Mathematica). Die oft komplizierten und zeitaufwendigen Berechnungen der tensoriellen Gleichungen der ART können mit Hilfe dieser Programme erleichtert werden. Diverse Anwendungen der Einstein- und Geodätengleichung werden in Maple implementiert, quasi analytische Berechnungen durchgeführt und entsprechende Lösungen berechnet und visualisiert. Der Schwerpunkt liegt hierbei sowohl auf dem Verständnis des Bildes des schwarzen Lochs in M87 als auch auf der Vermittlung spezieller Programmierkenntnisse.  

Im folgenden betrachten wir die Geodätengleichung in der Schwarzschildmetrik. Diese lässt sich 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 \quad, \] wobei die $\Gamma^\mu_{\nu \rho}$ die Christoffel Symbole zweiter Art und $\lambda$ ein affiner Parameter (z.B. die Eigenzeit) darstellen.  

>

 

Grundlegende Größen der Allgemeinen Relativitätstheorie 

Im folgenden werden einige grundlegende Größen der allgemeinen Relativitätstheorie am Beispiel der Schwarzschildmetrik eines nicht-rotierenden schwarzen Lochs (SLs) erläutert. 

Zunächst wird das "tensor"-Paket und weitere Zusatzpakete eingebunden: 

>	restart: with( tensor ): with ( plots ): with( plottools ):

 

Definition der kovarianten Raumzeit-Metrik eines schwarzen Lochs der Masse M (wir werden die Masse des SLs später auf 6.5*10^(9) Sonnenmassen setzen) in Schwarzschildkoordinaten: 

>	coord := [t, r, theta, phi]: g_compts := array(symmetric,sparse, 1..4, 1..4): g_compts[1,1] := 1-2*M/r: g_compts[2,2] := -1/g_compts[1,1]: g_compts[3,3] := -r^2:    g_compts[4,4] := -r^2*sin(theta)^2: g := create( [-1,-1], eval(g_compts));

 

(2.1)

 

Kontravariante Raumzeit-Metrik: 

>	ginv := invert( g, 'detg' ):

 

Erste partielle Ableitung der Metrik: 

>	D1g := d1metric ( g, coord ):

 

Chistoffel Symbole (kontravariante Form): 

>	Cf1 := Christoffel1 ( D1g ):

 

Chistoffel Symbole (erster Index kontravariant, zweiter und dritter kontravariant): 

>	Cf2:= Christoffel2( ginv, Cf1 ):

 

Angabe einzelner Komponenten eines Tensors: 

>	get_compts(Cf2)[2,3,3]:

 

Riemann Tensor: 

>	D2g := d2metric ( D1g, coord ): RMN := Riemann( ginv, D2g, Cf1 ):

 

Ricci Tensor: 

>	RICCI := Ricci( ginv, RMN ):

 

Ricci Skalar: 

>	RS := Ricciscalar( ginv, RICCI ):

 

Einsteintensor: 

>	Estn := Einstein(g, RICCI, RS):

 

Berechnung des infinitesimalen Weglängenelements ds²: 

>	dx:=create([1], array([dt,dr,dtheta,dphi])): ds2:=get_compts(prod(dx,lower(g,dx,1),[1,1])): ds2:=collect( simplify(ds2), [dt,dr,dtheta,dphi]):

 

Hoch- und herunterheben der raumzeitlichen Indices: 

>	raise(ginv,Cf1,1): lower(g,Cf2,1):

 

Beispiel: Riemann Tensor mit allen Indices oben: 

>	RMNinv:=raise(ginv,RMN,1,2,3,4):

 

Das Produkt zweier Tensoren: 

>	Cf1Cf2:=prod(Cf1,Cf2):

 

Contraction zweier Indices: 

>	contract(Cf1Cf2,[4,6]):

 

Kontrahiert man das quadratische Produkt der Riemann Tensors vollständig, so erhält man den Kretschmann-Skalar: 

>	prod(RMN,RMNinv,[1,1],[2,2],[3,3],[4,4]):

 

Beispiel 1: Bewegung eines Probekörpers um ein schwarzes Loch in der Ebene der Akkretionsscheibe 

Basierend auf den im ersten Teil berechneten Größen wird im folgenden eine Klassifizierung der Umlaufbahnen eines Probekörpers um ein schwarzes Loch vorgestell. Im folgenden wird die Geodätengleichung in vorgegebener Schwarzschild Raumzeit betrachtet. Die Geodätengleichung beschreibt wie sich ein Probekörper (Masse m -> 0) im Raum bewegt und sagt voraus, dass diese Bewegung sich stehts entlang der kürzesten Kurve, in der durch die Metrik beschriebenen gekrümmten Raumzeit, vollzieht.Die Geodätengleichung ist ein System aus vier gekoppelten Differentialgleichungen, welche von einem affinen Parameter lambda (später die Eigenzeit tau des in das SL einfallenden Probekörpers) funktional abhängt: 

>	eqns:=geodesic_eqns( coord, lambda, Cf2 );

 

(3.1)

 

Wir lassen nur ebene Bewegungen zu (theta=Pi/2, dtheta=0) und setzen M auf den experimentell gefundenen Wert (6.5*10^9 Sonnenmassen; entspricht 6.5*10^(9)/1.4766 km). Der Ereignishorizont (Schwarzschild Radius)  RS = 2M  des SLs in M87 befindet sich somit bei ungefähr 8.8*10^9 km. 

>

setM:=6.5*10^9*1.4766: eq1:=subs({r=r(lambda),sin(theta)=1,cos(theta)=0,diff(theta(lambda),lambda)=0,theta=Pi/2,M=setM},eqns[1]): eq2:=subs({r=r(lambda),sin(theta)=1,cos(theta)=0,diff(theta(lambda),lambda)=0,theta=Pi/2,M=setM},eqns[2]): eq3:=subs({r=r(lambda),sin(theta)=1,cos(theta)=0,diff(theta(lambda),lambda)=0,theta=Pi/2,M=setM},eqns[3]): eq4:=subs({r=r(lambda),sin(theta)=1,cos(theta)=0,diff(theta(lambda),lambda)=0,theta=Pi/2,M=setM},eqns[4]): eq1:=subs({r(lambda)(lambda)=r(lambda)},eq1); eq2:=subs({r(lambda)(lambda)=r(lambda)},eq2); eq3:=subs({r(lambda)(lambda)=r(lambda)},eq3); eq4:=subs({r(lambda)(lambda)=r(lambda)},eq4);

 

 

 

 

	(3.2)

 

In Abhängigkeit von den Anfangswerten können unterschiedliche Bahnen der Bewegung entstehen. Wir wählen zunächst als Beispiel die Anfangswerte einer geschlossenen Bahn:  

Zur Zeit t=0 sei der Körper bei einem Radius von r=4.5*(Schwarzschildradius), die Anfangsgeschwindigkeit des Körpers sei in radialer Richtung dr=0 und in phi-Richtung dphi=.5.2*10^(-12). Wir beschreiben die Bewegung aus der Sichtweise eines im Unendlichen ruhenden Beobachters. 

>	r0:=4.5*2*setM: t0:=0: phi0:=0: theta0:=Pi/2: dtheta0:=0: dr0:=0: dphi0:=5.2*10^(-12): dt0:=evalf(solve(subs({theta=theta0,dphi=dphi0,dr=dr0,dtheta=dtheta0,M=setM,r=r0},ds2)=1,dt)[1]): dtheta0:=0:

 

In der Literatur wird die Bewegung eines Probekörpers 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 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) charakteriesiert. 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 zweite Gleichung dieses Systems von Differentialgleichungen in die folgenden Gleichungen umschreiben läßt: \[ \begin{eqnarray} \hbox{1. Gleichung:}& \frac{d}{d\lambda} \left[ \left( 1 - \frac{2 M}{r} \right) \frac{dt}{d\lambda} \right] = 0 \quad \rightarrow \quad & \left( 1 - \frac{2 M}{r} \right) \frac{dt}{d\lambda} = E = \hbox{const} \\ \hbox{2. Gleichung:}& \frac{d}{d\lambda} \left( r^2 \hbox{sin}^2(\theta) \frac{d\phi}{d\lambda} \right) = 0 \quad \rightarrow \quad & r^2 \hbox{sin}^2(\theta) \frac{d\phi}{d\lambda} = l = \hbox{const} \quad , \end{eqnarray} \] wobei die während der Bewegungen erhaltenen Größen E (Teilchenenergie pro Masse) und l (Drehimpuls pro Masse m) sich aus der Definition des Viererimpulses $ p_\mu = m u_\mu $ ergeben (Beachte: Die im folgenden benutzten einzelnen Rechenschritte gelten speziell nur för die Schwarzschildmetrik und desweiteren ist $ \lambda=\tau $): \[ \begin{eqnarray} && 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_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} \] 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, 4. Geodätengleichung. So definieren z.B. die Literaturangaben 1.-3. das effektive Potentail wie folgt: \[ \begin{eqnarray} \hbox{4. 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{4. 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} \] Die folgende 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 obigen gewählten Anfangswerten: 

>

plotrmin:=0.28*r0: plotrmax:=3.0*r0: setl:=r0^2*dphi0; setE:=dt0*evalf((1-2*setM/r0)); VeffFb:=(r,M,l)->-M/r+l^2/(2*r^2)-M*l^2/r^3; VeffRez:=(r,M,l)->sqrt(1-2*M/r)*sqrt(1+l^2/r^2); Pot:=plot(VeffRez(r,setM,setl),r=plotrmin..plotrmax,title="Effektives Potential V(r) (Definition nach Referenz 4)"): Pot1:=plot(VeffFb(r,setM,setl),r=plotrmin..plotrmax,title="Effektives Potential V(r) (Definition nach Referenz 1-3)"): B:=pointplot({[r0, VeffRez(r0,setM,setl)]}, symbol=solidcircle,symbolsize=20,color=blue): B1:=pointplot({[r0, VeffFb(r0,setM,setl)]}, symbol=solidcircle,symbolsize=20,color=blue): display(Matrix(1,2,[display(Pot,B),display(Pot1,B1)]));

 

Numerisches Lösen der Geodätengleichung: 

>	Loes:=dsolve({eq1,eq2,eq4,t(0)=t0,r(0)=r0,D(r)(0)=dr0,D(t)(0)=dt0,phi(0)=phi0,D(phi)(0)=dphi0},{r(lambda),t(lambda),phi(lambda)},type=numeric,output=listprocedure);

 

(3.3)

 

Grafische Veranschaulichung der Lösung: 

>

lend:=5*10^12: Plot1:=odeplot(Loes,[lambda,t(lambda)],0..lend,numpoints=200,color=blue,thickness=2,title="Koordinatenzeit t vs affiner Parameter lambda"): Plot2:=odeplot(Loes,[lambda,r(lambda)],0..lend,numpoints=200,color=blue,thickness=2,title="radius vs affiner Parameter lambda"): Plot3:=odeplot(Loes,[r(lambda),t(lambda)],0..lend,numpoints=700,color=blue,thickness=2,title="Koordinatenzeit t vs radius"): display(Matrix(1,3,[Plot1,Plot2,Plot3]));

 

>

frames:=100: BH:=display(disk([0,0],2*setM,color=black)): for i from 0 by 1 to frames do Koerper[i]:=display(disk([rhs(Loes(i*lend/frames)[4])*cos(rhs(Loes(i*lend/frames)[2])),rhs(Loes(i*lend/frames)[4])*sin(rhs(Loes(i*lend/frames)[2]))],2*10^9,color=blue)): TrajK1[i]:=listplot([seq([rhs(Loes(j*lend/frames)[4])*cos(rhs(Loes(j*lend/frames)[2])),rhs(Loes(j*lend/frames)[4])*sin(rhs(Loes(j*lend/frames)[2]))], j = 0 .. i)],color=blue,linestyle=dot): Ani[i]:=display({Koerper[i],TrajK1[i],BH}); od:

 

>	display([seq(Ani[i],i=0..frames)],insequence=true,scaling=constrained):

 

Animation der Bewegung im effektiven Potential: 

>	frames:=100: Pot:=plot(VeffRez(r,setM,setl),r=plotrmin..plotrmax): for i from 0 by 1 to frames do Koerper[i]:=pointplot({[rhs(Loes(i*lend/frames)[4]), setE]}, symbol=solidcircle,symbolsize=20,color=blue): Ani1[i]:=display({Pot,Koerper[i]}); od:

 

>	display([seq(Ani1[i],i=0..frames)],insequence=true):

 

>	Animat1:=display([seq(Ani[i],i=0..frames)],insequence=true,scaling=constrained): Animat2:=display([seq(Ani1[i],i=0..frames)],insequence=true,title="Bewegung im effektives Potential V(r) (Def. nach Ref. 4)"): display(Array([Animat1,Animat2]));

 

EW_Sonne:=evalf(arctan(R_Sonne/A_ESonne)):
convert(evalf(round((i*lend/frames/(c*60*60))*1000)*0.001),string),"Stunden"])] 

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

Im folgenden werden wir die innerste stabile kreisförmige Bahnbewegung (Innermost Stable Circular Orbit (ISCO)) eines Probekörpers berechnen. Die ISCO Trajektorie hat die Eigenschaft: V'=0 und V''=0. Der Drehimpuls l muss somit einen Wert annehmen, dass das effektive Potential einen Sattelpunkt bei einem speziellen Radiuswert annimmt: \[ \begin{eqnarray} \frac{d V}{dr}=0 \,\,,\,\,\, \frac{d^2 V}{dr^2}=0 \quad \underbrace{\Rightarrow}_{\hbox{ISCO}} r=6 M \,\,,\,\,\, l = 2 \sqrt{3} M \end{eqnarray} \] Im vorigen Abschnitt wurde eine Beispielanimation einer Bahnbewegung gezeigt, wobei der  Drehimpuls l des Probekörpers auf den Wert festgelegt wurde. Das effektive Potential V(r,M,l) war somit nur noch von r abhängig, da als Masse die des SLs in M87 gewählt wurde. Die folgende Animation zeigt wie sich die Form des effektiven Potentials bei Variation des Parameters l von l=6 bis l=1 verändert. 

(4.1)

 

>	animate(VeffFb(r,setM,l),r=plotrmin..plotrmax,l=2*10^10..5*10^10,view=-0.4..0.2,title="Effektives Potential V(r) (Referenz 1-3), l=[2,6]",numpoints=500);

 

>	animate(VeffRez(r,setM,l),r=plotrmin..plotrmax,l=2*10^10..5*10^10,view=0.6..1.1,title="Effektives Potential V(r) (Referenz 4), l=[2,6]",numpoints=500);

 

Ist der Wert des Drehimpulses l oberhalb einer gewissen Grenze, so sind stabile Bahnen möglich, da ein Minimum im Potential vorhanden ist. Kreisförmige Bahnbewegungen sind dadurch charakterisiert, dass der Wert des Radiuses 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. Es muss somit gelten: 

>	Diff(V(r,M,l),r)=0; ExtremaV:=diff(VeffFb(r,M,l),r)=0;

 

 

	(4.2)

 

Lösen wir diese Gleichung nach r auf, so erhalten wir die folgenden Lösungen, wobei der erste Wert dem Minimum und der zweite dem Maximum entspricht: 

>	rextr:=solve(ExtremaV,r);

 

(4.3)

 

Darstellen des Radiuswertes im Minimum (grün) und Maximum (rot) des effektiven Potentials 

>	plot({subs(M=setM,rextr[1]),subs(M=setM,rextr[2])},l=2*10^10..5*10^10,title="Radiuswert im Minimum (grün) und Maximum (rot)");

 

Der Drehimpuls l für eine kreisförmige Bahnbewegung ergibt sich duch Auflösen von V'=0 nach l: 

>	lextr:=solve(ExtremaV,l);

 

(4.4)

 

Die beiden Lösungen entsprechen rechts und links rotierenden Probekörpern; Darstellung einer Lösung: 

>	plot(subs(M=setM,lextr[2]),r=plotrmin..plotrmax,title="Drehimpuls l einer kreisförmigen Bahn als Funktion des Radius");

 

Die innerste stabile Kreisbahn hat die Eigenschaft: V'=0 und V''=0. 

>	ExtremaV2:=diff(VeffFb(r,M,l),r,r)=0; lISCO:=solve(simplify(subs(r=rextr[1],ExtremaV2)),l); evalf(subs(M=setM,lISCO[1])); rISCO:=solve(simplify(subs(l=lextr[2],ExtremaV2)),r); evalf(subs(M=setM,rISCO));

 

 

 

 

 

	(4.5)

 

Definition der Anfangswerte der innersten stabilen Kreisbahn und zweier weiterer stabilen kreisförmigen Bahnbewegungen. 

>

t0:=0: phi0:=0: theta0:=Pi/2: dtheta0:=0: dr0:=0: setlA:=2*sqrt(3)*setM; r0A:=subs({M=setM,l=setlA},solve(simplify(diff(VeffFb(r,M,l),r)=0),r)[1]); dphi0A:=setlA/r0A^2: dt0A:=evalf(solve(subs({theta=theta0,dphi=dphi0A,dr=dr0,dtheta=dtheta0,M=setM,r=r0A},ds2)=1,dt)[1]): setEA:=dt0A*evalf((1-2/r0A)); print("------------------------------------------"); setlB:=1.02*setlA; r0B:=evalf(subs({M=setM,l=setlB},solve(simplify(diff(VeffFb(r,M,l),r)=0),r)[1])); dphi0B:=setlB/r0B^2: dt0B:=evalf(solve(subs({theta=theta0,dphi=dphi0B,dr=dr0,dtheta=dtheta0,M=setM,r=r0B},ds2)=1,dt)[1]): setEB:=dt0B*evalf((1-2/r0B)); print("------------------------------------------"); setlC:=1.10*setlA; r0C:=evalf(subs({M=setM,l=setlC},solve(simplify(diff(VeffFb(r,M,l),r)=0),r)[1])); dphi0C:=setlC/r0C^2: dt0C:=evalf(solve(subs({theta=theta0,dphi=dphi0C,dr=dr0,dtheta=dtheta0,M=setM,r=r0C},ds2)=1,dt)[1]): setEC:=dt0C*evalf((1-2/r0C));

 

 

 

 

 

 

 

 

 

 

 

	(4.6)

 

Struktur des effektiven Potentials der innersten stabilen Kreisbahn (blauer Kurve) und zweier weiterer stabilen kreisförmigen Bahnbewegungen. 

>

plotrmin:=0.28*r0: plotrmax:=2*r0: PotA:=plot(VeffFb(r,setM,setlA),r=plotrmin..plotrmax,color=blue): PotB:=plot(VeffFb(r,setM,setlB),r=plotrmin..plotrmax,color=red): PotC:=plot(VeffFb(r,setM,setlC),r=plotrmin..plotrmax,color=green): PlotEA:=pointplot({[r0A, VeffFb(r0A,setM,setlA)]}, symbol=solidcircle,symbolsize=20,color=blue): PlotEB:=pointplot({[r0B, VeffFb(r0B,setM,setlB)]}, symbol=solidcircle,symbolsize=20,color=red): PlotEC:=pointplot({[r0C, VeffFb(r0C,setM,setlC)]}, symbol=solidcircle,symbolsize=20,color=green): display(PotA,PotB,PotC,PlotEA,PlotEB,PlotEC,title="Effektives Potential V(r) (Def. nach Ref. 1-3) für drei verschiedene Drehimpulse",view=[plotrmin..plotrmax,-0.007..-0.1]);

 

Numerische Lösung der drei Bahnbewegungen: 

>

LoesA:=dsolve({eq1,eq2,eq4,t(0)=t0,r(0)=r0A,D(r)(0)=dr0,D(t)(0)=dt0A,phi(0)=phi0,D(phi)(0)=dphi0A},{r(lambda),t(lambda),phi(lambda)},type=numeric,output=listprocedure): LoesB:=dsolve({eq1,eq2,eq4,t(0)=t0,r(0)=r0B,D(r)(0)=dr0,D(t)(0)=dt0B,phi(0)=phi0,D(phi)(0)=dphi0B},{r(lambda),t(lambda),phi(lambda)},type=numeric,output=listprocedure): LoesC:=dsolve({eq1,eq2,eq4,t(0)=t0,r(0)=r0C,D(r)(0)=dr0,D(t)(0)=dt0C,phi(0)=phi0,D(phi)(0)=dphi0C},{r(lambda),t(lambda),phi(lambda)},type=numeric,output=listprocedure):

 

Animation der Bewegungen: Der ISCO ist die blaue Kurve 

>

frames:=80: lend:=0.5*10^12: BH:=display(disk([0,0],2*setM,color=black,transparency=0.1)): for i from 0 by 1 to frames do KoerperA[i]:=display(disk([rhs(LoesA(i*lend/frames)[4])*cos(rhs(LoesA(i*lend/frames)[2])),rhs(LoesA(i*lend/frames)[4])*sin(rhs(LoesA(i*lend/frames)[2]))],2*10^9,color=blue)): KoerperB[i]:=display(disk([rhs(LoesB(i*lend/frames)[4])*cos(rhs(LoesB(i*lend/frames)[2])),rhs(LoesB(i*lend/frames)[4])*sin(rhs(LoesB(i*lend/frames)[2]))],2*10^9,color=red)): KoerperC[i]:=display(disk([rhs(LoesC(i*lend/frames)[4])*cos(rhs(LoesC(i*lend/frames)[2])),rhs(LoesC(i*lend/frames)[4])*sin(rhs(LoesC(i*lend/frames)[2]))],2*10^9,color=green)): TrajA[i]:=listplot([seq([rhs(LoesA(j*lend/frames)[4])*cos(rhs(LoesA(j*lend/frames)[2])),rhs(LoesA(j*lend/frames)[4])*sin(rhs(LoesA(j*lend/frames)[2]))], j = 0 .. i)],color=blue,linestyle=dot): TrajB[i]:=listplot([seq([rhs(LoesB(j*lend/frames)[4])*cos(rhs(LoesB(j*lend/frames)[2])),rhs(LoesB(j*lend/frames)[4])*sin(rhs(LoesB(j*lend/frames)[2]))], j = 0 .. i)],color=red,linestyle=dot): TrajC[i]:=listplot([seq([rhs(LoesC(j*lend/frames)[4])*cos(rhs(LoesC(j*lend/frames)[2])),rhs(LoesC(j*lend/frames)[4])*sin(rhs(LoesC(j*lend/frames)[2]))], j = 0 .. i)],color=green,linestyle=dot): Ani[i]:=display({KoerperA[i],KoerperB[i],KoerperC[i],TrajA[i],TrajB[i],TrajC[i],BH}); od:

 

>	display([seq(Ani[i],i=0..frames)],insequence=true,scaling=constrained,title="Kreisförmige Bahnbewegungen und ISCO");