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 27.05.2021

Erster Vorlesungsteil: Allgemeine Relativitätstheorie mit Python

Die Tolman-Oppenheimer-Volkoff (TOV) Gleichung

Innenraum-Lösung eines sphärisch symmetrischen, statischen Objektes (z.B. Erde, Neutronenstern)

Von der Einstein Gleichung zur TOV Gleichung

In dem vorigen Jupyter Notebook ( siehe TOV.ipny ) hatten wir die TOV-Gleichungen analytisch aus der Einstein-Gleichung unter Zuhilfenahme der kovarianten Erhaltung des Energieimpulses hergeleitet

$$ G^\mu{}\!_\nu = R^\mu{}\!_\nu - \frac{1}{2}g^\mu{}\!_\nu R = 8\pi \, T^\mu{}\!_\nu \quad, \quad \nabla\!_\mu T^\mu{}\!_\nu = 0 \quad . $$

Dabei setzen wir die Metrik und den Energie-Impuls-Tensor im Inneren des Objektes wie folgt an:

$$g_{\mu\nu}=\left( \begin{array}{ccc} e^{2\Phi(r)} & 0 & 0 & 0\\ 0& - \left( 1 - \frac{2 m(r)}{r} \right)^{-1}& 0&0 \\ 0& 0& -r^2& 0\\ 0& 0& 0& -r^2 \hbox{sin}^2(\theta)\\ \end{array} \right) \quad , \quad T^\mu{}\!_\nu=\left( \begin{array}{ccc} e(r) & 0 & 0 & 0\\ 0& -p(r)& 0&0 \\ 0& 0& -p(r)& 0\\ 0& 0& 0& -p(r)\\ \end{array} \right) \quad . $$

Durch Umschreiben und Kombination der Gleichungen gelangt man zu einem System von drei gekoppelten Differentialgleichungen erster Ordnung, die sogenannten Tolman-Oppenheimer-Volkoff (TOV) Gleichungen:

$$ \begin{equation} \frac{dm}{dr} ~=~ 4\pi r^2~e ~, \quad \frac{dp}{dr} ~=~ -(e+p) \, \frac{m+4 \pi r^3 p}{r \left(r -2m \right)} ~, \quad \frac{d\Phi}{dr} ~=~ \frac{m+4 \pi r^3 p}{r \left(r -2m \right)} \quad . \end{equation} $$

Aufgrund des Birkoff-Theorems muss die Innenraum-Lösung der Metrik in die äußere Schwarzschildmetrik am Sternrand stetig ineinander übergehen. Am Rand des betrachteten Objektes ($r=R$) gilt für die $g_{tt}$-Komponente der Metrik:

$$ g_{tt}(R)=e^{2\,\Phi(R)} = \left( 1 - \frac{M}{R} \right) = \left( 1 - \frac{m(R)}{R} \right) \quad , $$

wobei $M=m(R)$ die Gesamtmasse des Sterns ist.

Wir hatten zusätzlich die numerische Lösung dieses System von gekoppelten Differentialgleichungen erster Ordnung mittels eines polytropen Ansatzes der Zustandsgleichung der Form $p(e)=K\,e^{\gamma}$ berechnet und uns die Eigenschaften der Neutronensterne visualisiert.

In diesem Jupyter Notebook werden wir die folgenden drei zusätzlichen Teilaspekte der TOV-Gleichungen behandeln:

• Spezialfall: Gravitationsfeld einer Kugel aus inkompressibler Flüssigkeit (konstante Dichte)
• Visualisierung der raumzeitlichen Struktur eines kompakten Sterns mittels des eingebetteten Diagramms der räumlichen Hypersphäre der Mannigfaltigkeit
• Numerisches Lösen der TOV-Gleichungen mittels des Euler-Verfahrens

Spezialfall: Gravitationsfeld einer Kugel aus inkompressibler Flüssigkeit (konstante Dichte)

Die TOV-Gleichungen wurden im Jahre 1939 gefunden und gründen auf den folgenden zwei Arbeiten: Tolman, Richard C. "Static solutions of Einstein's field equations for spheres of fluid." Physical Review 55.4 (1939): 364. und Oppenheimer, J. Robert, and George M. Volkoff. "On massive neutron cores." Physical Review 55.4 (1939): 374.. Es ist beachtlich, dass bereits im Jahre 1916, lange bevor die TOV-Gleichungen publiziert wurden, Karl Schwarzschild schon einen wichtigen Spezialfall der TOV-Gleichungen analytisch berechnete. Herr Schwarzschild betrachtete einen sphärisch symmetrischen Körper mit konstanter Dichte (siehe Karl Schwarzschild, "Gravitationsfeld einer Kugel aus inkompressibler Flüssigkeit ", Sitzungsberichte der Königlich-Preussischen Akademie der Wissenschaften. Reimer, Berlin 1916, S:424-434) und berechnete die Eigenschaften dieses Körpers in der Einsteins allgemeiner Relativitätstheorie.

Im Folgenden werden wir auch diesen Spezialfall betrachten und die analytische Lösung mit unseren numerischen Berechnungen vergleichen. Wir starten hingegen nicht, wie Schwarzschild es damals, von der Einsteingleichung, sondern nehmen die TOV-Gleichungen als gegeben an.

import numpy as np
from sympy import *
init_printing()
from einsteinpy.symbolic import *

r, Pi = symbols('r, pi', positive = True, real = True)
Fphi = Function('\Phi')(r)
Fm = Function('m')(r)
e = Function('e', positive = True, real = True)(r)
p = Function('p', positive = True, real = True)(r)

TOV1=Eq(Fm.diff(r),4*Pi*r**2*e)
TOV2=Eq(p.diff(r),-(e+p)*(Fm + 4*Pi*r**3*p)/(r*(r-2*Fm)))
TOV3=Eq(Fphi.diff(r),(Fm + 4*Pi*r**3*p)/(r*(r-2*Fm)))

Wir substituieren in den TOV-Gleichungen einen konstanten Wert der Energiedichte: $e(r)=e_0=$const

econst = symbols('e_0', positive = True, real = True)
TOV1a=TOV1.subs(e,econst)
TOV2a=TOV2.subs(e,econst)
TOV3a=TOV3
TOV1a,TOV2a,TOV3a

$$\left ( \frac{d}{d r} m{\left (r \right )} = 4 e_{0} \pi r^{2}, \quad \frac{d}{d r} p{\left (r \right )} = \frac{\left(- e_{0} - p{\left (r \right )}\right) \left(4 \pi r^{3} p{\left (r \right )} + m{\left (r \right )}\right)}{r \left(r - 2 m{\left (r \right )}\right)}, \quad \frac{d}{d r} \Phi{\left (r \right )} = \frac{4 \pi r^{3} p{\left (r \right )} + m{\left (r \right )}}{r \left(r - 2 m{\left (r \right )}\right)}\right )$$

