Python（SymPy）を使って学ぶSchwarzschildブラックホール解

今回は Python（SymPy）を活用する具体例として、４次元時空のSchwarzschild（シュワルツシルト）解がアインシュタイン方程式の解になっていることをPython（SymPy）を使って確認してみたいと思います。シュワルツシルト解は静的で球対称な場合のアインシュタイン方程式の真空解として求められます。

一般に、アインシュタイン方程式は時空の計量テンソルを偏微分してクリストッフェル記号、それを複雑に組み合わせた曲率テンソルを計算しないといけません。今回計算するような静的で球対称な場合には、仮定する変数依存性が少ないので手計算でもなんとかなるかもしれませんが、計算過程で求めなければいけない量が多く、そういった面倒な計算にはコンピュータを使ってラクをするのが良いでしょう。PythonのSymPyは、そのようなことを可能にする数式処理ライブラリです。

次の記事の中ではPython（SymPy）の使い方のごく初歩的な部分を紹介しています：

pianofisica.hatenablog.com

本記事では、４次元時空のSchwarzschild解をPython（SymPy）を使って求めます。おまけとして、５次元の場合についても同様に計算してみます。（考察する対象を容易に拡張できるのがコンピュータを使うメリットです）６次元の場合についても答えだけ載せておきました。
　

４次元時空のリーマン幾何
４次元Schwarzschild解
- Einstein方程式
５次元Schwarzschild解
６次元Schwarzschild解

４次元時空のリーマン幾何

以下、同じ項のなかで同じ文字の添字が上下に現れた場合、その添字については０から３までの和をとるものと約束します。（この約束をアインシュタインの規約といいます）
　

座標系の設定

４次元時空の座標系（いま３次元空間部分については極座標系をとることにします）

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

はSymPyで

import sympy as sp

sp.var('t, r, theta, phi')

xi = sp.Matrix([ [t],[r],[theta],[phi] ])

と入力します。後で他の次元に拡張しやすくするために、時空の次元を文字Dとしておきましょう。

D = len(xi)

計量テンソルの設定

上で導入した座標系 $\{x^i\}$ に付随する時空の計量テンソル $g_{ij}$ は不変距離 $ds^2$ に対して

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

として定義されます。計量テンソルは２階対称テンソル場なので、しばしば対称行列で次のように表されます。

$\qquad\displaystyle{g_{ij}=\left(\begin{array}{cccc} g_{00} & g_{01} & g_{02} & g_{03} \\ g_{10} & g_{11} & g_{12} & g_{13} \\ g_{20} & g_{21} & g_{22} & g_{23} \\ g_{30} & g_{31} & g_{32} & g_{33} \end{array} \right)}$
　
いま計量テンソルとして次の形を仮定してみましょう。

$\qquad\displaystyle{ds^2=-f(r)dt^2 + \frac{1}{f(r)}dr^2 + r^2d\theta^2 + r^2\sin^2\theta d\phi^2 }$

すなわち、計量の非自明な成分は以下の通りです。

$\qquad\displaystyle{g_{tt}=-f(r)\qquad g_{rr}=\frac{1}{f(r)} \qquad g_{\theta\theta} = r^2 \qquad g_{\phi\phi} = r^2\sin^2\theta }$

ただし $f=f(r)$ は動径 $r$ の関数で、その関数形は後でアインシュタイン方程式から決めます。これらをSymPyで入力するには

f = sp.Function('f')(r)

g = sp.Matrix([
[-f, 0, 0, 0],
[0, 1/f, 0, 0],
[0, 0, r**2, 0],
[0, 0, 0, r**2*sp.sin(theta)**2] ])

とします。 $g^{ik}g_{kj}=\delta^i_j$ として導入される上付き添字の $g^{ij}$ は

ginv = g.inv()

で計算されます。
　

Christoffel記号

天下りですが、Levi-Civita接続を定めるChristoffel記号は以下の量です。

$\qquad\displaystyle{\Gamma^i_{jk}=\frac12g^{in}\left(\frac{\partial g_{nj}}{\partial x^k}+\frac{\partial g_{nk}}{\partial x^j}-\frac{\partial g_{jk}}{\partial x^n}\right)}$

以下の曲率なども同様ですが、ここでは幾何学的な量の詳細については割愛します。理解したいという方は、以下のような一般相対論の教科書を参照してください。

Christoffel記号はSymPyでは次のようにして計算できます。

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

例えば

Gamma(1,3,3)

として

$\qquad\displaystyle{\Gamma^1_{33}=-rf(r)\sin^2\theta}$

などと求まります。
　

Riemannの曲率テンソル

リーマン曲率テンソルは以下の量です。

$\qquad\displaystyle{R^i_{jkl}=\frac{\partial\Gamma^i_{lj}}{\partial x^k}-\frac{\partial\Gamma^i_{kj}}{\partial x^l}+\Gamma^n_{lj}\Gamma^i_{kn}-\Gamma^n_{kj}\Gamma^i_{ln}}$

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

