Python（SymPy）を使って学ぶアインシュタイン方程式の膨張宇宙解と暗黒エネルギー

今回は Python（SymPy）を活用する具体例として、膨張する宇宙モデルを記述する４次元アインシュタイン方程式の解（フリードマン・ルメートル・ロバートソン・ウォーカー計量、略称FLRW計量）をPython（SymPy）を使って確認してみたいと思います。とくに膨張宇宙を議論するうえでの基礎方程式となる『フリードマン方程式』を導出します。また、現在観測されている膨張宇宙を説明する暗黒エネルギーについても簡単に触れます。

アインシュタイン方程式の解のいくつかは、過去の記事で見ています。あわせて眺めてみてください：

pianofisica.hatenablog.com
pianofisica.hatenablog.com
pianofisica.hatenablog.com
　

４次元FLRW計量の具体形
物質場を含めた場合
フリードマン方程式
エネルギー保存則
加速膨張宇宙とダークエネルギー

４次元FLRW計量の具体形

本記事内で未定義の量や用語の詳細については、上でリンクを貼った過去の記事を参照してください。

真空中のEinstein方程式は

$\qquad\displaystyle{G_{ij}=0}$

で与えられます。ここで $G_{ij}$ はEinsteinテンソルです。FLRW計量は、時刻 $t$ の未知関数 $a(t)$ （スケール因子）によって

$\qquad\displaystyle{g_{ij}=\left(\begin{array}{cccc} -c^2& 0 & 0 & 0 \\ 0 & a(t)^2\frac{1}{1-kr^2} & 0 & 0 \\ 0 & 0 & a(t)^2r^2 & 0 \\ 0 & 0 & 0 & a(t)^2r^2\sin^2\theta \end{array} \right)}$

で与えられる解です。ただし、これらの成分は座標系

$\qquad\displaystyle{\xi^{i}=\left(\begin{array}{c} t \\ r \\ \theta \\ \phi \end{array} \right)}$

に関するものとします。つまり

$\qquad\displaystyle{ds^2=g_{ij}\,d\xi^{i}d\xi^{j} }$

です。ここでパラメタ $k$ は時空の曲率に関する量で、曲率が正（閉じている）・ゼロ（平坦）・負（開いている）の場合に応じてそれぞれ $k=+1,0,-1$ これらをSymPyに入力します。

%reset
import sympy as sp

sp.var('c, k, t, r, theta, phi')
xi = sp.Matrix([ [t],[r],[theta],[phi] ])
D = len(xi)

a = sp.Function('a')(t)
g = sp.Matrix([
[-c**2, 0, 0, 0],
[0, a**2/(1-k*r**2), 0, 0],
[0, 0, a**2 * r**2, 0],
[0, 0, 0, a**2 * r**2 * sp.sin(theta)**2] ])

ginv = g.inv()

def Gamma(k,i,j):
    gm = 0
    for n in range(D):
        gm = gm + sp.Rational(1,2)*ginv[k,n]*(g[n,i].diff(xi[j])\
                +g[n,j].diff(xi[i])-g[i,j].diff(xi[n]))
    return gm

def R(i,j,k,l):
    rm = Gamma(i,l,j).diff(xi[k]) - Gamma(i,k,j).diff(xi[l]) 
    for n in range(D):
        rm = rm + Gamma(n,l,j)*Gamma(i,k,n) - Gamma(n,k,j)*Gamma(i,l,n)
    return rm

def Ricc(j,k):
    ric = 0
    for n in range(D):
        ric = ric + R(n,j,n,k)
    return ric

RiccMat = sp.Matrix([
[Ricc(0,0), Ricc(0,1), Ricc(0,2), Ricc(0,3)],
[Ricc(1,0), Ricc(1,1), Ricc(1,2), Ricc(1,3)],
[Ricc(2,0), Ricc(2,1), Ricc(2,2), Ricc(2,3)],
[Ricc(3,0), Ricc(3,1), Ricc(3,2), Ricc(3,3)] ])

R = sp.simplify((ginv*RiccMat).trace())

GMat = sp.simplify(RiccMat - sp.Rational(1,2)*g*R)
GMat

より、Einsteinテンソルが

$\displaystyle{G_{ij}=\begin{pmatrix}\frac{3 \left(c^{2} k + \left(\frac{d}{d t} a{\left(t \right)}\right)^{2}\right)}{a^{2}{\left(t \right)}} & 0 & 0 & 0\\0 & \frac{c^{2} 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}}{c^{2} \left(k r^{2} - 1\right)} & 0 & 0\\0 & 0 & \frac{r^{2} \left(- c^{2} 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}\right)}{c^{2}} & 0\\0 & 0 & 0 & \frac{r^{2} \left(- c^{2} 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}\right) \sin^{2}{\left(\theta \right)}}{c^{2}}\end{pmatrix} }$

