ニュートン力学における万有引力の2体問題についておさらいしておく。
このサイトの他のセクションでは，
- 4元ベクトルを太字で \(\boldsymbol{u}\)，その成分を \(u^{\mu} = (u^0, u^1, u^2, u^3)\)
- 3次元ベクトルを矢印で \(\vec{v} = (v_x, v_y, v_z)\)
などと区別して表記しているが，この節については，4元ベクトルは現れず，3次元ベクトルのみなので，
- ここでは3次元ベクトルを太字で \(\boldsymbol{v} = (v_x, v_y, v_z)\) と表すことにする。
まずは万有引力の1体問題の運動方程式
万有引力の1体問題とは，万有引力を及ぼす重力源である質量 $M$ の天体が原点に静止しており，そのまわりを質量 $m$ のテスト天体（質量が十分小さいため，重力源の天体に力を及ぼして動かすことはないと仮定できる天体のこと）が運動している，とする状況設定の問題のこと。
この場合，原点からテスト天体までの位置ベクトルを $\boldsymbol{r}$ とすると，運動方程式は
$$m \frac{d^2 \boldsymbol{r}}{dt^2} = – \frac{G M m}{r^3} \boldsymbol{r}$$
あるいは両辺の $m$ をキャンセルさせて
$$\frac{d^2 \boldsymbol{r}}{dt^2} = – \frac{G M}{r^3} \boldsymbol{r}$$
2体問題の運動方程式
- 天体1（太陽を想定）：質量 \(m_1\)，位置ベクトル \(\boldsymbol{r}_1\)
- 天体2（惑星を想定）：質量 \(m_2\)，位置ベクトル \(\boldsymbol{r}_2\)
の2つの天体が互いに万有引力を及ぼし合って運動する。運動方程式は，
$$m_1 \frac{d^2\boldsymbol{r}_1}{dt^2} = – \frac{Gm_1 m_2 (\boldsymbol{r}_1 – \boldsymbol{r}_2)}{|\boldsymbol{r}_1 – \boldsymbol{r}_2|^3} \tag{1}$$
$$m_2 \frac{d^2\boldsymbol{r}_2}{dt^2} = – \frac{Gm_2 m_1 (\boldsymbol{r}_2 – \boldsymbol{r}_1)}{|\boldsymbol{r}_2 – \boldsymbol{r}_1|^3} \tag{2}$$
\(\mbox{(1)}+\mbox{(2)}\) より
$$\frac{d^2}{dt^2} \left( m_1 \boldsymbol{r}_1 + m_2 \boldsymbol{r}_2 \right) = \boldsymbol{0}$$
質量中心の位置ベクトル \(\boldsymbol{R}\) を
$$\boldsymbol{R} \equiv \frac{m_1 \boldsymbol{r}_1 + m_2 \boldsymbol{r}_2}{M},
\quad M \equiv m_1 + m_2$$
と定義すると，上式から
$$\frac{d^2\boldsymbol{R}}{dt^2} = \boldsymbol{0}$$
これから，一般性を失うことなく，質量中心は静止しているとして話を進める。
\(\displaystyle \frac{1}{m_2}\times\mbox{(2)}-\frac{1}{m_1}\times\mbox{(1)}\) より，相対位置ベクトル \(\boldsymbol{r} \equiv \boldsymbol{r}_2 – \boldsymbol{r}_1\) に対して
$$\frac{d^2\boldsymbol{r}}{dt^2} = – \frac{GM \boldsymbol{r}}{r^3}
\tag{3}$$
2つの天体の位置ベクトル \(\boldsymbol{r}_1, \boldsymbol{r}_2\) に関する連立微分方程式だった2体問題の運動方程式が \(\boldsymbol{r}\) に対する1組の運動方程式になることを，業界用語で
2体問題は1体問題に帰着する
という。
保存量（運動の定数）
角運動量保存
まず，\(\mbox{(3)}\) 式に \(\boldsymbol{r}\) を外積してやると
$$\boldsymbol{r}\times \frac{d^2\boldsymbol{r}}{dt^2} = -\frac{GM}{r^3} \boldsymbol{r}\times \boldsymbol{r} = \boldsymbol{0}$$
$$\therefore\ \boldsymbol{r}\times \frac{d^2\boldsymbol{r}}{dt^2} =
\frac{d}{dt}\left(\boldsymbol{r}\times \frac{d\boldsymbol{r}}{dt} \right) = \boldsymbol{0}$$
$$\therefore\ \boldsymbol{r}\times \frac{d\boldsymbol{r}}{dt} = \mbox{const.} \equiv \boldsymbol{\ell}$$
定ベクトル \(\boldsymbol{\ell}\) は単位質量あたりの角運動量に相当するベクトルであり，一定であるその向きを \(z\) 軸にとることができる。
$$\boldsymbol{\ell} = (0, 0, \ell)$$
さらに，
$$\boldsymbol{r}\cdot\boldsymbol{\ell} =
\boldsymbol{r}\cdot\left(\boldsymbol{r}\times \frac{d\boldsymbol{r}}{dt} \right) =
\frac{d\boldsymbol{r}}{dt}\cdot(\boldsymbol{r}\times\boldsymbol{r}) = 0,
\quad\therefore\ \boldsymbol{r} \perp \boldsymbol{\ell}$$
となり，一般性を失うことなく \(\boldsymbol{r}\) を\(xy\) 平面（赤道面）上にとることができて
$$\boldsymbol{r} = (x, y, 0) = (r\cos\phi, r\sin\phi, 0)$$
\(\ell\) を極座標であらわすと
\begin{eqnarray}
\ell = |\boldsymbol{\ell}| =x\frac{dy}{dt} -y\frac{dx}{dt} = r^2 \frac{d\phi}{dt},
\quad \therefore\ \frac{d\phi}{dt} = \frac{\ell}{r^2} \tag{A}
\end{eqnarray}
エネルギー保存
もう一つの保存量は，\(\mbox{(3)}\) 式に \(\displaystyle\frac{d\boldsymbol{r}}{dt}\) を内積してやると得られる。
\begin{eqnarray}
\frac{d\boldsymbol{r}}{dt}\cdot\frac{d^2\boldsymbol{r}}{dt^2} &=&
– \frac{GM}{r^3} \boldsymbol{r}\cdot\frac{d\boldsymbol{r}}{dt}\\
\frac{d}{dt}\left(\frac{1}{2} \frac{d\boldsymbol{r}}{dt}\cdot\frac{d\boldsymbol{r}}{dt}\right)&=&
– \frac{GM}{r^3} \frac{d}{dt}\left(\frac{1}{2} \boldsymbol{r}\cdot\boldsymbol{r}\right)\\
&=&- \frac{GM}{r^3} \frac{d}{dt}\left(\frac{1}{2} r^2\right) \\
&=& \frac{d}{dt} \left(\frac{GM}{r}\right)\\
\end{eqnarray}
$$\therefore\ \ \frac{d}{dt} \left(\frac{1}{2} \frac{d\boldsymbol{r}}{dt}\cdot\frac{d\boldsymbol{r}}{dt}-\frac{GM}{r} \right) = 0$$
$$\therefore\ \ \frac{1}{2} \frac{d\boldsymbol{r}}{dt}\cdot\frac{d\boldsymbol{r}}{dt}-\frac{GM}{r} =
\mbox{const.} \equiv \varepsilon$$
\(\varepsilon\) は単位質量あたりの力学的エネルギー（全エネルギー）に相当する定数であり，極座標であらわすと
\begin{eqnarray}
\varepsilon &=& \frac{1}{2} \left\{\left(\frac{dx}{dt}\right)^2 + \left(\frac{dy}{dt}\right)^2 \right\} – \frac{GM}{r} \\
&=& \frac{1}{2}\left\{\left(\frac{dr}{dt}\right)^2 +r^2\left(\frac{d\phi}{dt}\right)^2 \right\} – \frac{GM}{r}\\
&=& \frac{1}{2}\left\{\left(\frac{dr}{dt}\right)^2 +\frac{\ell^2}{r^2} \right\} – \frac{GM}{r} \\
\therefore\ \ \left(\frac{dr}{dt}\right)^2 &=& 2 \varepsilon + \frac{2 G M}{r} – \frac{\ell^2}{r^2} \tag{B}
\end{eqnarray}
運動が有界な場合
運動が有界な束縛軌道である場合，
$$ 0 < r_{\rm min} \leq r \leq r_{\rm max}$$
となる。$r = r_{\rm min}$ および $r = r_{\rm max}$ では $r$ が極値をとることから $\displaystyle \frac{dr}{dt} = 0$ となるから (B) 式から
\begin{eqnarray}
0 &=& 2 \varepsilon + \frac{2 G M}{r_{\rm min}} – \frac{\ell^2}{r_{\rm min}^2} \\
0 &=& 2 \varepsilon + \frac{2 G M}{r_{\rm max}} – \frac{\ell^2}{r_{\rm max}^2}
\end{eqnarray}
この連立方程式を $\varepsilon, \ell^2$ について解くと
\begin{eqnarray}
\varepsilon &=& -\frac{G M}{r_{\rm max} + r_{\rm min}} \\
\ell^2 &=& \frac{2 G M \,r_{\rm max}\, r_{\rm min}}{r_{\rm max} + r_{\rm min}}
\end{eqnarray}
これらを
$$r_{\rm max} \equiv a (1 + e), \quad r_{\rm min} \equiv a (1 -e)$$
すなわち
\begin{eqnarray}
a &\equiv& \frac{r_{\rm max} + r_{\rm min}}{2} \\
e &\equiv& \frac{r_{\rm max} -r_{\rm min}}{r_{\rm max} + r_{\rm min}}
\end{eqnarray}
で定義される $a, \, e$ を使って書き表すと
\begin{eqnarray}
\varepsilon &=& -\frac{G M}{2 a} \\
\ell^2 &=& G M a ( 1 -e^2)
\end{eqnarray}
となる。現段階では $a$ は $r_{\rm max} $ と $r_{\rm min} $ の平均値，$e$ はそれらの差からつくられる無次元量という意味合いしかないが，後々にわかるように，軌道が楕円であることが解けたあかつきには，$a$ は軌道長半径，$e$ は離心率と呼ばれるようになる。
軌道を決める式
$r$ を $\phi$ を通して $t$ の関数であるとすると
$$ \frac{dr}{dt} = \frac{dr}{d\phi}\,\frac{d\phi}{dt} = \frac{\ell}{r^2} \frac{dr}{d\phi}$$
であるから (B) 式は
\begin{eqnarray}
\left(\frac{\ell}{r^2} \frac{dr}{d\phi}\right)^2 &=& 2 \varepsilon + \frac{2 G M}{r} – \frac{\ell^2}{r^2} \\
\therefore\ \ \left(\frac{1}{r^2} \frac{dr}{d\phi}\right)^2 &=& \frac{2 \varepsilon}{\ell^2} + \frac{2 G M}{\ell^2 }\frac{1}{r} -\frac{1}{r^2} \\
&=& -\frac{1}{a^2 (1 -e^2)} + \frac{2}{a (1 -e^2)} \frac{1}{r} -\frac{1}{r^2}
\end{eqnarray}
\(r\) がことごとく分母に現れている状況を鑑み，\(\displaystyle s \equiv \frac{1}{r}\) という変数で書き直すと
\begin{eqnarray}
\left(\frac{ds}{d\phi}\right)^2 &=& -\frac{1}{a^2 (1 -e^2)} + \frac{2}{a (1 -e^2)} s -s^2 \\
&=& -\frac{1}{a^2 (1 -e^2)} -\left(s -\frac{1}{a (1 -e^2)} \right)^2 + \frac{1}{a^2 (1 -e^2)^2} \\
&=& \frac{1}{b^2} -\left(s -\frac{1}{a (1 -e^2)} \right)^2
\end{eqnarray}
ここで
$$\frac{1}{b} \equiv \frac{e}{a (1 -e^2)} $$
シュバルツシルト時空における粒子の軌道を決める式との類似性に刮目せよ。
解：楕円軌道
軌道を決める式は，あらためて変数を
$$u_0 \equiv s -\frac{1}{a (1 -e^2)}$$
と置き直すと，
$$\left(\frac{du_0}{d\phi}\right)^2 = \frac{1}{b^2} – u_0^2$$
この式は変数分離形にできて
$$\pm \frac{d(b u_0)}{\sqrt{1 – (b u_0)^2}} = d\phi$$
「補足」に書いたように，初期条件として \(\phi = 0\) のときに \(u_0\) したがって \(r\) が極値をとるとすれば，ただちに解は以下のように求まる。
\begin{eqnarray}
u_0 &=& s -\frac{1}{a (1 -e^2)} = \frac{\cos\phi}{b} = \frac{e \cos\phi}{a (1 -e^2)} \\
\therefore\ \ s &=& \frac{1}{r} =\frac{1 + e \cos\phi}{a (1 -e^2)}
\end{eqnarray}
$$\therefore \ r = \frac{a(1-e^2)}{1 + e\cos\phi}$$
となり，この軌道は，長半径 \(a\)，離心率 \(e\) の楕円であることが一目瞭然となる。
\(a\)，\(e\) を保存量（運動の定数）を使って書き直してやると，束縛軌道の場合に \( \varepsilon = -|\varepsilon|\) となることを使うと
$$a = \frac{GM}{2|\varepsilon|}, \quad e = \sqrt{1 – \frac{2 |\varepsilon| \ell^2}{(GM)^2}}$$
楕円軌道であることはわかったものの…
ということで，万有引力の2体問題は，1体問題に帰着して（束縛軌道の場合には）楕円軌道が解となることが解析的に示されたわけだが，我々が素朴に思うところの「解」とは少し趣が違う。
たとえば，一様重力場中の斜方投射では以下のように \(x\) 座標と \(y\) 座標が時間 \(t\) の陽関数としてあからさまに解けた。
\begin{eqnarray}
x(t) &=& x_0 + v_{0x} t \\
y(t) &=& y_0 + v_{0y} t – \frac{1}{2} g t^2
\end{eqnarray}
しかし，万有引力の2体問題の解である楕円軌道は，\(r\) が \(\phi\) の関数として
$$r = \frac{a(1-e^2)}{1 + e\cos\phi}$$
のように与えられているだけで，時間 \(t\) の陽関数としてあからさまに解けているわけではない。したがって，ある時刻 \(t\) のときの天体の位置 \( (x(t), y(t)) \) あるいは \( (r(t), \phi(t)) \) を知りたいという場合には別の工夫が必要になる。このあたりについては以下の Memo を参照。