弱重力場中の粒子の軌道の近似解：近点移動

テスト粒子の軌道を決める式は，一般には解析的な厳密解を求めることができない。ここでは，粒子の軌道のいたるところで重力場が弱いという近似のもと，粒子（天体，人工衛星等）の軌道を近似的に解く。

念のため，「軌道のいたるところで重力場が弱い」とは，中心天体からの距離 $r$ が，中心天体の質量 $M$ で決まる重力半径 $r_g \equiv \dfrac{2GM}{c^2}$ の十分外側，つまり $ r_g \ll r$ であるような領域を運動しているということ。言い換えると

$$\frac{r_g}{r} \ll 1, \ \ \frac{r_g}{a} \ll 1, \quad a \equiv \frac{r_{\rm min} + r_{\rm max}}{2}$$ ということですよ。

世にあまたある相対論の教科書では，それぞれの著者が趣向をこらして近点移動（太陽のまわりの軌道については近日点移動）を導出しているが，ともすれば初めから近点移動ありきとして，$\cos \gamma \phi$ に比例する解を頭から仮定して導出する例もある。ここでは，なるべくシステマティックに，測地線方程式の積分から得られた微分方程式をうまく変形してやると，自然と近点移動する軌道が解となることがわかるような導出法をまとめておく。

先に答えを書いておく。

ほぼ \(r_g\) （厳密に言えば $\frac{r_g}{a}$）の1次までの範囲で（もう少し細かく言えば，\(O(r_g^2)\) 以上の項および \(O(r_g e^2)\) の項は無視するという近似で）求めたテスト粒子の軌道は以下のように書ける。

$$ r =  \frac{a(1-e^2)}{1 + e \cos(\gamma\phi) } $$

ここで，

$$ \gamma = \sqrt{ 1 -\frac{3 r_g}{a(1-e^2)}} \simeq 1 -\frac{3 r_g}{2a(1-e^2)}$$

ほぼ楕円軌道だが重要な違いは $\gamma$ が $1$ ではないことで，これが近点移動を表す。

また，$a$ および $e$ はニュートン理論では楕円の「軌道長半径」および「離心率」に対応するが，一般相対論的な近似解では，閉じた楕円にならないので，厳密には「軌道長半径」とか「離心率」とかなどと呼ぶことはできない。

しかし，軌道が厳密には解けないとしても，運動が有界であれば $r$ が無限大になったりゼロになったりすることなく，原点のまわりを有限の範囲

$$ r_g < r_{\rm min} \le r \le r_{\rm max}$$

で，有界な束縛運動をするハズである。そこで，

\begin{eqnarray}
r_{\rm max} &\equiv& a\,(1+e) \\
r_{\rm min} &\equiv& a\,(1-e)
\end{eqnarray}

として変数 $a$ および $e$ を定義する。（繰り返すが，閉じた楕円軌道の場合に限り，$a$ および $e$ はまさに「軌道長半径」および「離心率」である。）

また，\(O(r_g e^2)\) の項を無視せずに，$r_g$ の1次の項を全て含んだ場合の解についても追記している。

シュバルツシルト時空中のテスト粒子の軌道を決める式

$\displaystyle s \equiv \frac{1}{r}$ とすると，「シュバルツシルト時空中の粒子（観測者）の運動」のページにまとめたように（適宜移項して）

\begin{eqnarray}
\left( \frac{ds}{d\phi} \right)^2  +s^2 -\frac{2GM}{\ell^2} s  &=& \frac{\epsilon^2 c^2 -c^2}{\ell^2} +r_g\, s^3 \tag{1}
\end{eqnarray}

これが粒子の軌道を決める式であった。

ここで，$\epsilon$ および $\ell$ は測地線方程式を解く際に得られた運動の定数であり，

$$ \frac{d x^0}{d\tau} = \frac{c\, dt}{d\tau} = \frac{\epsilon\, c}{1 -\frac{r_g}{r}}, \quad \frac{d x^3}{d\tau} = \frac{d\phi}{d\tau} = \frac{\ell}{r^2}$$

