球面座標でLaplace方程式を丁寧に解く

はじめに

3次元球面座標で定義された有界なスカラー関数 $U(r,\theta ,\varphi )$ がLaplace方程式 ${\mathrm{\nabla }}^{2}U=0$ を満たすとき、$U$ は一般に

で書くことができる[ 小話1 ]。ここで、$P_n^m(\cos \theta )$ は第一種Legendre陪関数である。これを $U$ の球面調和関数展開という。
球面調和関数展開は電磁気学で多重極子を考えるときなどに登場する重要な概念である。本記事ではできるだけ議論を端折らずに上記の一般解を導出し、関連する数学の話をまとめる。

球面調和関数

調和関数とは、Laplace方程式の解となるスカラー関数 $f$ のことである。

特に、これを3次元球面座標系 $(r,\theta ,\varphi )$ において (変数分離を用いて) 解いた解 $f(r,\theta ,\varphi )=Y_n^m$ を球面調和関数という[ 小話2 ]。

実数値関数の範囲での球面調和関数の概形 *[1]

By Twistar48 - Own work , CC BY-SA 4.0, Link

球面調和関数を表すにはLegendre倍関数が不可欠である。Legendre陪関数とは、以下のLegendreの陪微分方程式の解である。

これが解を持つのは $m\le n$ のときである。Legendreの陪微分方程式は2階の常微分方程式なので、二つの線型独立な解として $P_n^m(x)$ と $Q_n^m(x)$ を持つ。物理学に耐えうる「良い性質」を満たす方を第一種Legendre陪関数 $P_n^m(x)$ と呼び、他方の $x=\pm 1$ で発散するものを第二種Legendre陪関数 $Q_n^m$ と呼ぶ。第一種Legendre関数は人気者なのだが、第二種の項は物理的要請 (連続性、有界性といった境界条件) によってほとんど消されてしまう。今回は $U$ が有界な場合を考えているので、登場するのは第一種のみである。

一般解の導出

球面座標における3次元Laplace方程式

$$\begin{array}{rl}0& ={\mathrm{\nabla }}^{2}U\\ & =\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2\frac{\partial U}{\partial r}\right)+\frac{1}{r^2\sin \theta }\frac{\partial }{\partial \theta }\left(\sin \theta \frac{\partial U}{\partial \theta }\right)+\frac{1}{r^2{\sin }^2\theta }\frac{{\partial }^{2}U}{\partial {\varphi }^2}.\end{array}$$

の一般解の形が $U(r,\theta ,\varphi )=R(r)\cdot \mathrm{\Theta }(\theta )\cdot \mathrm{\Phi }(\varphi )$ で書けると仮定して変数分離法で解く[ 小話3 ]。
方程式に代入して

$$\frac{\mathrm{\Theta }\mathrm{\Phi }}{r^2}\frac{d}{dr}\left(r^2\frac{dR}{dr}\right)+\frac{\mathrm{\Phi }R}{r^2\sin \theta }\frac{d}{d\theta }\left(\sin \theta \frac{d\mathrm{\Theta }}{d\theta }\right)+\frac{R\mathrm{\Theta }}{r^2{\sin }^2\theta }\frac{d^2\mathrm{\Phi }}{d{\varphi }^2}=0.$$

この両辺を $R\mathrm{\Theta }\mathrm{\Phi }/r^2$ で割ると[ 小話4 ]

$$\frac{1}{\mathrm{\Theta }\sin \theta }\frac{d}{d\theta }(\sin \theta \frac{d\mathrm{\Theta }}{d\theta })+\frac{1}{\mathrm{\Phi }{\sin }^2\theta }\frac{d^2\mathrm{\Phi }}{d{\varphi }^2}=-\frac{1}{R}\frac{d}{dr}\left(r^2\frac{dR}{dr}\right).$$

である。左辺は $\theta ,\varphi$ の関数、右辺は $r$ のみの関数となっているから、この等式が成り立つためには両辺が定数でなければならない。分離定数 $C$ を導入して、