Zunächst betrachten wir uns die ersten beiden TOV-Gleichungen. Die Lösung der ersten Gleichung lässt sich durch einfache Integration berechnen, wobei wir als Randbedingung $m(r=0)=0$ setzen:

Loesm=dsolve(TOV1a,ics={Fm.subs(r,0):0})
Loesm

$$m{\left (r \right )} = \frac{4 e_{0} \pi r^{3}}{3}$$

Wir nehmen im Folgenden an, dass die betrachtete Kugel einen Radius $R$ habe, sodass die gesamte gravitative Masse $M$ der Kugel sich durch $M=m(R)$ angeben lässt:

R, M = symbols('R M', positive = True, real = True)
EqMasse=Eq(M,Loesm.rhs.subs(r,R))
EqMasse

$$M = \frac{4 R^{3} e_{0} \pi}{3}$$

Die gefundene Lösung für $m(r)$ substituieren wir in die zweite TOV-Gleichung:

TOV2b=TOV2a.subs(Fm,Loesm.rhs)
TOV2b

$$\frac{d}{d r} p{\left (r \right )} = \frac{\left(- e_{0} - p{\left (r \right )}\right) \left(\frac{4 e_{0} \pi r^{3}}{3} + 4 \pi r^{3} p{\left (r \right )}\right)}{r \left(- \frac{8 e_{0} \pi r^{3}}{3} + r\right)}$$

Mittels des Befehls "dsolve" können wir die analytische Lösung dieser Differentialgleichung bestimmen. Man erhält zwei mögliche Lösungen, die noch von dem Parameter $C_1$ abhängen, den wir durch unsere Randbedingungen festlegen können.

LoesTOV2b=dsolve(TOV2b)
LoesTOV2b

$$\left [ p{\left (r \right )} = - \frac{e_{0} \left(- 24 e_{0} \pi r^{2} e^{C_{1} e_{0}} + 2 \sqrt{\left(8 e_{0} \pi r^{2} - 3\right) e^{C_{1} e_{0}}} + 9 e^{C_{1} e_{0}} + 1\right)}{3 \left(- 8 e_{0} \pi r^{2} e^{C_{1} e_{0}} + 3 e^{C_{1} e_{0}} + 1\right)}, \quad p{\left (r \right )} = \frac{e_{0} \left(24 e_{0} \pi r^{2} e^{C_{1} e_{0}} + 2 \sqrt{\left(8 e_{0} \pi r^{2} - 3\right) e^{C_{1} e_{0}}} - 9 e^{C_{1} e_{0}} - 1\right)}{3 \left(- 8 e_{0} \pi r^{2} e^{C_{1} e_{0}} + 3 e^{C_{1} e_{0}} + 1\right)}\right ]$$

Als Randbedingung verwenden wir die Eigenschaft, dass der Druck $p$ am Sternrand $r=R$ gleich null sein muss ($p(R)=0$). Für die zweite (positive) Lösung ergibt sich somit die folgende Bedingung:

EqC1=Eq(LoesTOV2b[1].subs(r,R).rhs,0)
EqC1

$$\frac{e_{0} \left(24 R^{2} e_{0} \pi e^{C_{1} e_{0}} + 2 \sqrt{\left(8 R^{2} e_{0} \pi - 3\right) e^{C_{1} e_{0}}} - 9 e^{C_{1} e_{0}} - 1\right)}{3 \left(- 8 R^{2} e_{0} \pi e^{C_{1} e_{0}} + 3 e^{C_{1} e_{0}} + 1\right)} = 0$$

Um den noch unbekannten Parameter $C_1$ zu bestimmen, lösen wir diese Gleichung nach $C_1$ auf und erhalten die folgende Lösung:

LoesC1=solve(EqC1,Symbol('C1'))
LoesC1[0]

$$\frac{\log{\left (\frac{1}{9 \left(8 R^{2} e_{0} \pi - 3\right)} \right )}}{e_{0}}$$

Wir setzen diese Lösung für $C_1$ ein und erhalten unsere analytische Druckverteilung im Körper: $p(r,e_0,R)$

LoesP=LoesTOV2b[1].subs(Symbol('C1'),LoesC1[0])
LoesP

$$p{\left (r \right )} = \frac{e_{0} \left(\frac{8 e_{0} \pi r^{2}}{3 \left(8 R^{2} e_{0} \pi - 3\right)} + \frac{2 \sqrt{\frac{8 e_{0} \pi r^{2} - 3}{8 R^{2} e_{0} \pi - 3}}}{3} - 1 - \frac{1}{8 R^{2} e_{0} \pi - 3}\right)}{3 \left(- \frac{8 e_{0} \pi r^{2}}{9 \left(8 R^{2} e_{0} \pi - 3\right)} + 1 + \frac{1}{3 \left(8 R^{2} e_{0} \pi - 3\right)}\right)}$$

Im Folgenden stellen wir uns die analytischen Lösungen für drei unterschiedliche konstante Energiedichten $e_0$ dar, wobei wir den Radius des Körpers auf 10 Kilometer festlegen.

import matplotlib.pyplot as plt 
import matplotlib
import matplotlib.gridspec as gridspec
params = {
    'figure.figsize'    : [16,5],
#    'text.usetex'       : True,
    'axes.titlesize' : 14,
    'axes.labelsize' : 16,  
    'xtick.labelsize' : 14 ,
    'ytick.labelsize' : 14 
}
matplotlib.rcParams.update(params)
Druckl=lambdify([r,(econst,R)], LoesP.rhs.subs(Pi,np.pi))
Massel=lambdify([r,econst], Loesm.rhs.subs(Pi,np.pi))

sete0=[0.0005,0.0006,0.0007]
setR=10
pts = np.linspace(0, setR, 10001)
setymax=sete0[-1]+0.0001
fig = plt.figure(figsize=(16,5))
gs = gridspec.GridSpec(1, 3, width_ratios=[1,1,1], wspace=0.40)
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1])
ax3 = plt.subplot(gs[2])

ax1.set_ylabel(r"$\rm Druck \,\,p$")
ax1.set_xlabel(r"$\rm Radius \,\,r$")
ax2.set_ylabel(r"$\rm Energie \,\,e $")
ax2.set_xlabel(r"$\rm Radius \,\,r$")
ax3.set_ylabel(r"$\rm Masse \,\,m$")
ax3.set_xlabel(r"$\rm Radius \,\,r$")
ax2.set_ylim(0,setymax)

ax1.plot(pts, Druckl(pts,(sete0[0],setR)), linewidth=2.5, c="blue")
ax2.plot(pts, [sete0[0]]*len(pts), linewidth=2.5, label=sete0[0], c="blue")
ax2.axvline(x=setR, ymax=sete0[0]/setymax, linewidth=2.5, c="blue")
ax3.plot(pts, Massel(pts,sete0[0]),linewidth=2.5, c="blue")

