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 [1]:
%display text latex
set_nproc()
omit_function_args(True)

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

In [2]:
M = Manifold(4, 'M')

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

In [3]:
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 [4]:
# 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 [5]:
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[5]:

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 [6]:
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 [7]:
g[0,0] = -c^2*exp(2*nu)
g[1,1] = 1/(1-2*G*m/(r*c^2))
g.display()
Out[7]:

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

In [8]:
Ricci = g.ricci(); Ricci.display()
Out[8]:

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

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

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

In [13]:
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[13]:

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

In [14]:
Ric_scalar == Ricci.up(g, 1).trace(0, 1)
Out[14]:

Left side of the Einstein equations is

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

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

In [16]:
u = M.vector_field('u')
u[0] = exp(-nu)
u.display()
Out[16]:

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 [17]:
umuumu = g.contract(0,u,0).contract(0,u,0).function_chart(X)
umuumu == -c^2
Out[17]:

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

In [18]:
umuumu = g(u,u)
umuumu == -c^2
Out[18]:

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

In [19]:
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[19]:
In [20]:
Ttrace = (T.up(g, 0)).trace(0, 1)
Ttrace.display()
Out[20]:

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

In [21]:
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 [22]:
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:

In [23]:
from sage.geometry.manifolds.utilities import ExpressionNice
In [24]:
ExpressionNice(E0)
Out[24]:

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

In [25]:
E1 = (E[1,1] - (8*pi*G/c^4)*T[1,1]).expr() == 0
In [26]:
E1 = solve((E1*(c^4*r^3 - 2*G*c^2*r^2*m)/2).expand().simplify_full(), nu.diff(r))[0]
ExpressionNice(E1)
Out[26]:

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 [27]:
Tup = T.up(g,0).up(g,1)
Tup[:]
Out[27]:

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

In [28]:
nab = g.connection()
nab.display()
Out[28]:
In [29]:
co = nab(Tup)

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

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

Solve for $dp/dr$:

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

Finally, the three TOV equations:

In [32]:
ExpressionNice(E0), ExpressionNice(E1), ExpressionNice(E2)
Out[32]: