球面調和関数

当サイトではアフィリエイトプログラムを利用して商品を紹介しています。

球面調和関数

ラプラス方程式 $\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)

Matplotlib 球面調和関数(電子軌道)の図

コメント