中心力場（万有引力、クーロン力など）を扱う問題は球対称性を持つため、球座標系で記述すると計算が大幅に簡単になります。ラプラス方程式の球面解（球面調和関数）や水素原子の波動関数も球座標系の産物です。
本記事では、球座標系における微分演算子の表現を計量係数から体系的に導出し、物理応用を解説します。
本記事の内容
- 球座標系の定義と座標変換
- 計量係数と微小体積要素
- 勾配・発散・回転・ラプラシアンの導出
- ラプラス方程式の球面解
- 中心力問題への応用
- Pythonでの可視化
前提知識
球座標系の定義
球座標系 $(r, \theta, \varphi)$ と直交座標系 $(x, y, z)$ の関係は：
$$ \boxed{x = r\sin\theta\cos\varphi, \quad y = r\sin\theta\sin\varphi, \quad z = r\cos\theta} $$
ここで各座標の範囲は： – $r \geq 0$: 原点からの距離 – $0 \leq \theta \leq \pi$: 天頂角（$z$ 軸からの角度） – $0 \leq \varphi < 2\pi$: 方位角（$x$ 軸からの角度）
逆変換は：
$$ r = \sqrt{x^2 + y^2 + z^2}, \quad \theta = \arccos\frac{z}{r}, \quad \varphi = \arctan\frac{y}{x} $$
単位ベクトル
位置ベクトル $\bm{r} = r\sin\theta\cos\varphi\,\bm{e}_x + r\sin\theta\sin\varphi\,\bm{e}_y + r\cos\theta\,\bm{e}_z$ を各座標で微分し、正規化します。
$$ \bm{e}_r = \sin\theta\cos\varphi\,\bm{e}_x + \sin\theta\sin\varphi\,\bm{e}_y + \cos\theta\,\bm{e}_z $$
$$ \bm{e}_\theta = \cos\theta\cos\varphi\,\bm{e}_x + \cos\theta\sin\varphi\,\bm{e}_y – \sin\theta\,\bm{e}_z $$
$$ \bm{e}_\varphi = -\sin\varphi\,\bm{e}_x + \cos\varphi\,\bm{e}_y $$
これらは互いに直交し、$\bm{e}_r\times\bm{e}_\theta = \bm{e}_\varphi$ を満たします。
計量係数の導出
微小変位ベクトルを計算します：
$$ \frac{\partial\bm{r}}{\partial r} = \sin\theta\cos\varphi\,\bm{e}_x + \sin\theta\sin\varphi\,\bm{e}_y + \cos\theta\,\bm{e}_z = \bm{e}_r $$
$$ \frac{\partial\bm{r}}{\partial\theta} = r\cos\theta\cos\varphi\,\bm{e}_x + r\cos\theta\sin\varphi\,\bm{e}_y – r\sin\theta\,\bm{e}_z = r\,\bm{e}_\theta $$
$$ \frac{\partial\bm{r}}{\partial\varphi} = -r\sin\theta\sin\varphi\,\bm{e}_x + r\sin\theta\cos\varphi\,\bm{e}_y = r\sin\theta\,\bm{e}_\varphi $$
各偏微分の大きさが計量係数です：
$$ \boxed{h_r = 1, \quad h_\theta = r, \quad h_\varphi = r\sin\theta} $$
微小変位ベクトル：
$$ d\bm{r} = dr\,\bm{e}_r + r\,d\theta\,\bm{e}_\theta + r\sin\theta\,d\varphi\,\bm{e}_\varphi $$
微小体積要素：
$$ dV = h_r h_\theta h_\varphi\,dr\,d\theta\,d\varphi = r^2\sin\theta\,dr\,d\theta\,d\varphi $$
勾配（grad）
直交曲線座標系の勾配の一般公式に計量係数を代入します：
$$ \nabla f = \frac{1}{h_r}\frac{\partial f}{\partial r}\bm{e}_r + \frac{1}{h_\theta}\frac{\partial f}{\partial\theta}\bm{e}_\theta + \frac{1}{h_\varphi}\frac{\partial f}{\partial\varphi}\bm{e}_\varphi $$
$$ \boxed{\nabla f = \frac{\partial f}{\partial r}\bm{e}_r + \frac{1}{r}\frac{\partial f}{\partial\theta}\bm{e}_\theta + \frac{1}{r\sin\theta}\frac{\partial f}{\partial\varphi}\bm{e}_\varphi} $$
発散（div）
一般公式：
$$ \nabla\cdot\bm{A} = \frac{1}{h_r h_\theta h_\varphi}\left[\frac{\partial}{\partial r}(h_\theta h_\varphi A_r) + \frac{\partial}{\partial\theta}(h_r h_\varphi A_\theta) + \frac{\partial}{\partial\varphi}(h_r h_\theta A_\varphi)\right] $$
$h_r h_\theta h_\varphi = r^2\sin\theta$ を代入すると：
$$ \nabla\cdot\bm{A} = \frac{1}{r^2\sin\theta}\left[\frac{\partial}{\partial r}(r^2\sin\theta \cdot A_r) + \frac{\partial}{\partial\theta}(r\sin\theta \cdot A_\theta) + \frac{\partial}{\partial\varphi}(r \cdot A_\varphi)\right] $$
$\sin\theta$ は $r$ に依存しないので第1項から外せ、$r$ は $\theta$ に依存しないので第2項から外せます：
$$ \boxed{\nabla\cdot\bm{A} = \frac{1}{r^2}\frac{\partial}{\partial r}(r^2 A_r) + \frac{1}{r\sin\theta}\frac{\partial}{\partial\theta}(\sin\theta\, A_\theta) + \frac{1}{r\sin\theta}\frac{\partial A_\varphi}{\partial\varphi}} $$
回転（curl）
一般公式の行列式表現に計量係数を代入します：
$$ \nabla\times\bm{A} = \frac{1}{r^2\sin\theta}\begin{vmatrix} \bm{e}_r & r\,\bm{e}_\theta & r\sin\theta\,\bm{e}_\varphi \\ \frac{\partial}{\partial r} & \frac{\partial}{\partial\theta} & \frac{\partial}{\partial\varphi} \\ A_r & rA_\theta & r\sin\theta\, A_\varphi \end{vmatrix} $$
各成分を展開します：
$$ (\nabla\times\bm{A})_r = \frac{1}{r\sin\theta}\left[\frac{\partial}{\partial\theta}(\sin\theta\, A_\varphi) – \frac{\partial A_\theta}{\partial\varphi}\right] $$
$$ (\nabla\times\bm{A})_\theta = \frac{1}{r}\left[\frac{1}{\sin\theta}\frac{\partial A_r}{\partial\varphi} – \frac{\partial}{\partial r}(rA_\varphi)\right] $$
$$ (\nabla\times\bm{A})_\varphi = \frac{1}{r}\left[\frac{\partial}{\partial r}(rA_\theta) – \frac{\partial A_r}{\partial\theta}\right] $$
ラプラシアンの導出
$\nabla^2 f = \nabla\cdot(\nabla f)$ に勾配と発散の公式を代入します。
$\nabla f$ の各成分は $A_r = \frac{\partial f}{\partial r}$, $A_\theta = \frac{1}{r}\frac{\partial f}{\partial\theta}$, $A_\varphi = \frac{1}{r\sin\theta}\frac{\partial f}{\partial\varphi}$ です。
第1項（$r$ 成分）：
$$ \frac{1}{r^2}\frac{\partial}{\partial r}\left(r^2 \cdot \frac{\partial f}{\partial r}\right) = \frac{1}{r^2}\frac{\partial}{\partial r}\left(r^2\frac{\partial f}{\partial r}\right) $$
第2項（$\theta$ 成分）：
$$ \frac{1}{r\sin\theta}\frac{\partial}{\partial\theta}\left(\sin\theta \cdot \frac{1}{r}\frac{\partial f}{\partial\theta}\right) = \frac{1}{r^2\sin\theta}\frac{\partial}{\partial\theta}\left(\sin\theta\frac{\partial f}{\partial\theta}\right) $$
第3項（$\varphi$ 成分）：
$$ \frac{1}{r\sin\theta}\frac{\partial}{\partial\varphi}\left(\frac{1}{r\sin\theta}\frac{\partial f}{\partial\varphi}\right) = \frac{1}{r^2\sin^2\theta}\frac{\partial^2 f}{\partial\varphi^2} $$
したがって：
$$ \boxed{\nabla^2 f = \frac{1}{r^2}\frac{\partial}{\partial r}\left(r^2\frac{\partial f}{\partial r}\right) + \frac{1}{r^2\sin\theta}\frac{\partial}{\partial\theta}\left(\sin\theta\frac{\partial f}{\partial\theta}\right) + \frac{1}{r^2\sin^2\theta}\frac{\partial^2 f}{\partial\varphi^2}} $$
第1項を展開すると $\frac{\partial^2 f}{\partial r^2} + \frac{2}{r}\frac{\partial f}{\partial r}$ となります。
ラプラス方程式の球面解
球対称なラプラス方程式 $\nabla^2 f = 0$ を考えます。$f = f(r)$（角度に依存しない）とすると：
$$ \frac{1}{r^2}\frac{d}{dr}\left(r^2\frac{df}{dr}\right) = 0 $$
$$ r^2\frac{df}{dr} = C_1 \implies \frac{df}{dr} = \frac{C_1}{r^2} $$
$$ f(r) = -\frac{C_1}{r} + C_2 $$
これは点電荷のポテンシャル $V = \frac{q}{4\pi\epsilon_0 r}$ の形です。
変数分離と球面調和関数
$f(r,\theta,\varphi) = R(r)\Theta(\theta)\Phi(\varphi)$ と変数分離すると、角度部分の解として球面調和関数 $Y_l^m(\theta,\varphi)$ が現れ、動径部分は $R(r) = Ar^l + Br^{-(l+1)}$ となります。
一般解は：
$$ f(r,\theta,\varphi) = \sum_{l=0}^{\infty}\sum_{m=-l}^{l}\left(A_{lm}r^l + B_{lm}r^{-(l+1)}\right)Y_l^m(\theta,\varphi) $$
応用例: 点電荷のポテンシャルと電場
原点に点電荷 $q$ がある場合のポテンシャルは：
$$ V(r) = \frac{q}{4\pi\epsilon_0 r} $$
電場を球座標系の勾配で計算します：
$$ \bm{E} = -\nabla V = -\frac{\partial V}{\partial r}\bm{e}_r = \frac{q}{4\pi\epsilon_0 r^2}\bm{e}_r $$
発散を計算すると：
$$ \nabla\cdot\bm{E} = \frac{1}{r^2}\frac{\partial}{\partial r}\left(r^2 \cdot \frac{q}{4\pi\epsilon_0 r^2}\right) = \frac{1}{r^2}\frac{\partial}{\partial r}\left(\frac{q}{4\pi\epsilon_0}\right) = 0 \quad (r \neq 0) $$
原点以外では $\nabla\cdot\bm{E} = 0$ であり、ガウスの法則と整合します。
Pythonでの可視化
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure(figsize=(18, 5))
# --- (1) 球座標系の座標面 ---
ax1 = fig.add_subplot(131, projection='3d')
# r=const の面（球面）
u = np.linspace(0, 2*np.pi, 40)
v = np.linspace(0, np.pi, 20)
U, V = np.meshgrid(u, v)
for r_val in [0.5, 1.0, 1.5]:
X = r_val * np.sin(V) * np.cos(U)
Y = r_val * np.sin(V) * np.sin(U)
Z = r_val * np.cos(V)
ax1.plot_wireframe(X, Y, Z, alpha=0.2, color='blue', linewidth=0.5)
ax1.set_xlabel('$x$'); ax1.set_ylabel('$y$'); ax1.set_zlabel('$z$')
ax1.set_title('Spherical coordinate surfaces', fontsize=12)
# --- (2) 点電荷の電場（断面） ---
ax2 = fig.add_subplot(132)
r = np.linspace(0.3, 3.0, 12)
theta = np.linspace(0, 2*np.pi, 16, endpoint=False)
R, Th = np.meshgrid(r, theta)
X2 = R * np.cos(Th)
Y2 = R * np.sin(Th)
# E ∝ 1/r^2, 動径方向
Er = 1.0 / R**2
Ex = Er * np.cos(Th)
Ey = Er * np.sin(Th)
ax2.quiver(X2, Y2, Ex/np.max(Er), Ey/np.max(Er), color='blue', alpha=0.7)
ax2.plot(0, 0, 'ro', markersize=10, label='Point charge')
ax2.set_xlim(-3.5, 3.5); ax2.set_ylim(-3.5, 3.5)
ax2.set_aspect('equal')
ax2.set_title('Electric field of point charge', fontsize=12)
ax2.set_xlabel('$x$'); ax2.set_ylabel('$y$')
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)
# --- (3) ラプラス方程式の球面解 ---
ax3 = fig.add_subplot(133)
r_plot = np.linspace(0.5, 5.0, 200)
# 球対称解: f(r) = A/r + B
solutions = [
(1.0, 0.0, '$1/r$'),
(-1.0, 2.0, '$-1/r + 2$'),
(2.0, -1.0, '$2/r - 1$'),
]
for A, B, label in solutions:
f = A / r_plot + B
ax3.plot(r_plot, f, lw=2, label=label)
ax3.set_xlabel('$r$', fontsize=12)
ax3.set_ylabel('$f(r)$', fontsize=12)
ax3.set_title('Spherically symmetric solutions of $\\nabla^2 f = 0$', fontsize=12)
ax3.legend(fontsize=11)
ax3.grid(True, alpha=0.3)
ax3.set_ylim(-3, 5)
ax3.axhline(0, color='black', lw=0.5)
plt.tight_layout()
plt.savefig('spherical_coordinates.png', dpi=150, bbox_inches='tight')
plt.show()
球面調和関数の可視化も行います。
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import sph_harm
theta = np.linspace(0, np.pi, 100)
phi = np.linspace(0, 2*np.pi, 100)
Theta, Phi = np.meshgrid(theta, phi)
fig, axes = plt.subplots(2, 3, figsize=(16, 10),
subplot_kw={'projection': '3d'})
harmonics = [
(0, 0, '$Y_0^0$'), (1, 0, '$Y_1^0$'), (1, 1, '$\\mathrm{Re}(Y_1^1)$'),
(2, 0, '$Y_2^0$'), (2, 1, '$\\mathrm{Re}(Y_2^1)$'), (2, 2, '$\\mathrm{Re}(Y_2^2)$'),
]
for ax, (l, m, title) in zip(axes.flat, harmonics):
Y = sph_harm(m, l, Phi, Theta)
if m != 0:
Y = np.real(Y)
else:
Y = np.real(Y)
R = np.abs(Y)
X = R * np.sin(Theta) * np.cos(Phi)
Y_3d = R * np.sin(Theta) * np.sin(Phi)
Z = R * np.cos(Theta)
# 色で符号を表現
colors = np.real(sph_harm(m, l, Phi, Theta))
norm = plt.Normalize(vmin=colors.min(), vmax=colors.max())
ax.plot_surface(X, Y_3d, Z, facecolors=plt.cm.RdBu_r(norm(colors)),
alpha=0.8, linewidth=0)
ax.set_title(title, fontsize=14)
ax.set_xlim(-0.6, 0.6); ax.set_ylim(-0.6, 0.6); ax.set_zlim(-0.6, 0.6)
ax.set_xlabel('$x$'); ax.set_ylabel('$y$'); ax.set_zlabel('$z$')
plt.suptitle('Spherical harmonics $Y_l^m(\\theta, \\varphi)$', fontsize=16, y=1.02)
plt.tight_layout()
plt.savefig('spherical_harmonics.png', dpi=150, bbox_inches='tight')
plt.show()
まとめ
本記事では、球座標系 $(r, \theta, \varphi)$ での微分演算子を導出しました。
- 計量係数: $h_r = 1$, $h_\theta = r$, $h_\varphi = r\sin\theta$
- 勾配: $\nabla f = \frac{\partial f}{\partial r}\bm{e}_r + \frac{1}{r}\frac{\partial f}{\partial\theta}\bm{e}_\theta + \frac{1}{r\sin\theta}\frac{\partial f}{\partial\varphi}\bm{e}_\varphi$
- 発散: $\nabla\cdot\bm{A} = \frac{1}{r^2}\frac{\partial}{\partial r}(r^2 A_r) + \frac{1}{r\sin\theta}\frac{\partial}{\partial\theta}(\sin\theta\, A_\theta) + \frac{1}{r\sin\theta}\frac{\partial A_\varphi}{\partial\varphi}$
- ラプラシアン: $\frac{1}{r^2}\frac{\partial}{\partial r}\left(r^2\frac{\partial f}{\partial r}\right) + \frac{1}{r^2\sin\theta}\frac{\partial}{\partial\theta}\left(\sin\theta\frac{\partial f}{\partial\theta}\right) + \frac{1}{r^2\sin^2\theta}\frac{\partial^2 f}{\partial\varphi^2}$
- ラプラス方程式の球対称解: $f(r) = A/r + B$
- 球面調和関数 $Y_l^m(\theta,\varphi)$ が角度方向の解
次のステップとして、以下の記事も参考にしてください。