Return to EinsteinPy や ctensor でアインシュタイン方程式を計算する

EinsteinPy でフリードマン方程式を求める

EinsteinPy を使ってフリードマン・ルメートル・ロバートソン・ウォーカー (FLRW) 計量からアインシュタイン方程式$$G^{\mu}_{\ \ \nu} + \Lambda \delta^{\mu}_{\ \ \nu} = 8\pi G T^{\mu}_{\ \ \nu}$$を使ってフリードマン方程式を求める。

ライブラリの import

In [1]:
from sympy.abc import *
from sympy import *

from einsteinpy.symbolic import *

フリードマン・ルメートル・ロバートソン・ウォーカー計量

Wikipedia の記述にそって，FLRW 計量を以下のようにおく。

$$ds^2 = -dt^2 + a^2(t) \left[\frac{dr^2}{1-k r^2} + r^2 \left(d\theta^2 + \sin^2\theta\,d\phi^2 \right)\right] $$

なぜこのような計量になるかについては，以下のページも参照：

a $\equiv a(t)$, g $\equiv g_{\mu\nu}$ と定義。

In [2]:
# スケール因子 a(t) は時間のみの関数
a = Function('a')(t)

# diag() を使うときは
# .tolist() をつける
Metric = diag(-1, 
              a**2/(1-k*r**2), 
              a**2*r**2, 
              a**2*r**2 * sin(theta)**2).tolist()

g = MetricTensor(Metric, [t, r, theta, phi])

# 計量テンソルの全成分の表示は .tensor() をつける
g.tensor()
Out[2]:
$\displaystyle \left[\begin{matrix}-1 & 0 & 0 & 0\\0 & \frac{a^{2}{\left(t \right)}}{-k r^{2} + 1} & 0 & 0\\0 & 0 & r^{2} a^{2}{\left(t \right)} & 0\\0 & 0 & 0 & r^{2} a^{2}{\left(t \right)} \sin^{2}{\left(\theta \right)}\end{matrix}\right]$

アインシュタイン・テンソル ein

ein $\equiv \displaystyle G^{\mu}_{\ \ \nu} = R^{\mu}_{\ \ \nu} -\frac{1}{2} R \delta^{\mu}_{\ \ \nu} $ と定義。（.change_config('ul') で上付下付に）

In [3]:
ein = EinsteinTensor.from_metric(g).change_config('ul')

アインシュタインテンソルの全成分を表示

ein $ = G^{\mu}_{\ \ \nu}$ を4行4列の行列として表示。

In [4]:
# アインシュタインテンソルの全成分表示は .tensor() をつけて...
ein.tensor()
Out[4]:
$\displaystyle \left[\begin{matrix}\frac{3 \left(-k -\left(\frac{d}{d t} a{\left(t \right)}\right)^{2}\right)}{a^{2}{\left(t \right)}} & 0 & 0 & 0\\0 & \frac{-k -2 a{\left(t \right)} \frac{d^{2}}{d t^{2}} a{\left(t \right)} -\left(\frac{d}{d t} a{\left(t \right)}\right)^{2}}{a^{2}{\left(t \right)}} & 0 & 0\\0 & 0 & \frac{-k -2 a{\left(t \right)} \frac{d^{2}}{d t^{2}} a{\left(t \right)} -\left(\frac{d}{d t} a{\left(t \right)}\right)^{2}}{a^{2}{\left(t \right)}} & 0\\0 & 0 & 0 & \frac{-k -2 a{\left(t \right)} \frac{d^{2}}{d t^{2}} a{\left(t \right)} -\left(\frac{d}{d t} a{\left(t \right)}\right)^{2}}{a^{2}{\left(t \right)}}\end{matrix}\right]$

$G^{\mu}_{\ \ \nu}$ のゼロでない成分を個別に表示

アインシュタインテンソル $G^{\mu}_{\ \ \nu}$ のゼロでない成分を個別に表示する。

ein[0, 0] $= G^{0}_{\ \ 0}$
In [5]:
ein[0,0]
Out[5]:
$\displaystyle \frac{3 \left(-k -\left(\frac{d}{d t} a{\left(t \right)}\right)^{2}\right)}{a^{2}{\left(t \right)}}$
ein[1, 1] $= G^{1}_{\ \ 1}$
In [6]:
ein[1,1]
Out[6]:
$\displaystyle \frac{-k -2 a{\left(t \right)} \frac{d^{2}}{d t^{2}} a{\left(t \right)} -\left(\frac{d}{d t} a{\left(t \right)}\right)^{2}}{a^{2}{\left(t \right)}}$
ein[2, 2] $= G^{2}_{\ \ 2}$
In [7]:
ein[2,2]
Out[7]:
$\displaystyle \frac{-k -2 a{\left(t \right)} \frac{d^{2}}{d t^{2}} a{\left(t \right)} -\left(\frac{d}{d t} a{\left(t \right)}\right)^{2}}{a^{2}{\left(t \right)}}$
ein[3, 3] $= G^{3}_{\ \ 3}$
In [8]:
ein[3,3]
Out[8]:
$\displaystyle \frac{-k -2 a{\left(t \right)} \frac{d^{2}}{d t^{2}} a{\left(t \right)} -\left(\frac{d}{d t} a{\left(t \right)}\right)^{2}}{a^{2}{\left(t \right)}}$

念のため，$G^{2}_{\ \ 2}$ と $G^{3}_{\ \ 3}$ が等しいことを確認します。

In [9]:
ein[2, 2] == ein[3, 3]
Out[9]:
True

完全流体のエネルギー運動量テンソル T(mu, nu, P)

\begin{eqnarray}
T^{\mu}_{\ \ \nu} &=& (\rho + P) u^{\mu} u_{\nu} + P \delta^{\mu}_{\ \ \nu} \\
u^{\mu} &=& (1, 0, 0, 0) \\
u_{\nu} &=& g_{\nu\mu} u^{\mu} = (-1, 0, 0, 0)
\end{eqnarray}

クロネッカーのデルタは，$\delta^{\mu}_{\ \ \nu} = $ KroneckerDelta(mu, nu)

  • rho $ \equiv \rho(t) $
  • P $\equiv P(t) $
  • uu $ \equiv u^{\mu}$
  • ud $ \equiv u_{\nu}$
  • T(mu, nu, P) $ \equiv T^{\mu}_{\ \ \nu}$

と定義。

In [10]:
# エネルギー密度 ρ
rho = Function('rho')(t)
# 圧力 P
P = Function('P')(t)

# u^{\mu}
uu = [1, 0, 0, 0]
# u_{\nu}
ud = [-1, 0, 0, 0]

def T(mu, nu, P):
    return (rho + P)*uu[mu]*ud[nu] + P*KroneckerDelta(mu, nu)

アインシュタイン方程式 EinEq(mu, nu, P)

\begin{eqnarray}
G^{\mu}_{\ \ \nu} +\Lambda \delta^{\mu}_{\ \ \nu} &=& 8\pi G T^{\mu}_{\ \ \nu} \\
\therefore\ \ G^{\mu}_{\ \ \nu} &=& 8\pi G T^{\mu}_{\ \ \nu} -\Lambda \delta^{\mu}_{\ \ \nu}
\end{eqnarray}

In [11]:
# 宇宙定数 Λ と sympy.core.function.Lambda がかち合うので
# あらためて変数として宣言
var('Lambda')
Out[11]:
$\displaystyle \Lambda$

EinEq(mu, nu, P): $ G^{\mu}_{\ \ \nu} = 8\pi G T^{\mu}_{\ \ \nu} -\Lambda \delta^{\mu}_{\ \ \nu}$ と定義。

In [12]:
def EinEq(mu, nu, P):
    return Eq(expand(ein[mu, nu]), 
              8*pi*G*T(mu, nu, P) - Lambda*KroneckerDelta(mu, nu))

ダスト物質の場合

$P = 0$

EinEq(0, 0, 0): $G^{0}_{\ \ 0} = 8\pi G T^{0}_{\ \ 0} -\Lambda \delta^{0}_{\ \ 0} $
In [13]:
EinEq(0,0,0)
Out[13]:
$\displaystyle -\frac{3 k}{a^{2}{\left(t \right)}} -\frac{3 \left(\frac{d}{d t} a{\left(t \right)}\right)^{2}}{a^{2}{\left(t \right)}} = -8 \pi G \rho{\left(t \right)} -\Lambda$

Eq(LHS, RHS) で表される方程式の両辺に定数 $-\tfrac{1}{3}$ をかけるには，

Eq(LHS, RHS)*(-1/3) ではなくて，

左辺 LHS，右辺 RHS それぞれに $-\tfrac{1}{3}$ をかけてからあらためて

Eq(LHS*(-1/3), RHS*(-1/3)) とする必要があります。

方程式 Eq() の左辺は Eq().lhs，右辺は Eq().rhs として参照できます。以下ではに $-\tfrac{1}{3}$ をかけます。

EinEq(0,0,0) * (-Rational(1,3)) すなわち

$\displaystyle -\frac{1}{3} G^{0}_{\ \ 0} = -\frac{8 \pi G}{3} + \frac{\Lambda}{3}$ を表示。

In [14]:
LHS = EinEq(0,0,0).lhs
RHS = EinEq(0,0,0).rhs

Eq(LHS * (-Rational(1,3)), 
   RHS * (-Rational(1,3)))
Out[14]:
$\displaystyle \frac{k}{a^{2}{\left(t \right)}} + \frac{\left(\frac{d}{d t} a{\left(t \right)}\right)^{2}}{a^{2}{\left(t \right)}} = \frac{8 \pi G \rho{\left(t \right)}}{3} + \frac{\Lambda}{3}$
フリードマン方程式

Friedmann $=$ EinEq(0,0,0) * (-Rational(1,3)) とおく。

In [15]:
# 上のセルの出力結果を変数 Friemann に代入し，
# Friedmann 方程式とする
Friedmann = _
EinEq(1, 1, 0): $G^{1}_{\ \ 1} = 8\pi G T^{1}_{\ \ 1} -\Lambda \delta^{1}_{\ \ 1} $
In [16]:
EinEq(1,1,0)
Out[16]:
$\displaystyle -\frac{k}{a^{2}{\left(t \right)}} -\frac{2 \frac{d^{2}}{d t^{2}} a{\left(t \right)}}{a{\left(t \right)}} -\frac{\left(\frac{d}{d t} a{\left(t \right)}\right)^{2}}{a^{2}{\left(t \right)}} = -\Lambda$

Friedmann すなわち EinEq(0,0,0) * (-Rational(1,3)) と EinEq(1,1,0) の両辺を足して $-2$ でわり，$k$ と $\displaystyle \left(\frac{\dot{a}}{a}\right)^2$ の項を消去し，$\displaystyle \frac{\ddot{a}}{a}$ の式にする。

In [17]:
# 上の式と Friedmann 方程式を足して -2 で割る
Eq((Friedmann.lhs + EinEq(1,1,0).lhs)*(-Rational(1,2)), 
   (Friedmann.rhs + EinEq(1,1,0).rhs)*(-Rational(1,2)))
Out[17]:
$\displaystyle \frac{\frac{d^{2}}{d t^{2}} a{\left(t \right)}}{a{\left(t \right)}} = -\frac{4 \pi G \rho{\left(t \right)}}{3} + \frac{\Lambda}{3}$
$\ddot{a}$ の式

上で得られた $\displaystyle \frac{\ddot{a}}{a}$ の式を
att とおく。

In [18]:
# 上のセルの出力結果を変数 att に代入し， 
# a の2階微分の方程式とする
att = _
エネルギー密度の式

Friedmann 方程式 Friedmann は $\displaystyle \left(\frac{\dot{a}}{a} \right)^2$ を含む式であり…

In [19]:
Friedmann
Out[19]:
$\displaystyle \frac{k}{a^{2}{\left(t \right)}} + \frac{\left(\frac{d}{d t} a{\left(t \right)}\right)^{2}}{a^{2}{\left(t \right)}} = \frac{8 \pi G \rho{\left(t \right)}}{3} + \frac{\Lambda}{3}$

att は $\displaystyle \frac{\ddot{a}}{a}$ についての方程式であり…

In [20]:
att
Out[20]:
$\displaystyle \frac{\frac{d^{2}}{d t^{2}} a{\left(t \right)}}{a{\left(t \right)}} = -\frac{4 \pi G \rho{\left(t \right)}}{3} + \frac{\Lambda}{3}$

Friedmann 方程式の両辺を $t$ で微分すると…

In [21]:
Eq(diff(Friedmann.lhs, t), 
   diff(Friedmann.rhs, t))
Out[21]:
$\displaystyle -\frac{2 k \frac{d}{d t} a{\left(t \right)}}{a^{3}{\left(t \right)}} + \frac{2 \frac{d}{d t} a{\left(t \right)} \frac{d^{2}}{d t^{2}} a{\left(t \right)}}{a^{2}{\left(t \right)}} -\frac{2 \left(\frac{d}{d t} a{\left(t \right)}\right)^{3}}{a^{3}{\left(t \right)}} = \frac{8 \pi G \frac{d}{d t} \rho{\left(t \right)}}{3}$

