衛星系研究会 平成25年8月7～9日 定山渓温泉 渓流荘（札幌）高軌道傾斜角を持つ三体系の階層安定性

斎藤 正也 (統計数理研究所)
Виктор В. Орлов (Санкт-Петербургский государственный университет)谷川 清隆 (国立天文台)

階層三体系
2つの交叉しない楕円軌道で近似できる三体系

m1
m0m2

外連星 (outer binary)

内連星 (inner binary)軌道要素
内連星: (m0,m1)外連星: ((m0+m1), m2)
位置ベクトル(ヤコビ座標系)

r1, r2m1r1

ωλ2 λ1 2

m0m2

♈
a1
近点距離

r2

q2Q2
軌道長半径 a2≡ (q2 + Q2)/2
離心率 e2≡ 1 - q2/a2

軌道要素
相互傾斜角 I
昇交点経度 Ω軌道長半径 a近点離心率 ef

近点引数 ωωΩ

軌道面1I
昇交点

真近点角 f遠点

運動方程式 (ヤコビ座標系)
運動方程式

dqj/dt = ∂H/∂pj, dpj /dt = - ∂H/∂pj, (j=1, 2)
ヤコビ系 (rj, pj) (j=1,2) と慣性系 (ρi, πi) (i=1,2,3)との対応

階層安定• 階層安定• 天体の順序の入れ替わりや、系の分裂が起きないこと• (階層)安定限界• 内外の連星が最も接近した階層安定な軌道配置q2• 一般に q2/a1の下限として記述される• ねらい– 観測された、またはシミュレーションで形成されa1

た、惑星系・連星系が長期間(母天体の寿命程度)に相続するかを軽量なモデルで網羅的に調べる
a1:内連星の軌道長半径q2:外連星の近点距離安定限界を策定する試み• ヒル安定性 (Golubev 1968, Marchal & Saari 1975; Marchal & Bozis1982; …): m0,m1の運動可能領域とm3の運動可能領域とが分離されている。制限三体問題(m3 = 0)の同様の概念の拡張。
• 数値計算 (Harrington 1977, Pilat-Lohinger & Dvorak 2003, Holman& Wiegert 1999): ヒル安定限界はあまり実用的でない。ヒル安定かつ脱出可能な領域は小さくない。ある決まった時間積分し、その間に脱出などの軌道不安定が発生するかどうか観察する。
• ランダムウォーク近似 (Orlov et al.,2008):• 共鳴重合 (Chirikov 1979, Wisdom 1980, Mydryk & Wu 2006)• 赤字は制限三体問題。一般三体で三次元の計算はすくない (ex.
Georgakarakos 2012, e1 = e2 = 0だが多数の質量比)

ランダムウォーク近似 (Orlov et al. 2008)• 1周ごとの軌道要素の変化は、近点通過時の離角Δλに依存し、についての平均は0である (摂動論)。
• しかし、毎周のΔλの実現はランダムと思ってよい(物理的直観)• 軌道要素は、ランダムウォークに従うと仮定することになる。
1Δλ2Δλ・・・• 10000周の間の軌道要素の変化が0.1%以下ならば階層安定と見なす

平均運動共鳴と共鳴重合• 平均運動共鳴: 公転周期の比が整数になる(k1 , k2を整数として、k1λ1 + k2 λ2 ≈ 0)場合。軌道長半径比上に離散的に出現する。
• 摂動ポテンシアルRの展開
R=∑ A (a, e, I ) cosψ , 臨界引数 ψ= k λ + k λ + k ϖ + k ϖ
kkk

1 123142k(ただし、k1 + k2 + k3 + k4 =0)• 共鳴の近傍で、このRの下での解を摂動展開のかたちで求めることはできないが、Rの中から永年項(すべてのk1 = k2 = 0の項)とひとつのk1, k2 ≠ 0である(k1, k2, k3, k4)を抽出したRkの下では、系は可積分となり、ψk = const.になる解 (周期軌道) が存在し、その解が線形安定である場合には、近傍に ψkが秤動する領域が広がる。

平均運動共鳴と共鳴重合
共鳴重合: カオスの発生モデル
Rkの下での運動

Rk’ の下での運動
安定領域

安定領域

Rk+Rk’ の下での運動

カオス的運動安定領域

安定領域

安定領域が縮小し、もとの安定領域の重なりの部分では運動がカオス的になる

平均運動共鳴と共鳴重合• m0とm1がcomparableなときは、sub-resonances (k1,k2を共通にする、異なるk3,k4に対応する共鳴の集合)の重なりが階層安定限界をきめる(Mydryk & Wu, 2006)。

