EinsteinPy を使ってフリードマン・ルメートル・ロバートソン・ウォーカー (FLRW) 計量からアインシュタイン方程式
$$G^{\mu}_{\ \ \nu} + \Lambda \delta^{\mu}_{\ \ \nu} = 8\pi G T^{\mu}_{\ \ \nu}$$
を使ってフリードマン方程式を求める。
必要なパッケージの import
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}$ と定義。
# スケール因子 a(t) は時間のみの関数
a = Function('a')(t)
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()
アインシュタイン・テンソル
ein $\equiv \displaystyle G^{\mu}_{\ \ \nu} = R^{\mu}_{\ \ \nu} – \frac{1}{2} R \delta^{\mu}_{\ \ \nu} $ と定義。（.change_config('ul') で上付下付に）
ein = EinsteinTensor.from_metric(g).change_config('ul')
ein $ = G^{\mu}_{\ \ \nu}$ を4行4列の行列として表示。
# アインシュタインテンソルの全成分表示は .tensor() をつけて...
ein.tensor()
アインシュタインテンソル $G^{\mu}_{\ \ \nu}$ のゼロでない成分を個別に表示する。
まず，ein[0, 0] $= G^{0}_{\ \ 0}$ を表示。
ein[0,0]
次に，ein[1, 1] $= G^{1}_{\ \ 1}$ を表示。
ein[1,1]
ein[2, 2] $= G^{2}_{\ \ 2}$ と，
ein[2,2]
ein[3, 3] $= G^{3}_{\ \ 3}$ も表示。
ein[3,3]
念のため，$G^{2}_{\ \ 2}$ と $G^{3}_{\ \ 3}$ が等しいことを確認します。
ein[2, 2] == ein[3, 3]
完全流体のエネルギー運動量テンソル
\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}$
と定義。
# エネルギー密度 ρ
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)
アインシュタイン方程式
$$G^{\mu}_{\ \ \nu} = 8\pi G T^{\mu}_{\ \ \nu} – \Lambda \delta^{\mu}_{\ \ \nu} $$
# 宇宙定数 Λ と sympy.core.function.Lambda がかち合うので
# あらためて変数として宣言
var('Lambda')
EinEq(mu, nu, P) $ \equiv G^{\mu}_{\ \ \nu} = 8\pi G T^{\mu}_{\ \ \nu} – \Lambda \delta^{\mu}_{\ \ \nu}$ と定義。
def EinEq(mu, nu, P):
return Eq(expand(ein[mu, nu]), 8*pi*G*T(mu, nu, P) - Lambda *KroneckerDelta(mu, nu))
ダスト物質の場合
$P = 0$
フリードマン方程式
$G^{0}_{\ \ 0} = 8\pi G T^{0}_{\ \ 0} – \Lambda \delta^{0}_{\ \ 0} $
EinEq(0, 0, 0) $= G^{0}_{\ \ 0}$ を表示。
EinEq(0,0,0)
方程式 Eq() の左辺は Eq().lhs，右辺は Eq().rhs として参照できます。以下では両辺を $-3$ で割ります。
EinEq(0,0,0) * (-Rational(1,3)) すなわち
$\displaystyle -\frac{1}{3} G^{0}_{\ \ 0} = -\frac{8 \pi G}{3} + \frac{\Lambda}{3}$ を表示。
# 両辺を -3 で割る
Eq(EinEq(0,0,0).lhs * (-Rational(1,3)),
EinEq(0,0,0).rhs * (-Rational(1,3)))
Friedmann $=$ EinEq(0,0,0) * (-Rational(1,3)) とおく。
# 上のセルの出力結果を変数 Friemann に代入し，
# Friedmann 方程式とする
Friedmann = _
$\ddot{a}$ の式
EinEq(1,1,0) すなわち $G^{1}_{\ \ 1} = 8\pi G T^{1}_{\ \ 1} – \Lambda \delta^{1}_{\ \ 1} $ を表示。
EinEq(1,1,0)
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}$ の式にする。
# 上の式と Friedmann 方程式を足して -2 で割る
Eq((Friedmann.lhs + EinEq(1,1,0).lhs)*(-Rational(1,2)),
(Friedmann.rhs + EinEq(1,1,0).rhs)*(-Rational(1,2)))
上で得られた $\displaystyle \frac{\ddot{a}}{a}$ の式を
att とおく。
# 上のセルの出力結果を変数 att に代入し，
# a の2階微分の方程式とする
att = _
エネルギー密度の式
Friedmann 方程式 Friedmann は $\displaystyle \left(\frac{\dot{a}}{a} \right)^2$ を含む式であり…
Friedmann
att は $\displaystyle \frac{\ddot{a}}{a}$ についての方程式であり…
att
Friedmann 方程式の両辺を $t$ で微分すると…
Eq(diff(Friedmann.lhs, t),
diff(Friedmann.rhs, t))
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}$
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)
ということで，以下の方程式が得られた。
\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$
フリードマン方程式
$G^{0}_{\ \ 0} = 8\pi G T^{0}_{\ \ 0} – \Lambda \delta^{0}_{\ \ 0} $
EinEq(0, 0, P) $= G^{0}_{\ \ 0}$ を表示。
EinEq(0,0,P)
EinEq(0,0,0) * (-Rational(1,3)) すなわち
$\displaystyle -\frac{1}{3} G^{0}_{\ \ 0} = -\frac{8 \pi G}{3} + \frac{\Lambda}{3}$ を表示。
# 両辺を -3 で割る
Eq(EinEq(0,0,P).lhs * (-Rational(1,3)),
EinEq(0,0,P).rhs * (-Rational(1,3)))
Friedmann $=$ EinEq(0,0,P) * (-Rational(1,3)) とおき，Friedmann 方程式として参照する。
# 上のセルの出力結果を変数 Friemann に代入し，
# Friedmann 方程式とする
Friedmann = _
$\ddot{a}$ の式
EinEq(1,1,P) すなわち $G^{1}_{\ \ 1} = 8\pi G T^{1}_{\ \ 1} – \Lambda \delta^{1}_{\ \ 1} $ を表示。
EinEq(1, 1, P)
Friedmann すなわち EinEq(0,0,P) * (-Rational(1,3)) と EinEq(1,1,P) の両辺を足して $-2$ でわり，$k$ と $\displaystyle \left(\frac{\dot{a}}{a}\right)^2$ の項を消去し，$\displaystyle \frac{\ddot{a}}{a}$ の式にする。
# 上の式と Friedmann 方程式を足して -2 で割る
Eq((Friedmann.lhs + EinEq(1,1,P).lhs)*(-Rational(1,2)),
(Friedmann.rhs + EinEq(1,1,P).rhs)*(-Rational(1,2)))
上で得られた $\displaystyle \frac{\ddot{a}}{a}$ の式を
att とおく。
# 上のセルの出力結果を変数 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}$
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)
ということで，以下の方程式が得られた。
\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}