SM_Tolman-Oppenheimer-Volkoff_equations.sagews

%auto
typeset_mode(True, display=False)

Tolman-Oppenheimer-Volkoff equations

This worksheet illustrates some features of SageManifolds (v0.8) on the derivation of the Tolman-Oppenheimer-Volkoff equations (spherically symmetric, stationary solution of general relativity).

We will calculate the Einstein equations $$R_{\mu\nu} - \frac{1}{2}Rg_{\mu\nu} = T_{\mu\nu}$$ for a corresponding spherically symmetric, stationary metric $g$. In the above, $R_{\mu\nu}$ is the Ricci tensor, $R=R^\mu_\mu$ is the Ricci scalar, and $T_{\mu\nu}$ is the energy-momentum tensor (left side of Einstein's equations describe the geometry of spacetime, and the right side the matter in the spacetime).

We first declare the spacetime Mas a general 4­dimensional manifold

M = Manifold(4, 'M', r'\mathcal{M}')
print M ; M

4-dimensional manifold 'M'

$\mathcal{M}$

with the standard spherical coordinates (X denotes the coordinate chart on M):

X.<t,r,th,ph> = M.chart(r't r:(0,+oo) th:(0,pi):\theta phi:(0,2*pi):\phi')

In order to define a general spherically symmetric, stationary metric one needs a few auxiliary functions of the radial coordinate $r$ - metric functions $\nu(r)$ and $\lambda(r)$, matter pressure $p(r)$ and energy density $\rho(r)$, as well as the mass $m(r)$ enclosed within the sphere of the radius $r$:

# metric functions:
nu = function("nu", r, latex_name=r"\nu")
lam = function("lambda", r, latex_name=r"\lambda")
# density and pressure:
rho = function("rho", r, latex_name=r"\rho")
p = function("P", r)
# mass enclosed in sphere of radius r:
m = function("m", r)

In general, such metric reads as follows:

g = M.lorentz_metric('g')
g[0,0] = -exp(2*nu)
g[1,1] = exp(2*lam)
g[2,2], g[3,3] = r^2, r^2*sin(th)^2
g.display()

$g = -e^{\left(2 \, \nu\left(r\right)\right)} \mathrm{d} t\otimes \mathrm{d} t + e^{\left(2 \, \lambda\left(r\right)\right)} \mathrm{d} r\otimes \mathrm{d} r + r^{2} \mathrm{d} {\theta}\otimes \mathrm{d} {\theta} + r^{2} \sin\left({\theta}\right)^{2} \mathrm{d} {\phi}\otimes \mathrm{d} {\phi}$

which works assuming that the physical constants G=c=1. Let's introduce G and c as variables to obtain the dimensional version of the equations:

var('G c pi'); assume(G>0); assume(c>0)

($G$, $c$, $\pi$)

From the Newtonian weak field limit considerations (Newtonian force far from the source) one may simplify the above expression and replace $\lambda(r)$ with $\frac{1-2Gm}{rc^2}$, as well as explicitly put $c^2$ in front of $g_{tt}$. Then

g[0,0] = -c^2*exp(2*nu)
g[1,1] = 1/(1-2*G*m/(r*c^2))
g.display()

$g = -c^{2} e^{\left(2 \, \nu\left(r\right)\right)} \mathrm{d} t\otimes \mathrm{d} t + \left( \frac{c^{2} r}{c^{2} r - 2 \, G m\left(r\right)} \right) \mathrm{d} r\otimes \mathrm{d} r + r^{2} \mathrm{d} {\theta}\otimes \mathrm{d} {\theta} + r^{2} \sin\left({\theta}\right)^{2} \mathrm{d} {\phi}\otimes \mathrm{d} {\phi}$

The Ricci tensor is a result of a method ricci() acting on the metric g:

Ricci = g.ricci(); Ricci.display()

$\mathrm{Ric}\left(g\right) = \left( \frac{c^{2} r^{2} e^{\left(2 \, \nu\left(r\right)\right)} \left(\frac{\partial\,\nu}{\partial r}\right)^{2} + c^{2} r^{2} e^{\left(2 \, \nu\left(r\right)\right)} \frac{\partial^2\,\nu}{\partial r^2} + 2 \, c^{2} r e^{\left(2 \, \nu\left(r\right)\right)} \frac{\partial\,\nu}{\partial r} - {\left(2 \, r e^{\left(2 \, \nu\left(r\right)\right)} m\left(r\right) \left(\frac{\partial\,\nu}{\partial r}\right)^{2} + 2 \, r e^{\left(2 \, \nu\left(r\right)\right)} m\left(r\right) \frac{\partial^2\,\nu}{\partial r^2} + {\left(r e^{\left(2 \, \nu\left(r\right)\right)} \frac{\partial\,m}{\partial r} + 3 \, e^{\left(2 \, \nu\left(r\right)\right)} m\left(r\right)\right)} \frac{\partial\,\nu}{\partial r}\right)} G}{r^{2}} \right) \mathrm{d} t\otimes \mathrm{d} t + \left( -\frac{c^{2} r^{3} \left(\frac{\partial\,\nu}{\partial r}\right)^{2} + c^{2} r^{3} \frac{\partial^2\,\nu}{\partial r^2} - {\left(2 \, r^{2} m\left(r\right) \left(\frac{\partial\,\nu}{\partial r}\right)^{2} + 2 \, r^{2} m\left(r\right) \frac{\partial^2\,\nu}{\partial r^2} + 2 \, r \frac{\partial\,m}{\partial r} + {\left(r^{2} \frac{\partial\,m}{\partial r} - r m\left(r\right)\right)} \frac{\partial\,\nu}{\partial r} - 2 \, m\left(r\right)\right)} G}{c^{2} r^{3} - 2 \, G r^{2} m\left(r\right)} \right) \mathrm{d} r\otimes \mathrm{d} r + \left( -\frac{c^{2} r^{2} \frac{\partial\,\nu}{\partial r} - {\left(2 \, r m\left(r\right) \frac{\partial\,\nu}{\partial r} + r \frac{\partial\,m}{\partial r} + m\left(r\right)\right)} G}{c^{2} r} \right) \mathrm{d} {\theta}\otimes \mathrm{d} {\theta} + \left( -\frac{c^{2} r^{2} \sin\left({\theta}\right)^{2} \frac{\partial\,\nu}{\partial r} - {\left(2 \, r m\left(r\right) \frac{\partial\,\nu}{\partial r} + r \frac{\partial\,m}{\partial r} + m\left(r\right)\right)} G \sin\left({\theta}\right)^{2}}{c^{2} r} \right) \mathrm{d} {\phi}\otimes \mathrm{d} {\phi}$

For example, the $R_{tt}$ component is

Ricci[0,0]

$\frac{c^{2} r^{2} e^{\left(2 \, \nu\left(r\right)\right)} \left(\frac{\partial\,\nu}{\partial r}\right)^{2} + c^{2} r^{2} e^{\left(2 \, \nu\left(r\right)\right)} \frac{\partial^2\,\nu}{\partial r^2} + 2 \, c^{2} r e^{\left(2 \, \nu\left(r\right)\right)} \frac{\partial\,\nu}{\partial r} - {\left(2 \, r e^{\left(2 \, \nu\left(r\right)\right)} m\left(r\right) \left(\frac{\partial\,\nu}{\partial r}\right)^{2} + 2 \, r e^{\left(2 \, \nu\left(r\right)\right)} m\left(r\right) \frac{\partial^2\,\nu}{\partial r^2} + {\left(r e^{\left(2 \, \nu\left(r\right)\right)} \frac{\partial\,m}{\partial r} + 3 \, e^{\left(2 \, \nu\left(r\right)\right)} m\left(r\right)\right)} \frac{\partial\,\nu}{\partial r}\right)} G}{r^{2}}$

Ricci[0,0].expand().collect(nu.diff(r)).collect(nu.diff(r,r)).collect(c^2*exp(2*nu)).collect_common_factors()

${\left(c^{2} \left(\frac{\partial\,\nu}{\partial r}\right)^{2} - \frac{2 \, G m\left(r\right) \left(\frac{\partial\,\nu}{\partial r}\right)^{2}}{r} + c^{2} \frac{\partial^2\,\nu}{\partial r^2} + \frac{2 \, c^{2} \frac{\partial\,\nu}{\partial r}}{r} - \frac{G \frac{\partial\,m}{\partial r} \frac{\partial\,\nu}{\partial r}}{r} - \frac{2 \, G m\left(r\right) \frac{\partial^2\,\nu}{\partial r^2}}{r} - \frac{3 \, G m\left(r\right) \frac{\partial\,\nu}{\partial r}}{r^{2}}\right)} e^{\left(2 \, \nu\left(r\right)\right)}$

Ricci[1,1].expand().collect_common_factors()

$-\frac{c^{2} r^{3} \left(\frac{\partial\,\nu}{\partial r}\right)^{2} - 2 \, G r^{2} m\left(r\right) \left(\frac{\partial\,\nu}{\partial r}\right)^{2} + c^{2} r^{3} \frac{\partial^2\,\nu}{\partial r^2} - G r^{2} \frac{\partial\,m}{\partial r} \frac{\partial\,\nu}{\partial r} - 2 \, G r^{2} m\left(r\right) \frac{\partial^2\,\nu}{\partial r^2} + G r m\left(r\right) \frac{\partial\,\nu}{\partial r} - 2 \, G r \frac{\partial\,m}{\partial r} + 2 \, G m\left(r\right)}{{\left(c^{2} r - 2 \, G m\left(r\right)\right)} r^{2}}$

Ricci[2,2].expand().collect(nu.diff(r)).collect(nu.diff(r,r)).collect(c^2*exp(2*nu))

$-r \frac{\partial\,\nu}{\partial r} + \frac{2 \, G m\left(r\right) \frac{\partial\,\nu}{\partial r}}{c^{2}} + \frac{G \frac{\partial\,m}{\partial r}}{c^{2}} + \frac{G m\left(r\right)}{c^{2} r}$

Ricci scalar is obtained by the ricci_scalar() method acting on g:

Ric_scalar = g.ricci_scalar()
(Ric_scalar.function_chart(X)).collect(nu.diff(r)).collect(nu.diff(r,r)).collect(c^2*exp(2*nu))

$-\frac{2 \, {\left(c^{2} r^{2} \left(\frac{\partial\,\nu}{\partial r}\right)^{2} + c^{2} r^{2} \frac{\partial^2\,\nu}{\partial r^2} + 2 \, c^{2} r \frac{\partial\,\nu}{\partial r} - {\left(2 \, r m\left(r\right) \left(\frac{\partial\,\nu}{\partial r}\right)^{2} + 2 \, r m\left(r\right) \frac{\partial^2\,\nu}{\partial r^2} + {\left(r \frac{\partial\,m}{\partial r} + 3 \, m\left(r\right)\right)} \frac{\partial\,\nu}{\partial r} + 2 \, \frac{\partial\,m}{\partial r}\right)} G\right)}}{c^{2} r^{2}}$

It is the trace of the Ricci tensor, $R = R_\mu^\mu$:

Ric_scalar == Ricci.up(g, 1).trace(0, 1)

$\mathrm{True}$

Left side of the Einstein equations is

E = Ricci - (Ric_scalar*g)/2; E.display()

$\mathrm{Ric}\left(g\right)-\mbox{unnamed metric} = \frac{2 \, G e^{\left(2 \, \nu\left(r\right)\right)} \frac{\partial\,m}{\partial r}}{r^{2}} \mathrm{d} t\otimes \mathrm{d} t + \left( \frac{2 \, {\left(c^{2} r^{2} \frac{\partial\,\nu}{\partial r} - {\left(2 \, r m\left(r\right) \frac{\partial\,\nu}{\partial r} + m\left(r\right)\right)} G\right)}}{c^{2} r^{3} - 2 \, G r^{2} m\left(r\right)} \right) \mathrm{d} r\otimes \mathrm{d} r + \left( \frac{c^{2} r^{3} \left(\frac{\partial\,\nu}{\partial r}\right)^{2} + c^{2} r^{3} \frac{\partial^2\,\nu}{\partial r^2} + c^{2} r^{2} \frac{\partial\,\nu}{\partial r} - {\left(2 \, r^{2} m\left(r\right) \left(\frac{\partial\,\nu}{\partial r}\right)^{2} + 2 \, r^{2} m\left(r\right) \frac{\partial^2\,\nu}{\partial r^2} + r \frac{\partial\,m}{\partial r} + {\left(r^{2} \frac{\partial\,m}{\partial r} + r m\left(r\right)\right)} \frac{\partial\,\nu}{\partial r} - m\left(r\right)\right)} G}{c^{2} r} \right) \mathrm{d} {\theta}\otimes \mathrm{d} {\theta} + \left( -\frac{{\left(2 \, r^{2} m\left(r\right) \left(\frac{\partial\,\nu}{\partial r}\right)^{2} + 2 \, r^{2} m\left(r\right) \frac{\partial^2\,\nu}{\partial r^2} + r \frac{\partial\,m}{\partial r} + {\left(r^{2} \frac{\partial\,m}{\partial r} + r m\left(r\right)\right)} \frac{\partial\,\nu}{\partial r} - m\left(r\right)\right)} G \sin\left({\theta}\right)^{2} - {\left(c^{2} r^{3} \left(\frac{\partial\,\nu}{\partial r}\right)^{2} + c^{2} r^{3} \frac{\partial^2\,\nu}{\partial r^2} + c^{2} r^{2} \frac{\partial\,\nu}{\partial r}\right)} \sin\left({\theta}\right)^{2}}{c^{2} r} \right) \mathrm{d} {\phi}\otimes \mathrm{d} {\phi}$

Now for the energy-momentum tensor, $T_{\mu\nu}$:

u = M.vector_field('u')
u[0] = exp(-nu)
u.display()

$u = e^{\left(-\nu\left(r\right)\right)} \frac{\partial}{\partial t }$

We can check if it is indeed the timelike 4-vector by checking $u_\mu u^\mu = -c^2$ by contracting it with the metric g using a method contract():

umuumu = g.contract(0,u,0).contract(0,u,0).function_chart(X)
umuumu == -c^2

$\mathrm{True}$

The product $u_\mu u^\mu$ can be also obtained in much a simpler way, by just invoking

umuumu = g(u,u)
umuumu == -c^2

$\mathrm{True}$

Let's now adopt $T_{\mu\nu}$ in perfect fluid form:

u_form = u.down(g)
T = (rho + p/c^2)*(u_form*u_form) + p*g
T.set_name('T')
print T
T.display()

field of symmetric bilinear forms 'T' on the 4-dimensional manifold 'M'

$T = c^{4} e^{\left(2 \, \nu\left(r\right)\right)} \rho\left(r\right) \mathrm{d} t\otimes \mathrm{d} t + \left( \frac{c^{2} r P\left(r\right)}{c^{2} r - 2 \, G m\left(r\right)} \right) \mathrm{d} r\otimes \mathrm{d} r + r^{2} P\left(r\right) \mathrm{d} {\theta}\otimes \mathrm{d} {\theta} + r^{2} P\left(r\right) \sin\left({\theta}\right)^{2} \mathrm{d} {\phi}\otimes \mathrm{d} {\phi}$

Ttrace = (T.up(g, 0)).trace(0, 1)
Ttrace.display()

$\begin{array}{llcl} & \mathcal{M} & \longrightarrow & \mathbb{R} \\ & \left(t, r, {\theta}, {\phi}\right) & \longmapsto & -c^{2} \rho\left(r\right) + 3 \, P\left(r\right) \end{array}$

Three components of the Einstein equations are as follows - the $E_{tt}$ one:

E0=(E[0,0] - (8*pi*G/c^4)*T[0,0]).expr() == 0

A small reorganization of the first equation, using the function solve() to solve for $dm/dr$:

E0 = solve((E0*(-r^2/(G*exp(2*nu))/2)).expand().simplify(), m.diff(r))[0]

Using SageManifolds ExpressionNice to display the derivatives in textbook form:

from sage.geometry.manifolds.utilities import ExpressionNice

ExpressionNice(E0)

$\frac{\partial\,m}{\partial r} = 4 \, \pi r^{2} \rho\left(r\right)$

Radial component of Einstein's equations, $E_{rr}$:

E1 = (E[1,1] - (8*pi*G/c^4)*T[1,1]).expr() == 0

E1 = solve((E1*(c^4*r^3 - 2*G*c^2*r^2*m)/2).expand().simplify_full(), nu.diff(r))[0]
ExpressionNice(E1)

$\frac{\partial\,\nu}{\partial r} = \frac{{\left(4 \, \pi r^{3} P\left(r\right) + c^{2} m\left(r\right)\right)} G}{c^{4} r^{2} - 2 \, G c^{2} r m\left(r\right)}$

For the third equation we use radial part of the energy-momentum conservation equation $\nabla_\mu T^{\mu\nu}$. First, to define the energy-momentum tensor $T^{\mu\nu}$ itself:

Tup = T.up(g,0).up(g,1)
Tup[:]

$\left(\begin{array}{rrrr} e^{\left(-2 \, \nu\left(r\right)\right)} \rho\left(r\right) & 0 & 0 & 0 \\ 0 & \frac{c^{2} r P\left(r\right) - 2 \, G P\left(r\right) m\left(r\right)}{c^{2} r} & 0 & 0 \\ 0 & 0 & \frac{P\left(r\right)}{r^{2}} & 0 \\ 0 & 0 & 0 & \frac{P\left(r\right)}{r^{2} \sin\left({\theta}\right)^{2}} \end{array}\right)$

Connection ${\tt nab}$ for the covariant derivative, and the printout of the non-vanishing Christoffel symbols:

nab = g.connection()
nab.display()

$\begin{array}{lcl} \Gamma_{ \phantom{\, t } \, t \, r }^{ \, t \phantom{\, t } \phantom{\, r } } & = & \frac{\partial\,\nu}{\partial r} \\ \Gamma_{ \phantom{\, t } \, r \, t }^{ \, t \phantom{\, r } \phantom{\, t } } & = & \frac{\partial\,\nu}{\partial r} \\ \Gamma_{ \phantom{\, r } \, t \, t }^{ \, r \phantom{\, t } \phantom{\, t } } & = & \frac{c^{2} r e^{\left(2 \, \nu\left(r\right)\right)} \frac{\partial\,\nu}{\partial r} - 2 \, G e^{\left(2 \, \nu\left(r\right)\right)} m\left(r\right) \frac{\partial\,\nu}{\partial r}}{r} \\ \Gamma_{ \phantom{\, r } \, r \, r }^{ \, r \phantom{\, r } \phantom{\, r } } & = & \frac{{\left(r \frac{\partial\,m}{\partial r} - m\left(r\right)\right)} G}{c^{2} r^{2} - 2 \, G r m\left(r\right)} \\ \Gamma_{ \phantom{\, r } \, {\theta} \, {\theta} }^{ \, r \phantom{\, {\theta} } \phantom{\, {\theta} } } & = & -\frac{c^{2} r - 2 \, G m\left(r\right)}{c^{2}} \\ \Gamma_{ \phantom{\, r } \, {\phi} \, {\phi} }^{ \, r \phantom{\, {\phi} } \phantom{\, {\phi} } } & = & -\frac{c^{2} r \sin\left({\theta}\right)^{2} - 2 \, G m\left(r\right) \sin\left({\theta}\right)^{2}}{c^{2}} \\ \Gamma_{ \phantom{\, {\theta} } \, r \, {\theta} }^{ \, {\theta} \phantom{\, r } \phantom{\, {\theta} } } & = & \frac{1}{r} \\ \Gamma_{ \phantom{\, {\theta} } \, {\theta} \, r }^{ \, {\theta} \phantom{\, {\theta} } \phantom{\, r } } & = & \frac{1}{r} \\ \Gamma_{ \phantom{\, {\theta} } \, {\phi} \, {\phi} }^{ \, {\theta} \phantom{\, {\phi} } \phantom{\, {\phi} } } & = & -\cos\left({\theta}\right) \sin\left({\theta}\right) \\ \Gamma_{ \phantom{\, {\phi} } \, r \, {\phi} }^{ \, {\phi} \phantom{\, r } \phantom{\, {\phi} } } & = & \frac{1}{r} \\ \Gamma_{ \phantom{\, {\phi} } \, {\theta} \, {\phi} }^{ \, {\phi} \phantom{\, {\theta} } \phantom{\, {\phi} } } & = & \frac{\cos\left({\theta}\right)}{\sin\left({\theta}\right)} \\ \Gamma_{ \phantom{\, {\phi} } \, {\phi} \, r }^{ \, {\phi} \phantom{\, {\phi} } \phantom{\, r } } & = & \frac{1}{r} \\ \Gamma_{ \phantom{\, {\phi} } \, {\phi} \, {\theta} }^{ \, {\phi} \phantom{\, {\phi} } \phantom{\, {\theta} } } & = & \frac{\cos\left({\theta}\right)}{\sin\left({\theta}\right)} \end{array}$

co = nab(Tup)

The following calculates the radial component of $\nabla_\mu T^{\mu\nu}$:

cosum = 0
# radial component of the covariant derivative:
for i in M.irange():
    cosum += co[i,1,i]
cosum

$\frac{c^{2} r \frac{\partial\,P}{\partial r} - 2 \, {\left(m\left(r\right) \frac{\partial\,P}{\partial r} + {\left(c^{2} m\left(r\right) \rho\left(r\right) + P\left(r\right) m\left(r\right)\right)} \frac{\partial\,\nu}{\partial r}\right)} G + {\left(c^{4} r \rho\left(r\right) + c^{2} r P\left(r\right)\right)} \frac{\partial\,\nu}{\partial r}}{c^{2} r}$

Solve for $dp/dr$:

E2 = solve(cosum.expr(), p.diff(r))[0]
ExpressionNice(E2)

$\frac{\partial\,P}{\partial r} = -{\left(c^{2} \rho\left(r\right) + P\left(r\right)\right)} \frac{\partial\,\nu}{\partial r}$

Finally, the three TOV equations:

ExpressionNice(E0), ExpressionNice(E1), ExpressionNice(E2)

($\frac{\partial\,m}{\partial r} = 4 \, \pi r^{2} \rho\left(r\right)$, $\frac{\partial\,\nu}{\partial r} = \frac{{\left(4 \, \pi r^{3} P\left(r\right) + c^{2} m\left(r\right)\right)} G}{c^{4} r^{2} - 2 \, G c^{2} r m\left(r\right)}$, $\frac{\partial\,P}{\partial r} = -{\left(c^{2} \rho\left(r\right) + P\left(r\right)\right)} \frac{\partial\,\nu}{\partial r}$)