運動が有界である場合

運動が有界であれば，動径座標 $r$ の値は，ある有限の範囲内にとどまる，つまり $r_{\rm min} \leq r \leq r_{\rm max}$。

\begin{eqnarray}
r_{\rm min} &\equiv& a(1 -e) \\
r_{\rm max} &\equiv& a(1 +e)
\end{eqnarray}

とすると，$s = \frac{1}{r}$ に対しては

$$ \frac{1}{a(1+e)} \equiv \frac{1}{r_{\rm max}} \leq s \leq \frac{1}{r_{\rm min}} \equiv \frac{1}{a(1 -e)}$$

極値をとる点で $\displaystyle \frac{ds}{d\phi} =0$ であるから

\begin{eqnarray}
\left( \frac{1}{a (1+e)}\right)^2 -2 \frac{GM}{\ell^2} \left( \frac{1}{a (1+e)}\right) &=& \frac{\epsilon^2 c^2 -c^2}{\ell^2} +r_g\, \left( \frac{1}{a (1+e)}\right)^3 \\
\left( \frac{1}{a (1-e)}\right)^2 -2 \frac{GM}{\ell^2} \left( \frac{1}{a (1-e)}\right)  &=& \frac{\epsilon^2 c^2 -c^2}{\ell^2} +r_g\, \left( \frac{1}{a (1-e)}\right)^3
\end{eqnarray}

この連立方程式は以下のように解くことができる。（線形近似ではなく，厳密に解いているのであるが，結果は $r_g$ の1次までの形になっているのも興味深い。）

\begin{eqnarray}
\frac{GM}{\ell^2} &=& \frac{1}{a (1 -e^2)} -\frac{(3 + e^2) r_g}{2 a^2 (1 -e^2)^2} \tag{2}\\
\frac{\epsilon^2 c^2 -c^2}{\ell^2} &=& -\frac{1}{a^2 (1 -e^2)} + \frac{2 r_g}{a^3 (1 -e^2)^2} \\
&=& \left(\frac{e}{a (1-e^2)}\right)^2 -\left(\frac{1}{a (1 -e^2)} \right)^2+ \frac{2 r_g}{a^3 (1 -e^2)^2} \tag{3}
\end{eqnarray}

$a, \, e$ は楕円軌道の場合には軌道長半径，離心率と呼ばれるが，ここでは $r_{\rm max}, \, r_{\rm min}$ から以下のようにして決められる定数であることだけを覚えておく。

\begin{eqnarray}
a &\equiv& \frac{1}{2} (r_{\rm max} + r_{\rm min}) \\
e &\equiv& \frac{r_{\rm max} -r_{\rm min}}{r_{\rm max} + r_{\rm min}}
\end{eqnarray}

(2) 式と (3) 式を (1) 式に代入すると，(1) 式の左辺および右辺は

\begin{eqnarray}
\mbox{左辺} &=& \left(\frac{ds}{d\phi} \right)^2 + s^2 -2\left(\frac{1}{a (1 -e^2)} -\frac{(3 + e^2) r_g}{2 a^2 (1 -e^2)^2} \right)\, s \\
&=& \left(\frac{ds}{d\phi} \right)^2 + \left(s -\frac{1}{a (1 -e^2)} \right)^2 -\left(\frac{1}{a (1 -e^2)} \right)^2 \\
&& +  \frac{(3 + e^2) r_g}{a^3 (1 -e^2)^3} \\
&& + \frac{(3 + e^2) r_g}{a^2 (1 -e^2)^2} \left(s -\frac{1}{a (1 -e^2)} \right)\\ \ \\
\mbox{右辺} &=&
\left(\frac{e}{a (1-e^2)}\right)^2 -\left(\frac{1}{a (1 -e^2)} \right)^2 \\
&& + \frac{2 r_g}{a^3 (1 -e^2)^2} \\
&& + r_g \left\{\frac{1}{a (1 -e^2)} + \left(s -\frac{1}{a (1 -e^2)} \right) \right\}^3 \\
&=& \left(\frac{e}{a (1-e^2)}\right)^2 -\left(\frac{1}{a (1 -e^2)} \right)^2 \\
&& + \frac{(3-2 e^2) r_g}{a^3 (1 -e^2)^3} \\
&& + \frac{3 r_g}{a^2 (1 -e^2)^2} \left(s -\frac{1}{a (1 -e^2)} \right) \\
&& + \frac{3 r_g}{a (1 -e^2)} \left(s -\frac{1}{a (1 -e^2)} \right)^2 \\
&& + r_g \left(s -\frac{1}{a (1 -e^2)} \right)^3
\end{eqnarray}