ax1.plot(pts, Druckl(pts,(sete0[1],setR)), linewidth=2.5, c="red")
ax2.plot(pts, [sete0[1]]*len(pts), linewidth=2.5, label=sete0[1], c="red")
ax2.axvline(x=setR, ymax=sete0[1]/setymax, linewidth=2.5, c="red")
ax3.plot(pts, Massel(pts,sete0[1]),linewidth=2.5, c="red")

ax1.plot(pts, Druckl(pts,(sete0[2],setR)), linewidth=2.5, c="green")
ax2.plot(pts, [sete0[2]]*len(pts), linewidth=2.5, label=sete0[2], c="green")
ax2.axvline(x=setR, ymax=sete0[2]/setymax, linewidth=2.5, c="green")
ax3.plot(pts, Massel(pts,sete0[2]),linewidth=2.5, c="green")
ax2.legend(loc='best',fontsize=16);

In der linken oberen Abbildung sind die drei Druckprofile der Sterne dargestellt. Man erkennt, dass sich bei einer äquidistanten Erhöhung der Energiedichte (siehe mittlere Abbildung) der zentrale Druck im Stern ($p_c=p(r=0)$) nicht-linear vergrößert. Die rechte Abbildung stellt die Massenfunktion $m(r)$ für die drei Sterne dar, wobei sich die Gesamtmasse $M$ der Sterne am Ende der Kurven, bei $M=m(R)=m(10)$ ablesen lässt.

Bei weiterer Erhöhung des Wertes der konstanten Energiedichte $e_0$ tritt eine interessante Eigenschaft auf. Um dies zu verdeutlichen, stellen wir den zentralen Druck im Stern ($p_c=p(r=0)$) als Funktion der Energiedichte $e_0$ dar:

sete0anf=0.0005
sete0end=0.00105
setR=10
epts = np.linspace(sete0anf, sete0end, 1000)
setymax=sete0[-1]+0.0001
fig = plt.figure(figsize=(15,5))
gs = gridspec.GridSpec(1, 2, width_ratios=[1,1], wspace=0.40)
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1])

ax1.set_ylabel(r"$\rm Zentraler \, Druck \,\,p_c=p(0)$")
ax1.set_xlabel(r"$\rm Energie \,\,e$")
ax2.set_ylabel(r"$\rm Gesamtmasse \,\,M$")
ax2.set_xlabel(r"$\rm Energie \,\,e$")

ax1.plot(epts, Druckl(0,(epts,setR)), linewidth=2.5, c="blue", label=r"$\rm R=10 \, km$")
ax2.plot(epts, Massel(setR,epts), linewidth=2.5, c="blue", label=r"$\rm R=10 \, km$")

ax1.legend(loc='best',fontsize=16);

Man erkennt in der oberen linken Abbildung, dass der zentrale Druck $p_c$ bei einem gewissen Wert der Energiedichte $e_0=e_{\rm crit}$ divergiert und gegen Unendlich strebt ($p_c(e_0 \rightarrow e_{\rm crit}) \rightarrow \infty $). Wir bestimmen in Python die möglichen kritischen Werte, die sich beim zentralen Wert der Druckfunktion $p(r=0,e_0,R)$ ergeben können (hier speziell für R=10):

Loesecrit=solve(LoesP.rhs.subs({(r,0),(R,setR),(Pi,np.pi)}), check=False)
Loesecrit

$$\left [ 0.0, \quad 0.00106103295394597, \quad 0.00119366207318922\right ]$$

solve(LoesP.rhs.subs({(r,0)}), econst, check=False)

$$\left [ 0, \quad \frac{1}{3 R^{2} \pi}, \quad \frac{3}{8 R^{2} \pi}\right ]$$

Der Wert $e_{\rm crit}=0.00106103295394597$ ist die kritische Energiedichte, bei der der Druck divergiert; allgemein liegt er bei $e_{\rm crit}=\frac{1}{3 \, \pi \, R^2}$. Würde man umgekehrt einen willkürlichen Wert der Energiedichte $e_0$ wählen, so würde man zu einem kritischen Radiuswert $R_{\rm crit}$ gelangen, unterhalb dessen der zentrale Druck divergieren würde ($R_{\rm crit}=\frac{9}{8} R_S = \frac{9}{8} \, 2 M= \frac{1}{\sqrt{3 \, \pi \, e_0}} $).

solve(LoesP.rhs.subs({(r,0)}), R, check=False)

$$\left [ 0, \quad - \frac{\sqrt{3}}{3 \sqrt{e_{0}} \sqrt{\pi}}, \quad \frac{\sqrt{3}}{3 \sqrt{e_{0}} \sqrt{\pi}}, \quad - \frac{\sqrt{6}}{4 \sqrt{e_{0}} \sqrt{\pi}}, \quad \frac{\sqrt{6}}{4 \sqrt{e_{0}} \sqrt{\pi}}\right ]$$

Daraus folgt, dass Körper die kleiner als $\frac{9}{8}$ des Schwarzschildradius sind, instabil sein müssen und man kann zeigen, dass diese Grenze eine absolute Gültigkeit hat und auch unabhängig von der Zustandsgleichung der Materie ist. Diese absolute Grenze der Stabilität wird oft als "Buchdahl Grenze" bezeichnet ( siehe Buchdahl, Hans A. "General relativistic fluid spheres." Physical Review 116.4 (1959): 1027. ). Die einzige Ausnahme sind die sogenannten "Gravastars", die auch kleiner als diese Grenze werden können. Diese besitzen jedoch einen Bereich im Stern, der eine negative Energiedichte aufweist.

Im Folgenden werden wir diese analytischen Resultate mit den numerischen Simulationen vergleichen. Hierzu lösen wir das System von Differentialgleichung bei konstanter Dichte numerisch und stellen uns die Eigenschaften der Sterne dar. Zunächst definieren wir wie im vorigen Jupyter Notebook (TOV.ipny) das System der Differentialgleichungen:

from scipy import integrate

y1, y2, y3 = symbols('y1, y2, y3')
TOV1c=TOV1a.subs([(Fm,y1),(p,y2),(Pi,np.pi)])
TOV2c=TOV2a.subs([(Fm,y1),(p,y2),(Pi,np.pi)])
TOV3c=TOV3.subs([(Fm,y1),(p,y2),(Fphi,y3),(Pi,np.pi)])

TOV1d=lambdify([r,(y1,y2,y3,econst)], TOV1c.rhs)
TOV2d=lambdify([r,(y1,y2,y3,econst)], TOV2c.rhs)
TOV3d=lambdify([r,(y1,y2,y3,econst)], TOV3c.rhs)

def DGLSysTOV(y, r, econst):
    y1, y2, y3 = y
    dy1 = TOV1d(r,(y1,y2,y3,econst))
    dy2 = TOV2d(r,(y1,y2,y3,econst))
    dy3 = TOV3d(r,(y1,y2,y3,econst))
    return np.array([dy1, dy2, dy3])

Wir wählen die Randbedingungen von drei konstanten Dichte Sternen ( $e_0=[0.0008, \, 0.0009, \, 0.00103]$ ), wobei der letzte Wert nur ein wenig kleiner als der kritische Wert ist ( $e_{\rm crit}=\frac{1}{3 \, \pi \, R^2}=\frac{1}{3 \, \pi \, 10^2} \approx 0.00106$ ).

points=1001
setR=10
sete0=[0.0008,0.0009,0.00103]
setp0=[float(LoesP.rhs.subs({(r,0),(R,setR),(Pi,np.pi),(econst,sete0[0])})),float(LoesP.rhs.subs({(r,0),(R,setR),(Pi,np.pi),(econst,sete0[1])})),float(LoesP.rhs.subs({(r,0),(R,setR),(Pi,np.pi),(econst,sete0[2])}))]
setM=[float(EqMasse.subs([(R,setR),(econst,sete0[0]),(Pi,np.pi)]).rhs),float(EqMasse.subs([(R,setR),(econst,sete0[1]),(Pi,np.pi)]).rhs),float(EqMasse.subs([(R,setR),(econst,sete0[2]),(Pi,np.pi)]).rhs)]
pts = np.linspace(10**(-14), setR, 1001)
ptsnum = np.linspace(10**(-14), setR, points)
ptsaussen = np.linspace(setR, 20, 300)

Randbedingungen0 = np.array([0, setp0[0], 0])
Randbedingungen1 = np.array([0, setp0[1], 0])
Randbedingungen2 = np.array([0, setp0[2], 0])

LoesTOV0 = integrate.odeint(DGLSysTOV, Randbedingungen0, ptsnum, args=(sete0[0],))
LoesTOV1 = integrate.odeint(DGLSysTOV, Randbedingungen1, ptsnum, args=(sete0[1],))
LoesTOV2 = integrate.odeint(DGLSysTOV, Randbedingungen2, ptsnum, args=(sete0[2],))

Man kann zeigen (siehe z.B. Fließbach-Buch, S: 243), dass die Metrik-Komponenten im Innenraum einer Kugel konstanter Dichte mit Radius $R$ wie folgt lauten:

$$ \begin{eqnarray} g_{tt} (r) &=& \, \frac{1}{4} \left( 3 \sqrt{1-\frac{2M}{R}} - \sqrt{1-\frac{2M \, r^2}{R^3}} \right)^2 \quad \forall \, r \leq R \\ g_{tt} (r) &=& 1-\frac{2M}{r} \quad \forall \, r > R \\ g_{rr} (r) &=& \left( 1-\frac{2M \, r^2}{R^3} \right)^{-1} \quad \forall \, r \leq R \\ g_{rr} (r) &=& \left( 1-\frac{2M}{r} \right)^{-1} \quad \forall \, r > R \\ \end{eqnarray} $$

In den folgenden Abbildungen veranschaulichen wir die analytischen Resultate mit den numerischen Simulationen:

gtt_analytisch0=(( 3*sqrt(1-2*setM[0]/setR) - sqrt(1-2*setM[0]*r**2/setR**3) )**2)/4 
grr_analytisch0=-1/(1-2*setM[0]*r**2/setR**3)
gtt_analytisch1=(( 3*sqrt(1-2*setM[1]/setR) - sqrt(1-2*setM[1]*r**2/setR**3) )**2)/4 
grr_analytisch1=-1/(1-2*setM[1]*r**2/setR**3)
gtt_analytisch2=(( 3*sqrt(1-2*setM[2]/setR) - sqrt(1-2*setM[2]*r**2/setR**3) )**2)/4 
grr_analytisch2=-1/(1-2*setM[2]*r**2/setR**3)

PhiSM=[ln((1-2*setM[0]/setR))/2,ln((1-2*setM[1]/setR))/2,ln((1-2*setM[2]/setR))/2]
shiftphi=[PhiSM[0]-LoesTOV0[-1, 2],PhiSM[1]-LoesTOV1[-1, 2],PhiSM[2]-LoesTOV2[-1, 2]]
j=0
LoesTOV0gtt=[]
LoesTOV1gtt=[]
LoesTOV2gtt=[]
while j<points:
    LoesTOV0gtt.append(exp(2*(LoesTOV0[j, 2]+shiftphi[0])))
    LoesTOV1gtt.append(exp(2*(LoesTOV1[j, 2]+shiftphi[1])))
    LoesTOV2gtt.append(exp(2*(LoesTOV2[j, 2]+shiftphi[2])))
    j=j+1
j=0
Loesgtt_analyt0=[]
Loesgtt_analyt1=[]
Loesgtt_analyt2=[]
Loesgrr_analyt0=[]
Loesgrr_analyt1=[]
Loesgrr_analyt2=[]
while j<1001:
    Loesgtt_analyt0.append(gtt_analytisch0.subs(r,pts[j]))
    Loesgtt_analyt1.append(gtt_analytisch1.subs(r,pts[j]))
    Loesgtt_analyt2.append(gtt_analytisch2.subs(r,pts[j]))
    Loesgrr_analyt0.append(grr_analytisch0.subs(r,pts[j]))
    Loesgrr_analyt1.append(grr_analytisch1.subs(r,pts[j]))
    Loesgrr_analyt2.append(grr_analytisch2.subs(r,pts[j]))
    j=j+1

fig = plt.figure(figsize=(15,11))
gs = gridspec.GridSpec(2, 2, width_ratios=[1,1], wspace=0.25, hspace=0.25)
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1])
ax3 = plt.subplot(gs[2])
ax4 = plt.subplot(gs[3])
setymax=sete0[-1]+0.0001

ax1.set_ylabel(r"$\rm Druck \,\,p$")
ax1.set_xlabel(r"$\rm Radius \,\,r$")
ax3.set_ylabel(r"$\rm Energie \,\,e $")
ax3.set_xlabel(r"$\rm Radius \,\,r$")

ax2.set_ylabel(r"$\rm g_{tt} $")
ax2.set_xlabel(r"$\rm Radius \,\,r$")
ax4.set_ylabel(r"$\rm g_{rr} $")
ax4.set_xlabel(r"$\rm Radius \,\,r$")
ax3.set_ylim(0,setymax)

ax1.plot(pts, Druckl(pts,(sete0[0],setR)), label='Analytisch', linewidth=2, c="blue")
ax1.plot(ptsnum, LoesTOV0[:, 1], linewidth=6, label='Numerisch', c="black", linestyle=":")
ax1.plot(pts, Druckl(pts,(sete0[1],setR)), linewidth=2, c="red")
ax1.plot(ptsnum, LoesTOV1[:, 1], linewidth=6, c="black", linestyle=":")
ax1.plot(pts, Druckl(pts,(sete0[2],setR)), linewidth=2, c="green")
ax1.plot(ptsnum, LoesTOV2[:, 1], linewidth=6, c="black", linestyle=":")

ax2.plot(ptsnum, LoesTOV0gtt, linewidth=6, c="black", linestyle=":", label='Numerisch')
ax2.plot(ptsnum, LoesTOV1gtt, linewidth=6, c="black", linestyle=":")
ax2.plot(ptsnum, LoesTOV2gtt, linewidth=6, c="black", linestyle=":")

ax2.plot(pts, Loesgtt_analyt0, linewidth=2, linestyle="-", c="blue", label='Innenraum')
ax2.plot(pts, Loesgtt_analyt1, linewidth=2, linestyle="-", c="red")
ax2.plot(pts, Loesgtt_analyt2, linewidth=2, linestyle="-", c="green")

ax2.plot(ptsaussen, (1-2*setM[0]/ptsaussen), linewidth=2.5, linestyle="dashdot", c="grey", label='Außenraum')
ax2.plot(ptsaussen, (1-2*setM[1]/ptsaussen), linewidth=2.5, linestyle="dashdot", c="grey")
ax2.plot(ptsaussen, (1-2*setM[2]/ptsaussen), linewidth=2.5, linestyle="dashdot", c="grey")

ax3.plot(pts, [sete0[0]]*len(pts), linewidth=2.5, label=sete0[0], c="blue")
ax3.axvline(x=setR, ymax=sete0[0]/setymax, linewidth=2.5, c="blue")
ax3.plot(pts, [sete0[1]]*len(pts), linewidth=2.5, label=sete0[1], c="red")
ax3.axvline(x=setR, ymax=sete0[1]/setymax, linewidth=2.5, c="red")
ax3.plot(pts, [sete0[2]]*len(pts), linewidth=2.5, label=sete0[2], c="green")
ax3.axvline(x=setR, ymax=sete0[2]/setymax, linewidth=2.5, c="green")

ax4.plot(ptsnum, -(1-2*LoesTOV0[:, 0]/ptsnum[:])**(-1), linewidth=6, c="black", linestyle=":", label='Numerisch')
ax4.plot(ptsnum, -(1-2*LoesTOV1[:, 0]/ptsnum[:])**(-1), linewidth=6, c="black", linestyle=":")
ax4.plot(ptsnum, -(1-2*LoesTOV2[:, 0]/ptsnum[:])**(-1), linewidth=6, c="black", linestyle=":")

ax4.plot(pts, Loesgrr_analyt0, linewidth=2, linestyle="-", c="blue", label='Innenraum')
ax4.plot(pts, Loesgrr_analyt1, linewidth=2, linestyle="-", c="red")
ax4.plot(pts, Loesgrr_analyt2, linewidth=2, linestyle="-", c="green")

ax4.plot(ptsaussen, -(1-2*setM[0]/ptsaussen[:])**(-1), linewidth=2.5, linestyle="dashdot", c="grey", label='Außenraum')
ax4.plot(ptsaussen, -(1-2*setM[1]/ptsaussen[:])**(-1), linewidth=2.5, linestyle="dashdot", c="grey")
ax4.plot(ptsaussen, -(1-2*setM[2]/ptsaussen[:])**(-1), linewidth=2.5, linestyle="dashdot", c="grey")

ax1.legend(loc='best',fontsize=16)
ax2.legend(loc='best',fontsize=16)
ax3.legend(loc='best',fontsize=16)
ax4.legend(loc='best',fontsize=16);

Die grünen Kurven kennzeichnen die Ergebnisse der Kugel mit einer konstanten Dichte, die nahe dem kritischen Wert $e_{\rm crit}$ ist. Der Druck im Zentrum des Sterns hat einen hohen Wert, die Metrik-Komponente $g_tt$ im Zentrum ist nahe null und die Metrik-Komponente $g_rr$ am Sternrand ist sehr klein.

Visualisierung der raumzeitlichen Struktur eines kompakten Sterns

Eingebettetes Diagramme der räumlichen Hypersphäre der Mannigfaltigkeit

Wir betrachten uns nun das eingebettete Diagramm der räumlichen Hypersphäre $\Sigma_t$ der Mannigfaltigkeit ${\cal M}$ die durch ein massives kompaktes Objekt (beschrieben durch die TOV-Gleichungen) gekrümmt ist. Aufgrund der Zeitunabhängigkeit der Metrik betrachtet man sich die Raumzeit zu einem festen, beliebigen Zeitpunkt, d.h. auf einer räumlichen Hyperfläche $\Sigma_t$ mit t=const. Durch die sphärische Symmetrie des Objektes sind alle Flächen mit $\theta= const$ gleichbedeutend, sodass wir den Raum in einem Schnitt durch die Äquatorialebene $\theta= \pi/2$ betrachten können. Das Weglängenelement schreibt sich dadurch wie folgt:

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

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

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))
ds2.subs({(dt,0),(dtheta,0),(theta,pi/2)})

$$- d\phi^{2} r^{2} - \frac{dr^{2}}{1 - \frac{2 m{\left (r \right )}}{r}}$$

Die durch das Weglängenelement beschriebene gekrümmte zweidimensionale Geometrie betten wir nun (zur Visualisierung) in einen dreidimensionalen euklidischen Raum ein, den wir durch zylindrische Koordinaten $(z, r, \phi)$ beschreiben. Das Weglängenelement in diesem flachen euklidischen Visualisierungsraum lautet:

$$ds^2=-dz^2 - dr^2 -r^2\, d\phi^2 \quad .$$

Wir setzen nun beide Weglängenelemente gleich und beschreiben dadurch die Abweichungen der gekrümmten Raumzeit von der flachen, euklidischen Geometrie durch die Koordinate $z=z(r)$, den sogenannten ’Lift’.

dz = symbols('dz')
Eq1=Eq(ds2.subs({(dt,0),(dtheta,0),(theta,pi/2)}),-dz**2 - dr**2 - r**2*dphi**2)
Eq1

$$- d\phi^{2} r^{2} - \frac{dr^{2}}{1 - \frac{2 m{\left (r \right )}}{r}} = - d\phi^{2} r^{2} - dr^{2} - dz^{2}$$

Diese Gleichung lösen wir nach dz auf

LoesTOV4=solve(Eq1,dz)
LoesTOV4

$$\left [ - \sqrt{2} dr \sqrt{\frac{m{\left (r \right )}}{r - 2 m{\left (r \right )}}}, \quad \sqrt{2} dr \sqrt{\frac{m{\left (r \right )}}{r - 2 m{\left (r \right )}}}\right ]$$

Wir betrachten uns die positive Lösung und erhalten somit eine neue Differentialgleichung für die Funktion $z(r)$