例えば

R(1,3,1,3)

として

$\qquad\displaystyle{R^1_{313}=-\frac{r}{2}f'(r)\sin^2\theta}$

などと求まります。
　

Ricciテンソル

リッチテンソルは

$\qquad\displaystyle{R_{jk}=R^n_{jnk}}$

です。

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)] ])
sp.simplify(RiccMat)

より

${R_{jk}=\left(\begin{array}{cccc}\frac{\left(r \frac{d^{2}}{d r^{2}} f{\left(r \right)} + 2 \frac{d}{d r} f{\left(r \right)}\right) f{\left(r \right)}}{2 r} & 0 & 0 & 0\\0 & \frac{- \frac{r \frac{d^{2}}{d r^{2}} f{\left(r \right)}}{2} - \frac{d}{d r} f{\left(r \right)}}{r f{\left(r \right)}} & 0 & 0\\0 & 0 & - r \frac{d}{d r} f{\left(r \right)} - f{\left(r \right)} + 1 & 0\\0 & 0 & 0 & \left(- r \frac{d}{d r} f{\left(r \right)} - f{\left(r \right)} + 1\right) \sin^{2}{\left(\theta \right)}\end{array}\right)}$

と計算されます。
　

Ricciスカラー

リッチスカラーは以下の量です。

$\qquad\displaystyle{R=g^{jk}R_{jk}}$

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

より

$\qquad\displaystyle{R=\frac{- r^{2} \frac{d^{2}}{d r^{2}} f{\left(r \right)} - 4 r \frac{d}{d r} f{\left(r \right)} - 2 f{\left(r \right)} + 2}{r^{2}}}$

と計算されます。

Einsteinテンソル

アインシュタインテンソルは次の量です。

$\qquad\displaystyle{G_{ij}=R_{ij}-\frac12Rg_{ij}}$

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

より

$\displaystyle{G_{jk}=\left(\begin{array}{cccc}\frac{\left(- r \frac{d}{d r} f{\left(r \right)} - f{\left(r \right)} + 1\right) f{\left(r \right)}}{r^{2}} & 0 & 0 & 0\\0 & \frac{r \frac{d}{d r} f{\left(r \right)} + f{\left(r \right)} - 1}{r^{2} f{\left(r \right)}} & 0 & 0\\0 & 0 & \frac{r \left(r \frac{d^{2}}{d r^{2}} f{\left(r \right)} + 2 \frac{d}{d r} f{\left(r \right)}\right)}{2} & 0\\0 & 0 & 0 & \frac{r \left(r \frac{d^{2}}{d r^{2}} f{\left(r \right)} + 2 \frac{d}{d r} f{\left(r \right)}\right) \sin^{2}{\left(\theta \right)}}{2}\end{array}\right)}$

と計算されます。
　

４次元Schwarzschild解

Einstein方程式

計量テンソル $g_{ij}$ を決定する基本方程式が、次のEinstein方程式

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

です。（いま物質場も宇宙定数も存在しない場合を考えています）よって未知関数 $f=f(r)$ について

$\qquad\displaystyle{rf''+2f'=0\\ -rf'-f+1=0}$

が要請されます。ここで第１式は第２式から従うので、第２式だけ考えれば良いことがわかります。これは１階微分方程式であり、変数分離の方法を使う方針で

