「シミュレーションで学ぶ相対論入門」を参考にブラックホールのシミュレーションを作っていきます。

測地線方程式
- 粒子の時間による時間発展方程式
- 観測者の時間による時間発展方程式

測地線方程式

シミュレーションのためにシュワルツシルトブラックホール周りの光子や質点の運動を求めたいので、それを記述する測地線方程式（時間発展方程式）を求めます。シュワルツシルト計量は以下の通り。

$\displaystyle{ (cd\tau)^2=-\left(1-\frac{r_s}{r}\right)(cdt)^2+\frac{dr^2}{1-\frac{r_s}{r}}+r^2(d\theta^2+\sin^2\theta d\phi^2) }$

$r_s$ はシュワルツシルト半径です。ラグランジアン $\mathcal{L}$ は

$\displaystyle{ \mathcal{L} =-c^2\left(1-\frac{r_s}{r}\right)\dot{\tau}^2+\frac{\dot{r}^2}{1-\frac{r_s}{r}}+r^2(\dot{\tau}^2+\sin^2\theta \dot{\phi}^2)=c^2\epsilon \tag{1} }$

ただし、ドット( $\cdot$ )は固有時 $\tau$ による微分を表しています。 $\epsilon$ は質点の場合1、光子の場合は0を表しています。

オイラー・ラグランジュ方程式から測地線方程式を求めます。 $t$ と $\phi$ は循環座標なので運動の定数（保存量）が存在し、それぞれ $\tilde{E},\ L$ と表すことにします。まずはそれらから考えると良いです。 $r$ はひとまず式(1)を用いて次のように表します。

$\displaystyle{ \begin{aligned} \dot{\tau}&=\frac{\tilde{E}}{c^2(1-r_s/r)} \\ \dot{\phi}&=\frac{L}{r^2\sin^2\theta}\\ \ddot{\theta}&=-\frac{2}{r}\dot{r}\dot{\theta}+\frac{\cos \theta}{r^4\sin^3\theta}L^2\\ \dot{r}^2&=\frac{\tilde{E^2}}{c^2}-r^2\left(1-\frac{r_s}{r}\right)\dot{\theta}^2 -\frac{1}{r^2}\left(1-\frac{r_s}{r}\right)\frac{L^2}{\sin^2\theta}-c^2\left(1-\frac{r_s}{r}\right)\epsilon \end{aligned} }$

ニュートン力学的には角運動量が保存しているので運動面を $\theta=\pi/2$ ( $xy$ 平面)に限定してもよいです。実際、 $\theta$ の微分方程式について、 $\theta=\pi/2$ , $\dot{\theta}=0$ を初期条件として解くと、 $\theta=\pi/2$ (一定)の解を得ます。しかし、シミュレーションでは一つの平面に固定しない3次元的な動きを表現したいのでここでは $\theta$ を固定しないで求めます。

また、固有時間 $\tau$ と観測者にとっての時間 $t$ での時間発展では運動のふるまいが異なることはみなさんご存知だと思います。ここで固有時間 $\tau$ は注目する物体に伴って移動する座標系で計測した時間を表し、 $t$ は遠方で静止している観測者の時間を表しています。簡便のため以下では出来るだけ $\tau$ を粒子の時間、 $t$ を観測者の時間と呼ぶことにします。他にも $t$ を座標時間と読んだりもします。シミュレーションでは両方の立場で運動を表現したいので両方の時間発展方程式を求めます。

光子の計算では固有時間は0になってしまうので、ドット( $\cdot$ )は、 $\tau$ ではなく適当なパラメータ $\lambda$ による微分とします。このパラメータ $\lambda$ を削除して積分すれば光子の軌跡を求めることができるが、時間発展を求めたいので固有時間による光子の運動は考えないことにします。これらの式を使ってルンゲクッタ法により運動を求めることもできますが、プログラムの仕様上 $r,\ \phi$ の2階微分も求めます。

粒子の時間 $\tau$ による時間発展方程式

質点の固有時間の時間発展は

$\displaystyle{ \begin{aligned} \ddot{\phi}&=-\frac{2\dot{r}L}{r^3\sin^2\theta}-\frac{2L\dot{\theta}\cos{\theta}}{r^2\sin^3\theta}\\ \ddot{\theta}&=-\frac{2\dot{r}\dot{\theta}}{r}+\frac{L^2\cos \theta}{r^4\sin^3\theta}\\ \ddot{r}&=\left(r-\frac{3}{2}r_s\right)\dot{\theta}^2-\frac{c^2r_s\epsilon}{2r^2}+\frac{L^2}{r^4\sin^2\theta}\left(r-\frac{3r_s}{2}\right) \end{aligned} }$

