古典力学における楕円運動

古典力学における楕円運動

　3次元空間における楕円運動を考える。典型的には太陽のような大きな恒星の周りを周回する惑星の運動が該当するだろう。古典力学において楕円運動のエネルギー（ハミルトニアン）は

\[ \small E = H(p,q) = \frac{p^2}{2m} – \frac{GmM}{r} \]

で与えられる。一方の物質（恒星）がもう一方の物質（惑星）と比較して非常に大きな質量があるものと仮定し、その物質（恒星）は原点から動かないものと仮定する。現実には、一方の質量が非常に大きいという仮定は特に必要なく、2つの物質の相対的な運動しか考えない問題では一般性を失うことなく片方の座標を原点と仮定してよいだろう。このときのもう一方の物質（惑星）の運動方程式を求めよう。相対性理論では、シュバルツシルトの解と言われる問題に相当するだろう。計算が長くなるので、相対性理論における楕円運動は別の記事で取り扱うことにする。

　単純化のため、楕円運動の回転方向は\(\small xy\)平面とし、\(\small z=0\)と仮定する。加速度を計算すれば

\[ \small \begin{align*} &-\frac{\partial H}{\partial x}=F_x = m\frac{d^2x}{dt^2} = -\frac{GmM}{r^2}\frac{x}{r} \\ &-\frac{\partial H}{\partial y}=F_y = m\frac{d^2y}{dt^2} = -\frac{GmM}{r^2}\frac{y}{r} \\ &-\frac{\partial H}{\partial z}=F_z = m\frac{d^2z}{dt^2} = -\frac{GmM}{r^2}\frac{z}{r}=0 \end{align*} \]

となる。座標と力の外積を計算すると

\[ \small q \times F = \left[ \begin{array}{c} yF_z-zF_y \\ zF_x-xF_z \\ xF_y-yF_x \end{array} \right] = \left[ \begin{array}{c} 0 \\ 0 \\ 0 \end{array} \right] \]

を得る。これは角運動量\(\small q\times p\)​の時間微分に相当するが、それが0であるということであり、角運動量が定数（保存量）であることが分かる。力と角運動量の外積を計算すると

\[ \small \begin{align*} F \times l \; & = \left[ \begin{array}{c} F_yl_z -F_z l_y \\ F_zl_x -F_x l_z \\ F_xl_y -F_y l_x \end{array} \right] = -\frac{GmM}{r^3}\left[ \begin{array}{c} yl_z-zl_y \\ zl_x-x l_z \\ xl_y-y l_x \end{array} \right] \\ & = -\frac{GmM}{r^3}\left[ \begin{array}{c} x(xp_x+yp_y+zp_z)-(x^2+y^2+z^2)p_x \\ y(xp_x+yp_y+zp_z)-(x^2+y^2+z^2)p_y \\ z(xp_x+yp_y+zp_z)-(x^2+y^2+z^2)p_z \end{array} \right] \end{align*} \]

であり、右辺は

\[ \small Gm^2M\frac{d}{dt}\left(\frac{q}{r} \right) = -Gm^2M\frac{qr\frac{dr}{dt}-\frac{dq}{dt}r^2}{r^3}, \quad q = x,y,z \]

に相当するが、角運動量は定数であったから両辺を積分できて

\[ \small p \times l = Gm^2M\left[ \begin{array}{c} \frac{x}{r} + \Lambda_x \\ \frac{y}{r} + \Lambda_y \\ \frac{z}{r} + \Lambda_z \end{array} \right] \]

を得る。ここで、\(\small \Lambda​\)は積分定数である。

　具体的に計算しておくと

\[ \small \left[ \begin{array}{c} \Lambda_x \\ \Lambda_y \\ \Lambda_z \end{array} \right] = p\times l – \frac{Gm^2M}{r}\left[ \begin{array}{c} x \\ y \\ z \end{array} \right] \]

であるから、

\[ \small \left[ \begin{array}{c} \Lambda_x \\ \Lambda_y \\ \Lambda_z \end{array} \right] = \left[ \begin{array}{c} x(p_x^2+p_y^2+p_z^2)- p_x(xp_x + yp_y + zp_z) \\ y(p_x^2+p_y^2+p_z^2)- p_y(xp_x + yp_y + zp_z) \\ z(p_x^2+p_y^2+p_z^2)- p_y(xp_x + yp_y + zp_z) \end{array} \right]-\frac{Gm^2M}{r}\left[ \begin{array}{c} x \\ y \\ z \end{array} \right] \]