$$ \frac{dz}{dr} ~=~ \sqrt{\frac{2 \, m(r)}{r - 2 \, m(r)}} . $$

Im Fall der äußeren Schwarzschildmetrik (siehe ...) erhielten wir eine ähnliche Form der Differentialgleichung

$$ \frac{dz}{dr} ~=~ \sqrt{\frac{2 \, M}{r - 2 \, M}} , $$

die wir dann über den Radius dr integrierten. Wir erhielten $z(r) = \sqrt{8 M \, \left( r - 2M \right)} $, was uns zum Raumzeit-Trichter eines schwarzen Lochs führte.

Im Folgenden betrachten wir (wie im Jupyter Notebook TOV.ipny ) eine polytrope Zustandsgleichung der Form $p(e)=K\,e^{\gamma}$ (hier speziell $p(e)=10\,e^{5/3}$) und lösen die TOV-Gleichungen mit zusätzlicher Einbeziehung der neuen Differentialgleichung für die Funktion $z(r)$ numerisch. Die Randbedingungen für $z(r)$ müssen dann (ähnlich wie bei der Metrik $g_{tt}$) stetig in die Außenraummetrik übergehen.

from scipy import integrate

setK=10
setgamma=5.0/3.0
Fz = Function('z', real = True)(r)

TOV1a=TOV1.subs(e,(p/setK)**(1/setgamma))
TOV2a=TOV2.subs(e,(p/setK)**(1/setgamma))
TOV3a=TOV3
TOV4a=Eq(Fz.diff(r),LoesTOV4[1]/dr)

y4 = symbols('y4')
TOV1b=TOV1a.subs([(Fm,y1),(p,y2),(Pi,np.pi)])
TOV2b=TOV2a.subs([(Fm,y1),(p,y2),(Pi,np.pi)])
TOV3b=TOV3a.subs([(Fm,y1),(p,y2),(Fphi,y3),(Pi,np.pi)])
TOV4b=TOV4a.subs([(Fz,y4),(Fm,y1),(Pi,np.pi)])

TOV1c=lambdify([r,(y1,y2,y3,y4)], TOV1b.rhs)
TOV2c=lambdify([r,(y1,y2,y3,y4)], TOV2b.rhs)
TOV3c=lambdify([r,(y1,y2,y3,y4)], TOV3b.rhs)
TOV4c=lambdify([r,(y1,y2,y3,y4)], TOV4b.rhs)

def DGLSysTOV(y, r):
    y1, y2, y3 , y4= y
    dy1 = TOV1c(r,(y1,y2,y3,y4))
    dy2 = TOV2c(r,(y1,y2,y3,y4))
    dy3 = TOV3c(r,(y1,y2,y3,y4))
    dy4 = TOV4c(r,(y1,y2,y3,y4))
    return np.array([dy1, dy2, dy3, dy4])

Wir berechnen nun die numerischen Lösungen für zwei Neutronensterne mit unterschiedlicher zentraler Energiedichte $e_c=0.0004$ und $e_c=0.002$ und stellen uns die Eigenschaften der Sterne im Innen- und Außenraum dar:

pts = np.linspace(10**(-14), 30, 5001)
Randbedingungen1 = np.array([0, setK*(0.0004)**setgamma, 0, 0])
Randbedingungen2 = np.array([0, setK*(0.002)**setgamma, 0, 0])

LoesTOV1 = integrate.odeint(DGLSysTOV, Randbedingungen1, pts)
LoesTOV2 = integrate.odeint(DGLSysTOV, Randbedingungen2, pts)

<string>:3: RuntimeWarning: invalid value encountered in double_scalars
<string>:3: RuntimeWarning: invalid value encountered in double_scalars

j=0
while LoesTOV1[j, 1]>0 and np.isnan(LoesTOV1[j, 1])==False:
    j=j+1
itend1=j-1
R1=pts[itend1]

M1=LoesTOV1[itend1, 0]

PhiSM=ln((1-2*M1/R1))/2
shiftphi=PhiSM-LoesTOV1[itend1, 2]
j=0
LoesTOVgtt1=[]
while j<itend1:
    LoesTOVgtt1.append(exp(2*(LoesTOV1[j, 2]+shiftphi)))
    j=j+1

shiftz1=np.sqrt(8)*np.sqrt(M1*(pts[itend1] - 2*M1))-LoesTOV1[itend1, 3]

j=0
while LoesTOV2[j, 1]>0 and np.isnan(LoesTOV2[j, 1])==False:
    j=j+1
itend2=j-1
R2=pts[itend2]

M2=LoesTOV2[itend2, 0]

PhiSM=ln((1-2*M2/R2))/2
shiftphi=PhiSM-LoesTOV2[itend2, 2]
j=0
LoesTOVgtt2=[]
while j<itend2:
    LoesTOVgtt2.append(exp(2*(LoesTOV2[j, 2]+shiftphi)))
    j=j+1

shiftz2=np.sqrt(8)*np.sqrt(M2*(pts[itend2] - 2*M2))-LoesTOV2[itend2, 3]

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

ax1.set_ylabel(r"$\rm Einbettung \,\,z$")
ax1.set_xlabel(r"$\rm Radius \,\,r$")
ax2.set_ylabel(r"$\rm Energie \,\,e $")
ax2.set_xlabel(r"$\rm Radius \,\,r$")
ax3.set_ylabel(r"$\rm g_{tt} $")
ax3.set_xlabel(r"$\rm Radius \,\,r$")
ax4.set_ylabel(r"$\rm g_{rr} $")
ax4.set_xlabel(r"$\rm Radius \,\,r$")

ax1.plot(pts[:itend1], LoesTOV1[:itend1, 3]+shiftz1, linewidth=2.5, label='z(r) Innenraum', c="blue")
ax1.plot(pts[itend1:], np.sqrt(8)*np.sqrt(M1*(pts[itend1:] - 2*M1)), linewidth=2.5, linestyle=":", label='$z(r)$ Außenraum (Schwarzschild)', c="red")

ax2.plot(pts[:itend1], (LoesTOV1[:itend1, 1]/setK)**(1/setgamma), linewidth=2.5, label='e(r), $e_c=0.0004$', c="blue")
ax3.plot(pts[:itend1], LoesTOVgtt1, linewidth=2.5, linestyle="-", label='$g_{tt}$ Innenraum (TOV)', c="blue")
ax3.plot(pts[itend1:], (1-2*M1/pts[itend1:]), linewidth=2.5, linestyle=":", label='$g_{tt}$ Außenraum', c="red")
ax4.plot(pts[:itend1], -(1-2*LoesTOV1[:itend1, 0]/pts[:itend1])**(-1), linewidth=2.5, linestyle="-", label='$g_{rr}$ Innenraum (TOV)', c="blue")
ax4.plot(pts[itend1:], -(1-2*M1/pts[itend1:])**(-1), linewidth=2.5, linestyle=":", label='$g_{rr}$ Außenraum', c="red")