と求まります。Einstein方程式から、未知関数 $a$ について

$\qquad\displaystyle{\left(\frac{\dot{a}}{a}\right)^2+\frac{kc^2}{a^2}=0\\ 2\frac{\ddot{a}}{a}+\left(\frac{\dot{a}}{a}\right)^2+\frac{kc^2}{a^2}=0}$

が要請されます。ここで $\dot{a}=da/dt$ は未知関数 $a$ の時刻 $t$ に関する導関数です。さて、第１式を使うと第２式は

$\qquad\displaystyle{2\frac{\ddot{a}}{a}=0}$

です。さらに第１式の両辺を時間微分すると

$\qquad\displaystyle{\frac{d}{dt}\left[\left(\frac{\dot{a}}{a}\right)^2+\frac{kc^2}{a^2}\right]=2\frac{\dot{a}}{a}\left[\frac{\ddot{a}}{a}-\left(\frac{\dot{a}}{a}\right)^2-\frac{kc^2}{a^2}\right]=\frac{\dot{a}}{a}\left(2\frac{\ddot{a}}{a}\right)=0}$

となって、第２式が自動的に導かれることがわかります。よって第１式だけ考えれば良いことがわかります。すなわち

$\qquad\displaystyle{\left(\frac{\dot{a}}{a}\right)^2+\frac{kc^2}{a^2}=0}$

です。分母の $a^2$ を払えば

$\qquad\displaystyle{\dot{a}^2+kc^2=0}$

です。時空の曲率が平坦な場合（ $k=0$ ）には $a={\rm const.}$ となって（スケール倍を除いて）Minkowski時空に一致します。曲率が負の場合（ $k=-1$ ）には

$\qquad\displaystyle{a=ct+a_0}$

であることがわかります。

物質場を含めた場合

完全流体の物質場による次のエネルギー・運動量テンソルを考えます：

$\qquad\displaystyle{\begin{aligned} T_{ij}&=P\,g_{ij}+(P\,c^2+\rho\,c^4)\,\delta_{i,0}\delta_{j,0}\\ &= P\left(\begin{array}{cccc} -c^2& 0 & 0 & 0 \\ 0 & a(t)^2\frac{1}{1-kr^2} & 0 & 0 \\ 0 & 0 & a(t)^2r^2 & 0 \\ 0 & 0 & 0 & a(t)^2r^2\sin^2\theta \end{array} \right)+(P\,c^2+\rho\,c^4)\left(\begin{array}{cccc} 1& 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{array} \right)\\ &= \left(\begin{array}{cccc} \rho\,c^4& 0 & 0 & 0 \\ 0 & P\,a(t)^2\frac{1}{1-kr^2} & 0 & 0 \\ 0 & 0 & P\,a(t)^2r^2 & 0 \\ 0 & 0 & 0 & P\,a(t)^2r^2\sin^2\theta \end{array} \right) \end{aligned}}$

ここで $P$ は圧力、 $\rho$ は密度です。圧力と密度を関係付ける式を状態方程式といいます。とくに状態方程式が $P=-\rho\,c^2$ となる場合、エネルギー・運動量テンソルは時空の計量に比例し、従って宇宙項（ $\Lambda$ ）と同一視できます。

物質場が存在するとき、アインシュタイン方程式は

$\qquad\displaystyle{G_{ij}+\Lambda\,g_{ij}=\frac{8\pi G}{c^4}T_{ij}}$

と変更されます。ここでは宇宙定数 $\Lambda$ も含めました。これらをSymPyに入力して計算します：

%reset
import sympy as sp

sp.var('c, k, t, r, theta, phi, P, rho, G, Lambda')
xi = sp.Matrix([ [t],[r],[theta],[phi] ])
D = len(xi)

a = sp.Function('a')(t)
P = sp.Function('P')(t)
rho = sp.Function('rho')(t)
g = sp.Matrix([
[-c**2, 0, 0, 0],
[0, a**2/(1-k*r**2), 0, 0],
[0, 0, a**2 * r**2, 0],
[0, 0, 0, a**2 * r**2 * sp.sin(theta)**2] ])
E00 = sp.Matrix([
[1, 0, 0, 0],
[0, 0, 0, 0],
[0, 0, 0, 0],
[0, 0, 0, 0] ])
T = P*g + (P*c**2 + rho*c**4)*E00

ginv = g.inv()

def Gamma(k,i,j):
    gm = 0
    for n in range(D):
        gm = gm + sp.Rational(1,2)*ginv[k,n]*(g[n,i].diff(xi[j])\
                +g[n,j].diff(xi[i])-g[i,j].diff(xi[n]))
    return gm

