Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
Download

📚 The CoCalc Library - books, templates and other resources

133958 views
License: OTHER
Tweet
Share
Share
%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μν−12Rgμν=TμνR_{\mu\nu} - \frac{1}{2}Rg_{\mu\nu} = T_{\mu\nu} for a corresponding spherically symmetric, stationary metric gg. In the above, RμνR_{\mu\nu} is the Ricci tensor, R=RμμR=R^\mu_\mu is the Ricci scalar, and Tμν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 rr - metric functions ν(r)\nu(r) and λ(r)\lambda(r), matter pressure p(r)p(r) and energy density ρ(r)\rho(r), as well as the mass m(r)m(r) enclosed within the sphere of the radius rr:

# 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 λ(r)\lambda(r) with 1−2Gmrc2\frac{1-2Gm}{rc^2}, as well as explicitly put c2c^2 in front of gttg_{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 RttR_{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μμ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μν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μuμ=−c2u_\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μuμ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μν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 EttE_{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/drdm/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, ErrE_{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 ∇μTμν\nabla_\mu T^{\mu\nu}. First, to define the energy-momentum tensor Tμν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 nab{\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 ∇μTμν\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/drdp/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}$)