相対論的レイトレーシング

これは、以前ロシア語で公開した記事の、わずかに洗練されたバージョンです。この記事では、一般相対性理論に完全に準拠した、物理的に正確なブラックホールをレンダリングする方法を示します。そのために、光線の運動方程式を導出し、GLSLでルンゲ＝クッタ積分器を実装し、最終的にそれらを組み合わせて、カメラから過去へ光線を追跡するフラグメントシェーダーを作成します。 準備：計量、接続、測地線 光は直線を進みます。しかし、本当に「直線」とは何を意味するのでしょうか？線形方程式で記述されるもの？いいえ、そうではありません。極座標系を取ると、直線の方程式は線形とは程遠いものになります。そして、空間自体が平坦でない場合（このトーラスのように）、直線が線形になるような座標系の選択肢はありません。 ある点での曲線に沿ったベクトルを取り、その曲線をたどって移動させるとします。それはまだ曲線に沿った方向を指すでしょうか？もし任意の点でそうであれば、そのような曲線は測地線と呼ばれ、これは任意の多様体上で直線に最も近いものです。 4次元多様体上の測地線の方程式を導出するのは非常に骨の折れる作業です。幸いなことに、手作業で行う必要はありません。機械に骨の折れる作業をやらせましょう。PythonパッケージであるSymPyを使用します。 ```python from sympy import * from sympy.diffgeom import * import numpy as np ``` sympy.diffgeomモジュールはやや不便で扱いにくいですが、アトラスで記述された多様体上の偏微分方程式を導出するという仕事は reasonably well 行います。私たちの場合はアトラスは1つのパッチしか持ちませんが、いくつかの座標系を持ちます。一つは球座標系（シュワルツシルト計量を扱う際に便利）、もう一つは our future camera 用のデカルト座標系です。座標シンボルを定義しましょう。 ```python t_, x_, y_, z_ = symbols("t x y z") rho_, theta_, phi_ = symbols("\\rho \\theta \\phi") ``` （私は末尾に`_`を付けるのを、終端シンボルのための慣習として使用しています） 球座標系からデカルト座標系への遷移の方程式を作成しましょう。 ```python x = rho_ * sin(theta_) * sin(phi_) y = rho_ * sin(theta_) * cos(phi_) z = rho_ * cos(theta_) print(latex(x_) + " = " + latex(x)) print(latex(y_) + " = " + latex(y)) print(latex(z_) + " = " + latex(z)) ``` $ \begin{array}l x = \rho \sin{\left(\phi \right)} \sin{\left(\theta \right)} \\ y = \rho \sin{\left(\theta \right)} \cos{\left(\phi \right)} \\ z = \rho \cos{\left(\theta \right)} \\ \end{array} $ 新しい4次元多様体`spacetime`を作成し、そのアトラスに`patch`を追加し、2つの座標系をその`patch`にアタッチします。 ```python spacetime = Manifold("spacetime", 4) patch = Patch("patch", spacetime) relation_dict = { ("spherical", "cart"): [(t_, rho_, theta_, phi_), (x, y, z, t_)] } coord_sh = CoordSystem("spherical", patch, (t_, rho_, theta_, phi_), relation_dict) coord_cart = CoordSystem("cart", patch, (x_, y_, z_, t_), relation_dict) ``` symPy.diffgeomの不便な点の一つは、`CoordSystem`関数に渡された座標シンボルをそのまま使用できないことです。代わりに、そこから新しいシンボルを取得する必要があります。 ```python c_t, c_rho, c_theta, c_phi = coord_sh.base_scalars() print(",".join(map(latex, (c_t, c_rho, c_theta, c_phi)))) ``` $ \mathbf{t},\mathbf{\rho},\mathbf{\theta},\mathbf{\phi} $ ご覧の通り、シンボルは同じですが、新しいセットの変数のみがsympy.diffgeomで使用できます。 次に、シュワルツシルト計量を定義しましょう。 ```python Rs = symbols("r_s") d_t, d_rho, d_theta, d_phi = coord_sh.base_oneforms() TP = TensorProduct metric = ( + (1 - Rs/c_rho) * TP(d_t, d_t) - 1 / (1 - Rs / c_rho) * TP(d_rho, d_rho) - (c_rho**2) * TP(d_theta, d_theta) - (c_rho**2) * (sin(c_theta)**2) * TP(d_phi, d_phi) ) print(latex(metric)) ``` $ \left(- \frac{r_{s}}{\mathbf{\rho}} + 1\right) \operatorname{d}t \otimes \operatorname{d}t - \sin^{2}{\left(\mathbf{\theta} \right)} \mathbf{\rho}^{2} \operatorname{d}\phi \otimes \operatorname{d}\phi - \mathbf{\rho}^{2} \operatorname{d}\theta \otimes \operatorname{d}\theta - \frac{\operatorname{d}\rho \otimes \operatorname{d}\rho}{- \frac{r_{s}}{\mathbf{\rho}} + 1} $ ここで、テンソル積の1-形式を使用することを強制されるのは、もう一つの不便な点です。$g_{ij}$成分の2次元配列があればより簡単だと考えるかもしれませんが、残念ながらそうではありません。微分形式で進めます。いずれにせよ、計量をこの形式で得られれば、[レヴィ・チヴィタ接続](https://en.wikipedia.org/wiki/Levi-Civita_connection)（別名、2種クリストッフェル記号）の成分を計算するために使用できます。多様体上の接続は、ある点の接空間のベクトルが、異なる点の接空間に連続的に移動する際にどのように変換されるかを定義します。そして、レヴィ・チヴィタ接続は、そのような移動の下で計量テンソルを保存する特別な種類の接続です。これは、計量テンソルの共変微分がゼロであることを意味します（その成分は変化するかもしれませんが）。 symPy.diffgeomは、その計量テンソル保存条件に基づいて、クリストッフェル記号を導出する関数`metric_to_Christoffel_2nd`を提供しています。 ```python christoffel = simplify(metric_to_Christoffel_2nd(metric)) print("\\begin{cases}") # v v v intentional reordering here, to sort by lower indices first for (j,k,i) in np.ndindex(*np.shape(christoffel)): if not christoffel[i,j,k].is_zero: print(f"\\Gamma^{{{i}}}_{{{j}{k}}} =" + latex(christoffel[i,j,k]) + " \\") print("\\end{cases}") ``` $ \begin{cases} \Gamma^{1}_{00} =\frac{r_{s} \left(- r_{s} + \mathbf{\rho}\right)}{2 \mathbf{\rho}^{3}} \\ \Gamma^{0}_{01} =\frac{r_{s}}{2 \left(- r_{s} + \mathbf{\rho}\right) \mathbf{\rho}} \\ \Gamma^{0}_{10} =\frac{r_{s}}{2 \left(- r_{s} + \mathbf{\rho}\right) \mathbf{\rho}} \\ \Gamma^{1}_{11} =\frac{r_{s}}{2 \left(r_{s} - \mathbf{\rho}\right) \mathbf{\rho}} \\ \Gamma^{2}_{12} =\frac{1}{\mathbf{\rho}} \\ \Gamma^{3}_{13} =\frac{1}{\mathbf{\rho}} \\ \Gamma^{2}_{21} =\frac{1}{\mathbf{\rho}} \\ \Gamma^{1}_{22} =r_{s} - \mathbf{\rho} \\ \Gamma^{3}_{23} =\frac{1}{\tan{\left(\mathbf{\theta} \right)}} \\ \Gamma^{3}_{31} =\frac{1}{\mathbf{\rho}} \\ \Gamma^{3}_{32} =\frac{1}{\tan{\left(\mathbf{\theta} \right)}} \\ \Gamma^{1}_{33} =\left(r_{s} - \mathbf{\rho}\right) \sin^{2}{\left(\mathbf{\theta} \right)} \\ \Gamma^{2}_{33} =- \frac{\sin{\left(2 \mathbf{\theta} \right)}}{2} \\ \end{cases} $ 接続成分が手に入ったので、測地線方程式を導出できます。そして、それを手動で行う必要があります。幸いなことに、それほど難しくはありません。ある曲線がパラメータ方程式$x^\mu(t)$で定義されているとします。その場合、その接線ベクトルは$\\dot{x}^\mu(t) = \frac{d}{dt}x^\mu(t)$の形になります。測地線の定義によれば、このベクトルを自身に沿った共変微分はゼロになるはずです：$ \nabla_{\\dot{x}} \\dot{x} = 0 $ これをクリストッフェル記号で書き出すと： $ \nabla_{\\dot{x}} \\dot{x} = \left( \frac{\partial}{\partial{x^\mu}} \\dot{x}^\nu + \Gamma^{\nu}_{\kappa\mu}\\dot{x}^\kappa \right) \\dot{x}^\mu = \left( \frac{\partial}{\partial{x^\mu}} \frac{dx^\nu}{dt} \right) \frac{dx^\mu}{dt} + \Gamma^{\nu}_{\kappa\mu} \frac{dx^\kappa}{dt} \frac{dx^\mu}{dt} $ 最初の項で、因子の順序を入れ替え、連鎖律を使用しましょう。 $ \left( \frac{\partial}{\partial{x^\mu}} \frac{dx^\nu}{dt} \right) \frac{dx^\mu}{dt} = \frac{dx^\mu}{dt}\\frac{\partial}{\partial{x^\mu}} \left( \frac{dx^\nu}{dt} \right) = \frac{d}{dt}\left( \frac{dx^\nu}{dt} \right) = \frac{d^2 x^\nu}{dt^2} $ すると、方程式は次のようになります。 $ \frac{d^2 x^\nu}{dt^2} + \Gamma^{\nu}_{\kappa\mu} \frac{dx^\kappa}{dt} \frac{dx^\mu}{dt} = 0 $ または、単純に： $ \boxed{ \ddot{x}^\nu = -\Gamma^{\nu}_{