def R(i,j,k,l):
    rm = Gamma(i,l,j).diff(xi[k]) - Gamma(i,k,j).diff(xi[l]) 
    for n in range(D):
        rm = rm + Gamma(n,l,j)*Gamma(i,k,n) - Gamma(n,k,j)*Gamma(i,l,n)
    return rm

def Ricc(j,k):
    ric = 0
    for n in range(D):
        ric = ric + R(n,j,n,k)
    return ric

RiccMat = sp.Matrix([
[Ricc(0,0), Ricc(0,1), Ricc(0,2), Ricc(0,3)],
[Ricc(1,0), Ricc(1,1), Ricc(1,2), Ricc(1,3)],
[Ricc(2,0), Ricc(2,1), Ricc(2,2), Ricc(2,3)],
[Ricc(3,0), Ricc(3,1), Ricc(3,2), Ricc(3,3)] ])

R = sp.simplify((ginv*RiccMat).trace())

GMat = sp.simplify(RiccMat - sp.Rational(1,2)*g*R)
GMat + Lambda*g - (8*sp.pi*G/c**4)*T

すると、未知関数 $a$ に関する方程式系

$\qquad\displaystyle{\left(\frac{\dot{a}}{a}\right)^2+\frac{k c^2}{a^2}-\frac{\Lambda c^2}{3}=\frac{8 \pi G}{3} \rho\\ 2 \frac{\ddot{a}}{a}+\left(\frac{\dot{a}}{a}\right)^2+\frac{k c^2}{a^2}-\Lambda c^2=-\frac{8 \pi G}{c^2} P}$

を得ます。

フリードマン方程式

上で得た方程式系の第１式

$\qquad\displaystyle{\left(\frac{\dot{a}}{a}\right)^2+\frac{k c^2}{a^2}-\frac{\Lambda c^2}{3}=\frac{8 \pi G}{3} \rho}$

をとくにフリードマン方程式といい、これは膨張宇宙を議論する上での基礎となります。実は、真空のときと同様に上の方程式系の第２式は第１式から従います。このことを見るために、まず第１式と第２式を組み合わせて

$\qquad\displaystyle{\frac{\ddot{a}}{a}-\frac{\Lambda c^2}{3}=-\frac{4 \pi G}{3c^2} (c^2\rho+3P)\qquad (\ast)}$

としておきます。次に第１式を時間微分してみます。

$\qquad\displaystyle{\begin{aligned}0&=\frac{d}{dt}\left[\left(\frac{\dot{a}}{a}\right)^2+\frac{k c^2}{a^2}-\frac{\Lambda c^2}{3}-\frac{8 \pi G}{3} \rho\right]=2\frac{\dot{a}}{a}\left[\frac{\ddot{a}}{a}-\left(\frac{\dot{a}}{a}\right)^2-\frac{kc^2}{a^2}-\frac{4 \pi G}{3}\frac{a}{\dot{a}} \dot{\rho}\right]\\&=2\frac{\dot{a}}{a}\left[\frac{\ddot{a}}{a}-\frac{\Lambda c^2}{3}-\frac{8 \pi G}{3c^2}\,c^2 \rho-\frac{4 \pi G}{3c^2}\frac{a}{\dot{a}} c^2\dot{\rho}\right]\end{aligned}}$

ここで最後の等号にふたたび第１式を用いています。さて、これが（第１式と組み合わせた）第２式（式 $(\ast)$ ）と整合して自明となるためには

$\qquad\displaystyle{3\frac{\dot{a}}{a}(c^2\rho+P)+c^2\dot{\rho}=0}$

であれば良いことがわかります。実際、この式は物質場のエネルギー保存則から従います（このことは次の節で示します）。よって第２式は、第１式（とエネルギー保存則）から導かれることが示されました。

エネルギー保存則

上の議論で用いた、エネルギー保存則による関係式を示しましょう。それは曲がった時空中の連続の式

$\qquad\displaystyle{\nabla^kT_{k0}=0}$

を満たすというものです。共変微分を具体的に書き下すと

$\qquad\displaystyle{0=g^{kl}\nabla_kT_{l0}=g^{kl}\left(\partial_kT_{l0} -\Gamma^m_{kl}T_{m0}-\Gamma^m_{k0}T_{lm}\right)}$

となります。これを計算してみましょう：

eq = 0
for k in range(D):
    for l in range(D):
        for m in range(D):
            eq = eq + ginv[k,l]*( - Gamma(m,k,l)*T[m,0] - Gamma(m,k,0)*T[l,m] )
        eq = eq + ginv[k,l]*( T[l,0].diff(xi[k]) )
eq