ラグランジュ点の理論 — 制限三体問題の平衡解と宇宙ミッションへの応用

ジェイムズ・ウェッブ宇宙望遠鏡(JWST)は、地球から約150万km離れた「太陽-地球系の第2ラグランジュ点(L2)」と呼ばれる特別な場所に位置しています。地球の影に入ることなく、常に太陽と地球の反対側にとどまりながら深宇宙を観測しています。この位置が特別なのはなぜでしょうか? なぜ普通に地球を周回する軌道ではなく、150万kmも離れた場所が選ばれたのでしょうか?

答えはラグランジュ点にあります。太陽と地球という2つの大きな天体が作る重力場の中に、小さな物体(衛星)がほぼ静止していられる5つの特別な位置が存在します。これらの点では、太陽と地球の重力と遠心力がちょうど釣り合い、衛星は太陽と地球に対して一定の位置関係を保ったまま太陽のまわりを公転できます。

ラグランジュ点を理解すると、以下の応用が見えてきます。

  • 宇宙望遠鏡の配置: JWST(L2)、SOHO(L1)など、安定した観測環境を提供する位置の選定
  • 太陽観測: L1に配置した衛星は常に太陽を監視でき、太陽嵐の早期警報が可能
  • 深宇宙探査のゲートウェイ: L1やL2を中継点とした月・火星探査ミッションの計画
  • トロヤ群小惑星: 木星のL4/L5に集まる数千個の小惑星群の力学的な理解

本記事の内容

  • 三体問題の難しさと制限三体問題
  • 回転座標系での運動方程式
  • ヤコビ積分(エネルギー的な保存量)
  • 5つのラグランジュ点の導出
  • 安定性解析(L1-L3不安定、L4-L5条件付き安定)
  • ハロー軌道とリサージュ軌道
  • Pythonでゼロベロシティ曲面とラグランジュ点の可視化
  • JWST(L2)ミッションの実例

前提知識

この記事を読む前に、以下の記事を読んでおくと理解が深まります。

三体問題の難しさ

二体問題と三体問題の壁

二体問題(太陽と惑星、地球と衛星)はケプラーが解いた通り、解析的に完全に解ける問題です。軌道は楕円・放物線・双曲線のいずれかで、未来永劫にわたって正確に予測できます。

ところが、天体を3つにした途端、問題は劇的に難しくなります。一般の三体問題には解析的な一般解が存在しないことが、ポアンカレによって19世紀末に証明されました。3つの天体が互いに重力で引き合う運動は、初期条件のわずかな違いが指数関数的に増大するカオスの典型例です。

しかし、すべてが解析不可能というわけではありません。特定の仮定のもとでは、三体問題の重要な性質を解析的に調べることができます。その最も有名な例が制限三体問題です。

制限三体問題の仮定

制限三体問題(Restricted Three-Body Problem, R3BP)では、以下の仮定を置きます。

  1. 2つの主天体(太陽と地球、地球と月など)は互いの重力のみで運動し、ケプラー運動(二体問題の解)に従う
  2. 第3の天体(衛星)は質量が無視できるほど小さく、主天体の運動に影響を与えない
  3. 円制限三体問題(CR3BP): さらに主天体が円軌道で互いの周りを公転すると仮定

3番目の仮定(円軌道)を加えたものが円制限三体問題で、ラグランジュ点の解析に使われる標準的な枠組みです。地球と太陽の離心率は0.0167と非常に小さいので、この近似は十分に妥当です。

制限三体問題の枠組みが定まったところで、次に回転座標系での運動方程式を導出します。この座標系を使うことで、ラグランジュ点が「静止した平衡点」として自然に現れます。

回転座標系での運動方程式

座標系の設定

2つの主天体(質量 $m_1 > m_2$)が共通重心のまわりを角速度 $\omega$ で円運動している状況を考えます。この系と一緒に回転する座標系(回転座標系、synodic frame)を設定します。

回転座標系では

  • 原点: 2つの主天体の共通重心
  • $x$ 軸: $m_1$ から $m_2$ に向かう方向
  • $z$ 軸: 軌道面に垂直(角運動量方向)
  • $y$ 軸: $z \times x$ で完成する右手系

主天体間の距離を $L$ とし、これを長さの単位に取ります。質量パラメータを

$$ \mu = \frac{m_2}{m_1 + m_2} $$

と定義すると、$m_1$ は $(-\mu, 0, 0)$ に、$m_2$ は $(1-\mu, 0, 0)$ に位置します。回転の角速度は $\omega = 1$(時間の単位を公転周期/$2\pi$ に取る)です。

太陽-地球系では $\mu \approx 3.003 \times 10^{-6}$(太陽が圧倒的に重い)、地球-月系では $\mu \approx 0.01215$ です。

運動方程式の導出

回転座標系では、コリオリ力と遠心力が現れます。第3天体(衛星)の位置を $(x, y, z)$ とすると、運動方程式は

$$ \ddot{x} – 2\dot{y} = \frac{\partial U}{\partial x} $$

$$ \ddot{y} + 2\dot{x} = \frac{\partial U}{\partial y} $$

$$ \ddot{z} = \frac{\partial U}{\partial z} $$

ここで $U$ は有効ポテンシャル(遠心力ポテンシャルと重力ポテンシャルの和)です。

$$ U(x, y, z) = \frac{1}{2}(x^2 + y^2) + \frac{1-\mu}{r_1} + \frac{\mu}{r_2} $$

ここで

$$ r_1 = \sqrt{(x+\mu)^2 + y^2 + z^2}, \quad r_2 = \sqrt{(x-1+\mu)^2 + y^2 + z^2} $$

は衛星から $m_1$ と $m_2$ までの距離です。

運動方程式の左辺で $-2\dot{y}$ と $+2\dot{x}$ がコリオリ力(回転座標系で運動する物体に現れる見かけの力)、有効ポテンシャルの $\frac{1}{2}(x^2+y^2)$ 部分が遠心力に対応しています。

この運動方程式には重要な保存量が存在します。それがヤコビ積分です。

ヤコビ積分

ヤコビ定数の導出

運動方程式に $\dot{x}$、$\dot{y}$、$\dot{z}$ をそれぞれ掛けて足し合わせます。

$$ \dot{x}\ddot{x} + \dot{y}\ddot{y} + \dot{z}\ddot{z} = \dot{x}\frac{\partial U}{\partial x} + \dot{y}\frac{\partial U}{\partial y} + \dot{z}\frac{\partial U}{\partial z} $$

左辺のコリオリ項は打ち消し合います($-2\dot{y}\dot{x} + 2\dot{x}\dot{y} = 0$)。右辺は $U$ の全微分 $dU/dt$ です。したがって

$$ \frac{d}{dt}\left(\frac{1}{2}(\dot{x}^2+\dot{y}^2+\dot{z}^2)\right) = \frac{dU}{dt} $$

両辺を積分すると、回転座標系での速度の大きさ $v^2 = \dot{x}^2+\dot{y}^2+\dot{z}^2$ について

$$ v^2 = 2U(x,y,z) – C_J $$

ここで $C_J$ は積分定数で、ヤコビ定数(Jacobi constant)と呼ばれます。

$$ C_J = 2U(x,y,z) – v^2 $$

ヤコビ定数は円制限三体問題における唯一の積分定数です。二体問題ではエネルギーと角運動量がそれぞれ保存しましたが、回転座標系ではそれらが混ざり合って一つの保存量(ヤコビ定数)にまとまります。

ゼロベロシティ曲面

$v^2 \geq 0$ という物理的条件から、$2U(x,y,z) \geq C_J$ でなければなりません。したがって

$$ 2U(x,y,z) = C_J $$

を満たす曲面が、衛星の運動が許される領域と許されない領域の境界になります。これがゼロベロシティ曲面(Zero Velocity Surface, ZVS)です。

ゼロベロシティ曲面は「衛星がここを超えると速度の2乗がマイナスになってしまう(物理的にあり得ない)」境界です。つまり、ヤコビ定数の値によって衛星が行ける範囲が制限されます。これは制限三体問題の力学を理解する上で極めて重要な概念です。

ゼロベロシティ曲面の構造は $C_J$ の値に依存して劇的に変化し、その変化が起こる臨界的な $C_J$ の値がまさにラグランジュ点に対応します。では、ラグランジュ点を数学的に求めましょう。

5つのラグランジュ点の導出

平衡条件

ラグランジュ点は、回転座標系で衛星が静止できる点($\dot{x}=\dot{y}=\dot{z}=\ddot{x}=\ddot{y}=\ddot{z}=0$)です。運動方程式に代入すると

$$ \frac{\partial U}{\partial x} = 0, \quad \frac{\partial U}{\partial y} = 0, \quad \frac{\partial U}{\partial z} = 0 $$

$z$ 方向の条件から、$\partial U/\partial z = -z[(1-\mu)/r_1^3 + \mu/r_2^3] = 0$ となり、$z = 0$(軌道平面内)であることがわかります。ラグランジュ点はすべて軌道平面内に存在します。

共線ラグランジュ点(L1, L2, L3)

$y = 0$ とすると $\partial U/\partial y = 0$ は自動的に満たされます。$x$ 軸上で $\partial U/\partial x = 0$ の条件から

$$ x – \frac{(1-\mu)(x+\mu)}{|x+\mu|^3} – \frac{\mu(x-1+\mu)}{|x-1+\mu|^3} = 0 $$

この方程式は $x$ の区間によって3つの解を持ちます。

L1($m_1$ と $m_2$ の間): $m_2$ の「手前」に位置します。太陽-地球系では地球の内側約150万km。$m_1$ と $m_2$ の重力が互いに部分的にキャンセルし、遠心力と釣り合う点です。

L2($m_2$ の外側): $m_2$ の「向こう側」に位置します。太陽-地球系では地球の外側約150万km。JWSTが配置されている場所です。

L3($m_1$ の反対側): $m_2$ と反対側の $m_1$ の向こうに位置します。太陽-地球系では太陽の反対側約150万km。

$\mu \ll 1$ の場合、L1とL2の位置は $m_2$ からの距離で次のように近似できます。

$$ r_{L1} \approx r_{L2} \approx \left(\frac{\mu}{3}\right)^{1/3} L $$

太陽-地球系($\mu \approx 3 \times 10^{-6}$, $L = 1$ AU)では

$$ r_{L1} \approx r_{L2} \approx \left(\frac{3 \times 10^{-6}}{3}\right)^{1/3} \approx 0.01 \, \text{AU} \approx 1.5 \times 10^6 \, \text{km} $$

三角形ラグランジュ点(L4, L5)

$y \neq 0$ の場合、平衡条件は $r_1 = r_2 = 1$($m_1$ と $m_2$ と等距離)のときに満たされることが示せます。

$x$ 方向の平衡条件に $r_1 = r_2 = 1$ を代入すると、$x$ が直接求まります。$y$ 方向の条件を合わせると

$$ L4 = \left(\frac{1}{2}-\mu, \, \frac{\sqrt{3}}{2}, \, 0\right), \quad L5 = \left(\frac{1}{2}-\mu, \, -\frac{\sqrt{3}}{2}, \, 0\right) $$

L4とL5は、$m_1$ と $m_2$ を底辺とする正三角形の頂点に位置します。これがラグランジュが1772年に発見した「三角形解」です。驚くべきことに、L4とL5の位置は質量比 $\mu$ にほとんど依存しません。

