読者です 読者をやめる 読者になる 読者になる

倭算数理研究所

科学・数学・学習関連の記事を、「倭マン日記」とは別に書いていくのだ!

球対称な質量分布が作る重力ポテンシャル

3次元 極座標 重力 調和振動子

ちょっと所用で重力ポテンシャルを計算する必要ができたので、学部の演習問題でやった積分の計算を久し振りにやり直してみました。 どちらかというと電磁気学で静電ポテンシャルの計算としてやった記憶の方が強いですが、まぁ全く同じ計算でした。

質点と質量分布が作る重力ポテンシャル

まずは質点が作る重力ポテンシャル。 位置  { \textbf{r}' } にある質量  { m } の質点が位置  { \textbf{r} } に作る重力ポテンシャル  { \phi_{\textbf{r}'}(\textbf{r}) }

  { \displaystyle\begin{align*}
  \phi_{\textbf{r}'}(\textbf{r}) = -G \frac{m}{\left|\textbf{r} - \textbf{r}'\right|}
\end{align*}}

で与えられます(無限遠点を基準とする)。 ここで  { G }万有引力定数(ニュートン重力定数)です:

  { \displaystyle\begin{align*}
  G \fallingdotseq 6.67 \times 10^{-11} \quad [\textrm{m}^3 / (\textrm{kg} \cdot \textrm{s}^2)]
\end{align*}}

