Friedmann equations

This notebook demonstrates a few capabilities of SageMath in computations regarding cosmological spacetimes with Friedmann-Lemaître-Robertson-Walker (FLRW) metrics. The corresponding tools have been developed in the framework of the SageManifolds project.

NB: a version of SageMath at least equal to 8.2 is required to run this notebook:

In [1]:
version()
Out[1]:
'SageMath version 9.6, Release Date: 2022-05-15'

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

In [2]:
%display latex

We declare the spacetime $M$ as a 4-dimensional Lorentzian manifold:

In [3]:
M = Manifold(4, 'M', structure='Lorentzian')
print(M)
4-dimensional Lorentzian manifold M

We introduce the standard FLRW coordinates, via the method chart(), the argument of which is a string expressing the coordinates names, their ranges (the default is $(-\infty,+\infty)$) and their LaTeX symbols:

In [4]:
fr.<t,r,th,ph> = M.chart(r't r:[0,+oo) th:[0,pi]:\theta ph:[0,2*pi):\phi')
fr
Out[4]:
\(\displaystyle \left(M,(t, r, {\theta}, {\phi})\right)\)

Assuming units where that the speed of light is 1, let us define a few variables: Newton's constant $G$, the cosmological constant $\Lambda$, the spatial curvature constant $k$, the scale factor $a(t)$, the fluid proper density $\rho(t)$ and the fluid pressure $p(t)$:

In [5]:
G = var('G', domain='real')
Lambda = var('Lambda', domain='real')
k = var('k', domain='real')

a = M.scalar_field(function('a')(t), name='a')
rho = M.scalar_field(function('rho')(t), name='rho')
p = M.scalar_field(function('p')(t), name='p')

The FLRW metric is defined by its components in the manifold's default frame, i.e. the frame associated with the FLRW coordinates:

In [6]:
g = M.metric()
g[0,0] = -1
g[1,1] = a*a/(1 - k*r^2)
g[2,2] = a*a*r^2
g[3,3] = a*a*(r*sin(th))^2
g.display()
Out[6]:
\(\displaystyle g = -\mathrm{d} t\otimes \mathrm{d} t + \left( -\frac{a\left(t\right)^{2}}{k r^{2} - 1} \right) \mathrm{d} r\otimes \mathrm{d} r + r^{2} a\left(t\right)^{2} \mathrm{d} {\theta}\otimes \mathrm{d} {\theta} + r^{2} a\left(t\right)^{2} \sin\left({\theta}\right)^{2} \mathrm{d} {\phi}\otimes \mathrm{d} {\phi}\)

A matrix view of the metric components:

In [7]:
g[:]
Out[7]:
\(\displaystyle \left(\begin{array}{rrrr} -1 & 0 & 0 & 0 \\ 0 & -\frac{a\left(t\right)^{2}}{k r^{2} - 1} & 0 & 0 \\ 0 & 0 & r^{2} a\left(t\right)^{2} & 0 \\ 0 & 0 & 0 & r^{2} a\left(t\right)^{2} \sin\left({\theta}\right)^{2} \end{array}\right)\)

The Levi-Civita connection associated with the metric is computed:

In [8]:
nabla = g.connection()
g.christoffel_symbols_display()
Out[8]:
\(\displaystyle \begin{array}{lcl} \Gamma_{ \phantom{\, t} \, r \, r }^{ \, t \phantom{\, r} \phantom{\, r} } & = & -\frac{a\left(t\right) \frac{\partial\,a}{\partial t}}{k r^{2} - 1} \\ \Gamma_{ \phantom{\, t} \, {\theta} \, {\theta} }^{ \, t \phantom{\, {\theta}} \phantom{\, {\theta}} } & = & r^{2} a\left(t\right) \frac{\partial\,a}{\partial t} \\ \Gamma_{ \phantom{\, t} \, {\phi} \, {\phi} }^{ \, t \phantom{\, {\phi}} \phantom{\, {\phi}} } & = & r^{2} a\left(t\right) \sin\left({\theta}\right)^{2} \frac{\partial\,a}{\partial t} \\ \Gamma_{ \phantom{\, r} \, t \, r }^{ \, r \phantom{\, t} \phantom{\, r} } & = & \frac{\frac{\partial\,a}{\partial t}}{a\left(t\right)} \\ \Gamma_{ \phantom{\, r} \, r \, r }^{ \, r \phantom{\, r} \phantom{\, r} } & = & -\frac{k r}{k r^{2} - 1} \\ \Gamma_{ \phantom{\, r} \, {\theta} \, {\theta} }^{ \, r \phantom{\, {\theta}} \phantom{\, {\theta}} } & = & k r^{3} - r \\ \Gamma_{ \phantom{\, r} \, {\phi} \, {\phi} }^{ \, r \phantom{\, {\phi}} \phantom{\, {\phi}} } & = & {\left(k r^{3} - r\right)} \sin\left({\theta}\right)^{2} \\ \Gamma_{ \phantom{\, {\theta}} \, t \, {\theta} }^{ \, {\theta} \phantom{\, t} \phantom{\, {\theta}} } & = & \frac{\frac{\partial\,a}{\partial t}}{a\left(t\right)} \\ \Gamma_{ \phantom{\, {\theta}} \, r \, {\theta} }^{ \, {\theta} \phantom{\, r} \phantom{\, {\theta}} } & = & \frac{1}{r} \\ \Gamma_{ \phantom{\, {\theta}} \, {\phi} \, {\phi} }^{ \, {\theta} \phantom{\, {\phi}} \phantom{\, {\phi}} } & = & -\cos\left({\theta}\right) \sin\left({\theta}\right) \\ \Gamma_{ \phantom{\, {\phi}} \, t \, {\phi} }^{ \, {\phi} \phantom{\, t} \phantom{\, {\phi}} } & = & \frac{\frac{\partial\,a}{\partial t}}{a\left(t\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)} \end{array}\)

Ricci tensor:

In [9]:
Ricci = nabla.ricci()
Ricci.display()
Out[9]:
\(\displaystyle \mathrm{Ric}\left(g\right) = -\frac{3 \, \frac{\partial^2\,a}{\partial t ^ 2}}{a\left(t\right)} \mathrm{d} t\otimes \mathrm{d} t + \left( -\frac{2 \, \left(\frac{\partial\,a}{\partial t}\right)^{2} + a\left(t\right) \frac{\partial^2\,a}{\partial t ^ 2} + 2 \, k}{k r^{2} - 1} \right) \mathrm{d} r\otimes \mathrm{d} r + \left( 2 \, r^{2} \left(\frac{\partial\,a}{\partial t}\right)^{2} + r^{2} a\left(t\right) \frac{\partial^2\,a}{\partial t ^ 2} + 2 \, k r^{2} \right) \mathrm{d} {\theta}\otimes \mathrm{d} {\theta} + {\left(2 \, r^{2} \left(\frac{\partial\,a}{\partial t}\right)^{2} + r^{2} a\left(t\right) \frac{\partial^2\,a}{\partial t ^ 2} + 2 \, k r^{2}\right)} \sin\left({\theta}\right)^{2} \mathrm{d} {\phi}\otimes \mathrm{d} {\phi}\)
In [10]:
Ricci.display_comp()
Out[10]:
\(\displaystyle \begin{array}{lcl} \mathrm{Ric}\left(g\right)_{ \, t \, t }^{ \phantom{\, t}\phantom{\, t} } & = & -\frac{3 \, \frac{\partial^2\,a}{\partial t ^ 2}}{a\left(t\right)} \\ \mathrm{Ric}\left(g\right)_{ \, r \, r }^{ \phantom{\, r}\phantom{\, r} } & = & -\frac{2 \, \left(\frac{\partial\,a}{\partial t}\right)^{2} + a\left(t\right) \frac{\partial^2\,a}{\partial t ^ 2} + 2 \, k}{k r^{2} - 1} \\ \mathrm{Ric}\left(g\right)_{ \, {\theta} \, {\theta} }^{ \phantom{\, {\theta}}\phantom{\, {\theta}} } & = & 2 \, r^{2} \left(\frac{\partial\,a}{\partial t}\right)^{2} + r^{2} a\left(t\right) \frac{\partial^2\,a}{\partial t ^ 2} + 2 \, k r^{2} \\ \mathrm{Ric}\left(g\right)_{ \, {\phi} \, {\phi} }^{ \phantom{\, {\phi}}\phantom{\, {\phi}} } & = & {\left(2 \, r^{2} \left(\frac{\partial\,a}{\partial t}\right)^{2} + r^{2} a\left(t\right) \frac{\partial^2\,a}{\partial t ^ 2} + 2 \, k r^{2}\right)} \sin\left({\theta}\right)^{2} \end{array}\)