ax1.plot(pts[:itend2], LoesTOV2[:itend2, 3]+shiftz2, linewidth=2.5, c="green")
ax1.plot(pts[itend2:], np.sqrt(8)*np.sqrt(M2*(pts[itend2:] - 2*M2)), linewidth=2.5, linestyle=":", c="red")

ax2.plot(pts[:itend2], (LoesTOV2[:itend2, 1]/setK)**(1/setgamma), linewidth=2.5, label='e(r), $e_c=0.002$', c="green")
ax3.plot(pts[:itend2], LoesTOVgtt2, linewidth=2.5, linestyle="-", c="green")
ax3.plot(pts[itend2:], (1-2*M2/pts[itend2:]), linewidth=2.5, linestyle=":", c="red")
ax4.plot(pts[:itend2], -(1-2*LoesTOV2[:itend2, 0]/pts[:itend2])**(-1), linewidth=2.5, linestyle="-", c="green")
ax4.plot(pts[itend2:], -(1-2*M2/pts[itend2:])**(-1), linewidth=2.5, linestyle=":", c="red")




ax1.legend(loc='best',fontsize=16)
ax2.legend(loc='best',fontsize=16)
ax3.legend(loc='best',fontsize=16)
ax4.legend(loc='best',fontsize=16);

Die obere, linke Abbildung zeigt das eingebettete Diagramm der räumlichen Hypersphäre der Mannigfaltigkeit für die beiden Sterne. Wir visualisieren uns im Folgenden das Diagramm des grünen Sterns in einer dreidimensionalen Darstellung und vergleichen es mit dem "Raumzeit-Trichter" eines schwarzen Loches mit gleicher gravitativen Masse wie der Neutronenstern.

from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm

Rend=20
x,y = symbols('x,y')
zBett=sqrt(8*M2*(r-2*M2))
Rl = np.linspace(2.00001, Rend, 300)
PHIl = np.linspace(0, 2*np.pi, 300)
Rl, PHIl = np.meshgrid(Rl, PHIl)
X = Rl*np.cos(PHIl)
Y = Rl*np.sin(PHIl)
zlam = lambdify((x,y), zBett.subs({(r,sqrt(x**2+y**2))}))
ZF = np.array(list(np.real(zlam(X,Y))))

PHIl2 = np.linspace(0, 2*np.pi, 300)
Rl2, PHIl2a = np.meshgrid(pts[:itend2], PHIl2)
ZF2, PHIl2 = np.meshgrid(LoesTOV2[:itend2, 3]+shiftz2, PHIl2)
X2 = Rl2*np.cos(PHIl2)
Y2 = Rl2*np.sin(PHIl2)

Rl2a = np.linspace(pts[itend2], Rend, 300)
PHIl2a = np.linspace(0, 2*np.pi, 300)
Rl2a, PHIl2a = np.meshgrid(Rl2a, PHIl2a)
X2a = Rl2a*np.cos(PHIl2a)
Y2a = Rl2a*np.sin(PHIl2a)
ZF2a = np.array(list(np.real(zlam(X2a,Y2a))))

Rl2b = np.linspace(pts[itend2]-0.2, pts[itend2]+0.2, 30)
PHIl2b = np.linspace(0, 2*np.pi, 300)
Rl2b, PHIl2b = np.meshgrid(Rl2b, PHIl2b)
X2b = Rl2b*np.cos(PHIl2b)
Y2b = Rl2b*np.sin(PHIl2b)
ZF2b = np.array(list(np.real(zlam(X2b,Y2b))))

zend=15
fig = plt.figure(figsize=(14,6))
gs = gridspec.GridSpec(1, 2, width_ratios=[1,1], wspace=0.21)
ax1 = plt.subplot(gs[0],projection='3d')
ax2 = plt.subplot(gs[1],projection='3d')
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')
ax2.xaxis.pane.fill = False
ax2.yaxis.pane.fill = False
ax2.zaxis.pane.fill = False
ax2.xaxis.pane.set_edgecolor('white')
ax2.yaxis.pane.set_edgecolor('white')
ax2.zaxis.pane.set_edgecolor('white')
#ax1.grid(False)
#ax2.grid(False)
ax1.set_zlim(0, zend)
ax2.set_zlim(0, zend)

ax1.set_zlabel("Visualisierung z")
ax2.set_zlabel("Visualisierung z")
ax1.set_ylabel(r"$\rm y \,[km]$")
ax1.set_xlabel(r"$\rm x \,[km]$")
ax2.set_ylabel(r"$\rm y \,[km]$")
ax2.set_xlabel(r"$\rm x \,[km]$")

ax1.view_init(azim=35, elev=20)
ax2.view_init(azim=35, elev=20)

ax1.plot_surface(X2, Y2, ZF2, cmap=cm.nipy_spectral, linewidth=0, alpha=0.9, vmin=0, vmax=zend)
ax1.plot_surface(X2a, Y2a, ZF2a, cmap=cm.nipy_spectral, linewidth=0, alpha=0.9, vmin=0, vmax=zend)
ax1.plot_surface(X2b, Y2b, ZF2b, color="black")

ax2.plot_surface(X, Y, ZF, cmap=cm.nipy_spectral, linewidth=0, alpha=0.9, vmin=0, vmax=15);

Die linke Abbildung zeigt das eingebettete Diagramm der räumlichen Hypersphäre des Neutronensterns, wobei die schwarze Linie die Oberfläche des Sterns und somit die Grenze zwischen der Innenraum und Außenraum Metrik kennzeichnet. Die rechte Abbildung stellt die räumliche Hypersphäre eines schwarzen Lochs mit gleicher Masse dar.

Die Visualisierung kann man auch wieder mittels der interaktiven "Plotly Python Open Source Graphing Library" erstellen.

import plotly.graph_objects as go

zanf=min(LoesTOV2[:itend2, 3]+shiftz1)
Innenraum = go.Surface(x=X2, y=Y2, z=ZF2, colorscale="spectral_r", name="Innenraum", cmin=zanf, cmax=zend)
Aussenraum = go.Surface(x=X2a, y=Y2a, z=ZF2a, colorscale="spectral_r", name="Aussenraum", cmin=zanf, cmax=zend)
Surface = go.Surface(x=X2b, y=Y2b, z=ZF2b, colorscale="greys", name="Surface", showscale=False)

layout = go.Layout(autosize=True,
                  width=800, height=800,
                  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",
                    zaxis_range=[0,zend],
                    xaxis_title='x',
                    yaxis_title='y',
                    zaxis_title='f'),
                  legend=dict(
                    yanchor="top",
                    y=0.99,
                    xanchor="left",
                    x=0.01)
            )

data=[Innenraum,Aussenraum,Surface]
fig=go.Figure(data=data, layout=layout)

fig.show()