ということで，$\mbox{左辺} = \mbox{右辺}$ として，適宜移項してまとめると

\begin{eqnarray}
\left(\frac{ds}{d\phi} \right)^2 +\gamma^2 \left(s -\frac{1}{a (1 -e^2)} \right)^2 &=&
\gamma^2 \left(\frac{e}{a (1 -e^2)} \right)^2 \\
&& -\frac{e^2\, r_g}{a^2 (1 -e^2)^2} \left(s -\frac{1}{a (1 -e^2)} \right) \\
&& + r_g \left(s -\frac{1}{a (1 -e^2)} \right)^3 \\ \ \\
\gamma^2 &\equiv& \left(1 -\frac{3r_g}{a(1 -e^2)} \right)
\end{eqnarray}

ここまでは厳密な式である。このままでは解析に解けない。特に右辺の第3項（3乗の項）が曲者である。そこで，以下では何らかの方法で近似的に解くことになる。

軌道を決める式を2階微分方程式の形にして解く

このままでも近似的に解くことは可能であるが（別ページ参照），両辺を \(\phi\) で微分して

\begin{eqnarray}
2\frac{ds}{d\phi} \frac{d^2s}{d\phi^2}  + 2\gamma^2 \left( s – \frac{1}{a (1 -e^2)} \right) \frac{ds}{d\phi}
&=& -\frac{e^2\, r_g}{a^2 (1 -e^2)^2} \frac{ds}{d\phi} + 3 r_g \left( s – \frac{1}{a (1 -e^2)} \right)^2 \frac{ds}{d\phi}\\
\therefore\ \  \frac{d^2s}{d\phi^2} +\gamma^2 \left( s – \frac{1}{a (1 -e^2)} \right)
&=& -\frac{{\color{red}{e^2\, r_g}}}{2 a^2 (1 -e^2)^2} + \frac{3}{2}r_g \left( s – \frac{1}{a (1 -e^2)} \right)^2 \tag{4}
\end{eqnarray}

として，こちらを近似的に解く方法を紹介する。なんでもう1階微分して，この2階微分方程式にしたかというと，みてわかるように，全体を見わたすと懐かしい単振動の方程式に補正項がついた形になっていて，ひょっとしたら解きやすいかも… と思われるからである。

積分定数を決める初期条件

ただし，このままだと本来1階常微分方程式だから解は積分定数を1個持つはずが，もう1階微分して2階微分方程式となったので解は積分定数を2個持つことになってしまう。

微分の階数を人為的に上げたことで現れることになってしまった余分の積分定数も決めるための初期条件として，以下を採用する。

• $\phi = 0$ で $r$ は最小値 $r_{\rm min} = a(1-e)$ をとる。すなわち
  1. $\phi = 0$ で $\displaystyle \frac{ds}{d\phi} = 0$ および
  2. $\phi = 0$ で$\displaystyle s = \frac{1}{r_{\rm min}} = \frac{1}{a(1-e)}$

1階常微分方程式のままで解く場合は，積分定数1個を決める初期条件はひとつのみでよい。

$r_g$ のゼロ次解

天体の軌道は重力半径 $r_g$ の十分外側であるという状況では， \(\displaystyle 0 < \frac{r_g}{r} = r_g s \ll 1\)，あるいは同じことだが $\displaystyle \frac{r_g}{a} \ll 1$ としてよい。上の (4) 式で $r_g$ がかかっている項を無視した場合の解を $s_0$ と書くと，この近似では $\gamma = 1$ としてよいので，