$$\{\begin{array}{rlr}& -\frac{1}{R}\frac{d}{dr}\left(r^2\frac{dR}{dr}\right)=C.& \text{(R)}\\ & \frac{1}{\mathrm{\Theta }\sin \theta }\frac{d}{d\theta }\left(\sin \theta \frac{d\mathrm{\Theta }}{d\theta }\right)+\frac{1}{\mathrm{\Phi }{\sin }^2\theta }\frac{d^2\mathrm{\Phi }}{d{\varphi }^2}=C.& \text{(L)}\end{array}$$

を得る。

方程式(R) について

$$r^2\frac{d^{2}R}{d{r}^2}+2r\cdot \frac{dR}{dr}+CR=0.$$
である。まともに解くのは大変なので、解を推測して代入してみる。
式を見ると、「$r$ で2回微分して $r^2$ 倍した項」と「$r$ で1回微分して $r$ 倍した項」と「$r$ で0回微分して $r^0$ 倍した項」とが互いに打ち消しあっているから、これら3項の次数は等しいはずである。したがって、$R(r)$ は $r$ の $n$ 次関数 ($n$ は整数) である。そこで、$R(r)=A_1{r}^n$. ($A_1$ は定数) とおいて方程式(L) に代入すると、
$$\{n(n-1)+2n+C\}\cdot A_1{r}^n=0.$$
を得る。これが任意の $r\ge 0$ において成り立つためには、分離定数は
$$C=-(n^2-n+2n)=-n(n+1).$$
でなければならないとわかる[ 小話5 ]。改めて、方程式(L)
$$r^2\frac{d^{2}R}{d{r}^2}+2r\cdot \frac{dR}{dr}-n(n+1)R=0.$$
を解くと、この方程式は $R(r)=A_1{r}^n$. のほかに $R(r)=B_1{r}^{-(n+1)}$. ($B_1$ は定数) を解にもつ。$r^n$ と $r^{-(n+1)}$ は互いに線型独立な解だから、$R(r)$ の一般解はこれら二つの線型結合である[ 小話6 ]。
$$R(r)=A_1{r}^n+\frac{B_1}{r^{n+1}}.$$

方程式(L) について

$$\frac{1}{\mathrm{\Theta }\sin \theta }\frac{d}{d\theta }\left(\sin \theta \frac{d\mathrm{\Theta }}{d\theta }\right)+\frac{1}{\mathrm{\Phi }{\sin }^2\theta }\frac{d^2\mathrm{\Phi }}{d{\varphi }^2}=-n(n+1).$$
である。両辺を ${\sin }^2\theta$ 倍して、
$$-\frac{1}{\mathrm{\Phi }}\frac{d^2\mathrm{\Phi }}{d{\varphi }^2}=\frac{\sin \theta }{\mathrm{\Theta }}\frac{d}{d\theta }\left(\sin \theta \frac{d\mathrm{\Theta }}{d\theta }\right)+n(n+1){\sin }^2\theta .$$
この左辺は $\theta$ のみの関数、右辺は $\varphi$ のみの関数だから、この等式が成り立つためには両辺が定数でなければならない。分離定数 $D$ を導入して

$$\{\begin{array}{rlr}& -\frac{1}{\mathrm{\Phi }}\frac{d^2\mathrm{\Phi }}{d{\varphi }^2}=D.& \text{(LL)}\\ & \frac{\sin \theta }{\mathrm{\Theta }}\frac{d}{d\theta }\left(\sin \theta \frac{d\mathrm{\Theta }}{d\theta }\right)+n(n+1){\sin }^2\theta =D.& \text{(LR)}\end{array}$$
を得る。

---- 方程式(LL) について

$$\frac{d^2\mathrm{\Phi }}{d{\varphi }^2}=-D\mathrm{\Phi }.$$
の解は
$$\mathrm{\Phi }(\varphi )=A_2\cos \sqrt{D}\varphi +B_2\cos \sqrt{D}\varphi .$$
である。いま球面座標を考えているから、$\mathrm{\Phi }(\varphi +2\pi )=\mathrm{\Phi }(\varphi )$ が境界条件である。これを満たすためには、分離定数が $D=m^2,\sqrt{D}=m.$ ($m$ は $0$ 以上の整数) でなければならない。したがって、
$$\mathrm{\Phi }(\varphi )=A_2\cos m\varphi +B_2\cos m\varphi .$$
である。

---- 方程式(LR) について

$$\frac{\sin \theta }{\mathrm{\Theta }}\frac{d}{d\theta }\left(\sin \theta \frac{d\mathrm{\Theta }}{d\theta }\right)+n(n+1){\sin }^2\theta =m^2.$$
ここで、$x=\cos \theta$ とおき、$\mathrm{\Theta }(\theta )=\tilde{\mathrm{\Theta }}(\cos \theta )=\tilde{\mathrm{\Theta }}(x)$ と変数変換すると、
$$\{\begin{array}{r}\frac{d\mathrm{\Theta }}{d\theta }=\frac{dx}{d\theta }\frac{d\tilde{\mathrm{\Theta }}}{dx}=-\sin \theta \frac{d\tilde{\mathrm{\Theta }}}{dx}.\\ {\sin }^2\theta =1-{\cos }^2\theta =1-x^2.\end{array}$$
であるから、
$$\sin \theta \cdot (-\sin \theta )\frac{d}{dx}[\sin \theta \cdot (-\sin \theta )\frac{d\tilde{\mathrm{\Theta }}}{dx}]+\{n(n+1){\sin }^2\theta -m^2\}\tilde{\mathrm{\Theta }}=0.$$
整理すると
$$\frac{d}{dx}[(1-x^2)\frac{d\tilde{\mathrm{\Theta }}}{dx}]+[n(n+1)-\frac{m^2}{1-x^2}]\tilde{\mathrm{\Theta }}=0.$$
これはLegendreの陪微分方程式であり、$\mathrm{\Theta }(\theta )=\tilde{\mathrm{\Theta }}(\cos \theta )$ はLegendre陪関数
$$\mathrm{\Theta }(\theta )=P_n^m(\cos \theta ).$$
である。第二種Legendre陪関数 $Q_n^m(\cos \theta )$ は $\cos \theta =\pm 1$ のとき、すなわち北極 $(\theta =0)$ と南極 $(\theta =\varphi )$ において値が発散してしまうから、解として不適。

以上より、球面座標における3次元Laplace方程式の1つの解は

$$\begin{array}{rl}{Y}_n^m(r,\theta ,\varphi )& =R(r)\cdot \mathrm{\Theta }(\theta )\cdot \mathrm{\Phi }(\varphi )\\ & =\{(B_1{A}_2\cos m\varphi +B_1{B}_2\sin m\varphi )\frac{1}{r^{n+1}}+(A_1{A}_2\cos m\varphi +A_1{B}_2\sin m\varphi )r^n\}{P}_n^m(\cos \theta ).\end{array}$$

である。任意定数 $B_1{A}_2,B_1{B}_2,A_1{A}_2,A_1{B}_2$ は各 $n,m$ に対して境界条件を解くことによって定まるから、改めてそれぞれを $g_n^m,h_n^m,q_n^m,s_n^m$ とおくと、一般解はこれらの $n,m$ に関する和で表せる。

$$U(r,\theta ,\varphi )=\sum _{n=0}^{\mathrm{\infty }}\sum _{m=0}^n\{(g_n^m\cos m\varphi +h_n^m\sin m\varphi ){\left(\frac{1}{r}\right)}^{n+1}+(q_n^m\cos m\varphi +s_n^m\sin m\varphi )r^n\}{P}_n^m(\cos \theta ).$$

系が $z$ 軸対称 ($\varphi$ について回転対称) の場合

$U(r,\theta ,\varphi )$ の値は $\varphi$ に依存しないから、

$$\frac{dU}{d\varphi }=\frac{d\mathrm{\Phi }}{d\varphi }=0,\mathrm{\Phi }(\varphi )=\mathrm{const}.$$

である。変数分離して得られた式(LL) において

$$\frac{1}{\mathrm{\Phi }}\frac{d^2\mathrm{\Phi }}{d{\varphi }^2}=-m^2=0.$$

より、$m=0$ である。このとき、$\mathrm{\Theta }(\theta )$ に関するLegendreの陪微分方程式は

となる (ただし $x=\cos \theta$ とした)。これはLegendreの微分方程式と呼ばれ、その解をLegendre関数という。特に、$x=\pm 1$ において発散しない解 $P_n^0(x)=P_n(x)$ をLegendre多項式という[ 小話7 ]。ポテンシャルの一般解は

となる。これはポテンシャルの多重極展開の式 (単極子 + 双極子 + 四重極子 + 八重極子 + ...) にあたる。特に、$g_0^0$ がモノポールに、$q_0^0$ がポテンシャルの基準に、$g_1^0,q_1^0$ が双極子に、それぞれ対応する。

こんな展開いつ使うの？

地球の磁場は球面調和関数展開で表すことができる。
地球の中心を原点にとった球面座標で考える。$\theta ,\varphi$ がそれぞれ緯度、経度にあたる。地球の磁場 $\overrightarrow{B}(r,\theta ,\varphi )$ は、絶縁体中 (宇宙) で時間変化しない磁場とみなすことができる[ 小話8 ]。Maxwell方程式は

$$\mathrm{div}\overrightarrow{B}=0,\mathrm{rot}\overrightarrow{B}=\overrightarrow{0}.$$

となる。$\overrightarrow{B}(r,\theta ,\varphi )$ は渦なしの静ベクトル場だから、あるスカラー関数 $U(r,\theta ,\varphi )$ が存在して、

$$\overrightarrow{B}=-\mathrm{grad}U.$$

と表せる。このときの $U(r,\theta ,\varphi )$ を磁気ポテンシャルという。
磁場の発散はゼロであるから[ 小話9 ]、磁気ポテンシャル $U$ は

$$\mathrm{div}\overrightarrow{B}=-\mathrm{div}\mathrm{grad}U=-{\mathrm{\nabla }}^{2}U=0.$$

すなわち、Laplace方程式

$${\mathrm{\nabla }}^{2}U=0.$$

の解である。したがって、地磁気のポテンシャルは球面調和関数の和で表せる。このとき、$g_n^m,h_n^m$ は内部Gauss係数、$q_n^m,s_n^m$ は外部Gauss係数と呼ばれる。数学への多大な貢献で知られるガウスは、それまであまり調べられていなかった地磁気の「強さ」を研究した人物でもあった。地球は北極をS極、南極をN極とする棒磁石でおおむね近似できると言われているが、実際には地磁気の生成はそれほど単純ではないため、z軸対称な系ではない。

球面座標におけるその他の2階偏微分方程式

Helmholtz方程式

波動方程式や拡散方程式を球面座標で解こうとすると、以下のHelmholtz方程式が現れる。

これは $k=0$ でLaplace方程式になる。この一般解は

$$A(r,\theta ,\varphi )=\sum _k\sum _{l=0}^{\mathrm{\infty }}\sum _{m=-l}^l\{a_{lm}{j}_l(kr)+b_{lm}{n}_l(kr)\}\cdot Y_l^m(\theta ,\varphi ).$$

である。ここで、$\theta ,\varphi$ 方向の変動を展開する $Y_l^m(\theta ,\varphi )$ は球面調和関数である。$r$ 方向の変動を展開する $j_l(kr),n_l(kr)$ は球Bessel関数と呼ばれる。

球面座標におけるシュレーディンガー方程式

水素原子中の電子の波動関数を求めたいとき、以下のシュレーディンガー方程式を解くことになる。

これは $\psi (r,\theta ,\varphi )=R(r)\cdot Y(\theta ,\varphi )$ と変数分離して解くことができて、$\theta ,\varphi$ 方向の変動 $Y(\theta ,\varphi )$ は球面調和関数になる。$r$ 方向の変動はラゲールの陪多項式で表すことができる。

まとめ

3次元球面座標で定義されたスカラー関数 $U(r,\theta ,\varphi )$ がLaplace方程式 ${\mathrm{\nabla }}^{2}U=0$ を満たすとき、$U$ は球面調和関数の和で表すことができる。特に z軸対称 ($\varphi$ に関する回転対称) な場合には、解は多重極子 (単極子 + 双極子 + 四重極子 + ...) の和の形で表せる。Laplace方程式以外にも、2階偏微分方程式を球面座標で解くと$\theta ,\varphi$ 方向の変動は球面調和関数で展開できる場合が多い。

小話

あまり追いきれていない細かい話をメモとしてまとめた。真偽不明とまでは行かないものの、要出典の話が多い。

小話1

$r^n,1/r^{n+1}$ といったべき関数によるテイラー展開 (Laurent展開と言うべき？) のような側面と、$\cos m\varphi ,\sin m\varphi$ といった周波数成分によるフーリエ展開のような側面の両方を併せ持つ点が個人的な推しポイント。

小話2

球面調和関数の定義は他にも「$n$次元Laplace方程式の解のうち単位球面上に限った関数」というものもある。

小話3

今回のLaplace方程式だけでなく、波動方程式とかシュレーディンガー方程式とか、物理学で登場する偏微分方程式ではしばしばこの「一般解が変数分離できると仮定」して解く変数分離法が用いられる。これで実際うまく行っているから物理的には良いのだが、個人的には「変数分離された解を仮定して得られた解が、本当に『その微分方程式の解を全て表すことのできる一般解』である保証はないのではないか」という疑問をずっと拭えないでいる。まあ、物理で登場するような「性質のいい」偏微分方程式には常微分方程式と同様に解の存在と一意性が成り立つ[要出典]ので、積分定数の数が足りていれば大丈夫と考えて良いのだろうとは思うのだが...。今回のケースでは一般解の定数 $g_n^m,h_n^m,q_n^m,s_n^m$ が無限個ある ($n$ は $1$ から $\mathrm{\infty }$ まで動く) ので、どう理解すべきか悩ましい。

小話4

両辺を $R\mathrm{\Theta }\mathrm{\Phi }/r^2(=U/r^2)$ で割っているが、$U=0$ となる点においてゼロ除算をしていることにはならないのだろうか。どうやら $U$ の零点が高々可算個であれば、測度論を用いて「ほとんどいたるところで」大丈夫ということになるらしい[要出典]。

小話5

この段階では、厳密には $n$ は一般に実数の値を取ってよい。物理の文脈でLaplace方程式を解こうとすると、物理的に満たされるべき条件によって結局 $n$ が自然数に限られることが多い[要出典]。分離定数 $C$ は、実際には実数全体を動けるほど自由な定数ではなく、せいぜい $C\le 1/4$. である。

小話6

しれっと $n>0$ であるかのような解の書き方をしたが、$n<0$ でも同型の解になる。$n<0$ のとき、求めた解 $R(r)$ は

$$R(r)=\frac{A_1}{r^{-n}}+B_1{r}^{-n+1}.$$

となるが、改めて $l=-n+1(>0)$ と置き直すことで

$$R(r)=B_1{r}^l+\frac{A_1}{r^{l+1}}.$$

の形に表すことができる。

小話7

Legendre多項式 $P_n(x)$ は $n$ 次多項式であって

$$P_n(x)=\frac{1}{2^{n}n!}\frac{d^n}{d{x}^n}[(x^2-1{)}^n].$$

と表すことができる (ロドリゲスの公式)。これを用いて Legendre陪関数は

$$P_n^m(x)=(1-x^2{)}^{m/2}\frac{d^m{P}_n(x)}{d{x}^m}.$$

と表される。

小話8

厳密には地球の磁場も変化している。地球の磁場自体が数年のスケールでゆったりと変化していたり、地球の高層では磁気嵐が発生したりしている。しかし、それでは $\mathrm{rot}\overrightarrow{B}=0$ の仮定が満たされていないではないか！ということで、$\mathrm{rot}\overrightarrow{B}$ の効果を無視できる時間スケールを見積もってみよう。Maxwell方程式は

$$\mathrm{rot}\overrightarrow{B}={\mu }_0\overrightarrow{j}+{\mu }_0\frac{\partial \overrightarrow{D}}{\partial t}$$

であり、地球電磁気では

• Ohmの法則: $\overrightarrow{j}=\sigma \overrightarrow{E}$
  電気伝導度 $\sigma \approx {10}^{-15}[/\mathrm{\Omega }\mathrm{m}={\mathrm{Fs}}^{-1}{\mathrm{m}}^{-1}]$
  大気電流 $j=\sigma E\approx {10}^{-12}$[A/m]

• 物性: $\overrightarrow{D}=\epsilon \overrightarrow{E}$
  大気の誘電率 $\epsilon \approx {\epsilon }_0\approx {10}^{-11}$[F/m]

であるTY。これらを方程式に代入すると

$$\mathrm{rot}\overrightarrow{B}=(1+\frac{{\epsilon }_0}{\sigma }\frac{\partial }{\partial t}){\mu }_0\sigma \overrightarrow{E}$$

右辺の最後の ${\mu }_0\sigma \overrightarrow{E}$ の大きさは、真空の透磁率${\mu }_0\approx {10}^{-6}$[N/A] と、上で示した大気電流の大きさより、${\mu }_0\sigma E\approx {10}^{-18}$[N/m] くらいのオーダーになる。したがって、$\mathrm{rot}\overrightarrow{B}$ の効果が無視できるかどうかは $\frac{{\epsilon }_0}{\sigma }\frac{\partial }{\partial t}$ の大きさ次第であり、例えばこれが $O(1)$ となる

$$\partial t\approx \frac{\sigma }{{\epsilon }_0}\approx {10}^4\sec <1\mathrm{day}$$

くらいの時間スケールであれば、$\mathrm{rot}\overrightarrow{B}\approx 0$ とみなすことができる。すなわち、「◯月◯日のガウス係数はAで、その翌日のガウス係数はBだ」という議論は可能だが、「◯時◯分のガウス係数はコレで、その1時間後にはこう変化する」という主張は、$\mathrm{rot}\overrightarrow{B}=0$ の仮定が成り立たない可能性があり、危険な議論である。
近年のデータ社会の中で、物理的に意義のあるデータのサンプリング周期を検討することは割と重要な要素だと思うのだが、この点に触れている文献はあまり見られないので、ここに記した。

小話9

一般的な静電場とは異なり静磁場が球面調和関数で展開できるのは、単磁荷 (モノポール) が存在しないからだと言える。もし大きめのモノポールを作ってしまった方がおられるならば、地磁気に影響しないよう地球磁気圏の外側 (高度約 1,000,000 km 以上) へ移動して実験していただきたい。