である。\(\small \Lambda_x,\Lambda_y,\Lambda_z\)は定数であるから、\(\small |\Lambda|^2= \Lambda^2_x+\Lambda^2_y+\Lambda^2_z\)も定数であり、地道に計算すると

\[ \small \begin{align*} \Lambda^2\;&=G^2m^4M^2+\left(p_x^2+p_y^2+p_z^2- 2\frac{Gm^2M}{r}\right)(l_x^2+l_y^2+l_z^2) \\ & =G^2m^4M^2+2mEL^2 \end{align*} \]

を得る。ここで、\(\small E\)はエネルギーであり、\(\small L^2=l_x^2+l_y^2+l_z^2\)​は軌道角運動量の二乗である。座標と積分定数のベクトルの内積から角度を求めると

\[ \small \cos \omega = \frac{q_x\Lambda_x+q_y\Lambda_y+q_z\Lambda_z}{|\Lambda|r} = \frac{l_x^2+l_y^2+l_z^2}{|\Lambda|r}-\frac{Gm^2M}{|\Lambda|} \]

となり、軌道半径\(\small r\)を計算できて

\[ \small r = \frac{\frac{l_x^2+l_y^2+l_z^2}{Gm^2M}}{1+\frac{|\Lambda|}{Gm^2M}\cos \omega} = \frac{\frac{L^2}{Gm^2M}}{1+\sqrt{1+\frac{2EL^2}{G^2m^3M^2}}\cos\omega} \]

が得られる。この式は

\[ \small \begin{align*} &a= -\frac{1}{2}\frac{GmM}{E} \\ &\epsilon = \sqrt{1+\frac{2EL^2}{G^2m^3M^2}} \end{align*} \]

とおいて、

\[ \small r = \frac{a(1-\epsilon^2)}{1+\epsilon \cos \omega} \]

と変形することができる。\(\small a\)は長半径(楕円の長い方の半径)、\(\small \epsilon\)は離心率とそれぞれ言われる。式を反対にすると

\[ \small \begin{align*} &E = -\frac{1}{2}\frac{GmM}{a} \\ &L^2 = Gm^2Ma(1-\epsilon^2) \end{align*} \]

が成り立つ。

　運動方程式を特定するため、回転方向を固定して、\(\small z=0\)を仮定する。この場合における座標は

\[ \small \begin{align*} &x(t) = r(t)\cos(\omega (t)) \\ &y(t) = r(t)\sin(\omega (t)) \\ &z(t) = 0 \end{align*} \]

と表すことができる。円運動と異なり、\(\small r(t)\)は定数ではないことに注意する。この運動について、速度と加速度をそれぞれ求めよう。簡便のため、

\[ \small \dot{r} = \frac{dr(t)}{dt}, \quad \ddot{r}=\frac{d^2r(t)}{dt^2}, \quad \dot{\omega}=\frac{d\omega(t)}{dt}, \quad \ddot{\omega}=\frac{d^2\omega(t)}{dt^2} \]

と表す。計算すると

\[ \small \begin{align*} &\frac{dx}{dt} = \dot{r} \cos \omega(t) – r\dot{\omega} \sin \omega(t) \\ &\frac{dy}{dt} = \dot{r} \sin \omega(t) + r\dot{\omega} \cos \omega(t) \\ &\frac{d^2x}{dt^2} = (\ddot{r} – r\dot{\omega}^2) \cos \omega(t) – (r\ddot{\omega}+2\dot{r}\dot{\omega}) \sin \omega(t) \\ &\frac{d^2y}{dt^2} = (\ddot{r} – r\dot{\omega}^2) \sin \omega(t) + (r\ddot{\omega}+2\dot{r}\dot{\omega}) \cos \omega(t) \end{align*} \]

を得る。

\[ \small \begin{align*} &m\frac{d^2x}{dt^2} = -\frac{GmM}{r^2}\frac{x}{r} = -\frac{GmM}{r^2}\cos \omega(t) \\ &m\frac{d^2y}{dt^2} = -\frac{GmM}{r^2}\frac{y}{r} = -\frac{GmM}{r^2}\sin \omega(t) \end{align*} \]

であったから、

\[ \small \begin{align*} &m(\ddot{r} – r\dot{\omega}^2) = m\frac{d^2r}{dt^2} -mr \left(\frac{d\omega}{dt} \right)^2 = -\frac{GmM}{r^2} \\ &r\ddot{\omega}+2\dot{r}\dot{\omega} = r\frac{d^2\omega}{dt^2}+2\frac{dr}{dt}\frac{d\omega}{dt} =\frac{1}{r}\frac{d}{dt}\left(r^2 \frac{d\omega}{dt} \right) = 0 \end{align*} \]