$\qquad\displaystyle{\frac{f'}{1-f}=\frac{1}{r}}$

と変形すると $a$ を定数として

$\qquad\displaystyle{f=1-\frac{a}{r}}$

が解であることが容易にわかります。この解によって定まるEinstein方程式の真空解をSchwarzschild解といいます。その特徴は、計量テンソルが $r=a$ のところで特異性を示すことであり、この動径距離をシュワルツシルト半径といいます。

５次元Schwarzschild解

５次元の場合に、時間 $t$ と４次元空間の極座標 $(r,\theta,\phi,\psi)$ からなる座標系で

$\displaystyle{g_{ij}=\left(\begin{array}{ccccc} g_{tt} & g_{tr} & g_{t\theta} & g_{t\phi} & g_{t\psi} \\ g_{rt} & g_{rr} & g_{r\theta} & g_{r\phi} & g_{r\psi}\\ g_{\theta t} & g_{\theta r} & g_{\theta\theta} & g_{\theta\phi} & g_{\theta\psi} \\ g_{\phi t} & g_{\phi r} & g_{\phi\theta} & g_{\phi\phi} & g_{\phi\psi} \\ g_{\psi t} & g_{\psi r} & g_{\psi\theta} & g_{\psi\phi} & g_{\psi\psi} \end{array} \right)=\left(\begin{array}{ccccc} -f(r) & 0 & 0 & 0 & 0 \\ 0 & \frac{1}{f(r)} & 0 & 0 & 0 \\ 0 & 0 & r^2 & 0 & 0 \\ 0 & 0 & 0 & r^2\sin^2\theta &0 \\ 0 & 0 & 0 &0 & r^2\sin^2\theta\sin^2\phi \end{array} \right)}$

として４次元の場合と同様の計算をします：

%reset
import sympy as sp

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

f = sp.Function('f')(r)
g = sp.Matrix([
[-f, 0, 0, 0, 0],
[0, 1/f, 0, 0, 0],
[0, 0, r**2, 0, 0],
[0, 0, 0, r**2*sp.sin(theta)**2, 0],
[0, 0, 0, 0, r**2*sp.sin(theta)**2*sp.sin(phi)**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(0,4)],
[Ricc(1,0), Ricc(1,1), Ricc(1,2), Ricc(1,3), Ricc(1,4)],
[Ricc(2,0), Ricc(2,1), Ricc(2,2), Ricc(2,3), Ricc(2,4)],
[Ricc(3,0), Ricc(3,1), Ricc(3,2), Ricc(3,3), Ricc(3,4)],
[Ricc(4,0), Ricc(4,1), Ricc(4,2), Ricc(4,3), Ricc(4,4)] ])

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

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

より、Einsteinテンソル

$\displaystyle{G_{ij}=\left(\begin{array}{ccccc} -\frac{3f(rf'+2f-2)}{2r^2} & 0 & 0 & 0 & 0 \\ 0 & \frac{3(rf'+2f-2)}{2r^2f} & 0 & 0 &0 \\ 0 & 0 & \frac{r^2f''+4rf'+2f-2}{2} & 0 & 0 \\ 0 & 0 & 0 & \frac{r^2f''+4rf'+2f-2}{2}\sin^2\theta & 0 \\ 0 & 0 & 0 & 0 & \frac{r^2f''+4rf'+2f-2}{2}\sin^2\theta\sin^2\varphi \end{array} \right)}$

がわかります。Einstein方程式 $G_{ij}=0$ より

$\qquad\displaystyle{rf'+2f-2=0\\ r^2f''+4rf'+2f-2=0}$

が関数 $f$ に課されます。この場合についても、第２式は第１式から自動的に満たされて

$\qquad\displaystyle{\frac{f'}{1-f}=\frac{2}{r}}$

という $f$ に関する１階常微分方程式を得ます。その解は $a$ を定数として

$\qquad\displaystyle{f(r)=1-\frac{a}{r^2}}$

であることがわかります。この関数形によって定まる時空の計量

$\displaystyle{g_{ij}=\left(\begin{array}{ccccc} -\left(1-\frac{a}{r^2}\right) & 0 & 0 & 0 & 0 \\ 0 & \left(1-\frac{a}{r^2}\right)^{-1} & 0 & 0 & 0 \\ 0 & 0 & r^2 & 0 & 0 \\ 0 & 0 & 0 & r^2\sin^2\theta &0 \\ 0 & 0 & 0 &0 & r^2\sin^2\theta\sin^2\varphi \end{array} \right)}$

は５次元Schwarzschild解とよばれます。

６次元Schwarzschild解

６次元の場合に、時間 $t$ と５次元空間の極座標 $(r,\theta,\phi,\psi,\chi)$ からなる座標系で

$\displaystyle{g_{ij}=\left(\begin{array}{cccccc} -f(r) & 0 & 0 & 0 & 0 & 0 \\ 0 & \frac{1}{f(r)} & 0 & 0 & 0 & 0 \\ 0 & 0 & r^2 & 0 & 0 & 0 \\ 0 & 0 & 0 & r^2\sin^2\theta &0 & 0 \\ 0 & 0 & 0 &0 & r^2\sin^2\theta\sin^2\varphi & 0 \\ 0 & 0 & 0 &0 & 0 & r^2\sin^2\theta\sin^2\varphi\sin^2\psi \end{array} \right)}$

の仮定のもとで上と同様の解析をしてみると、Einstein方程式から $f$ に対して条件

$\qquad\displaystyle{\frac{f'}{1-f}=\frac{3}{r}}$

が課されます。その解は

$\qquad\displaystyle{f(r)=1-\frac{a}{r^3}}$

で与えられます。対応する計量が６次元Schwarzschild解です。

一般の n 次元の場合は、これまでの結果から予想できるでしょうか。

　
いかがだったでしょうか。今回は数式処理ライブラリSymPyを使って煩雑な計算をコンピュータに任せてみました。一度正しいコードを書いてしまえば、容易に他の次元に拡張できてしまう点などもコンピュータを使う利点です。アインシュタイン方程式の厳密解には、Schwarzschild解以外にもいくつか知られているものがあります。それらについてもコードを拡張して確認してみても良いかもしれません。

ちなみに、宇宙項を加えたドジッター解、電場を加えたライスナー・ノルドシュトロム解を以下の記事で扱っています。あわせて読んでみて下さい。

pianofisica.hatenablog.com

キーワード：Python、SymPy、Riemann幾何、一般相対論、Schwarzschild解