「Pythonプログラミング」（その他のトピック・重力レンズ）

Pythonプログラミング（重力レンズ）

このページでは、重力による光の湾曲のシミュレーションを試みる（ここはこれから）。

１．球対称な重力場：一般相対性理論からの結果

空間は極座標を考え、時空の座標を $x^0 = t,\ x^1=r,\ x^2=\theta,\ x^3=\phi$ で表すことにする。球対称な重力場について、アインシュタイン方程式を解いた結果はSchwarzschild解として知られており、計量（時空の「長さ」）$ds$ が $$ ds^2 = \left(1 - \frac{r_g}{r} \right) c^2 dt^2 - \frac{1}{\left( 1 - \frac{r_g}{r} \right)} dr^2 - r^2 d\theta^2 - r^2 \sin^2 \theta \phi^2 $$ で与えられる。ここで、$r_g = G M/c^2$ ($G$は重力定数、$M$は質量、$c$は光速）。対応する計量テンソルの成分は $$ \begin{eqnarray} g_{00} & = & \left(1 - \frac{r_g}{r} \right), \, g^{00} = \left(1 - \frac{r_g}{r} \right)^{-1}\\ g_{11} & = & - \left(1 - \frac{r_g}{r} \right)^{-1}, \, g^{11} = - \left(1 - \frac{r_g}{r} \right)\\ g_{22} & = & - r^2, \, g^{22} = - 1/r^2 \\ g_{33} & = & - r^2 \sin^2 \theta, \, g^{33} = - \frac{1}{r^2 \sin^2 \theta } \end{eqnarray} $$ で、これら以外の成分はすべて0である。

クリストッフェルの記号 $$ \Gamma^i_{k\ell} = \frac{1}{2} g^{im} \left(\frac{\partial g_{mk}}{\partial x^\ell} + \frac{\partial g_{m\ell}}{\partial x^k} - \frac{\partial g_{k\ell}}{\partial x^m} \right) $$ を使い、座標の微分を $$ u^i = \frac{dx^i}{ds} $$ と置くと、測地線の方程式（光や物体の時空での経路）は、座標についての二階の（非線形）微分方程式 $$ \frac{du^i}{ds} + \Gamma^i_{k\ell} u^k u^\ell = 0 $$ で表すことができる。

簡単のため、$\theta=\pi/2$となるような面内に限って考えることにすると、$t, r, \phi$、すなわち$x^0, x^1, x^3$ についてのみ考慮すれば十分である。 $i=0$について、測地線の方程式を具体的に書き下すと $$ \frac{d^2t}{ds^2} - \frac{r_g}{r^2} \left(1-\frac{r_g}{r}\right)^{-1} \frac{dt}{ds}\frac{dr}{ds} = 0 $$ となるが、これは $$ \frac{d}{ds} \left\{ \left(1 - \frac{r_g}{r} \right) \frac{dt}{ds} \right\} = 0 $$ と書き直せる。すなわち、$A$を定数として、保存則 $$ \left(1 - \frac{r_g}{r} \right) \frac{dt}{ds} = A $$ が成り立つ。

$i=3$についても同様に測地線方程式を書き下す $$ \frac{d^2\phi}{ds^2} + \frac{2}{r} \frac{dr}{ds}\frac{d\phi}{ds} = 0 $$ となるが、 $$ \frac{d}{ds} \left\{r^2 \frac{d\phi}{ds} \right\} = 0 $$ から、２つ目の保存則 $$ r^2 \frac{d\phi}{ds} = B $$ が得られる（ケプラーの面積速度一定の法則に対応）。

２．球対称な重力場の周りでの光の経路

次に、上で求めた測地線に沿った光の経路について考えてみる。

重力場中の光の経路はヌル測地線$ds=0$となること、およびこの場合の計量から、 $s$による微分をドット記号で表すことにすると、 $$ \left(1 - \frac{r_g}{r} \right) \dot{t}^2 - \frac{1}{\left( 1 - \frac{r_g}{r} \right)} \dot{r}^2 - r^2 \dot{\theta}^2 - r^2 \sin^2 \theta \dot{\phi}^2 = 0 $$ となる。ただし、ここで、$\dot{t} = c \frac{dt}{ds}$とした。これを変形して $$ \dot{r}^2 = \left( \left(1 - \frac{r_g}{r} \right) \dot{t} \right)^2 - \frac{1 - \frac{r_g}{r}}{r^2} \left( r^2 \dot{\phi}^2 \right) $$ 先に求めた保存則を思い出すと、この式は $$ \dot{r}^2 = (cA)^2 - \frac{1 - \frac{r_g}{r}}{r^2} B^2 $$ と書き直すことができる。両辺をさらに$s$で微分すると $$ 2 \dot{r} \ddot{r} = \left( \frac{2}{r^3} \dot{r} - \frac{3 r_g}{r^4} \dot{r} \right) B^2 $$ であるから、 $$ \ddot{r} = \left(\frac{1}{r^3} - \frac{3 r_g}{2 r^4} \right) B^2 $$ を得る。

こうして、定数$B$を決め、$r$と$\dot{r}$の初期値を与えることで$r(s)$が、$B$についての保存則から$\phi(s)$が決まることが分かった。

３．重力レンズのシミュレーション

では、上で求めた微分方程式によって、球対称な重力場の周りの光の経路を「レイトレース」してみよう。質量は原点に位置し、そのSchwartzchild半径 $r_g$ を変数 rg で指定する。光は $(x_0,y_0)$ をx軸に平行に出発するとする。そのとき、$\frac{dx(0)}{ds}=1$である。そのとき、その他の変数の初期条件は $$ \begin{eqnarray} r(0)&=&\sqrt{{x_0}^2 + {y_0}^2} \\ \phi(0) &=& \arctan(y_0/x_0) \\ \frac{dr(0)}{ds} &=& \frac{dx}{ds} \cos(\phi) \\ \frac{d\phi(0)}{ds} &=& - \frac{dx}{ds} \, \sin(\phi) / r \end{eqnarray} $$ となる。

そこから出発して、単純な差分法によって、経路長が$s=20$まで積分するコードの例を以下に示す。原点に半径$r_g$の円を描くようにした。

# coding: utf-8

import math
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as patches

rg=1.0

x0 = -10
y0 = 3.0

ds = 0.001
dxds = 1.0

r = math.sqrt(x0**2 + y0**2)
phi = math.atan2(y0,x0)
drds = dxds * math.cos(phi)
dphids = - dxds*math.sin(phi)/r
B = r**2 * dphids

R=[]
PHI=[]

s=0
while s<20:
    r += drds * ds    
    d2rds2 = (1/r**3 - 3*rg/(2*r**4))*B**2
    drds += d2rds2 * ds
    phi += dphids * ds    
    dphids = B/r**2
    s += ds
    PHI.append(phi)
    R.append(r)
    
N=len(R)
xs=[]
ys=[]
for i in range(N):
    xs.append( R[i]*math.cos(PHI[i]) )
    ys.append( R[i]*math.sin(PHI[i]) )

fig, ax = plt.subplots()    
ax.plot(xs,ys, color=(1.0,1.0,0.0), linewidth=2.0)
ax.set_xlabel('X')
ax.set_ylabel('Y')
circle = patches.Circle(xy=(0, 0), radius=rg, facecolor="pink", edgecolor="red",)
ax.add_patch(circle)
ax.set_xlim(-10,10)
ax.set_ylim(-10,10)
ax.grid(True)
plt.show()

　練習：解法の改良

上記のサンプルコードは経路が$r_g$に近づくにつれて計算誤差が増大したり、挙動が不安定になったりする。そこで、より精度の高いスキームを用いるように変更してみなさい。