という対応関係を発見することができる。これらは面積速度一定の法則（ケプラーの第２法則）として知られている方程式である。面積速度一定の法則はニュートン力学における角運動量保存の法則と等価である。軌道半径に関する加速度を通常の方法で計算すると

\[ \small \begin{align*} m\frac{d^2r}{dt^2} \; & = m\frac{d}{dt}\left(\frac{x\frac{dx}{dt}+y\frac{dy}{dt}}{r} \right) \\ &=\frac{m}{r}\left(x\frac{d^2x}{dt^2}+y\frac{d^2y}{dt^2}+\left(\frac{dx}{dt}\right)^2+\left(\frac{dy}{dt}\right)^2 -\frac{\left(x\frac{dx}{dt}+y\frac{dy}{dt}\right)^2}{r^2} \right) \\ &=\frac{1}{r}\left(-\frac{GmM}{r}+m\frac{(x^2+y^2)\left(\left(\frac{dx}{dt}\right)^2+\left(\frac{dy}{dt}\right)^2\right)-\left(x\frac{dx}{dt}+y\frac{dy}{dt}\right)^2}{r^2} \right) \\ & = \frac{1}{r}\left(-\frac{GmM}{r}+\frac{\left(xp_y-yp_x\right)^2}{mr^2} \right) \end{align*} \]

第二項は軌道角運動量の二乗\(\small L^2\)であったから

\[ \small m\frac{d^2r}{dt^2}-\frac{L^2}{mr^3} = -\frac{GmM}{r^2} \]

を得る。

\[ \small m\frac{d^2r}{dt^2} -mr \left(\frac{d\omega}{dt} \right)^2 = -\frac{GmM}{r^2} \]

であったから、２つの式を比較すると

\[ \small \frac{d\omega(t)}{dt} = \frac{\sqrt{L^2}}{mr^2} = \frac{|L|}{mr^2} \]

が成り立つ。

　したがって、最終的な運動方程式は

\[ \small \begin{align*} &r(t) = \frac{a(1-\epsilon^2)}{1+\epsilon \cos \omega(t)} \\ &x(t) = r(t)\cos \omega(t) \\ &y(t) = r(t)\sin \omega(t) \\ &z(t) = 0 \\ &\frac{d\omega}{dt} = \frac{|L|}{mr^2} \end{align*} \]

と得られる。角運動量は

\[ \small \begin{align*} &l_x = yp_z – zp_y = 0 \\ &l_y = zp_x – xp_z = 0 \\ &l_z = xp_y – yp_x = m\dot{\omega} r^2 \\ &L^2 = l_x^2+l_y^2+l_z^2 = m^2\dot{\omega}^2r^4 \end{align*} \]

であることに注意する。エネルギーは軌道半径の計算から

\[ \small E = \frac{p^2}{2m} – \frac{GmM}{r}=\frac{L^2}{2mr^2} – \frac{GmM}{r}=-\frac{1}{2}\frac{GmM}{a} \]

である。

　また、この楕円運動が描く楕円の面積は長半径\(\small a\)と短半径\(\small b = a\sqrt{1-\epsilon^2}\)から

\[ \small S = \pi a b = \pi a^2 \sqrt{1-\epsilon^2} \]

で与えられるが、面積速度が

\[ \small \frac{dS}{dt} = \frac{1}{2}r^2\dot{\omega} = \frac{|L|}{2m} \]

であるため、この楕円運動の周期（1周して元の座標に戻るまでの時間）\(\small T\)は

\[ \small T = \frac{S}{dS/dt} = \frac{2m\pi a^2\sqrt{1-\epsilon^2}}{\sqrt{Gm^2Ma(1-\epsilon^2)}}=\frac{2\pi}{\sqrt{GM}}a^{3/2} \]

と計算できる。これはケプラーの第３法則として知られている式であった。

任意の方向に回転する楕円運動

　等速円運動と同様に楕円運動についても任意方向の回転に拡張しおこう。回転座標変換行列

\[ \small \left[ \begin{array}{ccc} 1 & 0 & 0 \\ 0 & \cos \phi & \sin \phi \\ 0 & -\sin \phi & \cos \phi \end{array} \right] \left[ \begin{array}{ccc} \cos \theta & 0 & \sin \theta \\ 0 & 1 & 0 \\ -\sin \theta & 0 & \cos \theta \end{array} \right] \left[ \begin{array}{ccc} \cos \psi & \sin \psi & 0 \\ -\sin \psi & \cos \psi & 0 \\ 0 & 0 & 1 \end{array} \right] \]

