『重力の使命』であそぶ―(2)メスクリンの歳差運動

メスクリンの歳差運動

「もちろん、歳差はかなり急速なはずだ。赤道面が異常にふくらんでいるため、たとえこの惑星の質量の大部分が中心付近にあっても、太陽の引力に相当な手がかりを提供するからである。歳差の周期までは計算しなかったが、居住可能半球が二、三千年ごとに交替するため、その住民が高度な文明を建設する妨げになった、と考えたい人がいるなら、それに反対をとなえる気はない。」(『重力の使命』p.314)

「メスクリン創成期」のこの一節. ここがずっと気になっていた. 確かに扁平さは大きなモーメントを働かせる. しかし一方で自転角運動量の大きさはその軸の動きにくさにも寄与する. 「計算しなかった」と言っている以上は仕方ないが二, 三千年というのは妥当な値だろうか?
　これを計算するにあたり, もう一つメスクリンの特異な性質が絡んでくる. 軌道離心率の大きさだ. 軌道が円に近ければ(初歩的な力学の教科書に載っているという意味で)歳差運動の周期はすぐに計算できる. しかしメスクリンの軌道離心率は極めて大きい. 何せ遠日点では近日点より6倍も太陽から遠くなる.
　このことの重要性は地表の住民にとっても大きい. 近日点で夏を迎えるか冬を迎えるかによって太陽から受け取る放射エネルギーは大きく変わる.
　下の図に受け取るエネルギーの一日あたりの総量(カラーマップ)と一日に占める昼の長さの一年間での変化(モノクロ, 光度等の絶対値を考慮していない相対値)を示した. 縦軸は緯度(地理緯度), 横軸はメスクリン年. αは下で定義している座標だが, α=0のとき近日点で南半球が夏を迎えていると考えればよい. 0.5(近日点)で南半球の受け取る熱がきわめて多くなっていることがわかるだろう.

　実際歳差運動はどの程度の周期で起こるのか. 引き続き一様密度のメスクリンで計算を行う*1.
　扱う問題を改めてきっちり述べると, 重力下の束縛軌道にある剛体が自転軸の方向を変える周期はいくらになるか, ということになる. 剛体の運動方程式を一から求めて数値計算をすることにする.

角速度ベクトル

　『エターナル・フレイム』で習ったとおり, 4次元空間の回転は2つの単位四元数によって表せる. 3次元空間だとそれが1つになる.

　ではまず3次元空間の回転を四元数を使ってやっていく.
ここだけの記法として
・矢印付きは3-ベクトル：実数成分がゼロの四元数. クロス積・ドット積がある.
・太字は四元数.四元数どうしの積がある.
とする.
x(t)で位置を表す. 四元数qによってこれを回転させるときの角速度ベクトルを求める.

ここで,

とするとωは3-ベクトルで角速度になる.
なぜなら,

というように時間微分を表せるため.

　x(t)が剛体を構成する点であるとする. x0は剛体に固定された座標系での位置になる. Vをその領域とすると, 運動エネルギーの総和は,

という形で表せる.
ここでqを次のようにパラメーターα,β,γで表す.

4次元球面上を隈なく動くのですべてのqを表せることがわかる.
パラメーターの範囲としては

をとると一対一になる.(があとの計算でこれを考える必要はとくにない).

こう書くと見慣れたSO(3)の元であることがわかる.
角速度は

これは剛体に固定された系での角速度. 恒星系(原点を中心星にとる. 中心星は十分重く不動と仮定)では

慣性モーメント
　慣性モーメントは次のように求められる. なお惑星質量のほうに下添え字の○(マル)を, 恒星質量には＊(アスタリスク)を付けている.

扁平な回転楕円体として(引き続き)

をとると,慣性モーメントが求まる. z軸に平行な方向に||、垂直な方向に⊥を付けている.

剛体回転の運動エネルギーは

ハミルトニアンがH=K+(α,β,γの時間微分によらない項)という形なら,
共役運動量が求められる.

(この形を見ると自由回転ではpαとpγが保存することがわかる.)

重力四重極モーメント
　楕円体が外部に作る重力場を改めて考える. 3次元極座標を使うと十分遠方での重力場が求められる(あえてここで一から求めるようなことでもない気はするが確認がてら). ここだけ座標原点を剛体の中心にとる.
ポアソン方程式を解くと, 領域V内での密度分布がμ(r)のとき, 重力ポテンシャルは次のようになる

Vをz軸回転対称, その内部で密度は一様とすると,