Numerisches Lösen der TOV-Gleichungen mittels des Euler-Verfahrens

In diesem Unterpunkt befassen wir uns mit den Grundlagen der numerischen Lösung der TOV-Gleichungen. Mit Ausnahme der analytischen Lösung des "Gravitationsfeldes einer Kugel aus inkompressibler Flüssigkeit (konstante Dichte)" hatten wir die TOV-Gleichungen stets mit der im Modul "scipy" vordefinierten Funktion "odeint()" numerisch gelöst.

Im Folgenden werden wir die TOV-Gleichungen mittels des einfachen "Euler Verfahrens" lösen. In diesem Verfahren, welches man im Prinzip auf alle erste Ordnung Differentialgleichungssysteme anwenden kann, schreibt man die TOV-Gleichungen wie folgt um

$$ \begin{eqnarray} \frac{dm}{dr} ~=~ 4\pi r^2 e ~, \quad &\rightarrow& \quad dm ~=~ 4\pi r^2 e \, dr \\ \frac{dp}{dr} ~=~ -(e+p) \, \frac{m+4 \pi r^3 p}{r \left(r -2m \right)} \quad &\rightarrow& \quad dp ~=~ -(e+p) \, \frac{m+4 \pi r^3 p}{r \left(r -2m \right)} \, dr \\ \frac{d\Phi}{dr} ~=~ \frac{m+4 \pi r^3 p}{r \left(r -2m \right)} \quad &\rightarrow& \quad d\Phi ~=~ \frac{m+4 \pi r^3 p}{r \left(r -2m \right)} \, dr \quad , \end{eqnarray} $$

und geht nach folgendem Schema vor:

• Man definiert die Zustandsgleichung der Sternmaterie als eine Funktion e(p).
• Man startet im Sternzentrum und legt den Wert der zentralen Energiedichte und des zentralen Druckes fest. Da die TOV Gleichung bei $r=0$ singulär sind, wählt man hier einen sehr, sehr kleinen Wert (z.B. $r=10^{-14}$).
• Die TOV Gleichungen werden als Differenzengleichungen umgeschrieben und eine kleine Schrittweite $dr=\Delta r << 1$ wird festgelegt. In einer Schleife wird dann in jedem Radius-Schritt die Druck-, Massen- und $\Phi$-Änderungen berechnet ($\Delta p$, $\Delta m$ und $\Delta \Phi$) und die jeweiligen Größen beim nächsten Schritt um diesen Faktor erhöht bzw. verringert:
• Im Laufe der iterativen Lösung verringert sich der Druck ständig. Die Schleife wird so lange ausgeführt bis der Wert des Druckes gleich null bzw. negativ wird, da an der Sternoberfläche der Druck verschwindet.

Im Folgenden werden wir das Eulerverfahren in Python implementieren. Wir benutzen wieder die schon im vorigen Unterpunkt verwendete polytrope Zustandsgleichung der Form $p(e)=K\,e^{\gamma}$ (hier speziell $p(e)=10\,e^{5/3}$):

del(e,p,r,M,R,dt,dr,dphi)
# Definition der Zustandsgleichung
setK=10
setgamma=5/3
def eos(p):
    e=(p/setK)**(1/setgamma)
    return e

Wir implementieren nun das Euler-Verfahren und berechnen die TOV-Lösung eines Sterns mit zentraler Energiedichte $e_c=0.0004$. Als Radius-Schrittweite wählen wir $\Delta r=0.00001$ km und wir speichern uns die Energiedichten- und Druckprofile im Stern alle 10 Meter ( dr_plot=0.01 ).

M=0
Phi=0
r=10**(-14)
e=0.0004
p=setK*e**setgamma
dr=0.00001

dr_plot=0.01
r_plot=dr_plot

LoesE_r=[]
LoesE_e=[]
LoesE_p=[]
LoesE_r.append(r)
LoesE_e.append(e)
LoesE_p.append(p)

while p > 0:
    e=eos(p)                                        #Wert der Energiedichte bei momentanen Druck
    dM=4*np.pi*e*r*r*dr                             #Massenzunahme bei momentanem r und Schrittweite dr
    dp=-(p+e)*(M+4*np.pi*r*r*r*p)/(r*(r-2*M))*dr    #Druckzunahme bei momentanem r und Schrittweite dr (TOV-Gleichung)
    dPhi=(M + 4*np.pi*r*r*r*p)/(r*(r-2*M))*dr       #Metrikzunahme bei momentanem r und Schrittweite dr
    r=r+dr                                          #Momentaner Radius des Neutronensterns
    M=M+dM                                          #Momentane Masse des Neutronensterns innerhalb des Radius r
    p=p+dp                                          #Momentaner Druck des Neutronensterns beim Radius r
    Phi=Phi+dPhi                                    #Momentane Metrik-Komponente des Neutronensterns beim Radius r
    if r > r_plot:
        LoesE_r.append(r)
        LoesE_e.append(e)
        LoesE_p.append(p)
        r_plot=r_plot+dr_plot
R=r

In der folgenden Abbildung vergleichen wir die mittels des Euler Verfahrens berechneten Druck und Energiedichte Eigenschaften des Sterns (schwarze Punkte) mit den im vorigen Unterpunkt berechneten Eigenschaften (mittels odeint() berechnet, blaue Kurve).

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

ax1.set_ylabel(r"$\rm Druck \,\,p$")
ax1.set_xlabel(r"$\rm Radius \,\,r$")
ax2.set_ylabel(r"$\rm Energie \,\,e $")
ax2.set_xlabel(r"$\rm Radius \,\,r$")

ax1.plot(pts[:itend1], LoesTOV1[:itend1, 1], linewidth=2.5, label='Mit odeint()', c="blue")
ax1.plot(LoesE_r, LoesE_p, linewidth=6, c="black", linestyle=":", label='Euler Verfahren')

ax2.plot(pts[:itend1], (LoesTOV1[:itend1, 1]/setK)**(1/setgamma), linewidth=2.5, label='Mit odeint()', c="blue")
ax2.plot(LoesE_r, LoesE_e, linewidth=6, c="black", linestyle=":", label='Euler Verfahren')

ax1.legend(loc='best',fontsize=16)
ax2.legend(loc='best',fontsize=16);

Die beiden Abbildungen zeigen, dass man auch mittels des einfachen Euler Verfahrens eine gute Beschreibung der Sterneigenschaften erhält. Es ergeben sich, bei der gewählten Radius-Schrittweite, nur kleine Abweichungen zu den mittels "odeint" berechneten Werten.

Abweichungen der Gesamtmasse M: (Euler Verfahren, "odeint")

M, M1

$$\left ( 1.340626975352493, \quad 1.3406179288991404\right )$$

Abweichungen des Sternradius R: (Euler Verfahren, "odeint")

R, R1

$$\left ( 16.993719999525837, \quad 16.926000000000002\right )$$