木星のL4とL5にはトロヤ群小惑星と呼ばれる数千個の小惑星が集まっており、ラグランジュの理論的予測を壮大なスケールで実証しています。

5つのラグランジュ点の位置が求まりました。しかし、位置がわかっただけでは十分ではありません。衛星をそこに置いたとき、微小な摂動に対して元に戻るのか(安定)、それとも離れていくのか(不安定)を調べる必要があります。

安定性解析

線形化

ラグランジュ点 $(x_0, y_0)$ のまわりでの微小変位 $(\xi, \eta)$($x = x_0 + \xi$, $y = y_0 + \eta$)に対して、運動方程式を線形化します。

$$ \ddot{\xi} – 2\dot{\eta} = U_{xx}\xi + U_{xy}\eta $$

$$ \ddot{\eta} + 2\dot{\xi} = U_{xy}\xi + U_{yy}\eta $$

ここで $U_{xx}$, $U_{xy}$, $U_{yy}$ はラグランジュ点での有効ポテンシャルの2階偏微分です。

解を $\xi = A e^{\lambda t}$, $\eta = B e^{\lambda t}$ とおくと、特性方程式

$$ \lambda^4 + (4 – U_{xx} – U_{yy})\lambda^2 + U_{xx}U_{yy} – U_{xy}^2 = 0 $$

が得られます。これは $\lambda^2$ の二次方程式で、2つの $\lambda^2$ の根の符号が安定性を決定します。

共線ラグランジュ点(L1, L2, L3)の不安定性

共線点($y_0 = 0$)では $U_{xy} = 0$ で、$U_{xx} > 0$, $U_{yy} < 0$ となります。特性方程式の根は

$$ \lambda^2 = \frac{-(4 – U_{xx} – U_{yy}) \pm \sqrt{(4 – U_{xx} – U_{yy})^2 – 4(U_{xx}U_{yy})}}{2} $$

一方の $\lambda^2$ が正になるため、実数の正の $\lambda$ が存在し、共線ラグランジュ点は不安定です。微小な摂動が指数関数的に増大します。

不安定の時間スケール(e-folding time)はL1, L2で約23日(太陽-地球系)です。したがって、JWSTのような衛星はステーションキーピング噴射(年間数m/sのΔV)で位置を維持しています。

三角形ラグランジュ点(L4, L5)の条件付き安定性

三角形点では $U_{xy} \neq 0$ で、特性方程式のすべての $\lambda^2$ が負($\lambda$ が純虚数)になる条件は

$$ \mu < \mu_{\text{crit}} = \frac{1}{2}\left(1 - \sqrt{\frac{23}{27}}\right) \approx 0.03852 $$

この条件が満たされると、L4とL5は線形安定です。太陽-木星系($\mu \approx 0.000954$)や太陽-地球系($\mu \approx 3 \times 10^{-6}$)ではこの条件が満たされるため、トロヤ群小惑星が長期間安定に存在できます。

安定性の物理的な直感を述べると、L4/L5が「ポテンシャルの谷」ではなく「丘の上」にあるにもかかわらず安定なのは、コリオリ力のおかげです。微小な摂動を受けた物体はL4/L5のまわりをゆっくり回る「レトログレード運動」をし、ポテンシャルの丘から落ちることなく平衡点の周りにとどまります。これはコリオリ力がなければ起こりえない現象です。

安定性の理論が整ったところで、L1やL2付近の実用的な軌道(ハロー軌道)について見ていきましょう。

ハロー軌道とリサージュ軌道

L点付近の周期軌道

不安定なL1やL2は「そこに静止する」ことは(摂動なしでも長期的には)できませんが、その近傍には周期軌道が存在します。線形化した運動方程式の解のうち、振動解の部分がこの周期軌道に対応します。

リサージュ軌道(Lissajous orbit): 回転座標系のL点まわりでの面内振動と面外振動の周期が異なるため、衛星は一般にリサージュ図形のような軌跡を描きます。

