Tolman-Oppenheimer-Volkoff equations

This Jupyter/SageMath worksheet is relative to the lectures General relativity computations with SageManifolds given at the NewCompStar School 2016 (Coimbra, Portugal).

These computations are based on SageManifolds (v0.9)

Click here to download the worksheet file (ipynb format). To run it, you must start SageMath with the Jupyter notebook, with the command sage -n jupyter

This worksheet is divided in two parts:

  1. Deriving the TOV system from the Einstein equation
  2. Solving the TOV system to get stellar models

First we set up the notebook to display mathematical objects using LaTeX formatting:

In [1]:
%display latex

1. Deriving the TOV system from the Einstein equation

Spacetime

We declare the spacetime manifold $M$:

In [2]:
M = Manifold(4, 'M')
print(M)
4-dimensional differentiable manifold M
In [3]:
M?
In [4]:
M??

We declare the chart of spherical coordinates $(t,r,\theta,\phi)$:

In [5]:
M.chart?
In [6]:
X.<t,r,th,ph> = M.chart(r't r:(0,+oo) th:(0,pi):\theta ph:(0,2*pi):\phi')
X
Out[6]:

Metric tensor

The static and spherically symmetric metric ansatz, with the unknown functions $\nu(r)$ and $m(r)$:

In [7]:
g = M.lorentzian_metric('g')
nu = function('nu')
m = function('m')
g[0,0] = -exp(2*nu(r))
g[1,1] = 1/(1-2*m(r)/r)
g[2,2] = r^2
g[3,3] = (r*sin(th))^2
g.display()
Out[7]:

One can display the metric components as a list:

In [8]:
g.display_comp()
Out[8]:

By default, only the nonzero components are shown; to get all the components, set the option only_nonzero to False:

In [9]:
g.display_comp(only_nonzero=False)
Out[9]:

We can also display the metric components as a matrix, via the [] operator:

In [10]:
g[:]
Out[10]:

The [] operator can also be used to access to individual elements:

In [11]:
g[0,0]
Out[11]:

Einstein equation

Let us start by evaluating the Ricci tensor of $g$:

In [12]:
Ric = g.ricci()
print(Ric)
Field of symmetric bilinear forms Ric(g) on the 4-dimensional differentiable manifold M
In [13]:
Ric.display()
Out[13]:

The Ricci scalar is naturally obtained by the method ricci_scalar():

In [14]:
g.ricci_scalar?
In [15]:
g.ricci_scalar().display()
Out[15]:

As a check we can also compute it by taking the trace of the Ricci tensor with respect to $g$:

In [16]:
g.ricci_scalar() == g.inverse()['^{ab}']*Ric['_{ab}']
Out[16]:

The Einstein tensor $$ G_{ab} := R_{ab} - \frac{1}{2} R \, g_{ab},$$ or in index-free notation: $$ G := \mathrm{Ric}(g) - \frac{1}{2} r(g) \, g $$

In [17]:
G = Ric - 1/2*g.ricci_scalar() * g
G.set_name('G')
print(G)
Field of symmetric bilinear forms G on the 4-dimensional differentiable manifold M
In [18]:
G.display()
Out[18]:

The energy-momentum tensor

We consider a perfect fluid matter model. Let us first defined the fluid 4-velocity $u$:

In [19]:
u = M.vector_field('u')
u[0] = exp(-nu(r))
u.display()
Out[19]:
In [20]:
u[:]
Out[20]:
In [21]:
print(u.parent())
Free module X(M) of vector fields on the 4-dimensional differentiable manifold M

Let us check that $u$ is a normalized timelike vector, i.e. that $g_{ab} u^a u^b = -1$, or, in index-free notation, $g(u,u)=-1$:

In [22]:
g(u,u)
Out[22]:
In [23]:
print(g(u,u))
Scalar field g(u,u) on the 4-dimensional differentiable manifold M
In [24]:
g(u,u).display()
Out[24]:
In [25]:
g(u,u).parent()
Out[25]:
In [26]:
print(g(u,u).parent())
Algebra of differentiable scalar fields on the 4-dimensional differentiable manifold M
In [27]:
g(u,u) == -1
Out[27]:

To form the energy-momentum tensor, we need the 1-form $\underline{u}$ that is metric-dual to the vector $u$, i.e. $u_a = g_{ab} u^b$:

In [28]:
u_form = u.down(g)
print(u_form)
1-form on the 4-dimensional differentiable manifold M
In [29]:
u_form.display()
Out[29]:

The energy-momentum tensor is then $$ T_{ab} = (\rho + p) u_a u_b + p \, g_{ab},$$ or in index-free notation: $$ T = (\rho + p) \underline{u}\otimes\underline{u} + p \, g$$

Since the tensor product $\otimes$ is taken with the * operator, we write:

In [30]:
rho = function('rho')
p = function('p')
T = (rho(r)+p(r))* (u_form * u_form) + p(r) * g
T.set_name('T')
print(T)
Field of symmetric bilinear forms T on the 4-dimensional differentiable manifold M
In [31]:
T.display()
Out[31]:
In [32]:
T(u,u)
Out[32]:
In [33]:
print(T(u,u))
Scalar field T(u,u) on the 4-dimensional differentiable manifold M
In [34]:
T(u,u).display()
Out[34]:

The Einstein equation

The Einstein equation is $$ G = 8\pi T$$ We rewrite it as $E = 0$ with $E := G - 8\pi T$:

In [35]:
E = G - 8*pi*T
E.set_name('E')
print(E)
Field of symmetric bilinear forms E on the 4-dimensional differentiable manifold M
In [36]:
E.display()
Out[36]:
In [37]:
E.display_comp()
Out[37]:
In [38]:
E[0,0]
Out[38]:
In [39]:
EE0 = solve(E[0,0].expr()==0, diff(m(r),r))
EE0
Out[39]:

The notation D[0](m) for the first derivative of $m(r)$ is not very handy; this will be improved in SageMath 7.4; meanwhile, we may use the ExpressionNice utility to restore some textbook notation:

In [40]:
from sage.manifolds.utilities import ExpressionNice
EE0 = ExpressionNice(EE0[0])
EE0
Out[40]:
In [41]:
E[1,1]
Out[41]:
In [42]:
EE1 = solve(E[1,1].expr()==0, diff(nu(r),r))
EE1
Out[42]:
In [43]:
EE1 = ExpressionNice(EE1[0])
EE1
Out[43]:
In [44]:
E[3,3] == E[2,2]*sin(th)^2
Out[44]:
In [45]:
E[2,2]
Out[45]:

The energy-momentum conservation equation

The energy-momentum tensor must obey $$ \nabla_b T^b_{\ \, a} = 0$$ We first form the tensor $T^b_{\ \, a}$ by raising the first index of $T_{ab}$:

In [46]:
Tu = T.up(g, 0)
print(Tu)
Tensor field of type (1,1) on the 4-dimensional differentiable manifold M

We get the Levi-Civita connection $\nabla$ associated with the metric $g$:

In [47]:
nabla = g.connection()
print(nabla)
Levi-Civita connection nabla_g associated with the Lorentzian metric g on the 4-dimensional differentiable manifold M
In [48]:
nabla.display()
Out[48]:

We apply nabla to Tu to get the tensor $(\nabla T)^b_{\ \, ac} = \nabla_c T^b_{\ \, a}$ (MTW index convention):

In [49]:
dTu = nabla(Tu)
print(dTu)
Tensor field of type (1,2) on the 4-dimensional differentiable manifold M

The divergence $\nabla_b T^b_{\ \, a}$ is then computed as the trace of the tensor $(\nabla T)^b_{\ \, ac}$ on the first index ($b$, position 0) and last index ($c$, position 2):

In [50]:
divT = dTu.trace(0,2)
print(divT)
1-form on the 4-dimensional differentiable manifold M

We can also take the trace by using the index notation:

In [51]:
divT == dTu['^b_{ab}']
Out[51]:
In [52]:
print(divT)
1-form on the 4-dimensional differentiable manifold M
In [53]:
divT.display()
Out[53]:
In [54]:
divT[:]
Out[54]:

The only non trivially vanishing components is thus

In [55]:
divT[1]
Out[55]:

Hence the energy-momentum conservation equation $\nabla_b T^b_{\ \, a}=0$ reduces to

In [56]:
EE2 = solve(divT[1].expr()==0, diff(p(r),r))
EE2
Out[56]:
In [57]:
EE2 = ExpressionNice(EE2[0])
EE2
Out[57]:

The TOV system

Let us collect all the independent equations obtained so far:

In [58]:
EE0, EE1, EE2
Out[58]:

2. Solving the TOV system

In order to solve the TOV system, we need to specify a fluid equation of state. For simplicity, we select a polytropic one: $$ p = k \rho^2$$

In [59]:
var('k', domain='real')
p_eos(r) = k*rho(r)^2
p_eos(r)
Out[59]:

We substitute this expression for $p$ in the TOV equations:

In [60]:
EE1_rho = EE1.substitute_function(p, p_eos)
EE1_rho
Out[60]:
In [61]:
EE2_rho = EE2.substitute_function(p, p_eos)
EE2_rho
Out[61]:
In [62]:
EE2_rho = (EE2_rho / (2*k*rho(r))).simplify_full()
EE2_rho
Out[62]:

We subsitute the expression of equation EE1_rho for $\partial \nu/\partial r$, in order to get rid of any derivative in the right-hand sides:

In [63]:
EE2_rho = EE2_rho.subs({diff(nu(r),r): EE1_rho.rhs()}).simplify_full()
EE2_rho
Out[63]:

The system to solve for $(m(r), \nu(r), \rho(r))$ is thus

In [64]:
EE1_rho = ExpressionNice(EE1_rho)
EE2_rho = ExpressionNice(EE2_rho)
EE0, EE1_rho, EE2_rho
Out[64]:

Numerical resolution

Let us use a standard 4th-order Runge-Kutta method:

In [65]:
desolve_system_rk4?

We gather all equations in a list for the ease of manipulation:

In [66]:
eqs = [EE0, EE1_rho, EE2_rho]

To get a numerical solution, we have of course to specify some numerical value for the EOS constant $k$; let us choose $k=1/4$:

In [67]:
k0 = 1/4
rhs = [eq.rhs().subs(k=k0) for eq in eqs]
rhs
Out[67]:

The integration for $m(r)$ and $\rho(r)$ has to stop as soon as $\rho(r)$ become negative. An easy way to ensure this is to use the Heaviside function:

In [68]:
h(x) = (1+sign(x))/2
plot(h(x), (x,-4,4), thickness=2, aspect_ratio=1)
Out[68]:

and to multiply the r.h.s. of the $dm/dr$ and $d\rho/dr$ equations by $h(\rho)$:

In [69]:
rhs[0] = rhs[0] * h(rho(r))
rhs[2] = rhs[2] * h(rho(r))
rhs
Out[69]:

Let us add an extra equation, for the purpose of getting the star's radius as an output of the integration, via the equation $dR/dr = 1$, again multiplied by the heaviside function of $\rho$:

In [70]:
rhs.append(1 * h(rho(r)))
rhs
Out[70]:

For the purpose of the numerical integration via desolve_system_rk4, we have to replace the symbolic functions $m(r)$, $\nu(r)$ and $\rho(r)$ by some symbolic variables, $m_1$, $\nu_1$ and $\rho_1$, say:

In [71]:
var('m_1 nu_1 rho_1 r_1')
rhs = [y.subs({m(r): m_1, nu(r): nu_1, rho(r): rho_1}) for y in rhs]
rhs
Out[71]:

The integration parameters:

In [72]:
rho_c = 1
r_min = 1e-8
r_max = 1
np = 200
delta_r = (r_max - r_min) / (np-1)