Ricci scalar ($R^\mu_{\ \, \mu}$):

In [11]:
Ricci_scalar = g.ricci_scalar()
Ricci_scalar.display()
Out[11]:
\(\displaystyle \begin{array}{llcl} \mathrm{r}\left(g\right):& M & \longrightarrow & \mathbb{R} \\ & \left(t, r, {\theta}, {\phi}\right) & \longmapsto & \frac{6 \, {\left(\left(\frac{\partial\,a}{\partial t}\right)^{2} + a\left(t\right) \frac{\partial^2\,a}{\partial t ^ 2} + k\right)}}{a\left(t\right)^{2}} \end{array}\)

The fluid 4-velocity:

In [12]:
u = M.vector_field('u')
u[0] = 1
u.display()
Out[12]:
\(\displaystyle u = \frac{\partial}{\partial t }\)
In [13]:
g(u,u).expr()
Out[13]:
\(\displaystyle -1\)

u.dot(u) is equivalent to g(u,u):

In [14]:
u.dot(u).expr()
Out[14]:
\(\displaystyle -1\)

Perfect fluid energy-momentum tensor $T$:

In [15]:
u_form = u.down(g) # the 1-form associated to u by metric duality
T = (rho+p)*(u_form*u_form) + p*g
T.set_name('T')
print(T)
T.display()
Field of symmetric bilinear forms T on the 4-dimensional Lorentzian manifold M
Out[15]:
\(\displaystyle T = \rho\left(t\right) \mathrm{d} t\otimes \mathrm{d} t + \left( -\frac{a\left(t\right)^{2} p\left(t\right)}{k r^{2} - 1} \right) \mathrm{d} r\otimes \mathrm{d} r + r^{2} a\left(t\right)^{2} p\left(t\right) \mathrm{d} {\theta}\otimes \mathrm{d} {\theta} + r^{2} a\left(t\right)^{2} p\left(t\right) \sin\left({\theta}\right)^{2} \mathrm{d} {\phi}\otimes \mathrm{d} {\phi}\)

The trace of $T$ (we use index notation to denote the double contraction $g^{ab} T_{ab}$):

In [16]:
Ttrace = g.inverse()['^ab']*T['_ab']
Ttrace.display()
Out[16]:
\(\displaystyle \begin{array}{llcl} & M & \longrightarrow & \mathbb{R} \\ & \left(t, r, {\theta}, {\phi}\right) & \longmapsto & 3 \, p\left(t\right) - \rho\left(t\right) \end{array}\)

Einstein equation: $R_{\mu \nu} - {1 \over 2} R g_{\mu \nu} + \Lambda g_{\mu \nu} = {8 \pi G} T_{\mu \nu}$

In [17]:
E1 = Ricci - Ricci_scalar/2*g + Lambda*g - (8*pi*G)*T
print("First Friedmann equation:\n")
E1[0,0].expr().expand() == 0
First Friedmann equation:

Out[17]:
\(\displaystyle -8 \, \pi G \rho\left(t\right) - \Lambda + \frac{3 \, \frac{\partial}{\partial t}a\left(t\right)^{2}}{a\left(t\right)^{2}} + \frac{3 \, k}{a\left(t\right)^{2}} = 0\)

Trace-reversed version of the Einstein equation: $R_{\mu \nu} - \Lambda g_{\mu \nu} = {8 \pi G} \left(T_{\mu \nu} - {1 \over 2}T\,g_{\mu \nu}\right)$

In [18]:
E2 = Ricci - Lambda*g - (8*pi*G)*(T - Ttrace/2*g)
print("Second Friedmann equation:\n")
E2[0,0].expr().expand() == 0
Second Friedmann equation:

Out[18]:
\(\displaystyle -12 \, \pi G p\left(t\right) - 4 \, \pi G \rho\left(t\right) + \Lambda - \frac{3 \, \frac{\partial^{2}}{(\partial t)^{2}}a\left(t\right)}{a\left(t\right)} = 0\)