ハロー軌道(Halo orbit): 面内振動と面外振動の周期が一致する特殊な周期軌道です。回転座標系で見ると、L点のまわりのハロー(光輪)のような閉じた軌道を描きます。ハロー軌道は非線形効果によって存在し、線形解析では予測できません。

JWSTは太陽-地球L2のまわりのハロー軌道を周回しています。ハロー軌道のサイズは半径約100万kmで、6ヶ月の周期で一周します。L2に正確に位置するのではなく、そのまわりを大きく周回することで地球の影に入ることを避けています。

ステーションキーピング

L1/L2のハロー軌道は不安定であるため、定期的な軌道修正が必要です。典型的なステーションキーピングの $\Delta V$ は年間1〜5 m/sと非常に小さいです。これはGEO衛星の南北ステーションキーピング(年間40〜50 m/s)に比べて格段に少なく、L点ミッションの燃料効率が良い理由の一つです。

ハロー軌道の理論的背景を理解したところで、次にPythonでゼロベロシティ曲面とラグランジュ点を可視化してみましょう。

Pythonによる可視化

ラグランジュ点の数値計算

まず、5つのラグランジュ点の位置を数値的に求めます。

import numpy as np
from scipy.optimize import brentq, fsolve
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm

# --- 質量パラメータ ---
mu = 3.003e-6  # 太陽-地球系


def effective_potential(x, y, mu_val):
    """有効ポテンシャル U(x,y)"""
    r1 = np.sqrt((x + mu_val)**2 + y**2)
    r2 = np.sqrt((x - 1 + mu_val)**2 + y**2)
    return 0.5 * (x**2 + y**2) + (1 - mu_val) / r1 + mu_val / r2


def dUdx_on_axis(x, mu_val):
    """x軸上での ∂U/∂x"""
    r1 = abs(x + mu_val)
    r2 = abs(x - 1 + mu_val)
    return x - (1 - mu_val) * (x + mu_val) / r1**3 - mu_val * (x - 1 + mu_val) / r2**3


# --- 共線ラグランジュ点の計算 ---
# L1: m1とm2の間(1-mu_val より少し内側)
x_L1 = brentq(lambda x: dUdx_on_axis(x, mu), 0.5, 1 - mu - 1e-10)

# L2: m2の外側
x_L2 = brentq(lambda x: dUdx_on_axis(x, mu), 1 - mu + 1e-10, 1.5)

# L3: m1の反対側
x_L3 = brentq(lambda x: dUdx_on_axis(x, mu), -1.5, -mu - 1e-10)

# --- 三角形ラグランジュ点 ---
x_L4 = 0.5 - mu
y_L4 = np.sqrt(3) / 2

x_L5 = 0.5 - mu
y_L5 = -np.sqrt(3) / 2

print("=== ラグランジュ点の位置(無次元)===")
print(f"L1: x = {x_L1:.8f}")
print(f"L2: x = {x_L2:.8f}")
print(f"L3: x = {x_L3:.8f}")
print(f"L4: x = {x_L4:.8f}, y = {y_L4:.8f}")
print(f"L5: x = {x_L5:.8f}, y = {y_L5:.8f}")

# 物理的な距離に変換
AU = 1.496e8  # 1 AU [km]
print(f"\n=== 物理的な距離 ===")
print(f"L1: 地球から {(1-mu-x_L1)*AU:.0f} km 内側")
print(f"L2: 地球から {(x_L2-1+mu)*AU:.0f} km 外側")

ゼロベロシティ曲面の可視化

# --- ゼロベロシティ曲面の可視化 ---
fig, axes = plt.subplots(1, 2, figsize=(16, 7))

# 広域図
x_range = np.linspace(-1.5, 1.5, 800)
y_range = np.linspace(-1.5, 1.5, 800)
X, Y = np.meshgrid(x_range, y_range)

# 有効ポテンシャルの計算
U = np.zeros_like(X)
for i in range(X.shape[0]):
    for j in range(X.shape[1]):
        U[i, j] = effective_potential(X[i, j], Y[i, j], mu)

# 各ラグランジュ点でのヤコビ定数
CJ_L1 = 2 * effective_potential(x_L1, 0, mu)
CJ_L2 = 2 * effective_potential(x_L2, 0, mu)
CJ_L3 = 2 * effective_potential(x_L3, 0, mu)
CJ_L4 = 2 * effective_potential(x_L4, y_L4, mu)

# 等高線プロット
ax = axes[0]
levels = np.linspace(CJ_L4 * 0.98, CJ_L1 * 1.02, 30)
contour = ax.contourf(X, Y, 2 * U, levels=levels, cmap='RdYlBu_r', alpha=0.8)
ax.contour(X, Y, 2 * U, levels=[CJ_L1, CJ_L2, CJ_L3, CJ_L4],
           colors=['red', 'orange', 'yellow', 'lime'], linewidths=1.5)

# ラグランジュ点のプロット
L_points = [(x_L1, 0, 'L1'), (x_L2, 0, 'L2'), (x_L3, 0, 'L3'),
            (x_L4, y_L4, 'L4'), (x_L5, y_L5, 'L5')]
for xp, yp, label in L_points:
    ax.plot(xp, yp, 'k^', markersize=8)
    ax.annotate(label, (xp, yp), textcoords="offset points",
                xytext=(10, 5), fontsize=11, fontweight='bold')

# 主天体
ax.plot(-mu, 0, 'yo', markersize=12, label='Sun (m₁)')
ax.plot(1 - mu, 0, 'bo', markersize=8, label='Earth (m₂)')

ax.set_xlabel('x [L]')
ax.set_ylabel('y [L]')
ax.set_title('ゼロベロシティ曲面(広域図)')
ax.set_aspect('equal')
ax.legend(loc='lower right')

# L1, L2付近の拡大図
ax2 = axes[1]
x_range2 = np.linspace(0.95, 1.05, 500)
y_range2 = np.linspace(-0.05, 0.05, 500)
X2, Y2 = np.meshgrid(x_range2, y_range2)

U2 = np.zeros_like(X2)
for i in range(X2.shape[0]):
    for j in range(X2.shape[1]):
        U2[i, j] = effective_potential(X2[i, j], Y2[i, j], mu)

contour2 = ax2.contourf(X2, Y2, 2 * U2, levels=30, cmap='RdYlBu_r', alpha=0.8)
ax2.contour(X2, Y2, 2 * U2, levels=[CJ_L1, CJ_L2],
            colors=['red', 'orange'], linewidths=2)

for xp, yp, label in [(x_L1, 0, 'L1'), (x_L2, 0, 'L2')]:
    ax2.plot(xp, yp, 'k^', markersize=10)
    ax2.annotate(label, (xp, yp), textcoords="offset points",
                 xytext=(8, 5), fontsize=12, fontweight='bold')

ax2.plot(1 - mu, 0, 'bo', markersize=10, label='Earth')
ax2.set_xlabel('x [L]')
ax2.set_ylabel('y [L]')
ax2.set_title('L1-L2付近の拡大図')
ax2.set_aspect('equal')
ax2.legend()

plt.tight_layout()
plt.savefig('lagrange_points_zvs.png', dpi=150, bbox_inches='tight')
plt.show()

この図から、ラグランジュ点とゼロベロシティ曲面の構造が視覚的に理解できます。

  1. 広域図(左): 5つのラグランジュ点がすべて表示されています。L4とL5が太陽と地球から等距離にある正三角形の頂点に位置していることが確認できます。等高線(ゼロベロシティ曲面の断面)は、有効ポテンシャルの地形図として読めます。

  2. 拡大図(右): L1とL2が地球の両側に対称に近い位置にあることがわかります。L1とL2を通るゼロベロシティ曲面の等高線(赤と橙の線)は、地球の周辺で開口する「ネック」を形成しています。ヤコビ定数がL1の臨界値を超えると、このネックが開いて内側(地球側)と外側を行き来できるようになります。

