はじめに

下記の記事で一太陽年や一朔望月などの天体物理の基礎用語を調べた。

このうち地球が太陽の引力（潮汐力）によって地軸の向きが緩やかに変わること，すなわち歳差運動により北極星とされる星（現在はこぐま座のα星）が移り変わっていくことはもともと知っていた。しかし一太陽年は地球が太陽の周りを回って元の位置に戻る手前で定義されていること，すなわち一太陽年後に地球は元の位置には戻らないことを知って驚いた。地球の自転軸と太陽の向きの関係が変わらないだけで，地球の位置は少しずつ変わっていくのだ。この通りなら冬や夏の星座と呼ばれているものも2万6000年周期で入れ替わることになる・・・と思って調べてみたらその通りだった。いわゆる知り得た情報が知識として身に付いていなかったということだろう。

さて，今回はこの歳差周期2万6000年を高校レベルの数学で求めてみたい。

本記事では慣性モーメントや角運動量など剛体の物理を取り扱っており，内容自体は大学レベルだと思うが，あくまで使用している数学は高校レベルという意味である。

潮汐力を受ける部分の質量

地球は完全球体ではなく，自転による遠心力により赤道方向に僅かに膨らんだ回転楕円体となっている。地球は太陽および月などによって潮汐力を受けるが，地球が完全な球体ではないことにより，この潮汐力が地軸の傾きを小さくする方向へのトルクとして作用する。

簡単のため地球を完全球体とみなし，膨らんだ部分の余剰質量が細いリングとして赤道上に集約しているものと考える。リングの線密度 $\mu$ とおくと余剰質量の総和 $\Delta M$ は地球の半径 $r_p$ として

\Delta M = 2 \pi \cdot r_p \cdot \mu \tag{1}

となる。一方，地球の極半径 $r_p$，赤道半径 $r_e$ とおくと膨らんだ部分の余剰体積 $\Delta V$ は

\begin{align}
\Delta V &= \frac{4}{3}\pi \cdot r_p \cdot r_e^2 - \frac{4}{3}\pi \cdot r_p^3 \\
&= \frac{4}{3} \pi \cdot r_p \cdot ( r_e^2 - r_p^2)  \tag{2}
\end{align}

となる。膨らんだ部分は地球の中でも最外層の軽い部分であるから地殻密度 $\rho$ とおくと線密度 $\mu$ は

\begin{align}
\Delta M &= \rho \cdot \Delta V \\
2 \pi \cdot r_p \cdot \mu &= \frac{4}{3} \pi \cdot \rho \cdot r_p \cdot ( r_e^2 - r_p^2) \\
\mu &= \frac{2}{3} \cdot \rho \cdot (r_e^2 - r_p^2) \tag{3}
\end{align}

となる。

地球の受けるトルク

簡単のため傾いた地軸をZ軸にとる。地球と太陽の距離，すなわち地球の公転軌道半径 $R_e$ とおき，太陽はXZ平面上の座標 $(R_e \cos \varepsilon, 0, R_e \sin \varepsilon)$ に位置するものとする。ただし，太陽は十分遠方にあるものとし，地球上において受ける潮汐力は平行とする。

球対称の部分はトルクが相殺されるため無視してよい。このため非対称部分，すなわち赤道上に置いたリングに作用するトルクについてのみ考えればよい。

赤道リングの経度 $\phi$ において潮汐力を受ける部分の質量 $dm$ を示す。

dm = \mu \cdot r_p \cdot d\phi \tag{4}

太陽から受ける潮汐力 $df_s$ は地球の中心を通り，かつ太陽方向と鉛直な平面までの距離 $r_p \cos \phi \cos \varepsilon$ に比例する。

df_s = \frac{2GM_s}{R_e^3} \cdot dm \cdot r_p \cos \phi \cos \varepsilon \tag{5}

太陽から受けるトルク $d\tau_s$ は潮汐力 $df_s$ のうち赤道面と垂直な成分 $df_s \sin \varepsilon$ とY軸までの距離 $r_p\cos\phi$ の積である。

\begin{align}
d\tau_s &= df_s \cdot r_p \cos \phi \sin \varepsilon \\
&= \frac{GM_s}{R_e^3} \cdot dm \cdot r_p^2  \cos^2 \phi \sin 2\varepsilon \\
&= \frac{GM_s}{R_e^3} \cdot \mu \cdot r_p^3  \cos^2 \phi \sin 2\varepsilon \cdot d\phi \tag{6}
\end{align}

地球全体が太陽から受けるトルク $\tau_s$ は，これを積分したものなので

\begin{align}
\tau_s &= \int d\tau_s \\
&= \int_{-\pi}^{\pi} \frac{GM_s}{R_e^3} \cdot \mu \cdot r_p^3 \cos^2 \phi \sin 2 \varepsilon \cdot d\phi \\
&= \frac{GM_s}{R_e^3} \cdot \mu \cdot r_p^3 \sin 2 \varepsilon \int_{-\pi}^{\pi} \cos^2 \phi \cdot d\phi \\
&= \frac{GM_s}{R_e^3} \cdot \mu \cdot \pi \cdot r_p^3 \sin 2 \varepsilon \\
&= \frac{2GM_s}{3R_e^3} \cdot \rho \cdot \pi \cdot r_p^3 \cdot (r_e^2 - r_p^2) \sin 2 \varepsilon \tag{7}
\end{align}

となる。同様にして地球が月から受けるトルク $\tau_m$ を求めると

\tau_m = \frac{2GM_m}{3R_m^3} \cdot \rho \cdot \pi \cdot r_p^3 \cdot (r_e^2 - r_p^2) \sin 2\varepsilon \tag{8}

となる。表１より，それぞれのトルクを求めると

\left \{ \begin{aligned}
\tau_s &= 1.147 \times 10^{37} \\
\tau_m &= 2.421 \times 10^{38}
\end{aligned} \right. \tag{9}

となる。

潮汐力（トルク）に関しては太陽より月の影響のほうが２倍以上大きい。

表１　定数表
項目	記号	数値	単位
地球の赤道半径	$r_e$	$6,378,178$	[m]
地球の極半径	$r_p$	$6,356,752$	[m]
太陽質量	$M_s$	$1.989 \times 10^{30}$	[kg]
月質量	$M_m$	$7.346 \times 10^{22}$	[kg]
重力定数	$G$	$6.674 \times 10^{-11}$	[m³･kg⁻¹･s⁻²]
地球公転半径	$R_e$	$1.496 \times 10^{11}$	[m]
月公転半径	$R_m$	$3.884 \times 10^8$	[m]
地球の平均地殻密度	$\rho$	$2.7 \times 10^3$	[kg/m³]
地球の地軸の傾き	$\varepsilon$	$23.44$	[deg]

年間を通じた平均トルク

先ほど求めたのは（敢えて言明しなかったが）冬至のとき地球が受ける最大トルクである。春分・秋分のときは太陽が $(0,\pm R_e,0)$ の位置にあると考えると地球の受けるトルクはゼロとなる。一方，夏至のときは太陽が $(-R_e\cos(-\varepsilon),0,-R_e\sin(-\varepsilon))$ の位置にあると考えるとトルクは冬至のときと等しく最大値を取る。

最大値と最小値の間を直線で結ぶのかスムーズに正弦波状に繋げるべきかは不問として，年間を通じた平均トルクは夏至・冬至のとき（最大値）のちょうど半分になると考えられる。月も同様に考えると年間を通じた平均トルク $\bar{\tau}$ は

\bar{\tau} = \frac{1}{2} ( \tau_s + \tau_m ) \tag{10}

となると考えられる。

慣性モーメントを求める

簡単のため地球内部の密度を一様とすると自転の慣性モーメント $I_e$ は

I_e = \frac{2}{5} M_e \cdot r_p^2 = 9.653 \times 10^{37} \tag{11}

となるが，さすがに誤差が大きい。※実際の値に比べて２割ほど大きい。

表２　定数表
項目	記号	数値	単位
地球の質量	$M_e$	$5.9724 \times 10^{24}$	[kg]

地球の内部構造は詳しく分かっているので精密に求めてみよう。地球の内部密度は以下のように地球の中心に近づくほど高くなり，層の境界で不連続点を持つ。

参考文献[11]ではもっと細かく層が分かれているが，簡単のため６層にまとめた。また三次式でモデル化しているところ，境界値の間を単純に線形補間（一次近似）した。

まず半径 $r$，密度 $\rho(r)$ の球殻の慣性モーメント $I(r)$ を求める。

dI = \rho(r) \cdot r^2 \cos \theta \cdot d\theta \cdot d\phi \cdot (r \cos \theta)^2 \tag{12}

\begin{align}
I(r) &= \iint dI \\
&= \rho(r) \cdot r^4 \int_{-\pi/2}^{\pi/2} \cos^3 \theta \cdot d\theta \int_{-\pi}^{\pi} d\phi \\
&= \frac{8\pi}{3} \cdot \rho(r) \cdot r^4 \tag{13}
\end{align}

地球中心からの距離 $r_0 \le r \le r_1$ において密度関数 $\rho(r)$ の境界値が次のように与えられたとする。

\left\{ \begin{aligned}
\rho(r_0) &= \rho0 \\
\rho(r_1) &= \rho1 \\
\end{aligned} \right. \tag{14}

区間 $r_0 \le r \le r_1$ の間，密度は線形補間するものとする。

\begin{align}
\rho(r) &= \rho_0 + \frac{\rho_1 - \rho_0}{r_1 - r_0} \cdot (r - r_0) \\
&= \frac{\rho_0 \cdot r_1 - \rho_1 \cdot r_0}{r_1 - r_0} + \frac{\rho_1 - \rho_0}{r_1 - r_0} \cdot r \tag{15}
\end{align}

区間 $r_0 \le r \le r_1$ の慣性モーメント $I[r_0, r_1]$ を求めると

\begin{align}
I[r_0, r_1] &= \int_{r_0}^{r_1} I(r) \cdot dr \\
&= \int_{r_0}^{r_1} \frac{8\pi}{3} \cdot \rho(r) \cdot r^4 \cdot dr \\
&= \frac{8\pi}{3} \int_{r_0}^{r_1} \left( \frac{\rho_0 \cdot r_1 - \rho_1 \cdot r_0}{r_1 - r_0} + \frac{\rho_1 - \rho_0}{r_1 - r_0} \cdot r \right)  \cdot r^4 \cdot dr \\
&= \frac{8\pi}{3} \left[ \frac{\rho_0 \cdot r_1 - \rho_1 \cdot r_0}{r_1 - r_0} \cdot \frac{r^5}{5} + \frac{\rho_1 - \rho_0}{r_1 - r_0} \cdot \frac{r^6}{6} \right]_{r_0}^{r_1} \\
&= \frac{8\pi}{3} \left\{ \frac{\rho_0 \cdot r_1 - \rho_1 \cdot r_0}{r_1 - r_0} \cdot \frac{r_1^5 - r_0^5}{5} + \frac{\rho_1 - \rho_0}{r_1 - r_0} \cdot \frac{r_1^6 - r_0^6}{6} \right\} \tag{16}
\end{align}

となる。こうして各層ごとの慣性モーメントを求め，これらの総和により地球全体の慣性モーメントを求めたところ $I_e = 8.0359 \times 10^{37}$ を得られた。

表３　地球内部密度と慣性モーメント
名称	中心からの距離 [km]	密度 [g/cm³]	慣性モーメント [kg･m²]
内核	$0 \rightarrow 1,222$	$13.089 \rightarrow 12.764$	$5.8402 \times 10^{34}$
外核	$1,222 \rightarrow 3,480$	$12.166 \rightarrow\:\: 9.903$	$8.9096 \times 10^{36}$
下部マントル	$3,480 \rightarrow 5,701$	$5.566 \rightarrow\:\: 4.381$	$4.4296 \times 10^{37}$
上部マントル	$5,701 \rightarrow 6,347$	$3.992 \rightarrow\:\: 3.381$	$2.6248 \times 10^{36}$
地殻	$6,347 \rightarrow 6,368$	$2.900 \rightarrow\:\: 2.600$	$8.0520 \times 10^{35}$
海洋	$6,368 \rightarrow 6,371$	$1.020$	$4.2195 \times 10^{34}$
合計			$8.0359 \times 10^{37}$

同様にして地球の質量も計算できる。区間 $r_0 \le r \le r_1$ の質量 $M[r_0, r_1]$ は

M[r_0, r_1] = 4\pi \left\{ \frac{\rho_0 \cdot r_1 - \rho_1 \cdot r_0}{r_1 - r_0} \cdot \frac{r_1^3 - r_0^3}{3} + \frac{\rho_1 - \rho_0}{r_1 - r_0} \cdot \frac{r_1^4 - r_0^4}{4} \right\} \tag{17}

であるから，各層ごとの質量を求め，これらの総和により地球全体の質量を求めたところ $M_e = 5.9455 \times 10^{24} $ [kg] を得られた。これは表２に示す真値の 99.55% に当たる。

歳差周期を求める

地球の自転周期 $T_r = 2\pi / \omega_r$，自転角速度 $\omega_r$，慣性モーメント $I_e$ とおくと地球の自転による角運動量ベクトルの大きさ $L = I_e \cdot \omega_r$ となり，歳差運動による角速度 $\omega_p$ とおくと角運動量ベクトルの時間変化率が地球に加わるトルクと等しいので

\omega_p \cdot L \sin \varepsilon = \bar{\tau} \tag{18}

となる。

あくまで角運動量ベクトルの大きさ $|\vec{L}| = L$ は不変であり，作用するトルクによってベクトルの向きが変わるだけと理解されたい。

歳差運動の角速度 $\omega_p$ は

\omega_p = \frac{\bar{\tau}}{L \sin \varepsilon} = \frac{1}{2} \cdot \frac{\tau_s + \tau_m}{\omega_r \cdot I_e \sin \varepsilon} \tag{19}

より $\omega_p = 7.653 \times 10^{-12}$ となり，歳差運動の周期 $T_p = 2\pi / \omega_r = 8.210 \times 10^{11}$ が得られる。

\begin{aligned}
\frac{8.210 \times 10^{11}\textsf{［秒］}}{24 \times 60^2 \textsf{［秒/日］}} &= 9.503 \times 10^6 \textsf{［日］} \\
\frac{9.503 \times 10^6 \textsf{［日］}}{365.25 \textsf{［日/年］}} &= 26,017 \textsf{［年］} 
\end{aligned} \tag{20}

表４　定数表
項目	記号	数値	単位	備考
地球自転周期	$T_r$	$86,164$	[s]	23時間56分04秒

まとめ

以上より歳差周期 2万6017年を得ることができた。真値2万5772年と比べると素晴らしい精度（+1.0%）で得られている。

このうち赤道リングモデルについては，地球が受けるトルクを（奇跡的に）ほぼ正しく計算することができたようだ。

地球の慣性モーメントについては地球の中心からの距離と密度のみで計算していることから球対称であり，冒頭で述べたように地球は赤道方向に少し膨らんだ回転楕円体であることを考慮していない。これに加えて粗い線形近似モデルにも関わらず，こちらもまた驚異的な精度を得られている。

参考文献

[1] 国立天文台＞暦計算室＞暦Wiki＞惑星　太陽系天体の半径
[2] 地球内部物理学 - wikipedia
[3] 地球の地殻 - wikipedia
[4] 地球の構造 - wikipedia
[5] 重力定数 - wikipedia
[6] 地球 - wikipedia
[7] 太陽 - wikipedia
[8] 月 - wikipedia
[9] 歳差 - wikipedia
[10] 地圏－層構造 - 名古屋市科学館
[11] Preliminary reference Earth model, Physics of the Earth and Planetary Interiors, 25 (1981)