The numerical resolution, with the initial conditions $$(r_{\rm min}, m(r_{\rm min}), \nu(r_{\rm min}), \rho(r_{\rm min}), R(r_{\rm min})) = (r_{\rm min},0,0,\rho_{\rm c},r_{\rm min}) $$ set in the parameter ics:

In [73]:
sol = desolve_system_rk4(rhs, vars=(m_1, nu_1, rho_1, r_1), ivar=r, 
                         ics=[r_min, 0, 0, rho_c, r_min], 
                         end_points=(r_min, r_max), step=delta_r)

The solution is returned as a list, the first 10 elements of which being:

In [74]:
sol[:10]
Out[74]:

Each element is of the form $[r, m(r), \nu(r), \rho(r), R(r)]$. So to get the list of $(r, \rho(r))$ values, we write

In [75]:
rho_sol = [(s[0], s[3]) for s in sol]
rho_sol[:10]
Out[75]:

We may then use this list to have some plot of $\rho(r)$, thanks to the function line, which transforms a list into a graphical object:

In [76]:
graph = line(rho_sol, axes_labels=[r'$r$', r'$\rho$'], gridlines=True)
graph
Out[76]:

The solution for $m(r)$:

In [77]:
m_sol = [(s[0], s[1]) for s in sol]
m_sol[:10]
Out[77]:
In [78]:
graph = line(m_sol, axes_labels=[r'$r$', r'$m$'], gridlines=True)
graph
Out[78]:

The solution for $\nu(r)$ (has to be rescaled by adding a constant to ensure $\nu(+\infty) = 1$):

In [79]:
nu_sol = [(s[0], s[2]) for s in sol]
nu_sol[:10]
Out[79]:
In [80]:
graph = line(nu_sol, axes_labels=[r'$r$', r'$\tilde{\nu}$'], gridlines=True)
graph
Out[80]:

The solution for $R(r)$:

In [81]:
r_sol = [(s[0], s[4]) for s in sol]
line(r_sol)
Out[81]:

The total gravitational mass of the star is obtained via the last element (index: -1) of the list of $(r,m(r))$ values:

In [82]:
M_grav = m_sol[-1][1]
M_grav
Out[82]:

Similarly, the stellar radius is obtained through the last element of the list of $(r,R(r))$ values:

In [83]:
R = r_sol[-1][1]
R
Out[83]:

The star's compactness:

In [84]:
M_grav/R
Out[84]:

Sequence of stellar models

Let us perform a loop on the central density. First we set up a list of values for $\rho_c$:

In [85]:
rho_c_min = 0.01
rho_c_max = 3
n_conf = 40
rho_c_list = [rho_c_min + i * (rho_c_max-rho_c_min)/(n_conf-1) for i in range(n_conf)]
rho_c_list
Out[85]:

The loop:

In [86]:
M_list = list()
R_list = list()
for rho_c in rho_c_list:
    sol = desolve_system_rk4(rhs, vars=(m_1, nu_1, rho_1, r_1), ivar=r, 
                             ics=[r_min, 0, 0, rho_c, r_min], 
                             end_points=(r_min, r_max), step=delta_r)
    M_list.append( sol[-1][1] )
    R_list.append( sol[-1][4] )

The mass along the sequence:

In [87]:
M_list
Out[87]:

The radius along the sequence:

In [88]:
R_list
Out[88]:

To draw $M$ as a function of $\rho_{\rm c}$, we use the Python function zip to construct a list of $(\rho_{\rm c}, M)$ values:

In [89]:
zip(rho_c_list, M_list)
Out[89]:
In [90]:
graph = line(zip(rho_c_list, M_list), axes_labels=[r'$\rho_c$', r'$M$'], gridlines=True)
graph
Out[90]:

Similarly, we draw $M$ as a function of $R$:

In [91]:
graph = line(zip(R_list, M_list), axes_labels=[r'$R$', r'$M$'], gridlines=True)
graph
Out[91]:

and we save the plot in a pdf file to use it in our next publication ;-)

In [92]:
graph.save('plot_M_R.pdf')