\begin{eqnarray}
\frac{d^2 s_0}{d\phi^2}  +\left(s_0 -\frac{1}{a (1 -e^2)}  \right) &=& 0
\end{eqnarray}

これは大学に入って最初に習う微分方程式の一つであり，一般解は脊髄反射で以下のように解ける。
$$ s_0 – \frac{1}{a(1-e^2)} =  A \cos\phi + B \sin \phi$$
積分定数 \(A, B\) は以下のようにして決定する。

• \(\phi = 0 \) で \( \displaystyle \frac{ds_0}{d\phi} = 0 \) より，\(B = 0 \)。
• $\phi = 0$ で $\displaystyle s_0 = \frac{1}{a(1-e)}$ より，
  \(\displaystyle A = \frac{e}{a(1-e^2)} \) と求められる。

まとめ：\(r_g\) の項を無視したときの解を \(r_g\) のゼロ次の解ということで以下のように \(s_0\) と書く。
$$s_0 =\frac{1}{r} =  \frac{1 + e\cos\phi}{a(1-e^2)} $$

つまり，

$$r = \frac{a(1-e^2)}{1 + e\cos\phi}$$

ということで，確かに楕円になっている。

\(O(e^2\,r_g)\) の項を無視した場合の解

(4) 式右辺第2項に $s_0$ を代入して評価してやると

\begin{eqnarray}
\frac{3}{2} r_g \left(s -\frac{1}{a (1 -e^2)}\right)^2 &\simeq& \frac{3}{2} r_g \left(s_0 -\frac{1}{a (1 -e^2)}\right)^2 \\
&=& \frac{3}{2} r_g \left(\frac{e\cos\phi}{a(1-e^2)}\right)^2 \\
&=& \frac{3\, {\color{red}{e^2\, r_g}}\cos\phi}{2 a^2 (1-e^2)^2}
\end{eqnarray}

となり，(4) 式右辺第1項とおなじ $O({\color{red}{e^2\, r_g}})$ の項になる。楕円軌道の場合に離心率と呼ばれる $e$ は一般に $0 \leq e < 1$ であるから $O({\color{red}{e^2\, r_g}})$ の項を無視すると右辺はゼロとなり，

\begin{eqnarray}
\frac{d^2s}{d\phi^2} +\gamma^2
\left(s -\frac{1}{a (1 -e^2)} \right) &=& 0
\end{eqnarray}

これも大学に入って最初に習う微分方程式の一つであり，一般解は脊髄反射で以下のように解ける。

$$ s -\frac{1}{a(1-e^2)} = A \cos\gamma\phi + B \sin \gamma \phi$$

同じ初期条件を使うと，積分定数 $A, B$ はただちに求まり，

$$B = 0, \quad A = \frac{e}{a(1-e^2)}$$

最終的に

\begin{eqnarray}
\frac{1}{s} = r &=& \frac{a(1-e^2)}{1 + e\cos (\gamma \phi)} \\
\gamma &\equiv& \sqrt{1 -\frac{3 r_g}{a(1-e^2)}} \simeq 1 -\frac{3 r_g}{2a(1-e^2)}
\end{eqnarray}

$\phi = 0$ で $r = r_{\rm min}$ の最小値からスタートした軌道が再び $r = r_{\rm min}$ になるのは，$\phi = 2\pi$ ではなく，

$$\gamma \phi = 2 \pi, \quad \therefore\ \ \phi = \frac{2\pi}{\gamma} > 2\pi$$

のときである。したがって近点は移動する。

… とまとめてはみたものの，

$$\frac{3 r_g}{2a(1-e^2)} -\frac{3 r_g}{2a} = \frac{3e^2  r_g}{2a(1-e^2)} = O({\color{red}{e^2\, r_g}})$$

であるので，$O({\color{red}{e^2\, r_g}})$ の項は無視するという今回の近似の方針に従えば，

$$\gamma  \simeq 1 -\frac{3 r_g}{2a(1-e^2)}$$

としてもよいし，

$$\gamma  \simeq 1 -\frac{3 r_g}{2a}$$