を掛ければ良い。それぞれ、\(\small x,y,z\)軸周りの回転を表す。\(\small z\)軸周りは位相（初期座標）の値を変えるだけ（\(\small \omega (t)\)を\(\small \omega (t) + \psi\)と置けばよい）であるから無視して計算すると

\[ \small \begin{align*} &x(t) = r(t)\cos\theta \cos(\omega (t))\cos \phi + r(t)\sin \phi \sin(\omega (t)) \\ &y(t) = r(t)\cos\theta \cos(\omega (t))(-\sin \phi) + r(t)\cos \phi \sin(\omega (t)) \\ &z(t) = r(t)\sin\theta \cos(\omega (t))\ \end{align*} \]

のように変形できる。\(\small \phi=\theta=0\)としたときに元の式になるように座標軸を入れ替えている。\(\small \theta,\phi\)の値を変えることで任意の方向に傾いた楕円運動を表現できるということである。

　速度や運動量、角運動量をそれぞれ求めておこう。等速円運動と異なり、\(\small r,\omega\)が定数ではないことに注意する。速度は

\[ \small \begin{align*} &v_x =(\dot{r}(t)\cos(\omega (t)) -r(t)\dot{\omega}(t)\sin(\omega (t)))\cos\theta\cos \phi \\ & \qquad + (\dot{r}(t)\sin(\omega (t))+r(t)\dot{\omega}(t)\cos(\omega (t)))\sin \phi \\ &v_y = (\dot{r}(t) \cos(\omega (t))-r(t)\dot{\omega}(t)\sin(\omega (t)))\cos\theta(-\sin \phi) \\ & \qquad + (\dot{r}(t)\sin(\omega (t))+r(t)\dot{\omega}(t)\cos(\omega (t)))\cos \phi \\ &v_z = (\dot{r}(t)\cos(\omega (t))-r(t)\dot{\omega}(t)\sin(\omega (t)))\sin\theta \\ &v = \sqrt{v_x^2+v_y^2+v_z^2} = \sqrt{r^2(t) \dot{\omega}^2(t)+\dot{r}^2(t)} \end{align*} \]

である。運動量は\(\small p=mv\)で計算すればよいだろう。角運動量を計算すると

\[ \small \begin{align*} &l_x = yp_z – zp_y = m\dot{\omega}(t) r^2(t)\sin\theta\cos \phi \\ &l_y = zp_x – xp_z = m\dot{\omega}(t) r^2(t)\sin\theta\sin \phi \\ &l_z = xp_y – yp_x = m\dot{\omega}(t) r^2(t)\cos\theta \\ &L^2 = l_x^2+l_y^2+l_z^2 = m^2\dot{\omega}^2(t)r^4(t) \ \end{align*} \]

となる。楕円運動では角運動量は定数であるから

\[ \small r^2(t)\dot{\omega}(t) = \text{const.} \]

であることは前節で述べたとおりである。

残された問題

　これで楕円運動の運動方程式を一通り計算できた気になるかもしれないが、実を言うと\(\small \omega(t)\)やそこから派生して計算される\(\small r(t)\)の具体的な関数形というのを計算していないことに気づくだろう。等速円運動のように

\[ \small \omega(t) = \bar{\omega} t \]

と置くことはできない。一般に、楕円運動では焦点に近づくほど速度が速く、焦点から遠いほど速度が遅くなる傾向がある。この性質を満たすように\(\small \omega(t)\)の関数形を求める必要がある。

　\(\small \omega(t)\)は微分方程式

\[ \small \frac{d\omega(t)}{dt} = \frac{|L|}{mr^2(t)} \]

を満たすが、

\[ \small \begin{align*} &r(t) = \frac{a(1-\epsilon^2)}{1+\epsilon \cos \omega(t)} \\ &L^2 = Gm^2Ma(1-\epsilon^2) \ \end{align*} \]

であったから代入すると

\[ \small \frac{d\omega(t)}{dt} = \sqrt{GM}\frac{(1+\epsilon \cos \omega(t))^2}{a^{3/2}(1-\epsilon^2)^{3/2}} \]

が成り立つ。この微分方程式を解けばよいということであるが、この微分方程式の解を求めることができるだろうか。まあ、IT技術が発展した現代ではコンピュータで計算すればいいだろうという話ではあるけど・・・これはもう少し異なる方程式から特定する必要があり、ケプラー方程式という代数方程式を解くことで解を求めることができる。その際にベッセル関数を利用することはベッセル関数の記事で述べたとおりである。ケプラー方程式については次の稿で詳しく考察しよう。