Pythonで一般相対性理論：導入編

本記事では、Pythonを使って、一般相対性理論で必要とされる
　・計量 $g_{\mu \nu}$をはじめとする各種テンソルの計算
　・Einstein方程式の簡単な計算
を行う方法について導入する。

[文献]
一般相対論におけるテンソル解析や代数計算に関するコンピューティングについては、
　Heinicke, C., et.al.- Computer Algebra in Gravity
　Korolkova, A., et.al.- Tensor computations in computer algebra systems
などの手頃なReview論文があるので、こちらも適宜参照されたい。
また、SageMathを用いた同様の研究に、
　Gourgoulhon, E., et.al. - Symbolic tensor calculus on manifolds
などが知られる。

#GraviPy モジュールの導入
GraviPyはPython３上で動くテンソル計算用のモジュール。

代数計算の記号処理に特化したSymPyと協働することで、ストレスフリーなテンソル計算の環境をPython上で実現してくれる。

以下に従って、SymPy含むGraviPyモジュールダウンロードする
GraviPy, Tensor Calculus Package for General Relativity (Version 0.1.0) (2014)
(アクセス日：2019年11月26日)

実際にはpipを使えば良く、

$ pip install GraviPy

で問題なければインストールされる。SymPyはPython環境がver.３.７以上であれば同時にインストールされているはず。

##シュバルツシルト計量を計算してみた

GraviPyを使ってシュバルツシルト計量からシュバルツシルト解を導出する。
###シュバルツシルト計量の導入
計量の符号をMTWに従って$(+,-,-,-)$とし、線素を次のように決める。
$$ds^2=g_{\mu\nu}dx^\mu dx^\nu .$$
重力定数と光速度について、$G=c=1$をとれば、シュバルツシルト計量は次のように書かれる。

g_{\mu\nu}=\left[\begin{array}{cccc}
\Big( 1- \frac{2M}{r} \Big) & 0 & 0 & 0	\\
0 & -\Big( 1- \frac{2M}{r} \Big)^{-1} & 0 & 0	\\
0 & 0 & -r^2 & 0	\\
0 & 0 & 0 & -r^2sin^2\theta
\end{array}\right]. 

ここで座標系として四次元座標$(t, r,\theta,\phi)$をとり、$M$をブラックホール質量とした。シュバルツシルト計量は動径方向について$r=２M$において特異点を示し、これがいわゆる「事象の地平」(Event Horizon)として知られる。

###GraviPyを用いた導出の仕方
上記の$g_{\mu\nu}$からシュバルツシルト解を導出していこう。
GraviPyを用いて以下のプロセスに従って計算していけばよい。

1. 時空の定義（計量の決定）：$g_{\mu\nu}$
2. クリストッフェル記号の計算：$\Gamma {}^\mu{}_{\nu\rho}$
3. リッチテンソルの計算：$R_{\mu\nu}$
4. アインシュタインテンソルの計算：$G_{\mu\nu}$

時空の定義

はじめに、時空を定義しよう。四次元時空座標として、$$x=x^\mu=(t,r,\theta,\phi)$$をとる。GraviPyでは次のように定義する。

GR.py

#!/usr/bin/env python3
from gravipy import *
from gravipy import tensorial as ten  
from sympy   import *
import inspect 

# Coordinates (\ chi is the four - vector of coordinates )
t, r, theta, phi, M = symbols('t , r , theta , phi , M ')
x = ten.Coordinates('\chi',[t, r, theta, phi])

ここで、ten.Coordinates()の第一引数'\chi'は第二引数が四元ベクトルであることを指定している。また、ten.~としたのはこのCoordinates()が明示的にtensorialを参照してやらないと動かない為。

以降、出現するテンソルは軒並みten.~で参照する。

###シュバルツシルト計量の定義
次に計量メトリック$g_{\mu \nu}$を定義しよう。先に述べたシュバルツシルト計量を次のようの記述してやればよい。

GR.py

# 続き
# Metric tensor
Metric = diag((1 -2* M / r ) , -1/(1 -2* M / r ) , -r **2 , -r **2* sin( theta ) **2)
g = ten.MetricTensor('g', x , Metric )

ここで、g = ten.MetricTensor('g',x,Metric)の出力結果は

Metric
#　あるいは
g(ten.All,ten.All)

で表示することができ、次のようになる。

g_{\mu\nu}=
\displaystyle \left[\begin{matrix}- 2M/r + 1 & 0 & 0 & 0\\0 & \displaystyle- \frac{1}{- 2M/r + 1} & 0 & 0\\0 & 0 & - r^{2} & 0\\0 & 0 & 0 & - r^{2} \sin^{2}{\left(\theta \right)}\end{matrix}\right] 

