球面調和関数

球面調和関数

　ラプラス方程式 $\nabla^2u=0$ を球座標で $u=R(r)Y(\theta,\phi)$ の形に変数分離して解いたときの角度部分の関数を球面調和関数 (spherical harmonics)とよびます。
 
具体的には ルジャンドル多項式 $P_n^m$ を用いて
 \[Y_l^m(\theta, \phi)=(-1)^\frac{l+|m|}{2}\sqrt{\frac{(2l+1)(l-|m|)!}{4\pi(n+|m|)!}}P_l^m(\cos\theta)e^{im\phi}\tag{1}\]
のように表されます。$l$ は $0$ 以上の整数、$m$ は $-l\leq m\leq l$ を満たす整数です。
 
$l=0,\ 1,\ 2$ について書き下すと以下のようになります。
 \[\begin{align*}&Y_0^0(\theta, \phi)=\frac{1}{\sqrt{4\pi}}\tag{2}\\[6pt]&Y_1^0(\theta, \phi)={\sqrt{\frac{3}4\pi}}\cos\theta\tag{3}\\[6pt]&Y_1^{\pm 1}(\theta, \phi)=\mp{\sqrt\frac{3}{8\pi}}\sin\theta e^{\pm i\phi}\tag{4}\\[6pt]&Y_2^{0}(\theta, \phi)={\sqrt\frac{5}{16\pi}}(3\cos^2\theta-1)\tag{5}\\[6pt]&Y_2^{\pm 1}(\theta, \phi)=\mp{\sqrt\frac{15}{8\pi}}\sin\theta\cos\theta e^{\pm i\phi}\tag{6}\\[6pt]&Y_2^{\pm 2}(\theta, \phi)=\sqrt{\frac{15}{32\pi}}\sin^2\theta e^{\pm2i\phi}\tag{7}\end{align*}\]

中心力場における粒子の存在確率

量子力学では中心力場のシュレーディンガー方程式
 \[\left[-\frac{\hbar^2}{2m}\nabla^2+V(r)\right]\phi(\boldsymbol{r})=E\,\phi(\boldsymbol{r})\tag{8}\]
を満たす波動関数 $\varphi(r,\theta,\phi)=R(r)Y_l^m(\theta, \phi)$ の形で 球面調和関数 が登場します。このとき、$l$ は方位量子数、$m$ は磁気量子数とよばれます。
 
$|Y_l^m(\theta, \phi)|^2$ は粒子の存在確率の角度依存性を示しますが、$e^{mi\phi}$ の部分は消えるので、存在確率は $\theta$ のみに依存することがわかります。
 
球面調和関数の詳細については『物理のための応用数学』などを参照してください。
 

scipy.special.sph_harm()

scipy.special.sph_harm(m, n, theta, phi) は球面調和関数
 \[Y_n^m(\theta, \phi)=(-1)^\frac{n+|m|}{2}\sqrt{\frac{(2n+1)(n-|m|)!}{4\pi(n+|m|)!}}P_n^m(\cos\phi)e^{im\theta}\tag{5}\]
を計算します。戻り値は複素数です。球面調和関数の一般的な表式 (1) とは角度変数 $\theta$ と $\phi$ が逆になっていることに注意してください。(1) の表式に合わせるならば、(m, l, phi, theta) の順に引数を渡すことになります。
 
(2) にあるように、$Y_0^0$ が $\phi$ や $\theta$ によらず $1/\sqrt{4\pi}$ であることを確認してみましょう。

# PYTHON_SPHERICA_HARMONICS_01

import numpy as np

# scipy.specialから球面調和関数をインポート
from scipy.special import sph_harm

# φ,θ=0,pi/2,piについて球面調和関数Y(m=0,l=0)を計算
for k in range(4):
    y = sph_harm(0, 0, k*np.pi/2, k*np.pi/2)
    print(y)

# (0.28209479177387814+0j)
# (0.28209479177387814+0j)
# (0.28209479177387814+0j)
# (0.28209479177387814+0j)

中心力場において極方程式 $r=|Y_l^m(\theta, \phi)|$ は軌道の形を表します。
極方程式を球座標形式から直交座標形式に変換する orbital() 関数を定義して、$l=3,\ m=1$ の軌道を描画してみます。

# PYTHON_SPHERICAL_HARMONICS_02

import numpy as np
from scipy.special import sph_harm
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# 軌道関数
# m：磁気量子数
# l：方位量子数
# p：x軸からの偏角φ
# t：z軸からの偏角θ
def orbital(m, l, p, t):
    r = np.abs(sph_harm(m, l, p, t))
    x = r*np.sin(t)*np.cos(p)
    y = r*np.sin(t)*np.sin(p)
    z = r*np.cos(t)
    return x, y, z

# FigureとAxes
fig = plt.figure(figsize = (8, 6))
ax = fig.add_subplot(111, projection='3d')

# 軸ラベルを設定
ax.set_xlabel("x", size = 16)
ax.set_ylabel("y", size = 16)
ax.set_zlabel("z", size = 16)

# アングルを設定
ax.view_init(45, 45)

# theta,phiの格子点データ
t = np.linspace(0, np.pi, 257)
p = np.linspace(0, 2*np.pi, 257)

# 格子点の作成
tt, pp = np.meshgrid(t, p)

# l=3,m=1軌道
x, y, z = orbital(1, 3, pp, tt)

# 軌道を描画
ax.plot_surface(x, y, z, color="darkorange",
                edgecolor="black", linewidth=0.4)