5:1

a/a'
6:1

e'

吉田(1975)の脱出条件
Outer binaryの軌道半径が単調増加するための十分条件

この条件を初めて満たす時刻tescを「脱出時間」と呼ぶ。

研究に利用するツール• 数値積分 (多項式による補外法)• 脱出 (楕円軌道+双曲軌道)の十分条件• 脱出時間を初期値平面 (e2, q2) 上の分布を観察する• スパコンを利用した広範な初期条件の長期間の積分

平面(I = 0)の場合

脱出軌道の特徴 (q2/a1が小さい場合)

質量比 m0=0.5, m1= 0.5, m2 = 1/100 (Saito et al., 2012)(会合周期)• 2～4回の急激なe2の変化を経て、脱出に至る軌道の例。
• q2/a1を大きくすると、回数が増えるので、進化パターンはRW的になる。

安定軌道の特徴 (q2/a1が大きい場合)

質量比 m0=0.5, m1= 0.5, m2 = 1/100

q2 = 2.4

q2 = 2.6

q2 = 2.8• Outer binaryが近点(f2=0)にいるときの軌跡• q2について単調増加な臨界値 e2,crit(q2)が存在する• e2 < e2,crit(q2) ならば、軌道は準周期的• e2 > e2,crit(q2) ならば、酔歩的な変動を経てe2 > 1(脱出)に至る。
• この境界を数値的な安定限界と考える。

初期値面上の構造のまとめ

質量比 m0=0.5, m1= 0.5, m2 = 1/100「脱出時間」(単位:会合周期。最大値の256は103単位時間程度)

0.6

不安定

0.5

l(0)
l2(0)==00°l(0)l2(0)==3030°l2(0)==6060°l(0)l2(0)==9090°l(0)l2(0)==120
120°
l(0)
l2(0)==150
150°
l(0)
e = e(qst)

e2(0)

0.4
0.3
0.20.1

安定

01

1.5

2

2.5
3
q2(0)/a1(0)3.544.5他の質量比への適用

m2/(m0 + m1) = 1/5000, 1/50, 1/20, 1/10, 1/5, 1/2, 1 (Saito et al., 2013)

11/5

1/20

1/50

3次元(I > 0)の場合

質量比 m0 = m1 = m2 = 1

I=0, tesc < 106

1:9
1:8
1:7
1:6

I=10, tesc < 106

1:9
1:8
1:7
1:6

I=10, tesc < 106

1:9
1:8
1:7

𝜆1 − 8𝜆 + 𝜛1 + 6𝜛2t < 2.4x10^5, t > 2.4x10^51:6e2

I=45, tesc < 106

1:9
1:8
1:7
1:6

Iの増加にともない、平均運動共鳴に対応する棒状の不安定領域が拡大する

I=50, tesc < 106

1:91:81:7?

Iの増加にともない、平均運動共鳴に対応する棒状の不安定領域が拡大する

I=55, tesc < 106

1:9
1:8
1:7Iの増加にともない、平均運動共鳴に対応する棒状の不安定領域が拡大する

I=60, tesc < 106

1:91:8

1:71:6Iの増加にともない、平均運動共鳴に対応する棒状の不安定領域が拡大する

I=90, tesc < 106

1:9
1:8

1:7
1:6Iの増加にともない、平均運動共鳴に対応する棒状の不安定領域が拡大する

考察• Iを大きくすると、この領域は太くなる。一般にe,Iが大きくなるとMMRの(aの幅で測って)stability regionが大きくなる。その結果、resonanceoverlapが起こる領域が拡大し、不安定領域が拡大すると解釈できる。
• Mydruk & Wu(2006)を踏まえると、ひとつひとつの「棒」はひとつの周期比に対応するMMRのsubresonances同士の重なりであろう。m0≈ m1の制限三体問題よりさらに相互作用が強いので、共鳴重合はより強化されるはずである。
• Iをさらに大きくすると(I > 45)、異なる周期比のMMR同士でresonance overlapが起こり、不安定領域の本体と融合していくと解釈できる。
• 脱出に至るまでの時間の増大。数値計算の限界。

計算量
ひとつのI値に対する計算
N(q_2) = (5.0 – 3.0) / 0.005 = 400 , N(e_2) = (1.0 – 0.0)/ 0.02 = 50
N_orbit = 400x50 = 20,000
t_orbit = 10min (t_esc < 10^6の場合)
N_proc = 512
t_run = t_orbit x N_orbit / N_proc = 390.625min = 6.5h
(cf. t_run1 = t_orbit x N_orbit = 140 days)