次は質量分布が作る重力ポテンシャル。 質量分布が  { \rho(\textbf{r}) } で与えられているとき、この分布が作る重力ポテンシャル  { \Phi(\textbf{r}) } は、質点の重力ポテンシャルで質量を質量分布に置き換えて、分布を積分すれば求められます:

  { \displaystyle\begin{align*}
  \Phi(\textbf{r})
    &= -G \int \frac{\rho(\textbf{r}')}{\left|\textbf{r} - \textbf{r}'\right|} dV'
\end{align*}}

まぁ、重力ポテンシャルは重ね合わせができるよってくらいの話ですが*1。 さてこれを踏まえて、球対称な質量分布が作る重力ポテンシャルを計算してみましょう。

球対称な質量分布が作る重力ポテンシャル

(原点に関して)球対称な質量分布とは、分布関数が位置ベクトル  { \textbf{r} } の動径成分  { r } のみに依存する分布です:

  { \displaystyle\begin{align*}
  \rho(\textbf{r}) = \rho(r)
\end{align*}}

極座標の角度  { \theta,\,\varphi } には依存しない、と言った方が分かりやすいですかね。 このとき、前節の質量分布が作る重力ポテンシャル  { \Phi(\textbf{r}) } の式で、角度変数の積分が実行できます(3次元の体積要素は昔の記事で導いたことあったなぁ):

  { \displaystyle\begin{align*}
  \Phi(\textbf{r})
    &= - G \int \frac{\rho(\textbf{r}')}{\left|\textbf{r} - \textbf{r}'\right|} dV' \\
    &= - G \int_0^{2\pi} \int_0^\pi \int_0^\infty \frac{\rho(r')}{|\textbf{r} - \textbf{r}'|} r'^2
      \sin\theta' dr'd\theta' d\varphi' \\
    &= - G \int_0^{2\pi} \int_0^\pi \int_0^\infty \frac{\rho(r')}{\sqrt{r^2 - 2rr' \cos\theta' + r'^2}}
      r'^2 \sin\theta' dr'd\theta' d\varphi' \\
    &= - 2\pi G \int_0^\pi \int_0^\infty \frac{\rho(r')}{\sqrt{r^2 - 2rr' \cos\theta' + r'^2}} r'^2 \sin\theta' dr'd\theta' 
      \qquad(\varphi' \, \textrm{積分実行})\\
    &= - 2\pi G \int_0^\infty \rho(r') \int_0^\pi \frac{\sin\theta' d\theta'}{\sqrt{r^2 - 2rr' \cos\theta' + r'^2}} r'^2 dr'
\end{align*}}

ごちゃごちゃしてきたので  { \theta' } 積分だけ抜き出して実行しましょう:

  { \displaystyle\begin{align*}
  \int_0^\pi \frac{\sin\theta' d\theta'}{\sqrt{r^2 - 2rr' \cos\theta' + r'^2}}
    &= \int_{-1}^1 \frac{d(\cos\theta')}{\sqrt{r^2 - 2rr' \cos\theta' + r'^2}} \\
    &= \left[-\frac{1}{rr'}\sqrt{r^2 - 2rr' \cos\theta' + r'^2}\right]_{-1}^1 \\
    &= -\frac{1}{rr'}\sqrt{r^2 - 2rr' + r'^2} + \frac{1}{rr'}\sqrt{r^2 + 2rr' + r'^2} \\
    &= -\frac{1}{rr'}\left(|r-r'| - |r+r'|\right)  \qquad \left(\because \sqrt{A^2} = |A|\right) \\
    &= \begin{cases}
        \dfrac{2}{r} & (0 \leqq r' \leqq r) \\[2mm]
        \dfrac{2}{r'} & (r' > r)
      \end{cases}
\end{align*}}

最初の形からは想像できないくらい(天才にはできるのか?)綺麗な形にまとまりました。  { r } { r' } の大小関係で形が変わるのはちょっと面倒そうだけど。 さて、これを元の  { \Phi(\textbf{r}) } の式に戻して整理すると以下の様になります:

  { \displaystyle\begin{align*}
  \Phi(\textbf{r})
    &= -\dfrac{4\pi G}{r} \int_0^r \rho(r') r'^2 dr'
      - 4\pi G \int_r^\infty \rho(r') r' dr'
\end{align*}}

抜き出して計算した部分の表式を使うために、 { r' } 積分積分範囲を  { r } を境にして2つに分け、それぞれの範囲の表式を代入しています。 これで積分できるところはしてしまったので、これ以上計算をする( { r' } 積分を実行する)には質量分布の具体的な式が必要となります。

球体内に一様分布する質量分布

原点を中心とし、半径が  { R } の球体内に一様に質量が分布している系を考えましょう。 質量分布は以下で与えられます:

  { \displaystyle\begin{align*}
  \rho(r)
    &= \begin{cases}
        \rho & (0 \leqq r \leqq R) \\[1mm]
        0 & (r > R)
      \end{cases} \\[2mm]
    &= \rho\,\theta(R-r)
\end{align*}}

最後の式は階段関数(Heaviside 関数)  { \theta(x) } を使って書いてみただけで特に深い意味はありません(この  { \theta }極座標の角度とは関係ないよ)。 さて、この分布の式を前節の  { \Phi(\textbf{r}) } の式に代入して  { r' } 積分を実行していきましょう。 注意が必要なのは、今度は  { r } { R } の大小関係で場合分けする必要があるところ。

 { r > R } のとき
まずは  { r > R } のとき。 これは質量が分布している球体よりも外側にできる重力ポテンシャルを与えます:

  { \displaystyle\begin{align*}
  \Phi(\textbf{r})
    &= -\dfrac{4\pi G}{r} \int_0^R\rho r'^2 dr' \\
    &= - \frac{4\pi G}{r} \cdot \frac{\rho R^3}{3} \\
    &= - G\frac{M}{r} \qquad \left(M = \rho \cdot \tfrac{4}{3} \pi R^3\right)
\end{align*}}

最後の式の  { M } は球体内に分布している全質量です。 この結果を見ると、球体の外側では球の中心に全質量が集まったときにできる質点の重力ポテンシャルと同じ形をしていることが分かります。

 { r \leqq R } のとき
次は  { r \leqq R } のとき。 これは球体の内側の重力ポテンシャルを与えます:

  { \displaystyle\begin{align*}
  \Phi(\textbf{r})
    &= -\dfrac{4\pi G}{r} \int_0^r \rho r'^2 dr' - 4\pi G \int_r^R \rho r' dr'\\
    &= - \frac{4\pi G}{r} \cdot \frac{\rho r^3}{3} - 2 \pi G \rho \left(R^2 - r^2\right) \\
    &= \frac{2\pi G \rho}{3} r^2 - 2\pi G \rho R^2
\end{align*}}

これは  { r = R } で上記の球体外のポテンシャルと連続的につながっています。 第2項は単なる定数なのでまぁいいでしょう。 第1項は  { r^2 } に比例しているので、これは調和振動子のポテンシャルに等しいことが分かります。

球面上に一様分布する質量分布

次は原点を中心とし、半径が  { R } の球面上に一様に質量が分布している系を考えましょう。 質量分布は以下で与えられます:

  { \displaystyle\begin{align*}
  \rho(r) = \rho\,\delta(r-R)
\end{align*}}

 { \delta(x) }ディラックデルタ関数です*2。 では、この分布の式を  { \Phi(\textbf{r}) } の式に代入して  { r' } 積分を実行していきましょう。 やることは球体内の一様分布の場合と同じです(というかもう少し簡単)。

 { r > R } のとき
まずは  { r > R } の球殻外でのポテンシャル:

  { \displaystyle\begin{align*}
  \Phi(\textbf{r})
    &= -\dfrac{4\pi G}{r} \int_0^r\rho \delta(r'-R) r'^2 dr' \\
    &= -\dfrac{4\pi G \rho}{r} R^2 \\
    &= - G\frac{M}{r} \qquad \left(M = \rho \cdot 4 \pi R^2\right)
\end{align*}}

 { M } の定義は前節のものと異なっていますが、やはり球面上に分布している全質量です。 よって、今の場合も球殻外では球の中心に全質量が集まったときにできる質点の重力ポテンシャルと同じ形をしていることが分かります。

 { r \leqq R } のとき
次は  { r \leqq R } の球殻内の重力ポテンシャル:

  { \displaystyle\begin{align*}
  \Phi(\textbf{r})
    &= - 4\pi G \int_r^R \rho \delta(r'-R) r' dr'\\
    &= - 4 \pi G \rho R
\end{align*}}

これはやはり  { r = R } で上記の球殻外のポテンシャルと連続的につながっています。 また、 { r } に依存しなくなっているので、球殻内では重力が働かないことが分かります。

修正
最初の方のポテンシャルの表式で負符号が抜けてたので修正しました。

*1:質点の重力ポテンシャルは、某方程式の Green 関数だよ、とも言える。

*2:引数が0となる場所が積分範囲に入っていれば1、そうでなければ0となる(超)関数。