としてもよいことになるなぁ…  ということで，以下では \(O({\color{red}{e^2\, r_g}})\) の項を無視せずに，$r_g$ の1次の項を全て含んだ式を解いてみる。

\(O(e^2\,r_g)\) の項を無視しない場合の式

(4) 式右辺第2項を \(O(r_g e^2)\) の項を無視した場合の解を代入して評価してやると，$r_g$ の1次の項を全て含んだ式は以下のようになる。

\begin{eqnarray}
\frac{d^2s}{d\phi^2} +\gamma^2
\left(s -\frac{1}{a (1 -e^2)} \right)
&\simeq&
\frac{\color{red}{r_g e^2}}{2 a^2 (1-e^2)^2} \left( 3\cos^2(\gamma \phi) -1\right)
\end{eqnarray}

\(O(e^2\,r_g)\) の項を無視しない場合の解

${\color{red}{r_g e^2}}$ に比例する項は，2階常微分方程式の非同次項であるから，\(O({\color{red}{r_g e^2}}) \) の項を無視して求めた同次方程式の解に，非同次方程式の特殊解を加えたものが解となり，以下のように書ける。

\begin{eqnarray}
s= \frac{1}{r} &=& \frac{1 + e\cos (\gamma \phi)} {a(1-e^2)}
+ \frac{{\color{red}{r_g e^2}}}{2 a^2 (1-e^2)^2 } \frac{\sin^2(\gamma \phi)}{\gamma^2} \\
&\simeq& \frac{1 + e\cos (\gamma \phi)}{a(1-e^2)}
\left\{ 1 + \frac{{\color{red}{r_g e^2}}}{2 a (1-e^2)} \frac{\sin^2 (\gamma \phi)}{1+e\cos(\gamma\phi)}\right\}
\end{eqnarray}

非同次方程式の特殊解はロンスキアンを使った公式で求めることができる。授業でやりました。以下のページを参照：

• • 定数係数2階線形非同次方程式

いちいち手計算するのが面倒な人のために，Maxima や Python の SymPy で解く例も以下に示しておきます。

• • 参考：Maxima で非同次2階線形常微分方程式を解く
  • 参考：SymPy で非同次2階線形常微分方程式を解く

最終的に $r$ は

$$
r = \frac{a(1-e^2)} {1 + e\cos (\gamma \phi)} \left\{ 1 -\frac{{\color{red}{r_g e^2}}}{2 a (1-e^2)} \frac{\sin^2 (\gamma \phi)}{1+e\cos(\gamma\phi)}\right\}
$$

すでに上で定義されているが，あらためて $\gamma$ も書いておくと，

$$\gamma \equiv \sqrt{1 -\frac{3 r_g}{a(1-e^2)}} \simeq 1 -\frac{3 r_g}{2a(1-e^2)}$$

この場合であっても，$\phi = 0$ で $r = r_{\rm min}$ の最小値からスタートした軌道が再び $r = r_{\rm min}$ になるのは，$\phi = 2\pi$ ではなく，

$$\gamma \phi = 2 \pi, \quad \therefore\ \ \phi = \frac{2\pi}{\gamma} > 2\pi$$

のときである。したがって近点は移動することに変わりはない。

assume(gamma > 0)$

eq: 'diff(X, phi, 2) = - gamma**2 * X + A*(3*cos(gamma*phi)**2-1);

/* 微分方程式を解く */
ode2(eq, X, phi)$

/* 初期条件を課す */
ic2(%, phi=0, X=e/(a*(1-e**2)), 'diff(X, phi)=0);

from sympy.abc import *
from sympy import * 

X = Function('X')
var('phi')
var('gamma', positive = True)

eq = Eq(Derivative(X(phi), phi, 2), 
        - gamma**2 * X(phi) + A*(3*cos(gamma*phi)**2 - 1))
eq

dsolve(eq, X(phi), 
       ics = {X(phi).subs(phi, 0):e/(a*(1-e**2)), 
              diff(X(phi), phi).subs(phi, 0): 0}
      )