初期条件は $\vec{x}_0=(x_0,y_0,z_0)$ 、 $\dot{\vec{x}}_0=(v_{x0},v_{y0},v_{z0})$ を与えればよい。保存量 $L$ は $L=xv_{y0}-yv_{x0}$ から求めます。

観測者の時間 $t$ による時間発展方程式

続いて遠方の静止系による表式を求めます。遠方の静止形での時間 $t$ として、 $r,\ \theta,\ \phi$ の $t$ による2階微分を求めます。次の関係を利用して上記の式を変形していきます。

$\displaystyle{ \begin{aligned} \frac{d}{d\tau}&=\frac{dt}{d\tau}\frac{d}{dt}=\frac{\tilde{E}}{c^2\left(1-r_s/r\right)}\frac{d}{dt}\\ \frac{d^2}{d\tau^2}&=\frac{dt}{d\tau}\frac{d}{dt}\left(\frac{dt}{d\tau}\frac{d}{dt}\right)=\frac{\tilde{E}^2}{c^4\left(1-r_s/r\right)^2}\left[\frac{d^2}{dt^2}-\frac{r_s}{r^2(1-r_s/r)}\frac{dr}{dt}\frac{d}{dt}\right] \end{aligned} }$

$t$ による微分をプライム( $^\prime$ )で表すことにすると以下のようになります。

$\displaystyle{ \begin{aligned} \theta^{\prime\prime}&=\frac{3r_s-2r}{r^2(1-r_s/r)}r^{\prime}\theta ^{\prime}+c^2D^2\frac{(1-r_s/r)^2}{r^4}\frac{\cos \theta}{\sin^3 \theta}\\ \phi^{\prime}&=cD\frac{1-r_s/r}{r^2}\frac{1}{\sin^2\theta}\\ \phi^{\prime\prime}&=-cD\left(\frac{(2r-3r_s)r^{\prime}}{r^4\sin^2\theta}+\frac{2(1-r_s/r)\theta^{\prime}\cos\theta}{r^2\sin^3\theta}\right)\\ {r^{\prime}}^2&=c^2\left(1-\frac{r_s}{r}\right)^2-c^2\hat{E}\left(1-\frac{r_s}{r}\right)^3-c^2D^2\frac{(1-r_s/r)^3}{r^2\sin^2\theta}-r^2\left(1-\frac{r_s}{r}\right){\theta^{\prime}}^2\\ r^{\prime\prime}&=\frac{5r_s-2r}{2r^2(1-r_s/r)}{r^{\prime}}^2-\frac{c^2(3r_s-2r)}{2r^2}\left(1-\frac{r_s}{r}\right)-c^2\hat{E}\frac{(1-r_s/r)^3}{r} \end{aligned} }$

ただし

$\displaystyle{ D=cL/\tilde{E},\ \hat{E}=c^4\epsilon /\tilde{E}^2 %\label{eq:DE} }$

としています。初期条件は書籍による座標の設定では以下の図の通りで、

$\displaystyle{ \begin{aligned} r^{\prime}&=-v_\epsilon\left(1-\frac{r_s}{r}\right)\cos\Psi\\ \theta^{\prime}&=-\frac{v_\epsilon \sqrt{1-r_s/r}}{r}\sin\Psi\cos\Phi\\ \phi^{\prime}&=\frac{v_\epsilon\sqrt{1-r_s/r}}{r\sin\theta}\sin\Psi\sin\Phi \end{aligned} }$

と表されます。運動の定数 $D,\hat{E}$ は初期条件を使って、

$\displaystyle{ \begin{aligned} D&=\frac{r\sin\theta}{\sqrt{1-r_s/r}}\frac{v_\epsilon}{c}\sin\Psi\sin\Phi\\ \hat{E}&=\frac{1-{v_\epsilon}^2/c^2}{1-r_s/r} \end{aligned} }$

と表せます。ただし、光子の場合は $v_\epsilon=c,\ \hat{E}=0$ です。初期条件は $\vec{x}_0=(x_0,y_0,z_0)$ 、 $\vec{v}_\epsilon=(v_{\epsilon x},v_{\epsilon y},v_{\epsilon z})$ を与えればよいです。