δはrとr'のなす角.ここで1行目から2行目ではルジャンドル多項式の加法定理を使っている. 回転対称なので1次以上のルジャンドル陪多項式の項は消える.

　積分を計算する.ルジャンドル多項式は次数lの偶奇と関数の偶奇が一致するため, zについて対称な区間で積分するとlが奇数のとき消える.
l=0の項

l=2の項

l=2までの近似(つまり四重極までの展開)で重力ポテンシャルは

となる. 座標を恒星系に戻すと, 楕円体に固定された軸と, 重心から恒星への方向となす角Δに対して 重力エネルギーは

となる. 点重力源の場合と違って角度に依存する項があるのが重要(この項がないと剛体は自由運動をすることになり自転角運動量は変化しない).

　重心(Center of Mass)の運動エネルギーも極座標で表す.

これらから(ラグランジアンを経て)ハミルトニアンは

という形になる. ここでΔは恒星-惑星方向を向くベクトルと惑星固定のz'軸がなす角で

　さらに近似を使って解析的に扱うこともできるはずだが, ここから数値計算に進んで「実験」してしまう. シンプレクティック数値積分法を使おう.

シンプレクティック数値積分法
　ハミルトン形式で, 時間発展とはハミルトニアンを生成子とした無限小変換である. (q,p)→(q~,p~)の座標変換を微小時間τの時間発展のつもりでとる. 次の母関数Wによる正準変換として近似することになる.

すると座標変換は次のようになる.

　一般的には陰的で次のステップに進むために連立方程式を数値的に解かねばならないが, (幸運にも)この場合は陽的に計算できる. ちょっと式変形して並べ替えると次のようにチルダ付き変数(次のステップ)の値が決まっていく.

上から順に計算すれば左辺はチルダなし(前のステップ)座標のみからただちに計算できることに注意. 以下チルダ付き運動量から6つのチルダ付き座標が定まる.

　無次元化を行うために軌道半径, 軌道角運動量, 自転角運動量の基準を与える. 遠日点をt=0にとるのが都合がいい(公転1周分の変化量をとるにあたり遠日点での時間変化率が最も穏やか).

　Aは軌道長半径, Ω0は平均角速度, ω0は自転角速度を想定している.　(上の式を見ればわかる通りpα+pφも保存量になっている. これはz方向角運動量の和)

t=0を遠日点とするため初期値は次のようにとる.

そして無次元化は次のハット付き変数を使うことで行う.

さらに時間の基準としてTdayとTyearを導入する. 文字通り1日の長さと1年の長さ.

　時間ステップτについては, Tyearで割った値hをとる(h=τ/Tyear). γのみがTdayのスケールで変動する(TdayはTyearよりじゅうぶん短いと仮定)が, ほかの座標を求めるうえでは必要がない. グリニッジ天文台が何時なのか知らなくてもいい.
　これによって以下のように計算機に入れられる形で式を得る(次ステップを＋の上添え字付きで表した).

　歳差運動はαの変化である. つまり歳差運動を見るというのはαの変化を見ることに他ならない.

　上の式を使って時間に従ってαがどう変化するか見る. なお代入する数値に関して, 中心星の質量は0.70太陽質量としている. これははくちょう座61番星aの質量の観測値. βは地軸の傾きで28°. 1日の長さは17.75分, 1年の長さは1800地球日(設定に従っている)
　離心率はややはっきりしていないが, 「メスクリン創成期」の図や, 高熱にさらされるのが1年の4分の1におよぶとの記述により0.71(≒1/√2)を選んだ.

　グラフはαの初期値(遠日点)ごとに, φに対して変量Δαをプロットしたもの. 縦軸の単位に注意. 六十分法で1000倍されている.

　公転一周当たりの変化量に対して曲線フィッティングを行うときれいに定数+コサインの定数倍に乗る.

適用された"curve fit"の関数f(α)は

αは上で見たように公転一周の間にかなり大きく時間変化率が変わるが, このフィッティング関数にしたがって粗視化する*2. すなわち,

という微分方程式に従うとする. これはすぐに解けて,

ただし(0:π)区間の外では連続となるように定数部分を変えて接続する.
　これをプロットすると下のようなグラフになる*3.

　グラフを見ても分かるが, αが一周する(2π変化する)のにかかる時間は

およそ10万メスクリン年, つまり50万地球年であることがわかる.

　これは地球の歳差運動の周期*4より20倍長い. 少なくともUメスクリンに関して言えば, 近日点で夏を迎える極が二, 三千年であるとの見積もりは短すぎる.
　では, より中心に集中した密度分布だとどうなるだろうか. 角運動量は小さくなる. 潮汐力によるモーメントも小さくなる. 上のハミルトニアンを見ればわかる通り, ハミルトニアンのうち剛体運動部分と重力四重極子ポテンシャルはともに質量にはよらず慣性モーメントだけで書けている*5. おそらく慣性モーメントの縦横比がおおきく変わらない限り歳差運動周期もオーダーは同程度だろう.

軌道離心率
　軌道離心率は歳差運動にどのように影響しているだろうか. 軌道長半径, 質量, 慣性モーメント等は変えずに軌道離心率とそれに伴って変わる近日点距離, 軌道角運動量のみを変化させる. 上と同じフィッティング関数のp, qによってそれぞれの場合が特徴づけられ, 歳差運動周期が求められる.


　eに対して単調減少で, 離心率がゼロなら周期は約10倍の100万メスクリン年になる. eに関して逆S型の関数(下のグラフはセミログ)になっている様子は分かるが, これにうまくあてはめられる関数は今のところ見つけられてない.

　なおpとqの関係は次のようになっている(横がp,縦がq). 一見すると二次曲線的だがあてはめはうまくいかなかい.

　解析的な取り扱いが今後の課題として残る. 他の変数に関する依存性なども考慮して帰納的に調べていくことができるだろう.

まとめ
　一様密度の扁平な回転楕円体の, 重力束縛軌道下での歳差運動について調べた. とくに16木星質量の惑星が0.70太陽質量の恒星のまわりを1800地球日, 離心率0.71, 地軸の傾き28°で自転する惑星メスクリンを想定した.
　メスクリンが一様密度ならその歳差運動周期は10万年ほどだろう. 地球基準で言えば*6, 進化スケールの時間で常にその気候への影響にさらされるはずだが, メスクリンの生き物たちはこれにどう適応しているのだろうか.
　文明の興隆スケールの時間ではそう何度も経験していないかもしれない. メスクリン人たちは気候変動にどう対処していくだろうか.
　歳差運動はSFになる. ほかならぬ地球の気候変動にも関わっているという. エキセントリックプラネットでは特にその影響は大きく, 歳差運動周期に適応した生命たちを想像するとおもしろい. ロバート・J・ソウヤー『イリーガルエイリアン』のネタバレになるが, この作品ではアルファ・ケンタウリの三重連星にある架空の惑星の軌道運動の特異な周期が物語に大きく関わってくる. そういうイメージ.

*1:密度分布が中心に集中したときの影響はきちっと計算する必要がある. しかし極端なことを言えばδ関数的に分布していたらそもそも歳差運動は起こらないしそれはもはや剛体ではない. またいずれ.

*2:こんなことをせずともαが2π変化するまで計算を続ければいいように思われるかもしれない. しかし実はこの計算では公転一周を100万分等分しており自分のPCではminオーダーの時間がかかる.これより少ないとなぜか誤差が収束しないのだ. おそらくハミルトニアンの剛体項の特異性（ケプラー運動部分はかなり粗くしても妥当な値を出す）などに由来すると考えられるが原因を特定できていない. ただしこれより刻み幅を小さくしていくとそれに比例して誤差が小さくなっていく様子が観察される. これはすなわち一次解法であることを意味する. 二次以上に改善することもできるがそうすると今度は陰的に解かねばならなくなる. 　シンプレクティック数値積分法のよいところは, （陽的）ルンゲ・クッタ法のように誤差が蓄積せずハミルトン力学系を数値解析するにあたり信頼性が高いところだ. 実装が簡単という利点ももちろんあるが一周で打ち切るならややもったいない.

*3:自転の方向は順行, つまり自転と公転の角速度ベクトルのz成分が同じ向きになるとしている. しかしクレメントは逆行している可能性について言及している(作中には反映されていないが). その場合, 単にこのグラフは左右が入れ替わるだけで周期は変わらない. このことは一日の長さの符号を反転することで確かめられる.

*4:地球の歳差運動は太陽より月の影響が大きい

*5:対称コマの重力四重極モーメントが慣性モーメントで表せることをMcCulloughの公式と呼ぶらしい.参照 McCullough's Formula

*6:メスクリンはきわめて冷たい.『重力の使命』でメスクリン人の思考があまりにも人間的すぎる, という指摘があるが, それにもまして重大なのは同じ時間の中で交流できている点かもしれない. メスクリン人もやはり化学ベースで代謝を行っているので, 地球人よりゆっくりした時間を生きていると考えたほうが自然に思われる. 『竜の卵』のチーラと人間の比ほどではないが, コミュニケーションにおいて似たような状況が生まれるだろう.