若干、表式に差異はあるが、先ほどの定義が正確に反映されている。
例えば、動径方向の成分$g_{rr}=g_{11}$を参照する場合には

g(1,1)

とすればよく、次のように要素別で取得することができる。

g_{11}=- \frac{1}{\displaystyle- \frac{2M}{r} + 1} 

###クリストッフェルの計算
続いて、$g_{\mu \nu}$をもとにして、クリストッフェル記号$\Gamma^\mu{}_{\nu\rho}$を計算してみよう。

GR.py

# 続き
# Christoffel symbol
Ga = ten.Christoffel('Ga', g )

Gaが$g_{\mu\nu}$に対するクリストッフェル記号である。クリストッフェルは３階のテンソルであるから、$\mu=0$に制限してやると次のような2階のテンソルを得る。

\Gamma^{ 0}{}_{\nu\rho}=
\left[
\begin{matrix}
0 & \frac{M}{r^{2}} & 0 & 0 \\
\frac{M}{r^{2}} & 0 & 0 & 0 \\
0 & 0 & 0 & 0\\
0 & 0 & 0 & 0
\end{matrix}\right]

全成分を得るには

Ga(ten.All,ten.All,ten.All)

としてやればよく、

\displaystyle
\Gamma^\mu{}_{\nu\rho}=
 \left[\begin{matrix}\left[\begin{matrix}0 & \frac{M}{r^{2}} & 0 & 0\\\frac{M}{r^{2}} & 0 & 0 & 0\\0 & 0 & 0 & 0\\0 & 0 & 0 & 0\end{matrix}\right] & \left[\begin{matrix}- \frac{M}{r^{2}} & 0 & 0 & 0\\0 & \frac{M}{\left(2 M - r\right)^{2}} & 0 & 0\\0 & 0 & r & 0\\0 & 0 & 0 & r \sin^{2}{\left(\theta \right)}\end{matrix}\right] & \left[\begin{matrix}0 & 0 & 0 & 0\\0 & 0 & - r & 0\\0 & - r & 0 & 0\\0 & 0 & 0 & \frac{r^{2} \sin{\left(2 \theta \right)}}{2}\end{matrix}\right] & \left[\begin{matrix}0 & 0 & 0 & 0\\0 & 0 & 0 & - r \sin^{2}{\left(\theta \right)}\\0 & 0 & 0 & - \frac{r^{2} \sin{\left(2 \theta \right)}}{2}\\0 & - r \sin^{2}{\left(\theta \right)} & - \frac{r^{2} \sin{\left(2 \theta \right)}}{2} & 0\end{matrix}\right]\end{matrix}\right]

となる。3階のテンソルなので、形式上「2次元配列の1次元配列」（＝１次元配列の要素が2次元配列になっている形式）として得ることができる。

###リッチテンソルの計算
次に、リッチテンソルを計算する。（この名称は数学者Ricciに由来する。）

# Ricci tensor
Ri = ten.Ricci ('Ri ', g )
# Display all compon
Ri(ten.All,ten.All)

シュバルツシルト計量におけるRicciテンソルは全成分がゼロに落ちる。

R_{\mu\nu}=
\displaystyle \left[\begin{matrix}0 & 0 & 0 & 0\\0 & 0 & 0 & 0\\0 & 0 & 0 & 0\\0 & 0 & 0 & 0\end{matrix}\right]

もちろん、これは特殊なことであり、シュバルツシルト計量の表式で
$$r \longrightarrow r^{3.5}$$
などと置換してリッチテンソルを計算してみると、全く異なる値になるので実験してみるとよい。
###アインシュタインテンソルの計算
最後に、リッチテンソルを用いてアインシュタインテンソルを出力してみよう。

# Einstein tensor
G = ten.Einstein ('G', Ri )
G(ten.All,ten.All)

結果は、

G_{\mu\nu}=
\displaystyle \left[\begin{matrix}0 & 0 & 0 & 0\\0 & 0 & 0 & 0\\0 & 0 & 0 & 0\\0 & 0 & 0 & 0\end{matrix}\right]

となっている。この結果は、シュバルツシルト解の仮定である、真空条件
$$T_{\mu\nu}=0$$
にも適合しており、アインシュタイン方程式
$$G_{\mu\nu}=8\pi T_{\mu\nu}$$
を満たしている。

#参考文献
GraviPyのチュートリアルが公開されており、以下を参考にした。

　https://github.com/wojciechczaja/GraviPy

また、SymPyのサポートページは以下にある。

　[https://docs.sympy.org/dev/index.html]
(https://docs.sympy.org/dev/index.html)