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).
%display text latex
set_nproc()
omit_function_args(True)
We first declare the spacetime M
as a general 4-dimensional manifold,
M = Manifold(4, '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()
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)
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()
The Ricci tensor is a result of a method ricci()
acting on the metric g
:
Ricci = g.ricci(); Ricci.display()
For example, the $R_{tt}$ component is
Ricci[0,0]
Ricci[0,0].expand().collect(nu.diff(r)).collect(nu.diff(r,r)).collect(c^2*exp(2*nu)).collect_common_factors()
Ricci[1,1].expand().collect_common_factors()
Ricci[2,2].expand().collect(nu.diff(r)).collect(nu.diff(r,r)).collect(c^2*exp(2*nu))
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))
It is the trace of the Ricci tensor, $R = R_\mu^\mu$:
Ric_scalar == Ricci.up(g, 1).trace(0, 1)
Left side of the Einstein equations is
E = Ricci - (Ric_scalar*g)/2; E.display()
Now for the energy-momentum tensor, $T_{\mu\nu}$:
u = M.vector_field('u')
u[0] = exp(-nu)
u.display()
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
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
Let's now addopt $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'
Ttrace = (T.up(g, 0)).trace(0, 1)
Ttrace.display()
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)
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)
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[:]
Connection ${\tt nab}$ for the covariant derivative, and the printout of the non-vanishing Christoffel symbols:
nab = g.connection()
nab.display()
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
Solve for $dp/dr$:
E2 = solve(cosum.expr(), p.diff(r))[0]
ExpressionNice(E2)
Finally, the three TOV equations:
ExpressionNice(E0), ExpressionNice(E1), ExpressionNice(E2)