Friedmann 式とその微分，そして att 式を使って，$\rho$ が満たす方程式を求める。

$\displaystyle \Bigl(\displaystyle \frac{d}{dt}$ (Friedmann) + $\displaystyle 2 \frac{\dot{a}}{a} \times$ (Friedmann) -$\displaystyle 2 \frac{\dot{a}}{a} \times$ (att) $\displaystyle \Bigr) \times \frac{3}{8 \pi G}$

In [22]:
LHS = ((diff(Friedmann.lhs, t)  
        + 2*diff(a,t)/a * Friedmann.lhs  
        - 2*diff(a,t)/a * att.lhs)* 3/(8*pi*G)).expand()

RHS = ((diff(Friedmann.rhs, t)  
        + 2*diff(a,t)/a * Friedmann.rhs  
        - 2*diff(a,t)/a * att.rhs)* 3/(8*pi*G)).expand()

Eq(RHS, LHS)
Out[22]:
$\displaystyle \frac{d}{d t} \rho{\left(t \right)} + \frac{3 \rho{\left(t \right)} \frac{d}{d t} a{\left(t \right)}}{a{\left(t \right)}} = 0$
In [23]:
eqrho = _
In [24]:
display(Friedmann)
display(att)
display(eqrho)
$\displaystyle \frac{k}{a^{2}{\left(t \right)}} + \frac{\left(\frac{d}{d t} a{\left(t \right)}\right)^{2}}{a^{2}{\left(t \right)}} = \frac{8 \pi G \rho{\left(t \right)}}{3} + \frac{\Lambda}{3}$
$\displaystyle \frac{\frac{d^{2}}{d t^{2}} a{\left(t \right)}}{a{\left(t \right)}} = -\frac{4 \pi G \rho{\left(t \right)}}{3} + \frac{\Lambda}{3}$
$\displaystyle \frac{d}{d t} \rho{\left(t \right)} + \frac{3 \rho{\left(t \right)} \frac{d}{d t} a{\left(t \right)}}{a{\left(t \right)}} = 0$

ということで，以下の方程式が得られた。

\begin{eqnarray}
\left(\frac{\dot{a}}{a}\right)^2 + \frac{k}{a^2} &=& \frac{8\pi G}{3} \rho + \frac{\Lambda}{3}\\
\frac{\ddot{a}}{a} &=& -\frac{4\pi G}{3} \rho + \frac{\Lambda}{3} \\
\dot{\rho} + 3 \frac{\dot{a}}{a} \rho &=& 0
\end{eqnarray}

圧力がある物質の場合

$P \neq 0$

EinEq(0, 0, P): $G^{0}_{\ \ 0} = 8\pi G T^{0}_{\ \ 0} -\Lambda \delta^{0}_{\ \ 0} $
In [25]:
EinEq(0,0,P)
Out[25]:
$\displaystyle -\frac{3 k}{a^{2}{\left(t \right)}} -\frac{3 \left(\frac{d}{d t} a{\left(t \right)}\right)^{2}}{a^{2}{\left(t \right)}} = -8 \pi G \rho{\left(t \right)} -\Lambda$

EinEq(0,0,0) * (-Rational(1,3)) すなわち

$\displaystyle -\frac{1}{3} G^{0}_{\ \ 0} = -\frac{8 \pi G}{3} + \frac{\Lambda}{3}$ を表示。

In [26]:
# 両辺を -3 で割る
Eq(EinEq(0,0,P).lhs * (-Rational(1,3)), 
   EinEq(0,0,P).rhs * (-Rational(1,3)))
Out[26]:
$\displaystyle \frac{k}{a^{2}{\left(t \right)}} + \frac{\left(\frac{d}{d t} a{\left(t \right)}\right)^{2}}{a^{2}{\left(t \right)}} = \frac{8 \pi G \rho{\left(t \right)}}{3} + \frac{\Lambda}{3}$
フリードマン方程式

Friedmann $=$ EinEq(0,0,P) * (-Rational(1,3)) とおき，Friedmann 方程式として参照する。フリードマン方程式は $P=0$ の場合と同じ。

