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

In :
%display text latex
set_nproc()
omit_function_args(True)


We first declare the spacetime M as a general 4-dimensional manifold,

In :
M = Manifold(4, 'M')


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

In :
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$:

In :
# 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:

In :
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()

Out:

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:

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


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

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

Out:

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

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

Out:

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

In :
Ricci[0,0]

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

Out:
In :
Ricci[1,1].expand().collect_common_factors()

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

Out:

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

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

Out:

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

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

Out:

Left side of the Einstein equations is

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

Out:

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

In :
u = M.vector_field('u')
u = exp(-nu)
u.display()

Out:

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

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

Out:

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

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

Out:

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

In :
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'

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

Out:

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

In :
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$:

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


Using SageManifolds ExpressionNice to display the derivatives in textbook form:

In :
from sage.geometry.manifolds.utilities import ExpressionNice

In :
ExpressionNice(E0)

Out:

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

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

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

Out:

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:

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

Out:

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

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

Out:
In :
co = nab(Tup)


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

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

Out:

Solve for $dp/dr$:

In :
E2 = solve(cosum.expr(), p.diff(r))
ExpressionNice(E2)

Out:

Finally, the three TOV equations:

In :
ExpressionNice(E0), ExpressionNice(E1), ExpressionNice(E2)

Out: