地球の大きさが無視できなすぎる

問題

地球や惑星、月などの軌道計算をする時、通常は、これらの星の大きさは無視して、質点として扱う。真面目に考えると、これには、結構色々な近似が入ってる。

重力を受ける側として見ると、固体部分の弾性変形を無視しているし、更に、剛体として扱っても位置と姿勢の自由度が残るけど、質点近似では、位置の自由度だけ取り出してることになる。量子統計力学で、分子を扱う時、並進自由度、回転自由度、内部振動の自由度を独立に扱うのと似てなくもない。この3種類の自由度には、相互作用があるが、それを無視している。

一方、重力を作る側としては、質量分布が球対称であれば、星の外部の重力は、質量中心に全質量が集まってると考えて計算して問題ない。これはニュートンが証明したらしい。多分、プリンキピアのどこかに書いてあるのだろう。

素晴らしい定理に思えるが、地球の形状は、別に球でも何でもない。表面は凸凹してるし、大陸の分布とか眺めても、明らかに球対称でも軸対称でもない。そうすると、球対称近似には、誤差があるはずで、どれくらいあるか評価できないと、質点近似を行っていいのか判断できない。

歴史的には、そういう評価ができるようになる前に、軌道計算をやって、測定と合い、実用的に十分だからヨシとして進んできた。そういう適当さが重要な時もあるが、いつまでも18世紀気分で勉強するのも、よくなかろう。ここでは特に、地球が球対称じゃないことが、重力に、どう影響するか評価することを考える。

重力ポテンシャルの多重極展開

どの教科書にも書いてあると思うけど、地球の(剛体近似下での)質量分布$${\rho(\mathbf{y})}$$が与えられた時、地球が作る重力ポテンシャル$${V(\mathbf{x})}$$は

$$
V(\mathbf{x}) = -G \displaystyle \int_{\mathrm{Earth}} \dfrac{\rho(\mathbf{y})}{|\mathbf{x} - \mathbf{y}|} d \mathbf{y}
$$

で与えられる。しばしば暗黙の内に、この積分は、以下の式に置き換えられている。

$$
M = \displaystyle\int_{\mathrm{Earth}} \rho(\mathbf{y}) d \mathbf{y}
$$

$$
V(\mathbf{x}) = -\dfrac{G}{ |\mathbf{x} - \mathbf{y_{0}}| } M
$$

$$
\mathbf{y_{0}} = \dfrac{1}{M} \displaystyle \int_{\mathrm{Earth}} \mathbf{y} \rho(\mathbf{y}) d \mathbf{y}
$$

説明なしに見ると、数学的に飛躍があると思うので、導出手順を踏む。

$$
\dfrac{1}{|\mathbf{x} - \mathbf{y}|} = \dfrac{1}{\sqrt{|\mathbf{x}|^2 - 2 \mathbf{x} \cdot \mathbf{y} + |\mathbf{y}|^2}} = \dfrac{1}{|\mathbf{x}|} \dfrac{1}{\sqrt{1 - 2 \frac{|\mathbf{y}|}{|\mathbf{x}|} \cos \theta + (\frac{|\mathbf{y}|}{|\mathbf{x}|})^2}}
$$

という形になる。但し、$${\cos \theta = \dfrac{\mathbf{x} \cdot \mathbf{y}}{ | \mathbf{x} | | \mathbf{y} | }}$$

$${\mathbf{y}}$$は地球内部を走り、$${\mathbf{x}}$$は、地球外部にあるけど、地球から十分遠くにあれば、$${|\mathbf{y}| / |\mathbf{x}| < 1}$$となる。座標を、地球の質量中心$${\mathbf{y_{0}}}$$が原点に来るように取れば、$${|\mathbf{x}|}$$が、"地球半径"より大きければいい。

$${a = \cos \theta}$$,$${z=|\mathbf{y}| / |\mathbf{x}|}$$として、$${\dfrac{1}{\sqrt{1-2a z +z^2}}}$$を$${z=0}$$の周りでテイラー展開する。

手で計算してもいいが面倒なので、Wolfram Alphaに答えを聞くと

$$
\dfrac{1}{\sqrt{1-2a z +z^2}} = 1 + a z + ((-1 + 3 a^2) z^2)/2 + (a (-3 + 5 a^2) z^3)/2 + ((3 - 30 a^2 + 35 a^4) z^4)/8 + (a (15 - 70 a^2 + 63 a^4) z^5)/8 + O(z^6)
$$

この展開は、$${|a| \leq 1 , |z| \lt 1}$$で収束する(ほんとは証明すべきだが略)

 

以上から、

$$
V(\mathbf{x}) = - \dfrac{G}{|\mathbf{x}|} \displaystyle \int_{\mathrm{Earth}} \left( 1 + a z + \dfrac{3a^2 -1}{2} z^2 + \cdots \right) \rho(\mathbf{y}) d\mathbf{y}
$$

積分内の第一項と第二項だけ取ると

$$
V(\mathbf{x}) \approx -\dfrac{G}{ |\mathbf{x}| } \left( M + \dfrac{\mathbf{x} \cdot \mathbf{y}_{0}}{|\mathbf{x}|^2} M \right)
$$

既に、地球の質量中心$${\mathbf{y_{0}}}$$を原点に取ると決めたので、

$$
V(\mathbf{x}) \approx -\dfrac{GM}{ |\mathbf{x}| }
$$

となって、地球の質量$${M}$$が原点に集中していると仮定した時の重力ポテンシャルが得られる。この式は、$${|\mathbf{x}|}$$が"地球半径"より大きい時のみ正しいが、そういうわけで、地球がどんな形状をしてるかに関係なく、質点近似が出てくる。

誤差項は、テイラー展開の3項目以降で評価される。3項目まで入れると

$$
V(\mathbf{x}) \approx -\dfrac{GM}{ |\mathbf{x}| } - \dfrac{G}{|\mathbf{x}|} \displaystyle \int_{\mathrm{Earth}} \left( \dfrac{|\mathbf{y}|^2}{|\mathbf{x}|^2} \dfrac{3 (\cos \theta)^2 - 1}{2} \rho(\mathbf{y}) \right) d\mathbf{y}
$$

となる。2項目の積分を計算するためには、変数$${\mathbf{x}}$$と$${\mathbf{y}}$$をバラさないと難しい。球面調和関数の一般論を知ってれば楽できるけど、アホみたいに展開してみると

$$
|\mathbf{x}|^2|\mathbf{y}|^2 \dfrac{3 (\cos \theta)^2 -1}{2} = 3(x_1 x_2 y_1 y_2 + x_1 x_3 y_1 y_3 + x_2 x_3 y_2 y_3) - \dfrac{1}{2} (x_1^2(-2 y_1^2+y_2^2+y_3^2) + x_2^2(y_1^2-2 y_2^2+y_3^2) + x_3^2(y_1^2+y_2^2 - 2 y_3^2))
$$

のようになっている。項数が多くて複雑に見えるけど、以下の形の積分

$$
\displaystyle \int_{\mathrm{Earth}} y_{i} y_{j} \rho(\mathbf{y}) d\mathbf{y}
$$

を計算できれば十分ということになる。適当な一次結合を取ると、慣性モーメントが求まり、逆に慣性モーメントから、この積分を決定することもできる。

 

やや唐突ではあるけど、以下のように新しい量を定義する。

$$
Q_{ij} = \displaystyle \int_{\mathrm{Earth}} \dfrac{\rho(\mathbf{y})}{M} \left( y_{i} y_{j} - \dfrac{1}{3} |\mathbf{y}|^2 \delta_{ij} \right) d\mathbf{y}
$$

$$
K_{ij} = x_{i} x_{j} - \dfrac{1}{3} |\mathbf{x}|^2 \delta_{i}{j}
$$

これを使うと、以下のように整理される。

$$
V(\mathbf{x}) \approx -\dfrac{GM}{ |\mathbf{x}| } - \dfrac{3}{2} \dfrac{GM}{|\mathbf{x}|^3} \displaystyle \sum_{i,j} \dfrac{K_{i,j}}{|\mathbf{x}|^2} Q_{ij}
$$

 

地球の細かい凹凸のことは忘れて、自転軸周りで軸対称ということにする。自転軸方向を$${(0,0,1)}$$とする。この時は、$${y_1,y_2}$$の奇数次の項を含む積分は0になる。つまり、$${i \neq j}$$なら、

$$
\displaystyle \int_{\mathrm{Earth}} y_{i} y_{j} \rho(\mathbf{y}) d \mathbf{y} = 0
$$

となる。また

$$
\displaystyle \int \rho(\mathbf{y}) y_1^2 d\mathbf{y} = \displaystyle \int \rho(\mathbf{y}) y_2^2 d\mathbf{y}
$$

が成り立つ。$${y_2}$$を消去すると、以下の形になる。

$$
\dfrac{3}{2} \displaystyle \sum_{i,j} \dfrac{K_{i,j}}{|\mathbf{x}|^2} Q_{ij} = \dfrac{1}{2} (x_1^2+x_2^2-2 x_3^2) \displaystyle \int_{\mathrm{Earth}} (-y_1^2+y_3^2) \rho(\mathbf{y}) d\mathbf{y}
$$

再度、

$$
\displaystyle \int \rho(\mathbf{y}) y_1^2 d\mathbf{y} = \displaystyle \int \rho(\mathbf{y}) y_2^2 d\mathbf{y}
$$

を使うと、

$$
V(\mathbf{x}) \approx -\dfrac{GM}{ |\mathbf{x}| } - \dfrac{GM}{4|\mathbf{x}|^3} \dfrac{ (3 K_{33}) }{|\mathbf{x}|^2} (3 Q_{33})
$$

と整理される。他の項と高次の項で軸対称近似からのズレを評価できる(軸対称近似で消えない高次の項もあるが、ここでは省略)

 

最初に、地球の密度が一様として、$${Q_{33}}$$のオーダーを見積もる。赤道半径を$${a}$$、極半径を$${b}$$とすると、大学一年生レベルの積分計算で

$$
3 Q_{33} = \dfrac{2}{5} (b^2 - a^2)
$$

を得る。

 

改めて、

$$
V(\mathbf{x}) \approx -\dfrac{GM}{ |\mathbf{x}| } - \dfrac{GM}{4|\mathbf{x}|^3} \dfrac{ (3 K_{33}) }{|\mathbf{x}|^2} (3 Q_{33})
$$

に戻る。この式で例えば、月に働く重力ポテンシャルが、どのくらい影響を受けるか調べる。

$${b=a(1-f)}$$によって、扁平率$${f}$$を定義すると、$${f \approx 1/300}$$で

$$
3 Q_{33} = 0.4 a^2(f^2 - 2f) \approx (-0.0026622) a^2
$$

月までの距離は、およそ地球半径の60倍なので、$${|\mathbf{x}|=60a}$$の時を考える。

$$
V(\mathbf{x}) \approx -\dfrac{GM}{ 60 a } + (1.85 \times 10^{-7}) \dfrac{GM}{60 a} \dfrac{ (3 K_{33}) }{|\mathbf{x}|^2}
$$

となる。$${\dfrac{ (3 K_{33}) }{|\mathbf{x}|^2}}$$は最大で$${2}$$、最小で$${-1}$$であることに注意。実際は、白道面の傾斜から計算できるが省略。

球対称近似と軸対称近似の場合で、誤差は、$${10^{-7}}$$程度の水準で、これは、殆ど無視できるとして差し支えないだろう。

地球表面付近だと、重力ポテンシャルに1/1000くらいの影響は与えるので、これは測定できない水準ではない。が、赤道では、自転による遠心力の寄与が、重力の3.5/1000くらいあって、それよりは小さく、通常は問題にならない。

地球密度の二相近似

ところで、20世紀に明らかになったところでは、地球の密度は、それほど一様でない。地球の平均密度は、約$${5.51 \mathrm{g/cm^3}}$$であるが、表面に近いほど軽く、中心に沈むほど重いそうである。検索した所、地殻の密度が$${2.8 \mathrm{g/cm^3}}$$で、マントルの密度は$${3.4 \mathrm{g/cm^3}}$$、外核と内核は、$${10 \sim 13 \mathrm{g/cm^3}}$$というような数値が出てきた。理由は知らないけど、核と、それ以外で、大きな密度の差がある。

この結果を踏まえて、地球密度の単純な二相モデルを考えて、$${Q_{33}}$$を見積もることを考える。

$$
y_1 = a r \sin \phi \cos \lambda
$$

$$
y_2 = a r \sin \phi \sin \lambda
$$

$$
y_3 = b r \cos \phi
$$

という球面座標系を導入して、$${0 \leq r \leq 1}$$とする。地球密度の分布を

$$
\rho(\mathbf{y}) = \rho_{0} c(r)
$$

という形だとする。$${\rho_0}$$は平均密度で、

$$
c(r) = c_0 (0 \leq r \leq r_0)
$$

$$
c(r) = c_1 (r_0 \lt r \leq 1)
$$

と仮定する。$${c_{0}}$$は、核の密度と平均密度の比なので、大雑把には2くらい。$${c_1}$$は地殻やマントルの密度と平均密度の比なので、大体、0.6くらい。$${r_0}$$は核半径と地球半径の比で、核半径が3500kmくらいらしいから、0.55くらい。但し

$$
\displaystyle \int_{\mathrm{Earth}} \rho(\mathbf{y}) d \mathbf{y} = M
$$

という条件を満たす必要がある。計算すると

$$
(1-r_{0}^3)c_1 + r_{0}^3 c_0 = 1
$$

という当たり前の条件になる。$${r_{0}=0.55 , c_{0}=2.0 , c_{1}=0.6}$$を代入して左辺を計算すると、0.833となる。若干小さい。二相近似だから、こんなもんか

とりあえず、$${Q_{33}}$$の計算を行うと、

$$
3 Q_{33} = \dfrac{2}{5} (b^2 - a^2)(r_{0}^{5} c_{0} + (1-r_{0}^5) c_{1})
$$

仮に、特に根拠のない数値$${r_{0}=0.57 , c_{0}=2.5 , c_{1}=0.7}$$を代入してみると

$$
(1-r_{0}^3)c_1 + r_{0}^3 c_0 \approx 0.9985472
$$

$$
r_{0}^{5} c_{0} + (1-r_{0}^5) c_{1} \approx 0.793624
$$

などとなり、地球密度の非一様性の影響により、$${Q_{33}}$$が少し小さくなることが予想される。

 

地球内部の密度分布を詳細に決定するのは、現在の人類の知識では難しい。代わりに、地球密度の非一様性を考慮して

$$
{3 Q_{33} = -J_{2} a^2}
$$

によって、新しい無次元量$${J_{2}}$$を導入し、測定値から$${J_2}$$を決める。

 

平衡状態の液体の自由表面が等ポテンシャル面になることは知られている。

地球の海が平衡状態にあるかは疑問があるが、そう仮定すると、極と赤道上の海面が等ポテンシャルになる必要がある。極半径と赤道半径の差は、遠心力で生じていると考えられているから、重力ポテンシャルに遠心力ポテンシャルを足す必要がある。

すると、$${\omega}$$を自転速度として

$$
-\dfrac{GM}{b} + \dfrac{2GM a^2 J_2}{4 b^3} = -\dfrac{GM}{a} - \dfrac{GM a^2 J_2}{4 a^3} - \dfrac{1}{2}\omega^2 a^2
$$

$${b=a(1-f)}$$で$f$が小さければ、$${\dfrac{1}{b} \approx \dfrac{1}{a}(1+f)}$$なので

$$
-\dfrac{GM}{a}(1+f) + \dfrac{2 G M J_2}{4 a} (1+3f) \approx -\dfrac{GM}{a} - \dfrac{GM J_2}{4 a} - \dfrac{1}{2}\omega^2 a^2
$$

近似重力加速度$${g_{0} = \dfrac{GM}{a^2}}$$を使って整理すると、以下の関係式が得られる。

$$
g_{0}(2 f - \dfrac{3}{2} J_{2}) \approx \omega^2 a
$$

 

また、重力ポテンシャルから、極と赤道の重力加速度を計算する(赤道上では遠心力が働くことを忘れない)と

$$
g_{pole} = \dfrac{GM}{b^2} - \dfrac{3GMa^2 J_{2}}{2 b^4} \approx g_{0}(1+2f - \dfrac{3}{2} J_{2})
$$

$$
g_{eq} = \dfrac{GM}{a^2} + \dfrac{3GMJ_{2}}{4a^3} - \omega a^2 = g_{0}(1 -2 f + \dfrac{9}{4} J_{2})
$$

重力加速度の差は

$${g_{pole} - g_{eq} = g_{0}(4 f - \dfrac{15}{4} J_{2})}$$

 

$${4f - \dfrac{15}{4} J_{2} = A}$$

$${2f - \dfrac{3}{2} J_{2} = B}$$

とすると

$${J_{2} = -\dfrac{4}{3}(A - 2B)}$$

$${f = -A + \dfrac{5}{2} B}$$

と解ける。

$${a = 6378 \mathrm{km} , g_{0} = 9.8 \mathrm{m/s^2} , \omega = 2 \pi /86400 \mathrm{rad/s},  g_{pole}=9.832 \mathrm{m/s^2} , g_{eq}=9.780 \mathrm{m/s^2}}$$として計算すると

$${J_{2} = 0.00210}$$

$${1/f = 303.17}$$

扁平率は少しずれてるが、定義値$${1/f = 298.25722210}$$に近い。地球密度の一様性を仮定した時の$${J_{2}}$$は(計算で出した方の扁平率を使うと)

$${-0.4(f^2 - 2f) = 0.0026344}$$

で、上の$${J_{2}}$$は、その0.8倍くらいになってる。これは、地球密度の二相近似で予想された数値と概ね合っている。

 

以上で、球対称性よりは弱い軸対称性を仮定して、測定値から$${Q_{33}}$$の値を決定できた。

けれど当然、地球は、別に軸対称じゃないので、測定値から$${Q_{ij}}$$を全部、更には、もっと高次の項までも決定するのが理想的であり、地球各地の重力データを利用すれば、おそらくfittingで決定できるだろう。ここまでで意外と長くなったので、またの機会。