In [27]:
# 上のセルの出力結果を変数 Friemann に代入し，
# Friedmann 方程式とする
Friedmann = _
EinEq(1, 1, P): $G^{1}_{\ \ 1} = 8\pi G T^{1}_{\ \ 1} -\Lambda \delta^{1}_{\ \ 1} $
In [28]:
EinEq(1, 1, P)
Out[28]:
$\displaystyle -\frac{k}{a^{2}{\left(t \right)}} -\frac{2 \frac{d^{2}}{d t^{2}} a{\left(t \right)}}{a{\left(t \right)}} -\frac{\left(\frac{d}{d t} a{\left(t \right)}\right)^{2}}{a^{2}{\left(t \right)}} = 8 \pi G P{\left(t \right)} -\Lambda$

この式と Friedmann の両辺を足して $-2$ でわり，$k$ と $\displaystyle \left(\frac{\dot{a}}{a}\right)^2$ の項を消去し，$\displaystyle \frac{\ddot{a}}{a}$ の式にする。

In [29]:
# 上の式と Friedmann 方程式を足して -2 で割る
Eq((Friedmann.lhs + EinEq(1,1,P).lhs)*(-Rational(1,2)), 
   (Friedmann.rhs + EinEq(1,1,P).rhs)*(-Rational(1,2)))
Out[29]:
$\displaystyle \frac{\frac{d^{2}}{d t^{2}} a{\left(t \right)}}{a{\left(t \right)}} = -4 \pi G P{\left(t \right)} -\frac{4 \pi G \rho{\left(t \right)}}{3} + \frac{\Lambda}{3}$
$\ddot{a}$ の式

上で得られた $\displaystyle \frac{\ddot{a}}{a}$ の式を
att とおく。

In [30]:
# 上のセルの出力結果を変数 att に代入し， 
# a の2階微分の方程式とする
att = _
エネルギー密度の式

Friedmann 式とその微分，そして att 式を使って，$\rho$ が満たす方程式を求める。

$\displaystyle \Bigl(\displaystyle \frac{d}{dt}$ (Friedmann) + $\displaystyle 2 \frac{\dot{a}}{a} \times$ (Friedmann) -$\displaystyle 2 \frac{\dot{a}}{a} \times$ (att) $\displaystyle \Bigr) \times \frac{3}{8 \pi G}$

In [31]:
LHS = ((diff(Friedmann.lhs, t)  
        + 2*diff(a,t)/a * Friedmann.lhs  
        - 2*diff(a,t)/a * att.lhs)* 3/(8*pi*G)).expand()

RHS = ((diff(Friedmann.rhs, t)  
        + 2*diff(a,t)/a * Friedmann.rhs  
        - 2*diff(a,t)/a * att.rhs)* 3/(8*pi*G)).expand()

Eq(RHS, LHS)
Out[31]:
$\displaystyle \frac{3 P{\left(t \right)} \frac{d}{d t} a{\left(t \right)}}{a{\left(t \right)}} + \frac{d}{d t} \rho{\left(t \right)} + \frac{3 \rho{\left(t \right)} \frac{d}{d t} a{\left(t \right)}}{a{\left(t \right)}} = 0$
In [32]:
eqrho = _
In [33]:
display(Friedmann)
display(att)
display(eqrho)
$\displaystyle \frac{k}{a^{2}{\left(t \right)}} + \frac{\left(\frac{d}{d t} a{\left(t \right)}\right)^{2}}{a^{2}{\left(t \right)}} = \frac{8 \pi G \rho{\left(t \right)}}{3} + \frac{\Lambda}{3}$
$\displaystyle \frac{\frac{d^{2}}{d t^{2}} a{\left(t \right)}}{a{\left(t \right)}} = -4 \pi G P{\left(t \right)} -\frac{4 \pi G \rho{\left(t \right)}}{3} + \frac{\Lambda}{3}$
$\displaystyle \frac{3 P{\left(t \right)} \frac{d}{d t} a{\left(t \right)}}{a{\left(t \right)}} + \frac{d}{d t} \rho{\left(t \right)} + \frac{3 \rho{\left(t \right)} \frac{d}{d t} a{\left(t \right)}}{a{\left(t \right)}} = 0$

ということで，$P \neq 0$ の場合に以下の方程式が得られた。

\begin{eqnarray}
\left(\frac{\dot{a}}{a}\right)^2 + \frac{k}{a^2} &=& \frac{8\pi G}{3} \rho + \frac{\Lambda}{3}\\
\frac{\ddot{a}}{a} &=& -\frac{4\pi G}{3} (\rho + 3 P) + \frac{\Lambda}{3} \\
\dot{\rho} + 3 \frac{\dot{a}}{a} (\rho + P) &=& 0
\end{eqnarray}