CR3BPでの軌道のシミュレーション

L点付近の軌道の振る舞いを見るため、CR3BPの運動方程式を数値積分してリサージュ軌道を計算します。

from scipy.integrate import solve_ivp

def cr3bp_eom(t, state, mu_val):
    """円制限三体問題の運動方程式"""
    x, y, z, vx, vy, vz = state

    r1 = np.sqrt((x + mu_val)**2 + y**2 + z**2)
    r2 = np.sqrt((x - 1 + mu_val)**2 + y**2 + z**2)

    ax = 2*vy + x - (1-mu_val)*(x+mu_val)/r1**3 - mu_val*(x-1+mu_val)/r2**3
    ay = -2*vx + y - (1-mu_val)*y/r1**3 - mu_val*y/r2**3
    az = -(1-mu_val)*z/r1**3 - mu_val*z/r2**3

    return [vx, vy, vz, ax, ay, az]

# L2付近のリサージュ軌道の初期条件
# L2点からわずかにずれた位置・速度
x0 = x_L2 + 0.001
y0 = 0.0
z0 = 0.001
vx0 = 0.0
vy0 = -0.01
vz0 = 0.0

state0 = [x0, y0, z0, vx0, vy0, vz0]

# 無次元時間で10単位(約10/(2π)公転周期)
t_span = (0, 10)
t_eval = np.linspace(0, 10, 10000)

sol = solve_ivp(lambda t, s: cr3bp_eom(t, s, mu), t_span, state0,
                method='DOP853', t_eval=t_eval,
                rtol=1e-12, atol=1e-14)

# ヤコビ定数の確認
CJ_traj = np.zeros(len(sol.t))
for k in range(len(sol.t)):
    x_k, y_k, z_k = sol.y[0, k], sol.y[1, k], sol.y[2, k]
    vx_k, vy_k, vz_k = sol.y[3, k], sol.y[4, k], sol.y[5, k]
    U_k = effective_potential(x_k, y_k, mu)
    # z成分を含む有効ポテンシャルの修正
    r1_k = np.sqrt((x_k+mu)**2 + y_k**2 + z_k**2)
    r2_k = np.sqrt((x_k-1+mu)**2 + y_k**2 + z_k**2)
    U_k_full = 0.5*(x_k**2 + y_k**2) + (1-mu)/r1_k + mu/r2_k
    v2_k = vx_k**2 + vy_k**2 + vz_k**2
    CJ_traj[k] = 2*U_k_full - v2_k
fig = plt.figure(figsize=(14, 5))

# xy平面の軌跡
ax1 = fig.add_subplot(131)
ax1.plot(sol.y[0], sol.y[1], 'b-', linewidth=0.3, alpha=0.7)
ax1.plot(x_L2, 0, 'r^', markersize=10, label='L2')
ax1.plot(1-mu, 0, 'go', markersize=8, label='Earth')
ax1.set_xlabel('x [L]')
ax1.set_ylabel('y [L]')
ax1.set_title('L2付近の軌道(xy平面)')
ax1.set_aspect('equal')
ax1.legend()
ax1.grid(True, alpha=0.3)

# xz平面の軌跡
ax2 = fig.add_subplot(132)
ax2.plot(sol.y[0], sol.y[2], 'b-', linewidth=0.3, alpha=0.7)
ax2.plot(x_L2, 0, 'r^', markersize=10, label='L2')
ax2.set_xlabel('x [L]')
ax2.set_ylabel('z [L]')
ax2.set_title('L2付近の軌道(xz平面)')
ax2.set_aspect('equal')
ax2.legend()
ax2.grid(True, alpha=0.3)

# ヤコビ定数の保存
ax3 = fig.add_subplot(133)
ax3.plot(sol.t, (CJ_traj - CJ_traj[0]) / CJ_traj[0], 'b-', linewidth=0.5)
ax3.set_xlabel('Time [non-dim]')
ax3.set_ylabel('ΔC_J / C_J')
ax3.set_title('ヤコビ定数の保存(相対誤差)')
ax3.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('l2_orbit_simulation.png', dpi=150, bbox_inches='tight')
plt.show()

この結果から、L2付近の軌道の特徴が読み取れます。

  1. xy平面(左図): L2付近で複雑な軌跡を描いています。初期条件によっては、L2のまわりを周回するリサージュ軌道になることもあれば、L2から離れていく不安定軌道になることもあります。これはL2の不安定性を反映しています。

  2. xz平面(中央図): 面外方向($z$)の振動が確認できます。面内振動と面外振動の周期が一般に異なるため、リサージュ図形のパターンが現れます。

  3. ヤコビ定数の保存(右図): ヤコビ定数の相対誤差が $10^{-12}$ 以下に保たれており、数値積分の精度が十分であることが確認できます。CR3BPでのヤコビ定数は理論的に厳密な保存量なので、これを監視することで積分の品質を検証できます。

JWST(L2)ミッションの実例

なぜL2が選ばれたか

JWSTが太陽-地球系のL2に配置された理由は主に3つあります。

熱環境の安定性: L2では太陽、地球、月がすべて同じ方向にあるため、大型のサンシールドで一度に遮蔽できます。赤外線望遠鏡であるJWSTは検出器を-233°C以下に冷却する必要があり、熱源を一方向にまとめられるL2は理想的な場所です。

連続観測: 地球周回軌道では地球の影に入ったり出たりするため、観測が中断されます。L2のハロー軌道では常に太陽光が当たり(かつサンシールドで遮蔽し)、連続的な観測が可能です。

通信: L2は地球から約150万km離れていますが、常に地球と同じ側にあるため、通信リンクが安定しています。

JWSTの軌道

JWSTは太陽-地球L2のまわりの大型ハロー軌道を周回しています。軌道のサイズは約100万km×50万kmで、周期は約6ヶ月です。ステーションキーピングの $\Delta V$ は年間約2 m/sと非常に小さく、搭載燃料でミッション寿命10年(目標20年以上)の運用が可能です。

実際のJWSTの軌道投入は、2022年1月24日にL2ハロー軌道への最終軌道修正噴射(MCC-2)が行われ、成功しました。打ち上げロケット(アリアン5)の精度が非常に高かったため、計画より少ない燃料消費で軌道投入でき、ミッション寿命が当初の10年から20年以上に延長される見込みとなりました。

まとめ

本記事では、ラグランジュ点の理論的な導出から実際のミッションへの応用まで解説しました。

  • 制限三体問題: 2つの主天体が円軌道で公転する系に質量ゼロの第3天体を置く枠組み。回転座標系で記述すると平衡点が自然に現れる
  • ヤコビ積分: CR3BPの唯一の保存量。ゼロベロシティ曲面を通じて衛星の運動可能領域を制限する
  • 5つのラグランジュ点: 共線点L1-L3($x$軸上)と三角形点L4-L5(正三角形の頂点)。L1/L2は主天体から $(\mu/3)^{1/3}L$ の距離
  • 安定性: L1-L3は不安定(指数的に発散)。L4-L5は $\mu < 0.0385$ で線形安定(コリオリ力による安定化)
  • ハロー軌道: L1/L2付近の周期軌道。不安定だがステーションキーピングで維持可能(年間1〜5 m/s)
  • JWST: 太陽-地球L2のハロー軌道に配置。熱環境、連続観測、通信の観点から最適な位置

次のステップとして、以下の記事も参考にしてください。